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.