summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/matrixOperations/logm/cbal.c177
-rw-r--r--src/matrixOperations/logm/comqr3.c246
-rw-r--r--src/matrixOperations/logm/corth.c112
-rw-r--r--src/matrixOperations/logm/cortr.c72
-rw-r--r--src/matrixOperations/logm/logm_internal.h37
-rw-r--r--src/matrixOperations/logm/wbdiag.c306
6 files changed, 0 insertions, 950 deletions
diff --git a/src/matrixOperations/logm/cbal.c b/src/matrixOperations/logm/cbal.c
deleted file mode 100644
index ff866cd1..00000000
--- a/src/matrixOperations/logm/cbal.c
+++ /dev/null
@@ -1,177 +0,0 @@
-/*
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008 - INRIA - Arnaud TORSET
- *
- * This file must be used under the terms of the CeCILL.
- * This source file is licensed as described in the file COPYING, which
- * you should have received as part of this distribution. The terms
- * are also available at
- * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
- *
- */
-
-
- /*
- This is a transcription of the fortran subroutine cbal.f
- */
-
-#include "multiplication.h"
-#include "logm_internal.h"
-#include "abs.h"
-
-/******************/
-/* Sous-Fonctions */
-/******************/
-static int boucle_100(doubleComplex* in, int size, int l, int* j, int* m, int *iexc){
- int jj=0, i=0;
- for (jj=1;jj<=l;jj++){
- *j=l+1-jj;
- for(i=1;i<=l;i++){
- if(i!=(*j)) {
- if ( (zreals(in[((*j)-1)*size+i-1])!=0) || (zimags(in[((*j)-1)*size+i-1])!=0) ) break;
- }
- }
- if (i!=(l+1)) break;
- else {
- *m=l;
- *iexc=1;
- return 20;
- }
- }
- return 140;
-}
-
-
-static int boucle_140(doubleComplex* in, int size, int k, int l, int* j, int* m, int *iexc){
- int i=0;
- for ((*j)=k;(*j)<=l;(*j)++){
- for(i=k;i<=l;i++){
- if(i!=(*j)) {
- if ( (zreals(in[(i-1)*size+(*j)-1])!=0) || (zimags(in[(i-1)*size+(*j)-1])!=0) ) break;
- }
- }
- if (i!=(l+1)) break;
- else{
- *m=k;
- *iexc=2;
- return 20;
- }
- }
- return 170;
-}
-
-
-/***********************/
-/* Fonction principale */
-/***********************/
-
-void cbal(doubleComplex* in, int n, int* low, int* high, double* scale){
- const double radix=2;
- double b2=0, c=0, r=0, g=0, f=0, s=0;
- int k=0, l=0, continuer=1, j=0, i=0, iexc=0, noconv=0;
- int result, m;
- doubleComplex fCpx;
-
- b2=radix*radix;
- k=1;
- l=n;
-
- result= boucle_100(in,n,l,&j,&m, &iexc);
- while(continuer){
- if(result==20){
- scale[m-1] = j;
- if (j != m) {
- for(i=1; i<=l;i++){
- fCpx = in[(i-1)*n+j-1];
- in[(i-1)*n+j-1] = in[(i-1)*n+m-1];
- in[(i-1)*n+m-1] = fCpx;
- }
- for(i=k; i<=m; i++){
- fCpx = in[(j-1)*n+i-1];
- in[(j-1)*n+i-1] = in[(m-1)*n+i-1];
- in[(m-1)*n+i-1] = fCpx;
- }
- }
-
- if (iexc==1){
- if (l==1) {
- continuer=0;
- break;
- }
- else {
- l--;
- result = boucle_100(in,n,l,&j,&m, &iexc);
- }
- }
- else if (iexc==2){
- k++;
- result = boucle_140(in,n,k,l,&j,&m, &iexc);
- }
-
-
- }
-
- if (result == 140){
- result = boucle_140(in,n,k,l,&j,&m, &iexc);
- }
-
- if (result == 170){
- continuer=0;
- }
- }
-
- for (i=k;i<=l;i++){
- scale[i-1]=1;
- }
-
- do{
- noconv=0; /* noconv=false */
- for (i=k;i<=l;i++){
- c=0;
- r=0;
-
-
- for (j=k;j<=l;j++){
- if(j!=i){
- c += ( dabss(zreals(in[(i-1)*n+j-1])) + dabss(zimags(in[(i-1)*n+j-1])) );
- r += ( dabss(zreals(in[(j-1)*n+i-1])) + dabss(zimags(in[(j-1)*n+i-1])) );
- }
- }
-
- /* :::::::::: guard against zero c or r due to underflow :::::::::: */
- if ( (r!=0) && (c!=0) ){
- g=r/radix;
- f=1;
- s=c+r;
-
- while (c < g){
- f = f * radix;
- c = c * b2;
- }
-
- g = r * radix;
- while(c >= g) {
- f = f / radix;
- c = c / b2;
-
- }
- /* :::::::::: now balance :::::::::: */
- if (((c + r) / f) < 0.950 * s) {
- g = 1 / f;
- scale[i-1] = scale[i-1] * f;
- noconv = 1;
-
- for (j = k; j <= n; j++)
- in[(j-1)*n+i-1] = zmuls(in[(j-1)*n+i-1],DoubleComplex(g,0));
-
- for (j = 1; j <= l; j++)
- in[(i-1)*n+j-1] = zmuls(in[(i-1)*n+j-1],DoubleComplex(f,0));
- }
- }
- }
- }while(noconv); /* while (noconv == true) */
-
- *low = k;
- *high = l;
-}
-
diff --git a/src/matrixOperations/logm/comqr3.c b/src/matrixOperations/logm/comqr3.c
deleted file mode 100644
index 1e8fb063..00000000
--- a/src/matrixOperations/logm/comqr3.c
+++ /dev/null
@@ -1,246 +0,0 @@
-/*
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008 - INRIA - Arnaud TORSET
- *
- * This file must be used under the terms of the CeCILL.
- * This source file is licensed as described in the file COPYING, which
- * you should have received as part of this distribution. The terms
- * are also available at
- * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
- *
- */
-
- /*
- This is a transcription of the fortran subroutine comqr3.f
- */
-
-
-
-#include "logm_internal.h"
-#include "sqrt.h"
-#include "abs.h"
-#include "addition.h"
-#include "subtraction.h"
-#include "multiplication.h"
-#include "division.h"
-#include "pythag.h"
-#include "conj.h"
-#include "min.h"
-
-
-void comqr3(doubleComplex* in, int n, doubleComplex* w, int low, int high, int job, doubleComplex* z){
- /* Variables */
- int i=0, j=0, l=0, ll=0, ip1=0;
- int ierr=0, iend=0, jx=0, jy=0, en=0, itn=0, its=0, enm1=0, lp1=0;
- int boucle_220=0, boucle_240=0;
- double norm;
- doubleComplex s,t,x,y,zz;
-
-
- ierr=0;
- jx=job/10;
- jy=job%10;
- /*......... create real subdiagonal elements ..........*/
- iend = high-low-1;
-
- if (iend>=0){
- l=low+1;
-
- for(i=l;i<=high;i++){
- ll=min(i+1,high);
- if (zimags(in[(i-2)*n+i-1])!=0){
- norm = dpythags(zreals(in[(i-2)*n+i-1]),zimags(in[(i-2)*n+i-1]));
- y=zrdivs(in[(i-2)*n+i-1],DoubleComplex(norm,0));
- in[(i-2)*n+i-1]=DoubleComplex(norm,0);
-
- for(j=i;j<=n;j++){
- s=DoubleComplex(0, zimags(zmuls(y,zconjs(in[(j-1)*n+i-1]))));
- in[(j-1)*n+i-1] = DoubleComplex(zreals(zmuls(zconjs(y),in[(j-1)*n+i-1])),zimags(s));
- }
-
- for(j=1;j<=ll;j++){
- s=DoubleComplex(0, zimags(zmuls(zconjs(y),in[(i-1)*n+j-1])));
- in[(i-1)*n+j-1] = DoubleComplex(zreals(zmuls(zconjs(y),in[(i-1)*n+j-1])),zimags(s));
- }
-
- if (jy!=0){
- for (j=low;j<=high;j++){
- s=DoubleComplex(0, zimags(zmuls(y,z[(i-1)*n+j-1])));
- z[(i-1)*n+j-1] = DoubleComplex(zreals(zmuls(y,z[(i-1)*n+j-1])),zimags(s));
- }
- }
- }
- }
- }
-
- /*.......... store roots isolated by cbal ..........*/
- for(i=1;i<=n;i++){
- if( (i<low) || (i>high) ){
- w[i-1]=in[(i-1)*n+i-1];
- }
- }
-
- en = high;
- t=DoubleComplex(0,0);
- itn = 30*n;
- boucle_220=1;
-
- /*.......... search for next eigenvalue ..........*/
- while (boucle_220 && (en>=low)){
- its=0;
- enm1= en-1;
- /* .......... look for single small sub-diagonal element
- for l=en step -1 until low do -- ..........*/
- boucle_240=1;
- while(boucle_240){
-
- for(ll=low; ll<=en;ll++){
- l = en + low - ll;
- if (l!=low){
- x = DoubleComplex(dabss(zreals(in[(l-2)*n+l-2])) + dabss(zimags(in[(l-2)*n+l-2])) + dabss(zreals(in[(l-1)*n+l-1])) +dabss(zimags(in[(l-1)*n+l-1])),0);
- y = DoubleComplex(zreals(x) + dabss(zreals(in[(l-2)*n+l-1])),0);
- if (zreals(x) == zreals(y)) break;
- }
- }
- /*.......... form shift ..........*/
- if (l==en) {
- boucle_240=0;
- break;
- }
- if (itn==0) {
- boucle_240=0;
- break;
- }
- if ( (its == 10) || (its == 20) ){
- /*.......... form exceptional shift ..........*/
- s = DoubleComplex(dabss( zreals(in[(enm1-1)*n+en-1]) ) + dabss( zreals(in[(en-3)*n+enm1-1]) ),0);
- }
- else {
- s = in[(en-1)*n+en-1];
-
- x = zmuls(in[(en-1)*n+enm1-1],DoubleComplex(zreals(in[(enm1-1)*n+en-1]),0));
-
- if ( (zreals(x)!=0) || (zimags(x)!=0) ){
- y=zrdivs(zdiffs(in[(enm1-1)*n+enm1-1],s),DoubleComplex(2,0));
-
- zz=zsqrts(zadds(zmuls(y,y),x));
-
- if ( (zreals(y)*zreals(zz) + zimags(y)*zimags(zz)) < 0){
- zz = DoubleComplex(-zreals(zz),-zimags(zz));
- }
- zz=zrdivs(x,zadds(y,zz));
- s = zdiffs(s,zz);
- }
- }
-
- for(i=low;i<=en;i++){
- in[(i-1)*n+i-1]=zdiffs(in[(i-1)*n+i-1],s);
- }
-
- t = zadds(t,s);
- its++;
- itn--;
- /*.......... reduce to triangle (rows) ..........*/
- lp1=l+1;
-
- for(i=lp1;i<=en;i++){
- s=DoubleComplex(zreals(in[(i-2)*n+i-1]),zimags(s));
- norm = dpythags(dpythags(zreals(in[(i-2)*n+i-2]),zimags(in[(i-2)*n+i-2])),zreals(s));
- x=zrdivs(in[(i-2)*n+i-2],DoubleComplex(norm,0));
- w[i-2]=x;
- in[(i-2)*n+i-2]=DoubleComplex(norm, 0);
- in[(i-2)*n+i-1]=DoubleComplex(0, zreals(s)/norm);
-
- for (j=i;j<=n;j++){
- y = in[(j-1)*n+i-2];
- zz = in[(j-1)*n+i-1];
- in[(j-1)*n+i-2] = zadds(zmuls(zconjs(x),y),zmuls(DoubleComplex(zimags(in[(i-2)*n+i-1]),0),zz));
- in[(j-1)*n+i-1] = zdiffs(zmuls(x,zz),zmuls(DoubleComplex(zimags(in[(i-2)*n+i-1]),0),y));
- }
- }
- /* FIXME s make sign problem with the exemple */
- s=DoubleComplex(zreals(s),zimags(in[(en-1)*n+en-1]));
- if( zimags(s)!=0){
-
- norm = dpythags(zreals(in[(en-1)*n+en-1]),zimags(s));
- s = DoubleComplex( zreals(in[(en-1)*n+en-1])/norm , zimags(s)/norm);
-
- in[(en-1)*n+en-1]=DoubleComplex( norm, 0);
- if (en!=n){
-
- ip1=en+1;
- for(j=ip1;j<=n;j++){
- y=in[(j-1)*n+en-1];
- in[(j-1)*n+en-1]=zmuls(zconjs(s),y);
- }
- }
- }
-
- /*.......... inverse operation (columns) ..........*/
- for (j=lp1;j<=en;j++){
- x=w[j-2];
- for(i=1;i<=j;i++){
- y=DoubleComplex(zreals(in[(j-2)*n+i-1]),0);
- zz=in[(j-1)*n+i-1];
- if(i!=j){
-
- y=DoubleComplex(zreals(y), zimags(in[(j-2)*n+i-1]));
-
- in[(j-2)*n+i-1] = DoubleComplex(zreals(in[(j-2)*n+i-1]),
- zreals(x)*zimags(y) + zimags(x)*zreals(y) + zimags(in[(j-2)*n+j-1])*zimags(zz));
- }
-
- in[(j-2)*n+i-1] = DoubleComplex(
- zreals(x)*zreals(y) - zimags(x)*zimags(y) + zimags(in[(j-2)*n+j-1])*zreals(zz),
- zimags(in[(j-2)*n+i-1]));
-
- in[(j-1)*n+i-1] = zdiffs(zmuls(zconjs(x),zz),zmuls(DoubleComplex(zimags(in[(j-2)*n+j-1]),0),y));
- }
-
-
-
- if(jy!=0){
- for(i=low;i<=high;i++){
- y=z[(j-2)*n+i-1];
- zz=z[(j-1)*n+i-1];
- z[(j-2)*n+i-1] = zadds(zmuls(x,y),zmuls(DoubleComplex(zimags(in[(j-2)*n+j-1]),0),zz));
-
- z[(j-1)*n+i-1] = zdiffs(zmuls(zconjs(x),zz),zmuls(DoubleComplex(zimags(in[(j-2)*n+j-1]),0),y));
- }
- }
- }
-
- if (zimags(s)!=0){
- for (i=1;i<=en;i++){
- y=in[(en-1)*n+i-1];
- in[(en-1)*n+i-1]=zmuls(s,y);
- }
-
- if(jy!=0){
- for(i=low;i<=high;i++){
- y=z[(en-1)*n+i-1];
- z[(en-1)*n+i-1]=zmuls(s,y);
- }
- }
- }
-
- }
-
-
- if( (l!=en) && (itn==0) ){
- boucle_220=0;
- break;
- }
- in[(en-1)*n+en-1]=zadds(in[(en-1)*n+en-1],t);
- w[en-1] = in[(en-1)*n+en-1];
- en = enm1;
-
- }
-
- ierr=en;
-}
-
-
-
-
-
diff --git a/src/matrixOperations/logm/corth.c b/src/matrixOperations/logm/corth.c
deleted file mode 100644
index a98edefc..00000000
--- a/src/matrixOperations/logm/corth.c
+++ /dev/null
@@ -1,112 +0,0 @@
-/*
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008 - INRIA - Arnaud TORSET
- *
- * This file must be used under the terms of the CeCILL.
- * This source file is licensed as described in the file COPYING, which
- * you should have received as part of this distribution. The terms
- * are also available at
- * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
- *
- */
-
-/*
- This is a transcription of the fortran subroutine corth.f
- */
-
-
-#include "logm_internal.h"
-#include "abs.h"
-#include "addition.h"
-#include "subtraction.h"
-#include "multiplication.h"
-#include "division.h"
-#include "conj.h"
-
-void corth(doubleComplex* in , int n, int low,int high, doubleComplex* ort){
- int la=0, kp1=0, m=0, mp=0, i=0, j=0, ii=0, jj=0;
- double h=0, scale=0, g=0, f=0;
- doubleComplex fCpx;
-
- la = high - 1 ;
- kp1 = low + 1 ;
- if (la >= kp1) {
-
- for (m = kp1; m <= la; m++){
- h = 0;
- ort[m-1]=DoubleComplex(0,0);
- scale = 0;
- /* :::::::::: scale column (algol tol then not needed) :::::::::: */
- for (i = m; i <= high; i++)
- scale = scale + dabss(zreals(in[(m-2)*n+i-1])) + dabss(zimags(in[(m-2)*n+i-1]));
-
-
-
- if (scale != 0) {
- mp = m + high;
- /* :::::::::: for i=igh step -1 until m do -- :::::::::: */
- for (ii = m; ii <= high; ii++){
- i = mp - ii;
- ort[i-1] = zrdivs(in[(m-2)*n+i-1],DoubleComplex(scale,0));
- h = h + zreals(ort[i-1]) * zreals(ort[i-1]) + zimags(ort[i-1]) * zimags(ort[i-1]);
- }
-
-
- g = dsqrts(h);
-
- f = dsqrts(zreals(ort[m-1]) * zreals(ort[m-1]) + zimags(ort[m-1]) * zimags(ort[m-1]));
-
- if (f == 0){
- ort[m-1]=DoubleComplex(g,zimags(ort[m-1]));
- in[(m-2)*n+m-1]=DoubleComplex(scale,zreals(in[(m-2)*n+m-1]) );
- }
- else {
- h = h + f * g;
- g = g / f;
- ort[m-1]=zmuls(DoubleComplex(1+g,0),ort[m-1]);
- }
-
- /* :::::::::: form (i-(u*ut)/h) * a :::::::::: */
- for (j=m;j<=n;j++){
-
- fCpx = DoubleComplex(0,0);
-
- /* :::::::::: for i=igh step -1 until m do -- :::::::::: */
- for (ii=m;ii<=high;ii++){
- i = mp - ii;
- fCpx = zadds(fCpx, zmuls(zconjs(ort[i-1]),in[(j-1)*n+i-1]));
- }
-
- fCpx=zrdivs(fCpx,DoubleComplex(h,0));
-
- for (i=m;i<=high;i++){
- in[(j-1)*n+i-1] = zdiffs(in[(j-1)*n+i-1],zmuls(zconjs(fCpx),ort[i-1]));
- }
- }
-
- /* :::::::::: form (i-(u*ut)/h)*a*(i-(u*ut)/h) :::::::::: */
- for(i=1;i<=high;i++){
- fCpx = DoubleComplex(0,0);
- /* :::::::::: for j=igh step -1 until m do -- :::::::::: */
- for (jj=m;jj<=high;jj++){
- j = mp - jj;
- fCpx = zadds(fCpx, zmuls(ort[j-1],in[(j-1)*n+i-1]));
- }
-
- fCpx=zrdivs(fCpx,DoubleComplex(h,0));
-
- for (j=m;j<=high;j++){
- in[(j-1)*n+i-1] = zadds(in[(j-1)*n+i-1], zmuls(fCpx,DoubleComplex(-zreals(ort[j-1]),zimags(ort[j-1]))));
- }
- }
-
- ort[m-1] = zmuls(DoubleComplex(scale,0),ort[m-1]);
- in[(m-2)*n+m-1] = zmuls(DoubleComplex(-g,0),in[(m-2)*n+m-1]);
-
-
- }
- }
-
-
- }
-}
diff --git a/src/matrixOperations/logm/cortr.c b/src/matrixOperations/logm/cortr.c
deleted file mode 100644
index 609afe86..00000000
--- a/src/matrixOperations/logm/cortr.c
+++ /dev/null
@@ -1,72 +0,0 @@
-/*
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008 - INRIA - Arnaud TORSET
- *
- * This file must be used under the terms of the CeCILL.
- * This source file is licensed as described in the file COPYING, which
- * you should have received as part of this distribution. The terms
- * are also available at
- * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
- *
- */
-
- /*
- This is a transcription of the fortran subroutine cortr.f
- */
-
-#include "logm_internal.h"
-#include "conj.h"
-#include "addition.h"
-#include "multiplication.h"
-#include "division.h"
-
-
-
-void cortr(doubleComplex* in, int n, int low, int high, doubleComplex* ort, doubleComplex* z){
- int i=0, j=0, k=0, ii=0, iend=0, ip1=0;
- double norm;
- doubleComplex s;
-
-
- /*.......... initialize eigenvector matrix ..........*/
- for (i=1;i<=n;i++){
- for (j=1;j<=n;j++){
- if (i==j) z[(j-1)*n+i-1]=DoubleComplex(1,0);
- else z[(j-1)*n+i-1]=DoubleComplex(0,0);
- }
- }
-
- /*.......... form the matrix of accumulated transformations
- from the information left by corth ..........*/
- iend = high -low -1;
- if (iend >= 0){
- for (ii=1;ii<=iend;ii++){
- i=high-ii;
- if ( ( !((zreals(ort[i-1])==0) && (zimags(ort[i-1])==0) ))
- && ( !((zreals(in[(i-2)*n+i-1])==0) && (zimags(in[(i-2)*n+i-1])==0) ) ) ){
- /* .......... norm below is negative of h formed in corth .........*/
- norm = zreals(zmuls(zconjs(in[(i-2)*n+i-1]),ort[i-1]));
- if (norm != 0 ){
- ip1=i+1;
- for(k=ip1;k<=high;k++){
- ort[k-1] = DoubleComplex( zreals(in[(i-2)*n+k-1]), zimags(in[(i-2)*n+k-1]) );
- }
-
- for (j=i;j<=high;j++){
- s=DoubleComplex(0,0);
- for(k=i;k<=high;k++){
- s=zadds(s, zmuls(zconjs(ort[k-1]),z[(j-1)*n+k-1]));
- }
-
- s = zrdivs(s,DoubleComplex(norm,0));
-
- for(k=i;k<=high;k++){
- z[(j-1)*n+k-1]=zadds(z[(j-1)*n+k-1], zmuls(s,ort[k-1]));
- }
- }
- }
- }
- }
- }
-}
-
diff --git a/src/matrixOperations/logm/logm_internal.h b/src/matrixOperations/logm/logm_internal.h
deleted file mode 100644
index e0263cfc..00000000
--- a/src/matrixOperations/logm/logm_internal.h
+++ /dev/null
@@ -1,37 +0,0 @@
-/*
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008 - INRIA - Arnaud TORSET
- *
- * This file must be used under the terms of the CeCILL.
- * This source file is licensed as described in the file COPYING, which
- * you should have received as part of this distribution. The terms
- * are also available at
- * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
- *
- */
-
-
-
-#ifndef __LOGM_INTERNAL_H__
-#define __LOGM_INTERNAL_H__
-
-#include "doubleComplex.h"
-
-void wbdiag(doubleComplex* in, int size, doubleComplex* out);
-
-void cbal(doubleComplex* in, int n, int* low, int* high, double* scale);
-
-void corth(doubleComplex* in , int n, int low,int high, doubleComplex* ort);
-
-void cortr(doubleComplex* in, int n, int low, int high, doubleComplex* ort, doubleComplex* z);
-
-void comqr3(doubleComplex* in, int n, doubleComplex* w, int low, int high, int job, doubleComplex* z);
-
-void wexchn();
-void dad(doubleComplex* in, int size_1, int l11_1, int l22m1, int l11_2, int size_2, int one, int n);
-void balbak();
-void dset();
-void dad();
-void wshrsl();
-doubleComplex ddot();
-#endif /*__LOGM_INTERNAL_H__*/
diff --git a/src/matrixOperations/logm/wbdiag.c b/src/matrixOperations/logm/wbdiag.c
deleted file mode 100644
index da6fb1d6..00000000
--- a/src/matrixOperations/logm/wbdiag.c
+++ /dev/null
@@ -1,306 +0,0 @@
-/*
- * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
- * Copyright (C) 2008 - INRIA - Arnaud TORSET
- *
- * This file must be used under the terms of the CeCILL.
- * This source file is licensed as described in the file COPYING, which
- * you should have received as part of this distribution. The terms
- * are also available at
- * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
- *
- */
-
-#include "logm_internal.h"
-#include "pow.h"
-#include "abs.h"
-#include "addition.h"
-#include "subtraction.h"
-#include <malloc.h>
-#include <stdio.h>
-
-void wbdiag( doubleComplex* in, int size, doubleComplex* out){
- double norm=0, rav, cav,c,d;
- int i=0, j=0, ok=0, l1=0, l11pi=0,k,l, *bs, fail;
- int l11=0, l22=0, l22m1=0, da22=0, da11=0, one=1, mone, km1;
- int job=0; /* It seems job is always zero for a call with a complex matrix. It must be confirmed.
- If it's true, we could remove some parts of the code that use j<>0 */
- int low=0, high=0;
- double* scale;
- doubleComplex *eigenvalue, *x, *y, temp;
- int boucle_50, boucle_60, boucle_110, boucle_420, boucle_460, boucle_520, boucle_560;
-
-
- scale = malloc((uint)size*sizeof(double));
- eigenvalue = malloc((uint)size*sizeof(doubleComplex));
- x = malloc((uint)size*sizeof(doubleComplex));
- y = malloc((uint)size*sizeof(doubleComplex));
- bs = malloc((uint)size*sizeof(int));
-
-
- /* Compute l1 norm of in */
- for (i=0;i<size*size;i++){
- norm += (zreals(in[i])+zimags(in[i]));
- }
-
-
- /* convert a to upper hessenberg form */
- cbal(in, size, &low, &high, scale);
- corth(in, size, 1, size, eigenvalue);
- cortr(in, size, 1, size, eigenvalue, x);
-
-
- /* convert a to upper triangular form by qr method */
- comqr3(in, size, eigenvalue, 1, size, 11, x);
-
- /* reduce a to block diagonal form
-
- segment a into 4 matrices: a11, a 1 x 1 block
- whose (1,1)-element is at a(l11,l11)) a22, a 1 x 1
- block whose (1,1)-element is at a(l22,l22)) a12,
- a 1 x 1 block whose (1,1)-element is at a(l11,l22))
- and a21, a 1 x 1 block = 0 whose (1,1)-
- element is at a(l22,l11).
-
-
-
- this loop uses l11 as loop index and splits off a block
- starting at a(l11,l11). */
-
- l11=1;
-
- while(l11<=size){/* debut boucle_40 */
- l22=l11;
- /* this loop uses da11 as loop variable and attempts to split
- c off a block of size da11 starting at a(l11,l11) */
- while (boucle_50){
- if (l22!=l11) boucle_60=1;
- else {
- da11 = 1;
- l22 = l11 + 1;
- l22m1 = l22 - 1;
- boucle_60=0;
- }
- if (boucle_60==1){
- /* compute the average of the eigenvalues in a11 */
- rav = 0;
- cav = 0;
- for (i=l11;i<l22m1;i++){
- rav = rav + zreals(eigenvalue[i]);
- cav = cav + dabss(zimags(eigenvalue[i]));
- }
-
- rav = rav/da11;
- cav = cav/da11;
- /* loop on eigenvalues of a22 to find the one closest to the av */
- d = dpows(rav-zreals(eigenvalue[l22-1]),2) + dpows(cav-zimags(eigenvalue[l22-1]),2);
- k = l22;
- l = l22 + 1;
-
- if(l<size){
- c = dpows(rav-zreals(eigenvalue[l-1]),2) + dpows(cav-zimags(eigenvalue[l-1]),2);
- if (c < d){
- k = l;
- d = c;
- }
- l = l + 1;
- }
-
- /* loop to move the eigenvalue just located
- into first position of block a22.
- the block we're moving to add to a11 is a 1 x 1 */
- boucle_110=1;
- while (boucle_110){
- if (k==l22) boucle_110=0;
- else{
- km1 = k - 1;
- wexchn();
- temp = eigenvalue[k-1];
- eigenvalue[k-1] = eigenvalue[km1-1];
- eigenvalue[km1-1] = temp;
- k = km1;
-
- if (k<=l22) boucle_110=0;
- }
- }
- da11 = da11 + 1;
- l22 = l11 + da11;
- l22m1 = l22 - 1;
- }
- if (l22>=size) boucle_50=0;
- else {
- /* attempt to split off a block of size da11. */
- da22 = size - l22 + 1;
-
- /* save a12 in its transpose form in block a21. */
- for (j=l11;j<l22m1;j++){
- for (i=l22;i<size;i++){
- in[(i-1)*size+j]=in[(j-1)*size+i];
- }
- }
- /* convert a11 to lower quasi-triangular and multiply it by -1 and
- c a12 appropriately (for solving -a11*p+p*a22=a12).
-
- c write(6,'(''da11='',i2,''da22='',i2)') da11,da22
- c write(6,'(''a'')')
- c call wmdsp(ar,ai,n,n,n,10,1,80,6,cw,iw)
- */
-
- dad(in, size, l11, l22m1, l11, size, one, 0);
- dad(in, size, l11, l22m1, l11, l22m1, mone, 1);
-
- /* solve -a11*p + p*a22 = a12. */
- wshrsl();
- if (ok) boucle_50=0;
- else {
- /* change a11 back to upper quasi-triangular. */
- dad(in, size, l11, l22m1, l11, l22m1, one, 1);
- dad(in, size, l11, l22m1, l11, l22m1, mone, 0);
-
- /* move saved a12 back into its correct position. */
- for (j=l11;j<l22m1;j++){
- for(i=l22;i<size;i++){
- in[(j-1)*size+i-1] = in[(i-1)*size+j-1];
- in[(i-1)*size+j-1] = DoubleComplex(0,0);
- }
-
- }
- }
-
- }
- }/*boucle_50*/
- /* change solution to p to proper form */
- if (l22<=size){
- dad(in,size,l11,l22m1,l11,size,one,0);
- dad(in,size,l11,l22m1,l11,l22m1,mone,1);
- }
- bs[l11-1]=da11;
- j=da11-1;
- if (j!=0){
- for(i=1;i<=j;i++){
- l11pi=l1+i;
- bs[l11pi-1]=-(da11-i);
- }
- l11=l22;
- }/* boucle_40 */
-
-
- fail=1;
- /*
- set transformations matrices as required
- if (job == 3) return;
-
- compute inverse tranformation
- if (job ==1){
- for (i=0;i<size;i++){
- for(j=0;j<size;j++){
- y[i*size+j-1]=x[j*size+i-1];
- }
- }
-
- l22=1;
- boucle_420=1;
- while (boucle_420){
- l11=l22;
- l22=l11+bs[l11-1];
- if (l22>size) {
- boucle_420=0;
- break;
- }
- l22m1=l22-1;
-
- for (i=l11;i<=l22m1;i++){
- for(j=1;j<=size;j++){
- y[i*size+j-1]=zadds(zdiffs(y[i*size+j-1],ddot()),ddot());
- }
- }
- }
-
- if (high<>low){
- for (j=low;j<=high;j++){
- temp=1/scale[j-1];
- for(i=1;i<=size;i++){
- y[i*size+j-1]=y[i*size+j-1]*temp;
- }
- }
- }
-
- for (ii=1;ii<+size;i++){
- i=ii;
- if ( (i<low) || (i>high) ){
- if (i<low) i=low-ii;
- k=scale[i-1];
- if (k<>i){
- for (j=1;j<=size;j++){
- temp=y[j*size+i-1];
- y[j*size+i-1] = y[j*size+k-1];
- y[j*size+k-1] = temp;
- }
- }
- }
- }
- }
- */
- if (job!=2){
- /* Compute right transformation */
- l22=1;
- boucle_460=1;
- while(boucle_460){
- l11=l22;
- l22=l11+bs[l11-1];
- if (l22>size) {
- boucle_420=0;
- break;
- }
-
- for (j=l22;j<=size;j++){
- for(i=1;i<=size;i++){
- x[i*size+j-1]=zadds(zdiffs(x[i*size+j-1],ddot()),ddot());
- }
- }
- }
- balbak();
- }
-
- /* extract non orthogonal transformation from in */
- for (j=1;j<size;j++){
- dset();
- }
-
- dset();
- l22=1;
- boucle_520=1;
- while(boucle_520){
- l11=l22;
- if(l11>size) {
- boucle_520=0;
- break;
- }
- l22=l11+bs[l11-1];
- for (j=l22;j<=size;j++){
- for (i=1;i<=size;i++){
- y[i*size+j-1]=zdiffs(zadds(y[i*size+j-1],ddot()),ddot());
- }
- }
- }
- }
-
- /* set zeros in the matrix in */
- l11=1;
- boucle_560=1;
- while (boucle_560){
- l22=l11+bs[l11-1];
- if (l22>size){
- boucle_560=0;
- break;
- }
- l22m1=l22-1;
-
- for(j=l11;j<=l22m1;j++){
- dset();
- dset();
- }
- l11=l22;
- }
- out=in;
-
-}
generated by cgit (git 2.34.1) at 2025-03-26 13:58:15 +0000