--- a
+++ b/EigenLib3.mod
@@ -0,0 +1,869 @@
+IMPLEMENTATION MODULE EigenLib3;
+
+ (*------------------------------------------------------------------------*)
+ (* Routinen zu Berechnung der Eigensysteme komplexer genereller Matrizen *)
+ (* Routines to calculate the eigensystems of complex general matrices *)
+ (*------------------------------------------------------------------------*)
+ (* Letzt Bearbeitung: *)
+ (* *)
+ (* 07.08.93, MRi: Durchsicht *)
+ (* 06.02.15, MRi: Ausgetestete Routinen aus EisPack herausgenommen und *)
+ (* in EigenLib3.mod integriert (CSortEwEv,ComEig,ComHess *)
+ (* CHessLR,CHessInvIter) *)
+ (* 15.07.15, MRi: CSortEwEV,RHessQR,CHessLR auf ISO COMPLEX umgestellt *)
+ (* 18.07.15, MRi: Umstellung vom ComEig auf ISO-COMPLEX *)
+ (* 01.04.16, MRi: Aufraeumen wg. Compilerwarnungen, kleine Code- *)
+ (* veraenderung in CHessInvIter bei der Bestimmung von *)
+ (* Norm & VNorm *)
+ (* 20.09.17, MRi: Umstellung vom ComEig auf offene Felder *)
+ (* 21.09.17, MRi: Umstellen von ComHess und CHessLR auf ISO-COMPLX *)
+ (* 23.09.17, MRi: Umstellen von CHessInvIter auf ISO-COMPLX (neu) *)
+ (* 28.09.17, MRi: CBalance und CBalBak aufgenommen. In CBalance den *)
+ (* Parameter iFehl ergaenzt *)
+ (* 19.10.17, MRi: In ComEig die Transformationsschleifen von Komponenten- *)
+ (* arithmetik zur expliziter komplexer Arithmetik umge- *)
+ (* stellt, ebenfalls die Berechnung von H *)
+ (* 22.10.17, MRi: In ComEig MaxIter fix auf 70 begrenzt *)
+ (*------------------------------------------------------------------------*)
+ (* Offene Punkte *)
+ (* *)
+ (* - Umstellen auf offene Felder *)
+ (* - ComEig auf komplexe Arithmethikg durchsehen und vereinfachen *)
+ (* - Paramer "select" f��r CHessInvIter und ComEigEvEw einfuehren ? *)
+ (* - CHessInvIter laeuft fuer einen der beiden Eigenwert (1.0,1.0) aus *)
+ (* Num. Math. 16, 181-204 (1970) - BSP 2 auf einen Fehler - warum ??? *)
+ (* - CBalance liefert fuer "Test unbekannter Herkunft" unsinnige Werte *)
+ (*------------------------------------------------------------------------*)
+ (* Testprogramm ist TstCRM.mod *)
+ (*------------------------------------------------------------------------*)
+ (* Implementation : Michael Riedl *)
+ (* Licence : GNU Lesser General Public License (LGPL) *)
+ (*------------------------------------------------------------------------*)
+
+ (* $Id: EigenLib3.mod,v 1.4 2017/10/22 17:52:02 mriedl Exp mriedl $ *)
+
+FROM SYSTEM IMPORT TSIZE,ADR;
+FROM Storage IMPORT ALLOCATE,DEALLOCATE;
+FROM Errors IMPORT Fehler,Fehlerflag,ErrOut,FatalError;
+ IMPORT Errors;
+FROM Deklera IMPORT MaxDim,VEKTOR,CARDVEKTOR,CVEKTOR,CMATRIX;
+FROM LongMath IMPORT sqrt;
+ IMPORT LongComplexMath;
+FROM LongComplexMath IMPORT zero,one,scalarMult;
+FROM LMathLib IMPORT MachEps,MachEpsR4;
+FROM EigenLibAux IMPORT CNormMat;
+ IMPORT TIO;
+
+CONST debug = FALSE;
+
+(*================= Hilfsroutinen ==========================================*)
+
+PROCEDURE CEucNorm(VAR X : ARRAY OF LONGCOMPLEX; (* unchanged *)
+ n : INTEGER;
+ VAR Norm : LONGREAL);
+
+ (*--------------------------------------------------------*)
+ (* Berechnet die euclidische Norm des komplexen Vektors X *)
+ (* Modula-2 Transkription der LAPACK- Routine dznrm2 *)
+ (*--------------------------------------------------------*)
+
+ VAR scale,ssq,temp,t1 : LONGREAL;
+ i : INTEGER;
+BEGIN
+ IF (n < 1) THEN
+ Norm := 0.0;
+ ELSE
+ scale := 0.0;
+ ssq := 1.0;
+ FOR i:=0 TO n-1 DO
+ IF (RE(X[i]) # 0.0)THEN
+ temp := ABS(RE(X[i]));
+ IF (scale < temp) THEN
+ t1 := scale / temp;
+ ssq := 1.0 + ssq*t1*t1;
+ scale := temp;
+ ELSE
+ t1 := temp / scale;
+ ssq :=ssq + t1*t1;
+ END;
+ END;
+ IF (IM(X[i]) # 0.0) THEN
+ temp := ABS(IM(X[i]));
+ IF (scale < temp) THEN
+ t1 := scale / temp;
+ ssq := 1.0 + ssq*t1*t1;
+ scale := temp;
+ ELSE
+ t1 := temp/scale;
+ ssq := ssq+t1*t1;
+ END;
+ END;
+ END; (* FOR i *)
+ Norm := scale*sqrt(ssq);
+ END; (* IF *)
+END CEucNorm;
+
+PROCEDURE CNormMat2( dim : CARDINAL;
+ VAR CMat : CMATRIX);
+
+ VAR Norm : LONGREAL;
+ i,j : CARDINAL;
+BEGIN
+ FOR i:=1 TO dim DO
+ CEucNorm(CMat[i],dim,Norm);
+ IF (Norm # 0.0) THEN
+ Norm := 1.0 / Norm;
+ END;
+ FOR j:=1 TO dim DO
+ CMat[i,j] := scalarMult(Norm,CMat[i,j]);
+ END;
+ END;
+END CNormMat2;
+
+(*================= Ende Hilfsroutinen =====================================*)
+
+PROCEDURE CSortEwEv(VAR EW : CVEKTOR; (* Eigenwerte *)
+ VAR EV : CMATRIX; (* Eigenvektoren *)
+ dim : CARDINAL;
+ ord : INTEGER);
+
+ VAR i,l,iminmax : CARDINAL;
+ MinMax : LONGREAL;
+ Zw : LONGCOMPLEX;
+ AbsEW : VEKTOR;
+BEGIN
+ FOR i:=1 TO dim DO AbsEW[i]:=LongComplexMath.abs(EW[i]); END;
+ FOR l:=1 TO dim-1 DO
+ MinMax:=AbsEW[l]; iminmax:=l;
+ FOR i:=l+1 TO dim DO
+ IF (ord = 1) THEN (* Berechne Maximum von AbsEW[dim] *)
+ IF (AbsEW[i] > MinMax) THEN MinMax:=AbsEW[i]; iminmax:=i; END;
+ ELSE (* Berechne Minimum von AbsEW[dim] *)
+ IF (AbsEW[i] < MinMax) THEN MinMax:=AbsEW[i]; iminmax:=i; END;
+ END; (* IF Ordnung *)
+ END; (* FOR i *)
+ IF (iminmax <> l) THEN (* Vertausche EWi mit EWimax *)
+ Zw:=EW[l]; EW[l]:=EW[iminmax]; EW[iminmax]:=Zw;
+ AbsEW[iminmax]:=AbsEW[l];
+ FOR i:=1 TO dim DO (* Vertausche Eigenvektoren entsprechend *)
+ Zw:=EV[l,i]; EV[l,i]:=EV[iminmax,i]; EV[iminmax,i]:=Zw;
+ END; (* i *)
+ END; (* IF *)
+ END; (* l *)
+END CSortEwEv;
+
+PROCEDURE CBalance(VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
+ dim : CARDINAL;
+ VAR low,high : CARDINAL;
+ VAR Scale : ARRAY OF LONGREAL;
+ VAR iFehl : INTEGER);
+
+ VAR j,k,l,N : INTEGER;
+
+ PROCEDURE exc(m : INTEGER); (* lokal zu CBalance *)
+
+ VAR t : LONGCOMPLEX;
+ i : INTEGER;
+ BEGIN
+ Scale[m]:=LFLOAT(j);
+ IF (j # m) THEN
+ FOR i:=0 TO k DO t:=A[i,j]; A[i,j]:=A[i,m]; A[i,m]:=t; END;
+ FOR i:=l TO N-1 DO t:=A[j,i]; A[j,i]:=A[m,i]; A[m,i]:=t; END;
+ END;
+ END exc;
+
+ CONST radix = 2.0; (* Basis der Flie\3kommadarstellung *)
+ sqrdx = radix*radix;
+
+ VAR s,r,g,f,c : LONGREAL;
+ i : INTEGER;
+ Konvergenz,cont,flag,found : BOOLEAN;
+BEGIN
+ N:=VAL(INTEGER,dim);
+ iFehl:=0;
+
+ l := 0;
+ k := N-1;
+ REPEAT
+ (* Search for rows isolating an eigenvalue and push them down *)
+ j := k;
+ REPEAT
+ i := 0;
+ REPEAT
+ flag := (i # j) AND (A[j,i] # zero);
+ INC(i);
+ UNTIL flag OR (i > k);
+ found := NOT flag;
+ IF found THEN
+ exc(k);
+ DEC(k)
+ END;
+ DEC(j);
+ UNTIL found OR (j < 0);
+ UNTIL (NOT found) OR (k < 0);
+
+ IF (k < 0) THEN k := 0; END;
+ IF (k = 0) THEN
+ Errors.ErrOut("high = 0 (CBalance)");
+ low :=1;
+ high:=VAL(CARDINAL,k+1);
+ iFehl:=-1;
+ RETURN;
+ END;
+
+ REPEAT
+ (* Search for columns isolating an eigenvalue and push them left *)
+ j := l;
+ REPEAT
+ i := l;
+ REPEAT
+ flag := (i # j) AND (A[i,j] # zero);
+ INC(i);
+ UNTIL flag OR (i > k);
+ found := NOT flag;
+ IF found THEN
+ exc(l);
+ INC(l);
+ END;
+ INC(j);
+ UNTIL found OR (j > k);
+ UNTIL (NOT found);
+
+ low :=VAL(CARDINAL,l+1);
+ high:=VAL(CARDINAL,k+1);
+
+ (* Balanciere die Submatrix A[i,l..k] *)
+ FOR i:=l TO k DO Scale[i]:=1.0; END;
+ REPEAT
+ Konvergenz:=TRUE; (* Wird getestet *)
+ FOR i:=l TO k DO
+ c:=0.0; r:=0.0;
+ FOR j:=l TO k DO
+ IF (j # i) THEN
+ c:=c + ABS(RE(A[j,i])) + ABS(IM(A[j,i]));
+ r:=r + ABS(RE(A[i,j])) + ABS(IM(A[i,j])); (* ??? *)
+ END;
+ END;
+ IF (c # 0.0) AND (r # 0.0) THEN
+ g := r / radix;
+ f := 1.0;
+ s := c + r;
+ WHILE (c < g) DO f:=f * radix; c:=c * sqrdx; END;
+ g := r*radix;
+ WHILE (c >= g) DO f:=f / radix; c:=c / sqrdx; END;
+ IF (((c+r) / f) < 0.95*s) THEN
+ Konvergenz:=FALSE;
+ g := 1.0 / f;
+ Scale[i]:=Scale[i] * f;
+ FOR j:=l TO N-1 DO A[i,j]:=scalarMult(g,A[i,j]); END;
+ FOR j:=0 TO k DO A[j,i]:=scalarMult(f,A[j,i]); END;
+ END;
+ END; (* IF *)
+ END; (* FOR i *)
+ UNTIL Konvergenz;
+END CBalance;
+
+PROCEDURE CBalBak( dim : CARDINAL; (* Ordnung der Matrix Z *)
+ low,high : CARDINAL;
+ VAR Scale : ARRAY OF LONGREAL;
+ m : CARDINAL;
+ VAR Z : ARRAY OF ARRAY OF LONGCOMPLEX);
+
+ VAR i,j,k : CARDINAL;
+ S : LONGREAL;
+ Zw : LONGCOMPLEX;
+
+BEGIN
+ IF (m = 0) THEN RETURN END;
+
+ DEC(low); DEC(high);
+
+ IF (low # high) THEN
+ FOR i:=low TO high DO
+ S:=Scale[i]; (* Linksvektoren durch S:=1.0 / Scale[i] ! *)
+ FOR j:=0 TO m-1 DO
+ Z[j,i]:=scalarMult(S,Z[j,i]);
+ END;
+ END;
+ END;
+
+ (* for i=low-1 step -1 until 1,high+1 step 1 until dim do *)
+
+ IF (low > 0) THEN
+ FOR i:=low-1 TO 0 BY -1 DO
+ k := TRUNC(Scale[i]);
+ IF (k # i) THEN
+ FOR j:=0 TO m-1 DO Zw:=Z[j,i]; Z[j,i]:=Z[j,k]; Z[j,k]:=Zw; END;
+ END;
+ END;
+ END;
+ FOR i:=high+1 TO dim-1 DO
+ k := TRUNC(Scale[i]);
+ IF (k # i) THEN
+ FOR j:=0 TO m-1 DO Zw:=Z[j,i]; Z[j,i]:=Z[j,k]; Z[j,k]:=Zw; END;
+ END;
+ END;
+END CBalBak;
+
+PROCEDURE ComHess(VAR A : CMATRIX;
+ dim : CARDINAL;
+ k,l : CARDINAL;
+ VAR Index : CARDVEKTOR);
+
+ VAR i,j,m : CARDINAL;
+ X,Y : LONGCOMPLEX;
+BEGIN
+ Index[k]:=k; Index[l]:=l;
+ FOR m:=k+1 TO l-1 DO
+ i:=m; X:=CMPLX(0.0,0.0);
+ FOR j:=m TO l DO
+ IF ABS(RE(A[j,m-1]))+ABS(IM(A[j,m-1])) > ABS(RE(X))+ABS(IM(X)) THEN
+ X:=A[j,m-1]; i:=j;
+ END;
+ END;
+ Index[m]:=i;
+ IF (i # m) THEN (* interchange rows and columns of arrays a *)
+ FOR j:=m-1 TO dim DO Y:=A[i,j]; A[i,j]:=A[m,j]; A[m,j]:=Y; END;
+ FOR j:=1 TO l DO Y:=A[j,i]; A[j,i]:=A[j,m]; A[j,m]:=Y; END;
+ END;
+ IF (RE(X) # 0.0) OR (IM(X) # 0.0) THEN
+ FOR i:=m+1 TO l DO
+ Y:=A[i,m-1];
+ IF (RE(Y) # 0.0) OR (IM(Y) # 0.0) THEN
+ Y:=Y / X; A[i,m-1]:=Y;
+ FOR j:=m TO dim DO A[i,j]:=A[i,j] - Y*A[m,j]; END;
+ FOR j:=1 TO l DO A[j,m]:=A[j,m] + Y*A[j,i]; END;
+ END; (* IF Y *)
+ END; (* FOR i *)
+ END; (* IF X *)
+ END; (* FOR m *)
+END ComHess;
+
+PROCEDURE ComHessBak(VAR EV : CMATRIX;
+ VAR A : CMATRIX;
+ Index : CARDVEKTOR;
+ dim : CARDINAL;
+ EvAnz : CARDINAL;
+ k,l : CARDINAL);
+
+ VAR i,j,m : CARDINAL;
+ X : LONGCOMPLEX;
+BEGIN
+ IF (l > dim) THEN
+ Errors.FatalError("upp > dim (ComHessBak)");
+ END;
+ INC(k);
+ FOR m:=l-1 TO k BY -1 DO
+ FOR i:=m+1 TO l DO
+ X:=A[i,m-1];
+ IF (RE(X) # 0.0) OR (IM(X) # 0.0) THEN
+ FOR j:=1 TO EvAnz DO
+ EV[j,i]:=EV[j,i] + X*EV[j,m];
+ END;
+ END;
+ END; (* FOR i *)
+ i:=Index[m];
+ IF (i # m) THEN
+ FOR j:=1 TO EvAnz DO X:=EV[j,i]; EV[j,i]:=EV[j,m]; EV[j,m]:=X; END;
+ END;
+ END; (* FOR m *)
+END ComHessBak;
+
+PROCEDURE CHessLR(VAR H : CMATRIX; (* komplexe Hessenbergmatrix *)
+ VAR EW : CVEKTOR; (* komplexe Eigenwerte von H *)
+ dim : CARDINAL); (* Dimension von H *)
+
+ CONST MaxIter = 30;
+
+ VAR n,i,j,l,m,ll,mm,Iter : CARDINAL;
+ S,T,X,Y,Z : LONGCOMPLEX;
+ x,y,tmp : LONGREAL;
+BEGIN
+ X:=zero; Y:=zero; Z:=zero;
+ T:=zero; n:=dim;
+ WHILE (n > 0) DO
+ Iter:=0;
+ LOOP
+ LOOP (* Suche kleine Subdiagonalelemente *)
+ FOR l:=n TO 2 BY -1 DO
+ x:=ABS(RE(H[l ,l-1])) + ABS(IM(H[l ,l-1]));
+ y:=ABS(RE(H[l-1,l-1])) + ABS(IM(H[l-1,l-1])) +
+ ABS(RE(H[l ,l ])) + ABS(IM(H[l ,l ]));
+ IF (x <= MachEps*y) THEN ll:=l; EXIT; END;
+ END;
+ ll:=1; EXIT;
+ END;
+ l:=ll;
+ IF (l = n) THEN EXIT END; (* Wurzel gefunden ! *)
+
+ IF (Iter = MaxIter) THEN
+ Fehler:=TRUE; Fehlerflag:='MaxIter ueberschritten (CHessLR)';
+ RETURN;
+ END;
+ IF ((Iter MOD 10) = 0) AND (n > 2) THEN
+ S:=CMPLX(ABS(RE(H[n,n-1])) + ABS(RE(H[n-1,n-2])),
+ ABS(IM(H[n,n-1])) + ABS(IM(H[n-1,n-2])));
+ ELSE
+ S := H[n,n];
+ X := H[n-1,n]*H[n,n-1];
+ IF (RE(X) # 0.0) OR (IM(X) # 0.0) THEN
+ Y:=scalarMult(0.5,(H[n-1,n-1] - S));
+ Z:=CMPLX(
+ RE(Y)*RE(Y) - IM(Y)*IM(Y) + RE(X),
+ 2.0*RE(Y)*IM(Y) + IM(X)
+ );
+ Z:=LongComplexMath.sqrt(Z);
+ IF ((RE(Y)*RE(Z) + IM(Y)*IM(Z)) < 0.0) THEN
+ Z:= - Z;
+ END;
+ X:=X / (Y + Z);
+ S:=S - X;
+ END;
+ END; (* IF Iter *)
+
+ FOR i:=1 TO n DO H[i,i]:=H[i,i] - S; END;
+
+ (* Suche zwei aufeinander folgende kleine Subdiagonalelemente *)
+
+ X := CMPLX(ABS(RE(H[n-1,n-1])) + ABS(IM(H[n-1,n-1])),IM(X));
+ Y := CMPLX(ABS(RE(H[n ,n-1])) + ABS(IM(H[n ,n-1])),IM(Y));
+ Z := CMPLX(ABS(RE(H[n ,n ])) + ABS(IM(H[n ,n ])),IM(Z));
+
+ T:=T + S;
+ j:=l+1;
+ INC(Iter); (* !!! *)
+ LOOP
+ FOR m:=n-1 TO j BY -1 DO
+ Y := CMPLX(ABS(RE(H[m ,m-1])) + ABS(IM(H[m ,m-1])),RE(Y));
+ tmp := RE(X);
+ X := CMPLX(ABS(RE(H[m-1,m-1])) + ABS(IM(H[m-1,m-1])),RE(Z));
+ Z := CMPLX(tmp,IM(Z));
+ IF (RE(Y) < MachEps*(RE(Z) + RE(X) + IM(X))*(RE(Z)/IM(Y))) THEN
+ mm:=m; EXIT;
+ END;
+ END;
+ mm:=l; EXIT;
+ END;
+ m:=mm;
+ FOR i:=m+1 TO n DO
+ X:=H[i-1,i-1]; Y:=H[i,i-1];
+ IF ABS(RE(X)) + ABS(IM(X)) < ABS(RE(Y)) + ABS(IM(Y)) THEN
+ FOR j:=i-1 TO n DO Z:=H[i-1,j]; H[i-1,j]:=H[i,j]; H[i,j]:=Z; END;
+ Z := X / Y; EW[i]:=CMPLX(+1.0,IM(EW[i]));
+ ELSE
+ Z := Y / X; EW[i]:=CMPLX(-1.0,IM(EW[i]));
+ END;
+ H[i,i-1]:=Z;
+ FOR j:=i TO n DO H[i,j]:=H[i,j] - Z*H[i-1,j]; END;
+ END; (* FOR i *)
+ (* H = R*L *)
+ FOR j:=m+1 TO n DO
+ X:=H[j,j-1]; H[j,j-1]:=zero;
+ IF (RE(EW[j]) > 0.0) THEN
+ FOR i:=l TO j DO Z:=H[i,j-1]; H[i,j-1]:=H[i,j]; H[i,j]:=Z; END;
+ END;
+ FOR i:=l TO j DO H[i,j-1]:=H[i,j-1] + X*H[i,j]; END;
+ END; (* FOR j *)
+ END; (* LOOP *)
+ EW[n]:=H[n,n] + T;
+ DEC(n);
+ END; (* WHILE n > 0 *)
+END CHessLR;
+
+PROCEDURE CHessInvIter(VAR EV : CMATRIX; (* Eigenvektoren *)
+ VAR EW : CVEKTOR; (* Eigenwerte von A *)
+ VAR H : CMATRIX; (* Hessenbergmatrix *)
+ dim : CARDINAL; (* Dimension von H *)
+ VAR IErr : ARRAY OF CARDINAL;
+ VAR iFehl : INTEGER);
+
+ CONST MinEps = 100.0*MachEps;
+
+ VAR X,Y,eps3,Norm,NormV,GrowTo,SqrtUk : LONGREAL;
+ i,j,k,m,uk,Iter,ukk : CARDINAL;
+ Zw,Lambda : LONGCOMPLEX;
+ B : POINTER TO CMATRIX;
+ EVk : POINTER TO CVEKTOR;
+ gefunden : BOOLEAN;
+BEGIN
+ NEW(B);
+ IF (B = NIL) THEN
+ FatalError("kein Freispeicher vorhanden (EigenLib3.CHessInvIter) !");
+ END;
+ eps3:=0.0; GrowTo:=0.0; (* Wg. Compilerwarnungen *)
+
+ iFehl:=0; uk:=0; Fehler:=FALSE;
+ FOR k:=1 TO dim DO
+ EVk := ADR(EV[k]); (* Zeiger auf k-ten Eigenvektor *)
+ IF (uk < k) THEN
+ LOOP (* split *)
+ FOR uk:=k TO dim-1 DO
+ IF (RE(H[uk+1,uk]) = 0.0) AND (IM(H[uk+1,uk]) = 0.0) THEN
+ ukk:=uk;
+ EXIT;
+ END;
+ END;
+ ukk:=dim;
+ EXIT;
+ END;
+ uk:=ukk;
+ Norm:=0.0; m:=1;
+ FOR i:=1 TO uk DO (* norm of leading uk X uk matrix *)
+ X:=0.0;
+ FOR j:=m TO uk DO X:=X + LongComplexMath.abs(H[i,j]); END;
+ IF (X > Norm) THEN Norm:=X; END;
+ m:=i;
+ END; (* FOR *)
+
+ IF (Norm = 0.0) THEN Norm:=1.0; END;
+ eps3:=MachEps*Norm;
+ SqrtUk := sqrt(VAL(LONGREAL,uk));
+ GrowTo := 0.1 / SqrtUk
+ END; (* IF uk *)
+
+ Lambda := EW[k];
+ (* perturb eigenvalue if it is close to any previous eigenvalue *)
+ IF (k > 1) THEN
+ LOOP (* test *)
+ i:=k-1;
+ IF (i = 0) THEN EXIT; END;
+ REPEAT
+ gefunden:= (ABS(RE(EW[i]) - RE(Lambda)) < eps3) AND
+ (ABS(IM(EW[i]) - IM(Lambda)) < eps3);
+ DEC(i);
+ UNTIL gefunden OR (i = 0);
+ IF gefunden THEN
+ Lambda:=CMPLX(RE(Lambda) + eps3,IM(Lambda)); (* GOTO test: *)
+ ELSE
+ EXIT;
+ END;
+ END; (* LOOP *)
+ EW[k] := CMPLX(RE(Lambda),IM(EW[k]));;
+ END; (* IF *)
+
+ (* generate B = A -lambda x I, the matrix to be reduced, *)
+ (* and initial vector u + iv (EVk) *)
+ m:=1;
+ FOR i:=1 TO uk DO
+ FOR j:=m TO uk DO
+ B^[i,j] := H[i,j];
+ END;
+ B^[i,i]:=B^[i,i] - Lambda;
+ m:=i;
+ EVk^[i]:=CMPLX(eps3,0.0);
+ END; (* FOR *)
+
+ IF (uk # 1) THEN
+ (* reduce B to LU form *)
+ FOR i:=2 TO uk DO
+ m := i - 1;
+ IF (LongComplexMath.abs(B^[i,m]) >
+ LongComplexMath.abs(B^[m,m]))
+ THEN
+ FOR j:=m TO uk DO
+ Zw:=B^[i,j]; B^[i,j]:=B^[m,j]; B^[m,j]:=Zw;
+ END;
+ END; (* IF *)
+
+ IF (ABS(RE(B^[m,m])) <= MinEps) AND
+ (ABS(IM(B^[m,m])) <= MinEps) THEN
+ B^[m,m]:=CMPLX(eps3,IM(B^[m,m]));
+ END;
+ Zw := B^[i,m] / B^[m,m];
+ IF ((RE(Zw) # 0.0) OR (IM(Zw) # 0.0)) THEN
+ FOR j:=i TO uk DO
+ B^[i,j]:=B^[i,j] - Zw*B^[m,j];
+ END;
+ END; (* IF Zw = 0.0 *)
+ END; (* FOR i *)
+ END; (* IF uk # 1 *)
+ IF (RE(B^[uk,uk]) = 0.0) AND (IM(B^[uk,uk]) = 0.0) THEN
+ B^[uk,uk]:=CMPLX(eps3,0.0);
+ END;
+ Iter:=0;
+
+ LOOP (* back substitution *)
+ FOR i:=uk TO 1 BY -1 DO
+ Zw:=CMPLX(RE(EVk^[i]),0.0);
+ IF (i # uk) THEN
+ FOR j:=i+1 TO uk DO
+ Zw:=Zw - B^[i,j]*EVk^[j];
+ END; (* FOR *)
+ END; (* IF i # uk *)
+ EVk^[i] := Zw / B^[i,i];
+ END; (* FOR i *)
+
+ INC(Iter);
+
+ Norm:=0.0; NormV:=Norm; j:=1;
+ FOR i:=1 TO uk DO
+ X:= LongComplexMath.abs(EVk^[i]);
+ IF (NormV < X) THEN NormV:=X; j:=i; END;
+ Norm:=Norm + X;
+ END;
+
+ IF (Iter > 1) AND (Norm > GrowTo) THEN (* Akzeptabler Eigenvektor *)
+ Zw:=EVk^[j];
+ FOR i:=1 TO uk DO EVk^[i]:=EVk^[i] / Zw; END;
+ FOR i:=uk+1 TO dim DO EVk^[i]:=CMPLX(0.0,0.0); END;
+ j := uk + 1;
+ EXIT; (* !!! *)
+ END; (* IF *)
+
+ IF (Iter >= uk) THEN
+ IErr[iFehl]:=k; INC(iFehl); Fehler:=TRUE;
+ Fehlerflag:=
+ "nicht akzeptabler Eigenvektor (EigenLib3.CHessInvIter) !";
+ Errors.ErrOut(Fehlerflag);
+ Errors.WriteString(" uk = ");
+ Errors.WriteCard(uk);
+ Errors.WriteLn;
+ FOR i:=1 TO dim DO EVk^[i]:=CMPLX(0.0,0.0); END;
+ EXIT;
+ END;
+ (* Neuer Testvektor *)
+ X := sqrt(VAL(LONGREAL,dim));
+ Y := eps3 / (X + 1.0);
+ EVk^[1] := CMPLX(eps3,IM(EVk^[1]));
+ FOR i:=2 TO dim DO EVk^[i]:=CMPLX(Y,IM(EVk^[i])); END;
+ EVk^[j]:=EVk^[j] - CMPLX(eps3,0.0);
+ END; (* LOOP *)
+ END; (* FOR k *)
+ DISPOSE(B);
+END CHessInvIter;
+
+PROCEDURE ComEigEvEw(VAR A : CMATRIX;
+ dim : CARDINAL;
+ VAR EW : CVEKTOR;
+ VAR EV : CMATRIX;
+ Norm : CARDINAL;
+ ord : INTEGER;
+ VAR iFehl : INTEGER);
+
+ VAR i,j,u,low,upp : CARDINAL;
+ Index,IErr : CARDVEKTOR;
+ H : CMATRIX;
+
+BEGIN
+ low := 1; upp := dim;
+
+ FOR i:=1 TO dim DO Index[i]:=i; END;
+
+ ComHess(A,dim,low,upp,Index);
+
+ u:=1;
+ FOR i:=1 TO dim DO
+ FOR j:=1 TO u-1 DO H[i,j]:=CMPLX(0.0,0.0); END;
+ FOR j:=u TO dim DO H[i,j]:=A[i,j]; END;
+ IF (i > 2) THEN INC(u); END;
+ END;
+ CHessLR(H,EW,dim);
+
+ CHessInvIter(EV,EW,A,dim,IErr,iFehl);
+
+ ComHessBak(EV,A,Index,dim,dim,low,upp);
+
+ IF (Norm = 1) THEN
+ CNormMat(EV,dim,TRUE);
+ ELSE
+ CNormMat2(dim,EV);
+ END;
+
+ IF (ord # 0) THEN
+ IF (ord > 0) THEN ord:=1; ELSE ord:=-1; END;
+ CSortEwEv(EW,EV,dim,ord);
+ END;
+END ComEigEvEw;
+
+PROCEDURE ComEig(VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
+ VAR REV : ARRAY OF ARRAY OF LONGCOMPLEX; (* Rechts-EV *)
+ VAR LEV : ARRAY OF ARRAY OF LONGCOMPLEX; (* Links-EV *)
+ VAR EW : ARRAY OF LONGCOMPLEX; (* Eigenwerte *)
+ dim : CARDINAL;
+ genau : LONGREAL;
+ ord : INTEGER);
+
+ VAR i,j,k,m,Iter,MaxIter : CARDINAL;
+ eta,tse,tem,tep,tau,Max,Root1,Root2,Root : LONGREAL;
+ Hi,Hr,Hj,g,Te,Tee,Br,Bi,Er,Ei,Dr,Di,sig,isw : LONGREAL;
+ b,c,Ca,Cb,Ch,Cx,nC,e,s,Sa,Sb,Sh,Sx,d,De,nD : LONGREAL;
+ Cot2X,CotX,Cos2a,Sin2a,TanH : LONGREAL;
+ Aik,Aki,Aim,Ami,EVim,EVik,Zw,s1,s2,c1,c2,H : LONGCOMPLEX;
+ s1re,s1im,s2re,s2im : LONGREAL;
+ re,im : LONGREAL;
+ Mark : BOOLEAN;
+ En : POINTER TO ARRAY [0..MAX(INTEGER)-1] OF LONGREAL;
+BEGIN
+ ALLOCATE(En,dim*TSIZE(LONGREAL));
+ IF (genau > MachEpsR4) THEN genau:=MachEpsR4; END;
+ IF (genau < MachEps ) THEN genau:=MachEps; END;
+
+ MaxIter := 70; (* Sollte reichen, wie im Originalatikel *)
+
+ FOR i:=0 TO dim-1 DO (* Vorbelege Eigenvektormatritzen *)
+ FOR j:=0 TO dim-1 DO REV[i,j]:=zero; LEV[i,j]:=zero; END;
+ REV[i,i]:=one; LEV[i,i]:=one;
+ END;
+
+ Mark:=FALSE; Fehler:=FALSE; Iter:=0;
+ LOOP
+ INC(Iter);
+ IF Mark OR (Iter > MaxIter) THEN EXIT END; (* !!! *)
+ tau:=0.0;
+ FOR k:=0 TO dim-1 DO
+ tem:=0.0;
+ FOR i:=0 TO dim-1 DO
+ IF (i # k) THEN tem:=tem + ABS(RE(A[i,k])) + ABS(IM(A[i,k])); END;
+ END;
+ tau:=tau + tem;
+ En^[k]:=tem + ABS(RE(A[k,k])) + ABS(IM(A[k,k]));
+ END; (* FOR k *)
+
+ FOR k:=0 TO dim-2 DO (* Vertausche Zeilen und Spalten *)
+ Max:=En^[k]; i:=k;
+ FOR j:=k+1 TO dim-1 DO
+ IF (En^[j] > Max) THEN Max:=En^[j]; i:=j; END;
+ END;
+ IF (i # k) THEN
+ En^[i]:=En^[k];
+ FOR j:=0 TO dim-1 DO Zw:=A[k,j]; A[k,j]:=A[i,j]; A[i,j]:=Zw; END;
+ FOR j:=0 TO dim-1 DO
+ Zw:=A [j,k]; A [j,k]:=A [j,i]; A [j,i]:=Zw;
+ Zw:=REV[k,j]; REV[k,j]:=REV[i,j]; REV[i,j]:=Zw;
+ Zw:=LEV[k,j]; LEV[k,j]:=LEV[i,j]; LEV[i,j]:=Zw;
+ END;
+ END; (* IF i # k *)
+ END; (* FOR k *)
+
+ IF debug THEN
+ TIO.WrStr(" ComEig MaxIter,Iter, tau = ");
+ TIO.WrCard(Iter,4);
+ TIO.WrLngReal(tau,12,-5);
+ TIO.WrLn;
+ END;
+
+ IF (tau < genau) THEN EXIT END; (* !!! *)
+
+ Mark:=TRUE;
+ FOR k:=0 TO dim-1-1 DO
+ FOR m:=k+1 TO dim-1 DO
+ Hj:=0.0; g:=0.0; H:=zero;
+ FOR i:=0 TO dim-1 DO
+ IF (i # k) AND (i # m) THEN
+ H:=H + A[k,i]*A[m,i] - A[i,k]*A[i,m];
+ Te := RE(A[i,k])*RE(A[i,k]) + IM(A[i,k])*IM(A[i,k])
+ + RE(A[m,i])*RE(A[m,i]) + IM(A[m,i])*IM(A[m,i]);
+ Tee:= RE(A[i,m])*RE(A[i,m]) + IM(A[i,m])*IM(A[i,m])
+ + RE(A[k,i])*RE(A[k,i]) + IM(A[k,i])*IM(A[k,i]);
+ g := g + Te + Tee;
+ Hj := Hj - Te + Tee;
+ END; (* IF *)
+ END; (* FOR i *)
+ Hr:=RE(H); Hi:=IM(H);
+ Br:=RE(A[k,m]) + RE(A[m,k]); Bi:=IM(A[k,m]) + IM(A[m,k]);
+ Er:=RE(A[k,m]) - RE(A[m,k]); Ei:=IM(A[k,m]) - IM(A[m,k]);
+ Dr:=RE(A[k,k]) - RE(A[m,m]); Di:=IM(A[k,k]) - IM(A[m,m]);
+ Te := Br*Br + Ei*Ei + Dr*Dr;
+ Tee := Bi*Bi + Er*Er + Di*Di;
+ IF (Te >= Tee) THEN
+ isw:= 1.0; c:=Br; s:= Ei; d:=Dr; De:=Di; Root2:=sqrt(Te);
+ ELSE
+ isw:=-1.0; c:=Bi; s:=-Er; d:=Di; De:=Dr; Root2:=sqrt(Tee);
+ END;
+ Root1:= LongComplexMath.abs(CMPLX(s,c)); (* M.R. *)
+ Sa:=0.0;
+ IF (c >= 0.0) THEN Ca :=1.0; ELSE Ca :=-1.0; END;
+ IF (d >= 0.0) THEN sig:=1.0; ELSE sig:=-1.0; END;
+ IF (Root1 < MachEps) THEN
+ Sx:=0.0; Sa:=0.0; Cx:=1.0; Ca:=1.0;
+ IF (isw > 0.0) THEN e:=Er; b:=Bi; ELSE e:=Ei; b:=-Br; END;
+ nD:=d*d + De*De;
+ (* GOTO enter1; *)
+ ELSE
+ IF (ABS(s) > MachEps) THEN Ca := c / Root1; Sa := s / Root1; END;
+ Cot2X:= d / Root1;
+ CotX := Cot2X + (sig*sqrt(1.0 + Cot2X*Cot2X));
+ Sx := sig / sqrt(1.0 + CotX*CotX);
+ Cx := Sx*CotX;
+ eta := (Er*Br + Bi*Ei) / Root1;
+ tse := (Br*Bi - Er*Ei) / Root1;
+ Te := sig*( - Root1*De + tse*d) / Root2;
+ Tee := (d*De + Root1*tse) / Root2;
+ nD := Root2*Root2 + Tee*Tee;
+ Tee := Hj*Cx*Sx;
+ Cos2a:= Ca*Ca - Sa*Sa;
+ Sin2a:= 2.0*Ca*Sa;
+ tem := Hr*Cos2a + Hi*Sin2a;
+ tep := Hi*Cos2a - Hr*Sin2a;
+ Hr := Cx*Cx*Hr - Sx*Sx*tem - Ca*Tee;
+ Hi := Cx*Cx*Hi + Sx*Sx*tep - Sa*Tee;
+ b := isw*Te*Ca + eta*Sa;
+ e := eta*Ca - isw*Te*Sa;
+ END; (* IF (Root1 < MachEps) *)
+(*enter1:*) s := Hr - sig*Root2*e;
+ c := Hi - sig*Root2*b;
+ Root := LongComplexMath.abs(CMPLX(c,s)); (* M.R.*)
+ IF (Root < MachEps) THEN
+ Cb:=1.0; Ch:=1.0; Sb:=0.0; Sh:=0.0;
+ (* GOTO trans; *)
+ ELSE
+ Cb := - c / Root;
+ Sb := s / Root;
+ Tee := Cb*b - e*Sb;
+ nC := Tee*Tee;
+ TanH:= Root / (g + 2.0*(nC + nD));
+ Ch := 1.0 / sqrt(1.0 - TanH*TanH);
+ Sh := Ch*TanH;
+ END; (* IF *)
+(*trans:*) tem := Sx*Sh*(Sa*Cb - Sb*Ca);
+ im := - Sx*Sh*(Ca*Cb + Sa*Sb);
+ c1 := CMPLX(Cx*Ch - tem,im);
+ c2 := CMPLX(Cx*Ch + tem,im);
+ tep := Sx*Ch*Ca;
+ tem := Cx*Sh*Sb;
+ s1re := tep - tem;
+ s2re := - tep - tem;
+ tep := Sx*Ch*Sa;
+ tem := Cx*Sh*Cb;
+ s1im := tep + tem;
+ s2im := tep - tem;
+ s1 := CMPLX(s1re,s1im);
+ s2 := CMPLX(s2re,s2im);
+ tem := LongComplexMath.abs(s1); (* M.R.*)
+ tep := LongComplexMath.abs(s2);
+ IF (tem > MachEps) OR (tep > MachEps) THEN
+ Mark:=FALSE;
+ FOR i:=0 TO dim-1 DO (* Linksseitige Transformation *)
+ Aki:=A[k,i]; Ami:=A[m,i];
+ A[k,i] := c1*Aki + s1*Ami;
+ A[m,i] := s2*Aki + c2*Ami;
+ END;
+ FOR i:=0 TO dim-1 DO (* Linksvektoren *)
+ EVik:=LEV[k,i]; EVim:=LEV[m,i];
+ LEV[k,i] := c1*EVik + s1*EVim;
+ LEV[m,i] := s2*EVik + c2*EVim;
+ END;
+ FOR i:=0 TO dim-1 DO (* Rechtsseitige Transformation *)
+ Aik:=A[i,k]; Aim:=A[i,m];
+ A[i,k] := c2*Aik - s2*Aim;
+ A[i,m] := - s1*Aik + c1*Aim;
+ END;
+ FOR i:=0 TO dim-1 DO (* Rechtsvektoren *)
+ EVik:=REV[k,i]; EVim:=REV[m,i];
+ REV[k,i] := c2*EVik - s2*EVim;
+ REV[m,i] := - s1*EVik + c1*EVim;
+ END;
+ END; (* IF Transformation *)
+ END; (* FOR m *)
+ END; (* FOR k *)
+ (* tau2:=tau1; tau1:=tau; *)
+ END; (* LOOP *)
+ FOR i:=0 TO dim-1 DO EW[i]:=A[i,i]; END;
+ IF (Iter > MaxIter) THEN
+ Fehler:=TRUE; Fehlerflag:='MaxIter ueberschritten (ComEig).';
+ ErrOut(Fehlerflag);
+ END;
+ DEALLOCATE(En,dim*TSIZE(LONGREAL));
+END ComEig;
+
+END EigenLib3.