IMPLEMENTATION MODULE EigenLib2;
(*------------------------------------------------------------------------*)
(* Stellt Routinen fuer die Eigenwertberechung von unsymmetrischen *)
(* realwertigen Matrizen und symmetrischer hermitische Matrizen zur *)
(* Verfuegung. *)
(* Routines for the calculation of eigensystems of unsymmetric real *)
(* general and hermitian matrices *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung: *)
(* *)
(* 22.06.93, MRi: Prozeduren Balance und BalBak aus EISPACK.f nach M2 *)
(* "ubersetzt. *)
(* 17.07.93, MRi: Durchsicht *)
(* 24.03.15, MRi: Korrektur in Balance, BalBak, angleichen an Algol *)
(* 01.04.16, MRi: Elimination von Compilerwarnungen (7 auf 0) *)
(* MinGenau auf MachEps bezogen. *)
(* 15.10.16, MRi: Aufruf von QuadGl in HessQR angepasst ... *)
(* 18.08.17, MRi: Erste Versionen von RealVecHes und ComVecHes *)
(* 03.09.17, MRi: Erstellen der ersten uebersetzbaren HQR2 Pascal-Version *)
(* 10.09.17, MRi: Erstellen der ersten uebersetzbaren HQR2 Modula-Version *)
(* 14.09.17, MRi: Hinzufuegen von TfmReaHes und BakReaHes1 *)
(* 16.09.17, MRi: Erstellen von HQR aus alter RHessQrEwEv routine *)
(* 20.09.17, MRi: Zusammenfassen aller Routinen in neuer EigenLib2 *)
(* NormOne auf offene Felder umgestellt *)
(* 22.09.17, MRi: ElmHess, ElmHessTrans auf offene Felder umgestellt *)
(* 23.09.17, MRi: RHessQR eingefuegt und auf ISO COMPLEX umgestellet *)
(* 25.09.17, MRi: Eigen in EigenRG umbenannt und komplett ueberarbeitet *)
(* Balance und BalBak auf offene Felder umgestellt *)
(* 27.09.17, MRi: Tred2,EWEV,GivTri,GivTriBak nach EigenLib1 verschoben *)
(* 07.10.17, MRi: Erstellen der ersten Verion von HTriDi und HTriBak *)
(* 08.10.17, MRi: Umstellen von HTriBak auf ISO COMPLEX und offene *)
(* Felder (CHTriBak) *)
(* 09.10.17, MRi: Umstellen von HTriDi auf ISO COMPLEX und offene *)
(* Felder (CHTriDi) *)
(* Die komplexe Arithmetik ist so gut es geht implizit in *)
(* CHTriDi und CHTriBak. ZTransQL wurde aus TransQL *)
(* portiert um mit komplexen Eigenvektoren umgehen zu *)
(* koennen. *)
(* 10.10.17, MRi: EigenCHM2 eingefuegt, Fehler in CHTriDi korrigiert *)
(* 20.10.17, MRi: Erstellen der ersten Verison on CJacobi2D *)
(* 21.10.17, MRi: Heftige Anpassung von CJacobi2D *)
(* 22.10.17, MRi: Erstellen der ersten Version von CJacobi1D *)
(*------------------------------------------------------------------------*)
(* Offene Punkte *)
(* *)
(* - Alles auf offene Felder umstellen *)
(* - Prozedur Eigen ueberarbeiten - die iterative ermittlung der EV soll *)
(* eingebaut werden. Dabei auch "Verschieben" der Eigenwerte wie in *)
(* EigenLib3.CHessInvIter einbauen *)
(* - Auslagern der Hilfroutinen soweit moeglich/sinnvoll *)
(* - Komplexe Arithmetik konsequent einfuehren *)
(* - In CJacobi1D/2D das Abbruchkriterim und die Maximalzahl der *)
(* Iterationen noch einmal ueberdenken *)
(*------------------------------------------------------------------------*)
(* Testroutinen in TstRGM und TstCHM *)
(*------------------------------------------------------------------------*)
(* Hint to sources of some routines *)
(* *)
(* - Most routine originate from the entire Algol procedures of the *)
(* "Handbuch" (Wilkinson/Reinsch) *)
(* - Routine RealVecHes, ComVecHes, TfmReaHes und BakReaHes1 are *)
(* adoptions of Algol60 routine form NUMAL (Numerical Algol library). *)
(* - Improvement for ElmHes and ElmTrans are taken form the PMATH library *)
(* PMATH is provided by Jean Debord, l'Universite de Limoges, France *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: EigenLib2.mod,v 1.5 2017/10/22 17:52:11 mriedl Exp mriedl $ *)
FROM SYSTEM IMPORT TSIZE,ADR;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
FROM Deklera IMPORT MaxDim,VEKTOR,CVEKTOR,MATRIX,CARDVEKTOR,INTVEKTOR;
FROM Errors IMPORT Fehler,Fehlerflag,ErrOut;
FROM LongMath IMPORT sqrt;
FROM LMathLib IMPORT MachEps,MachEpsR4,Small,Pythag;
IMPORT LongComplexMath;
FROM CmplxMath IMPORT CDIV;
FROM F77func IMPORT DMAX,MIN0;
FROM MatLib IMPORT CopyMat,CSortVek;
FROM EigenLibAux IMPORT NormOne,NormTwo;
IMPORT TIO;
CONST MinGenau = 10.0*MachEps;
CONST debug = FALSE;
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
CONST zero = LongComplexMath.zero;
one = LongComplexMath.one;
(*
CABS = LongComplexMath.abs;
scalarMult = LongComplexMath.scalarMult;
conj = LongComplexMath.conj;
*)
PROCEDURE CABS(z : LONGCOMPLEX) : LONGREAL;
BEGIN
RETURN Pythag(RE(z),IM(z));
END CABS;
PROCEDURE scalarMult(a : LONGREAL;
z : LONGCOMPLEX) : LONGCOMPLEX;
BEGIN
RETURN CMPLX(a*RE(z),a*IM(z));
END scalarMult;
PROCEDURE conj(x : LONGCOMPLEX) : LONGCOMPLEX;
BEGIN
RETURN CMPLX(RE(x),-IM(x));
END conj;
(*================= Start Hilfsroutinen NUMAL =============================*)
PROCEDURE TamVec( l,u : INTEGER;
i : INTEGER;
VAR a : MATRIX;
VAR b : VEKTOR) : LONGREAL;
(*----------------------------------------------------------------*)
(* scalar product of the column vector given in array a[l:u,i:i] *)
(* and the vector given in array b[l:u]. *)
(*----------------------------------------------------------------*)
VAR k : INTEGER;
s : LONGREAL;
BEGIN
s:=0.0;
FOR k:=l TO u DO s:=s + a[k,i]*b[k]; END;
RETURN s;
END TamVec;
PROCEDURE VecVec( l,u : INTEGER;
shift : INTEGER;
VAR a : VEKTOR;
VAR b : VEKTOR) : LONGREAL;
(*----------------------------------------------------------------*)
(* scalar product of the vector given in array a[l:u] and the *)
(* vector given in array b[l+shift:u+shift] *)
(*----------------------------------------------------------------*)
VAR k : INTEGER;
s : LONGREAL;
BEGIN
s:=0.0;
FOR k:=l TO u DO s:=s + a[k]*b[shift + k]; END;
RETURN s;
END VecVec;
PROCEDURE MatVec( l,u : INTEGER;
i : INTEGER;
VAR a : MATRIX;
VAR b : VEKTOR) : LONGREAL;
(*----------------------------------------------------------------*)
(* scalar product of the row vector given in array a[i:i,l:u] and *)
(* the vector given in array b[l:u] *)
(*----------------------------------------------------------------*)
VAR k : INTEGER;
s : LONGREAL;
BEGIN
s:=0.0;
FOR k:=l TO u DO s:=s + a[i,k]*b[k]; END;
RETURN s;
END MatVec;
PROCEDURE MatMat( l,u : INTEGER;
i,j : INTEGER;
VAR a : MATRIX;
VAR b : MATRIX) : LONGREAL;
VAR k : INTEGER;
s : LONGREAL;
BEGIN
s:=0.0;
FOR k:=l TO u DO s:=s + a[i,k]*b[k,j]; END;
RETURN s;
END MatMat;
PROCEDURE ichrow( l,u : INTEGER;
i,j : INTEGER;
VAR a : MATRIX);
VAR r : LONGREAL;
k : INTEGER;
BEGIN
FOR k:=l TO u DO
r := a[i,k]; a[i,k] := a[j,k]; a[j,k] := r;
END;
END ichrow;
PROCEDURE ichcol( l,u : INTEGER;
i,j : INTEGER;
VAR a : MATRIX);
VAR r : LONGREAL;
k : INTEGER;
BEGIN
FOR k:=l TO u DO
r:= a[k,i]; a[k,i]:= a[k,j]; a[k,j]:= r;
END;
END ichcol;
(*============================== Sektion "Eigen" ==========================*)
PROCEDURE Balance(VAR A : ARRAY OF ARRAY OF LONGREAL;
dim : CARDINAL;
VAR low,high : CARDINAL;
VAR Scale : ARRAY OF LONGREAL);
(*----------------------------------------------------------------*)
(* Noch weiter Austesten *)
(*----------------------------------------------------------------*)
VAR i,j,k,l,n : CARDINAL;
PROCEDURE exc(m : CARDINAL); (* lokal zu Balance *)
VAR Zw : LONGREAL;
i : CARDINAL;
BEGIN
Scale[m]:=LFLOAT(j);
IF (j # m) THEN
FOR i:=0 TO k DO Zw:=A[i,j]; A[i,j]:=A[i,m]; A[i,m]:=Zw; END;
FOR i:=l TO n-1 DO Zw:=A[j,i]; A[j,i]:=A[m,i]; A[m,i]:=Zw; END;
END;
END exc;
CONST radix = 8.0; (* Basis der Flie\3kommadarstellung *)
sqrdx = radix*radix;
VAR s,r,g,f,c : LONGREAL;
Konvergenz,cont,null : BOOLEAN;
BEGIN
n:=VAL(INTEGER,dim);
l:=0; k:=n-1;
REPEAT
(* Suche Zeilen die einen Eigenwert isolieren *)
(* und schiebe diese nach unten. *)
LOOP (* outer *)
cont:=FALSE;
(* for j:=k step -1 until 1 do begin *)
FOR j:=k TO 0 BY -1 DO
LOOP (* inner *)
null:=FALSE;
IF (j > 0) THEN
FOR i:=0 TO j-1 DO
null:=(A[j,i] = 0.0);
IF NOT null THEN EXIT; (* exit inner loop *) END;
END;
END;
FOR i:=j+1 TO k DO
null:=(A[j,i] = 0.0);
IF NOT null THEN EXIT; (* exit inner loop *) END;
END;
EXIT; (* inner loop *)
END; (* inner LOOP *)
IF null THEN
exc(k); DEC(k); cont:=TRUE;
EXIT; (* outer loop *)
END;
END; (* FOR j *)
EXIT; (* exit outer loop *)
END; (* outer LOOP *)
UNTIL NOT cont;
REPEAT
(* Suche Spalten, die einen Eigenwert isolieren *)
(* und schiebe diese nach links. *)
LOOP (* outer *)
cont:=FALSE;
FOR j:=l TO k DO
LOOP (* inner *)
null:=FALSE;
IF (j > 0) THEN
FOR i:=l TO j-1 DO
null:=(A[i,j] = 0.0);
IF NOT null THEN EXIT; (* exit inner loop *) END;
END;
END;
FOR i:=j+1 TO k DO
null:=(A[i,j] = 0.0);
IF NOT null THEN EXIT; (* exit inner loop *) END;
END;
EXIT; (* exit inner loop *)
END; (* inner LOOP *)
IF null THEN
exc(l); INC(l); cont:=TRUE;
EXIT; (* exit outer loop *)
END;
END; (* FOR j *)
EXIT; (* exit outer loop *)
END; (* outer LOOP *)
UNTIL NOT cont;
low := VAL(CARDINAL,l+1);
high := VAL(CARDINAL,k+1);
FOR i:=l TO k DO Scale[i]:=1.0; END;
REPEAT (* Balanciere die Submatrix A[i,l..k] *)
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(A[j,i]); r:=r + ABS(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]:=A[i,j] * g; END;
FOR j:=0 TO k DO A[j,i]:=A[j,i] * f; END;
END;
END; (* IF *)
END; (* FOR i *)
UNTIL Konvergenz;
END Balance;
PROCEDURE BalBak( dim : CARDINAL; (* Ordnung der Matrix Z *)
low,high : CARDINAL; (* Von Balance ermittelt *)
VAR Scale : ARRAY OF LONGREAL; (* Von Balance ermittelt *)
m : CARDINAL;
VAR Z : ARRAY OF ARRAY OF LONGREAL); (* Eigenvektoren *)
VAR i,j,k : CARDINAL;
S,Zw : LONGREAL;
BEGIN
DEC(low); DEC(high);
IF (m = 0) THEN RETURN END;
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]:=Z[j,i]*S; END;
END;
END;
(* for i:=low-1 step -1 until 1,hi + 1 step 1 until n do begin *)
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 BalBak;
PROCEDURE GivHess(VAR A : MATRIX;
dim : CARDINAL;
genau : LONGREAL);
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
VAR i,j,k,p : CARDINAL;
Omega,cos,sin,Aij,Aik,Aki : LONGREAL;
BEGIN
FOR j:=1 TO dim-2 DO
p:=j+1;
FOR i:=j+2 TO dim DO
IF (A[i,j] # 0.0) THEN
Aij:=A[i,j];
IF (ABS(A[p,j]) < genau*ABS(Aij)) THEN
Omega:=-Aij;
cos:=0.0; sin:=1.0;
ELSE
IF (A[p,j] >= 0.0) THEN
Omega:=sqrt(sqr(A[p,j])+Aij*Aij)
ELSE
Omega:=-sqrt(sqr(A[p,j])+Aij*Aij);
END;
cos:=A[p,j]/Omega; sin:=-Aij/Omega;
END; (* IF *)
A[p,j]:=Omega;
A[i,j]:=0.0;
FOR k:=p TO dim DO
Aik:=A[i,k];
A[i,k]:=sin*A[p,k]+cos*Aik;
A[p,k]:=cos*A[p,k]-sin*Aik;
END; (* FOR k *)
FOR k:=1 TO dim DO
Aki:=A[k,i];
A[k,i]:=sin*A[k,p]+cos*Aki;
A[k,p]:=cos*A[k,p]-sin*Aki;
END; (* FOR k *)
END; (* IF Aij *)
END; (* FOR i *)
END; (* FOR j *)
END GivHess;
PROCEDURE ElmHess(VAR A : ARRAY OF ARRAY OF LONGREAL;
dim : CARDINAL;
k,l : CARDINAL; (* Low,High *)
VAR Index : ARRAY OF CARDINAL);
VAR i,j,m,la : CARDINAL;
x,y : LONGREAL;
BEGIN
DEC(k); DEC(l); (* Wg. offener Felder *)
la:=l-1;
FOR m:=k+1 TO la DO
i:=m; x:=0.0;
FOR j:=m TO l DO
IF (ABS(A[j,m-1]) > ABS(x)) THEN x:=A[j,m-1]; i:=j; END;
END;
Index[m]:=i; (* Fehler, war in FOR j - Schleife, 23.03.16 *)
IF (i # m) THEN (* Warung: Bedingung immer falsch !!! *)
(* interchange rows and columns of array A *)
FOR j:=m-1 TO dim-1 DO y:=A[i,j]; A[i,j]:=A[m,j]; A[m,j]:=y; END;
FOR j:=0 TO l DO y:=A[j,i]; A[j,i]:=A[j,m]; A[j,m]:=y; END;
END;
IF (x # 0.0) THEN
FOR i:=m+1 TO l DO
y:=A[i,m-1];
IF (y # 0.0) THEN
y:=y / x; A[i,m-1]:=y;
FOR j:=m TO dim-1 DO A[i,j]:=A[i,j] - y*A[m,j]; END;
FOR j:=0 TO l DO A[j,m]:=A[j,m] + y*A[j,i]; END;
END;
END; (* FOR i *)
END; (* IF *)
END; (* FOR m *)
END ElmHess;
PROCEDURE ElmHessTrans(VAR V : ARRAY OF ARRAY OF LONGREAL;
VAR H : ARRAY OF ARRAY OF LONGREAL;
VAR Index : ARRAY OF CARDINAL;
dim : CARDINAL;
low,upp : CARDINAL);
VAR i,j,k : CARDINAL;
BEGIN
DEC(low); DEC(upp); (* Wg. offener Felder *)
FOR i:=0 TO dim-1 DO
FOR j:=0 TO dim-1 DO V[i,j]:=0.0; END;
V[i,i]:=1.0;
END;
FOR i:=upp-1 TO low+1 BY -1 DO
j:=Index[i];
FOR k:=i+1 TO upp DO V[i,k]:=H[k,i-1]; END;
IF (i # j) THEN
FOR k:=i TO upp DO V[k,i]:=V[k,j]; V[k,j]:=0.0; END;
V[i,j]:=1.0;
END;
END;
END ElmHessTrans;
PROCEDURE RHessQR(VAR H : MATRIX; (* Hessenbergmatrix *)
VAR EW : ARRAY OF LONGCOMPLEX; (* Eigenwerte *)
dim : CARDINAL);
VAR nn,m,l,k,j,its,i,Mmin : CARDINAL;
z,y,x,w,v,u,t,s,r,q,p,ANorm,xp,xz : LONGREAL;
Konvergenz : BOOLEAN;
BEGIN
p:=0.0; q:=0.0; r:=0.0; (* Wg. Compilerwarnungen *)
ANorm := ABS(H[1,1]);
FOR i:=2 TO dim DO
FOR j:=i-1 TO dim DO ANorm := ANorm+ABS(H[i,j]); END;
END;
nn:=dim; t:=0.0;
WHILE (nn >= 1) DO
its := 0;
REPEAT (* UNTIL Konvergenz *)
Konvergenz:=TRUE; (* Wird getestet ! *)
l:=nn;
LOOP
IF (l = 1) THEN EXIT END;
s := ABS(H[l-1,l-1]) + ABS(H[l,l]);
IF s = 0.0 THEN s := ANorm; END;
IF ((ABS(H[l,l-1]) + s) = s) THEN EXIT END;
DEC(l);
END;
x := H[nn,nn];
IF (l = nn) THEN
DEC(nn);
EW[nn]:=CMPLX(x+t,0.0);
ELSE
y := H[nn-1,nn-1];
w := H[nn,nn-1]*H[nn-1,nn];
IF (l = nn-1) THEN
p := 0.5*(y-x);
q := p*p + w;
z := sqrt(ABS(q));
x := x+t;
IF (q >= 0.0) THEN
(* z := p + sign(z,p); *)
IF (p < 0.0) THEN z:= p - z; ELSE z:= p + z; END;
xz:=x+z;
IF (z # 0.0) THEN
EW[nn-1]:=CMPLX(x - w/z,0.0);
ELSE
EW[nn-1]:=CMPLX(xz,0.0);
END;
EW[nn-2]:=CMPLX(xz,0.0);
ELSE
xp:=x+p;
EW[nn-1]:=CMPLX(xp, z);
EW[nn-2]:=CMPLX(xp,-z);
END;
DEC(nn,2);
ELSE
IF its = 30 THEN
Fehler:=TRUE; Fehlerflag:='MaxIter ueberschritten !';
ErrOut(Fehlerflag); RETURN;
END;
IF (its = 10) OR (its = 20) THEN
t := t+x;
FOR i:=1 TO nn DO H[i,i]:=H[i,i] - x; END;
s := ABS(H[nn,nn-1])+ABS(H[nn-1,nn-2]);
x := 0.75*s;
y := x;
w := -0.4375*s*s;
END;
INC(its);
m := nn-1;
LOOP
DEC(m);
IF (m < l) THEN EXIT END;
z := H[m,m];
r := x-z;
s := y-z;
p := (r*s-w) / H[m+1,m] + H[m,m+1];
q := H[m+1,m+1]-z-r-s;
r := H[m+2,m+1];
s := ABS(p) + ABS(q) + ABS(r);
p := p / s;
q := q / s;
r := r / s;
IF (m = l) THEN EXIT END;
u := ABS(H[m,m-1])*(ABS(q) + ABS(r));
v := ABS(p)*(ABS(H[m-1,m-1]) + ABS(z) + ABS(H[m+1,m+1]));
IF ((u+v) = v) THEN EXIT END;
END;
FOR i:=m+2 TO nn DO
H[i,i-2] := 0.0;
IF (i # m+2) THEN H[i,i-3]:=0.0; END;
END;
FOR k:=m TO nn-1 DO
IF (k # m) THEN
p := H[k,k-1];
q := H[k+1,k-1];
IF (k # nn-1) THEN
r := H[k+2,k-1]
ELSE
r := 0.0;
END;
x := ABS(p) + ABS(q) + ABS(r);
IF (x # 0.0) THEN
p := p / x;
q := q / x;
r := r / x
END;
END; (* IF k *)
s := sqrt(p*p + q*q + r*r);
IF (p < 0.0) THEN s:= - s; END;
IF (s # 0.0) THEN
IF (k = m) THEN
IF (l # m) THEN H[k,k-1] := -H[k,k-1]; END;
ELSE
H[k,k-1] := -s*x;
END;
p := p + s;
x := p / s;
y := q / s;
z := r / s;
q := q / p;
r := r / p;
FOR j:=k TO nn DO
p := H[k,j]+q*H[k+1,j];
IF (k # nn-1) THEN
p := p+r*H[k+2,j];
H[k+2,j] := H[k+2,j]-p*z
END;
H[k+1,j] := H[k+1,j] - p*y;
H[k ,j] := H[k ,j] - p*x
END;
IF ((k+3) < nn) THEN Mmin:=k+3 ELSE Mmin:=nn END;
FOR i:=l TO Mmin DO
p := x*H[i,k]+y*H[i,k+1];
IF (k # nn-1) THEN
p:=p + z*H[i,k+2];
H[i,k+2]:=H[i,k+2] - p*r;
END;
H[i,k+1]:=H[i,k+1] - p*q;
H[i,k ]:=H[i,k ] - p;
END;
END; (* IF s *)
END; (* FOR k *)
Konvergenz:=FALSE;
END; (* ELSE (l # nn-1) *)
END; (* ELSE (l # nn) *)
UNTIL Konvergenz;
END; (* WHILE (nn >= 1) *)
END RHessQR;
PROCEDURE HQR(VAR H : MATRIX; (* Wird zerst"ort *)
VAR EW : ARRAY OF LONGCOMPLEX; (* Eigenwerte *)
n : CARDINAL; (* Dimension von H *)
low,upp : CARDINAL; (* Ausgabe von Balence *)
VAR cnt : ARRAY OF INTEGER);
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
CONST eps = 10.0*MachEps;
MaxIter = 60;
VAR i,j,k,l,m : CARDINAL;
its,en,na : CARDINAL;
p,q,r,s,t,w,x,y,z : LONGREAL;
Norm,ra,sa,a,b,tmp : LONGREAL;
ctmp,v : LONGCOMPLEX; (* Hilfsvariablen *)
Skip,notlast : BOOLEAN;
BEGIN (* hqr2 *)
FOR i:=1 TO low-1 DO EW[i-1]:=CMPLX(H[i,i],0.0); cnt[i]:=0; END;
FOR i:=upp+1 TO n DO EW[i-1]:=CMPLX(H[i,i],0.0); cnt[i]:=0; END;
p:=0.0; q:=0.0; r:=0.0; s:=0.0; z:=0.0; (* Wg. Compilerwarnungen *)
Norm := ABS(H[1,1]);
FOR i:=2 TO n DO
FOR j:=i-1 TO n DO Norm:=Norm + ABS(H[i,j]); END;
END;
en:=upp; t:=0.0;
WHILE (en >= low) DO (* Lauf Eigenwerte (nextw:) *)
its:=0; na:=en-1;
LOOP (* Iterationen (nextit:) *)
INC(its);
l:=en;
LOOP (* Suche kleines Subdiagonalelement *)
IF (l = low) THEN EXIT END;
s := ABS(H[l-1,l-1]) + ABS(H[l,l]);
IF s = 0.0 THEN s := Norm; END;
IF ((ABS(H[l,l-1]) + s) = s) THEN EXIT END;
DEC(l);
END;
x := H[en,en];
IF (l = en) THEN (* Eine Wurzel gefunden (onew) *)
EW[en-1] := CMPLX(x + t,0.0);
H[en,en] := x + t;
cnt[en-1] := its;
DEC(en);
EXIT; (* LOOP Iterationen *) (* !!! *)
END;
y:=H[na,na];
w:=H[en,na]*H[na,en];
IF (l = na) THEN (* Wurzelpaar gefunden (twow) *)
p := 0.5*(y - x);
q := p*p + w;
z := sqrt(ABS(q));
H[en,en] := x + t; x:=H[en,en];
H[na,na] := y + t;
cnt[en-1] := - VAL(INTEGER,its);
cnt[na-1] := its;
IF (q > 0.0) THEN (* Reeles Wurzelpaar *)
IF (p < 0.0) THEN z := p - z; ELSE z:= p + z; END;
EW[na-1] := CMPLX(x + z,0.0);
s := x - (w / z);
EW[en-1] := CMPLX(s,0.0);
x := H[en,na];
r := sqrt(x*x + z*z);
p := x / r;
q := z / r;
FOR j:=na TO n DO (* Zeilenmodifikation *)
z := H[na,j];
H[na,j] := q*z + p*H[en,j];
H[en,j] := q*H[en,j] - p*z;
END;
FOR i:=1 TO en DO (* Spaltenmodifikation *)
z := H[i,na];
H[i,na] := q*z + p*H[i,en];
H[i,en] := q*H[i,en] - p*z;
END;
ELSE (* Komplexes Wurzelpaar *)
EW[na-1] := CMPLX(x + p,z);
EW[en-1] := conj(EW[na-1]);
END; (* Komplexes Wurzelpaar *)
DEC(en,2);
EXIT; (* LOOP Iterationen *) (* !!! *)
END; (* IF *)
IF (its > MaxIter) THEN
Fehlerflag:=
"MaxIter ueberschritten (RHessQREwEv) Routine abgebrochen !";
Fehler:=TRUE; ErrOut(Fehlerflag);
RETURN;
END;
IF ((its MOD 10) = 0) AND (en > 2) THEN (* Berechne Ausnahmeschift *)
t:=t + x;
FOR i:=low TO en DO H[i,i]:=H[i,i] - x; END;
s := ABS(H[en,na]) + ABS(H[na,en-2]);
x := 0.75*s; y:=x;
w := -0.4375*s*s;
END;
m := en-1;
LOOP
DEC(m);
IF (m < l) THEN EXIT END;
z := H[m,m];
r := x - z;
s := y - z;
p := (r*s - w) / H[m+1,m] + H[m,m+1];
q := H[m+1,m+1] - z - r - s;
r := H[m+2,m+1];
s := ABS(p) + ABS(q) + ABS(r);
p := p / s;
q := q / s;
r := r / s;
IF (m = l) THEN EXIT END;
a := ABS(H[m,m-1])*(ABS(q) + ABS(r));
b := ABS(p)*(ABS(H[m-1,m-1]) + ABS(z) + ABS(H[m+1,m+1]));
IF ((a+b) = b) THEN EXIT END;
END;
FOR i:=m+2 TO en DO H[i,i-2]:=0.0; END;
FOR i:=m+3 TO en DO H[i,i-3]:=0.0; END;
FOR k:=m TO na DO (* Doppelter QR-Schritt, ver"andert die Zeilen *)
notlast:= k # na; (* [l..en] und Spalten [m ..en] von H. *)
Skip:=FALSE;
IF (k # m) THEN
p:=H[k,k-1];
q:=H[k+1,k-1];
IF notlast THEN r := H[k+2,k-1]; ELSE r:=0.0; END;
x := ABS(p) + ABS(q) + ABS(r);
IF (x = 0.0) THEN
Skip:=TRUE;
ELSE
p:=p / x; q:=q / x; r:=r / x;
END;
END;
IF NOT Skip THEN
s := sqrt(p*p + q*q + r*r);
IF (p < 0.0) THEN s:= -s; END;
IF (k # m) THEN
H[k,k-1]:= - s*x;
ELSIF (l # m) THEN
H[k,k-1] := - H[k,k-1];
END;
p := p + s;
x := p / s;
y := q / s;
z := r / s;
q := q / p;
r := r / p;
FOR j:=k TO n DO (* Zeilentransformation *)
p := H[k,j] + q*H[k+1,j];
IF notlast THEN
p := p + r*H[k+2,j];
H[k+2,j] := H[k+2,j] - p*z
END;
H[k+1,j] := H[k+1,j] - p*y;
H[k ,j] := H[k ,j] - p*x
END;
IF ((k+3) < en) THEN j:=k + 3; ELSE j:=en; END;
FOR i:=1 TO j DO (* Spaltentransformation *)
p := x*H[i,k] + y*H[i,k+1];
IF notlast THEN
p:=p + z*H[i,k+2];
H[i,k+2]:=H[i,k+2] - p*r;
END;
H[i,k+1]:=H[i,k+1] - p*q;
H[i,k ]:=H[i,k ] - p;
END;
END; (* IF Skip *)
END; (* FOR k *)
END; (* LOOP Iterationen *)
END; (* WHILE Eigenwerte *)
Norm:=0.0; k:=1;
FOR i:=1 TO n DO
FOR j:=k TO n DO Norm:=Norm + ABS(H[i,j]); END;
k:=i;
END;
FOR en:=n TO 1 BY -1 DO (* R"ucksubstitution *)
p:=RE(EW[en-1]);
q:=IM(EW[en-1]);
na := en - 1;
IF (q = 0.0) THEN (* Reeler Vektor *)
m:=en;
H[en,en]:=1.0;
FOR i:=na TO 1 BY -1 DO
w := H[i,i] - p;
r := H[i,en];
FOR j:=m TO na DO r:=r + H[i,j]*H[j,en]; END;
IF (IM(EW[i-1]) < 0.0) THEN
z := w; s := r;
ELSE
m:=i;
IF (IM(EW[i-1]) = 0.0) THEN
IF (w # 0.0) THEN
H[i,en] := - r / w;
ELSE
H[i,en] := - r / MachEps*Norm;
END;
ELSE (* L"ose |w x| |H[i,en | = |-r| *)
x := H[i,i+1]; (* |y z| |H[i+1,en| = |-s| *)
y := H[i+1,i];
q := RE(EW[i-1]) - p;
q := q*q + sqr(IM(EW[i-1]));
(* WITH EW[i] DO q := (re-p)*(re-p) + im*im; END; *)
t := (x*s - z*r) / q; H[i,en] := t;
IF ABS(x) > ABS(z) THEN
H[i+1,en] := (-r - w*t) / x;
ELSE
H[i+1,en] := (-s - y*t) / z;
END;
END; (* IF *)
END; (* ELSE (W[i].im >= 0.0) *)
END; (* FOR i *)
ELSIF (q < 0.0) THEN (* Komplexer Vektor von \lambda = p - \imag q *)
m:=na;
IF ABS(H[en,na]) > ABS(H[na,en]) THEN
H[na,na] := - (H[en,en] - p) / H[en,na];
H[na,en] := - q / H[en,na];
ELSE
ctmp := CMPLX(-H[na,en],0.0) / CMPLX(H[na,na]-p,q);
H[na,na]:=RE(ctmp); H[na,en]:=IM(ctmp);
H[en,na]:=1.0; H[en,en]:=0.0;
END;
FOR i:=na-1 TO 1 BY -1 DO
w := H[i,i] - p;
ra := H[i,en];
sa := 0.0;
FOR j:=m TO na DO
ra:=ra + H[i,j]*H[j,na];
sa:=sa + H[i,j]*H[j,en];
END;
IF (IM(EW[i-1]) < 0.0) THEN
z:=w; r:=ra; s:=sa;
ELSE
m:=i;
IF (IM(EW[i-1]) = 0.0) THEN
ctmp := CMPLX(-ra,-sa) / CMPLX(w,q);
H[i,na]:=RE(ctmp); H[i,en]:=IM(ctmp);
ELSE (* L"ose komplexe Gleichng *)
x := H[i,i+1]; y := H[i+1,1];
tmp := RE(EW[i-1]) - p;
v := CMPLX(tmp*tmp + sqr(IM(EW[i-1])) - q*q, 2.0*tmp*q);
IF (RE(v) = 0.0) AND (IM(v) = 0.0) THEN
tmp := MachEps*Norm*
(ABS(w) + ABS(q) + ABS(x) + ABS(y) + ABS(z));
v := CMPLX(tmp,0.0);
END;
ctmp := CMPLX(x*r-z*ra+q*sa ,x*s-z*sa-q*ra) / v;
H[i,na]:=RE(ctmp); H[i,en]:=IM(ctmp);
IF ABS(x) > (ABS(z) + ABS(q)) THEN
H[i+1,na] := ( - ra - w*H[i,na] + q*H[i,en]) / x;
H[i+1,en] := ( - sa - w*H[i,en] - q*H[i,na]) / x;
ELSE
ctmp:=CMPLX(-r-y*H[i,na] ,-s-y*H[i,en]) / CMPLX(z,q);
H[i+1,na]:=RE(ctmp); H[i+1,en]:=IM(ctmp);
END;
END; (* IF W[i].im > 0.0 *)
END; (* IF W[i].im >= 0.0 *)
END; (* FOR i *)
END; (* IF komplexer Vektor *)
END; (* FOR en *)
END HQR;
PROCEDURE HQR2( N : INTEGER;
Low : INTEGER;
High : INTEGER;
VAR H : MATRIX;
VAR Z : MATRIX;
VAR Wr : VEKTOR;
VAR Wi : VEKTOR;
VAR iErr : INTEGER);
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
CONST CmplxMath = TRUE;
VAR i,j,k,l,m,en,ii,jj,na,nn : INTEGER;
itn,its,mp2,enm2 : INTEGER;
p,q,r,s,t,w,x,y,zz : LONGREAL;
norm,tst1,tst2 : LONGREAL;
ra,sa,vi,vr : LONGREAL;
notlas,found,done : BOOLEAN;
ctmp : LONGCOMPLEX;
BEGIN
iErr:=0;
norm:=0.0;
k:=1;
(* store roots isolated by balanc and compute matrix norm *)
FOR i:=1 TO N DO
FOR j:=k TO N DO
norm:=norm + ABS(H[i,j]);
END; (* FOR *)
k:=i;
IF ((i < Low) OR (i > High)) THEN
Wr[i] := H[i,i];
Wi[i] := 0.0;
END; (* IF *)
END; (* FOR *)
p := 0.0;
q := 0.0;
r := 0.0;
s := 0.0;
zz := 0.0;
(*
m := 0;
l := 0;
*)
en:=High;
t:=0.0;
itn:=30*N + 2;
(* search for next eigenvalues *)
LOOP
IF (en < Low) THEN
(* all roots found. backsubstitute to find *)
(* vectors of upper triangular form *)
IF (norm <> 0.0) THEN
(* for en:=n step -1 until 1 do *)
FOR nn:=1 TO N DO
en := N + 1 - nn;
p := Wr[en];
q := Wi[en];
na := en - 1;
IF (q < 0.0) THEN
(* complex vector *)
m := na;
(* last vector component chosen imaginary so that *)
(* eigenvector matrix is triangular *)
IF (ABS(H[en,na]) <= ABS(H[na,en])) THEN
IF CmplxMath THEN
ctmp:=CMPLX(0.0,-H[na,en]) / CMPLX(H[na,na]-p,q);
H[na,na]:=RE(ctmp); H[na,en]:=IM(ctmp);
ELSE
CDIV(0.0,-H[na,en],H[na,na]-p,q,H[na,na],H[na,en]);
END;
ELSE
H[na,na] := q / H[en,na];
H[na,en] := - (H[en,en] - p) / H[en,na];
END; (* IF *)
H[en,na]:=0.0;
H[en,en]:=1.0;
enm2:=na - 1;
IF (enm2 <> 0) THEN
(* for i:=en-2 step -1 until 1 do *)
FOR ii:=1 TO enm2 DO
i:=na - ii;
w:=H[i,i] - p;
ra:=0.0;
sa:=0.0;
FOR j:=m TO en DO
ra:=ra + H[i,j]*H[j,na];
sa:=sa + H[i,j]*H[j,en];
END; (* FOR *)
IF (Wi[i] >= 0.0) THEN
m:=i;
IF (Wi[i] <> 0.0) THEN
(* solve complex linear system: *)
(* *)
(* |w+i*q x| |H[i,na] +i*H[i,en] | |-ra+i*sa| *)
(* | | | | = | | *)
(* |y z+i*q| |H[i+1,na]+i*H[i+1,en]| |-r+i*s | *)
x := H[i,i+1];
y := H[i+1,i];
vr := sqr(Wr[i] - p) + sqr(Wi[i]) - q*q;
vi := 2.0*q*(Wr[i] - p);
IF ((vr = 0.0) AND (vi = 0.0)) THEN
tst1:=norm*
(ABS(w) + ABS(q) + ABS(x) + ABS(y) + ABS(zz));
vr:=tst1;
REPEAT;
vr:=0.01*vr;
tst2:=tst1 + vr;
UNTIL (tst2 <= tst1);
END; (* IF *)
IF CmplxMath THEN
ctmp := CMPLX(x*r-zz*ra+q*sa,x*s-zz*sa-q*ra) /
CMPLX(vr,vi);
H[i,na] := RE(ctmp); H[i,en]:=IM(ctmp);
ELSE
CDIV(x*r-zz*ra+q*sa,x*s-zz*sa-q*ra,vr,vi,
H[i,na],H[i,en]);
END;
IF (ABS(x) <= ABS(zz)+ABS(q)) THEN
IF CmplxMath THEN
ctmp := CMPLX(-r-y*H[i,na],-s-y*H[i,en]) /
CMPLX(zz,q);
H[i+1,na] := RE(ctmp); H[i+1,en]:=IM(ctmp);
ELSE
CDIV(-r-y*H[i,na],-s-y*H[i,en],zz,q,
H[i+1,na],H[i+1,en]);
END;
ELSE
H[i+1,na] := (-ra - w*H[i,na] + q*H[i,en]) / x;
H[i+1,en] := (-sa - w*H[i,en] - q*H[i,na]) / x;
END; (* IF *)
ELSE
IF CmplxMath THEN
ctmp:=CMPLX(-ra,-sa) / CMPLX(w,q);
H[i,na]:=RE(ctmp); H[i,en]:=IM(ctmp);
ELSE
CDIV(-ra,-sa,w,q,H[i,na],H[i,en]);
END;
END; (* IF *)
(* overflow control *)
t:=DMAX(ABS(H[i,na]),ABS(H[i,en]));
IF (t <> 0.0) THEN
tst1:=t;
tst2:=tst1 + 1.0/tst1;
IF (tst2 <= tst1) THEN
FOR j:=i TO en DO
H[j,na]:=H[j,na] / t;
H[j,en]:=H[j,en] / t;
END; (* FOR *)
END; (* IF *)
END; (* IF *)
ELSE
zz := w;
r := ra;
s := sa;
END; (* IF *)
END; (* FOR *)
END; (* IF *)
ELSIF (q = 0.0) THEN
(* real vector *)
m:=en;
H[en,en]:=1.0;
IF (na <> 0) THEN
(* for i:=en-1 step -1 until 1 do *)
FOR ii:=1 TO na DO
i:=en - ii;
w:=H[i,i] - p;
r:=0.0;
FOR j:=m TO en DO
r:=r + H[i,j]*H[j,en];
END; (* FOR *)
IF (Wi[i] >= 0.0) THEN
m:=i;
IF (Wi[i] <> 0.0) THEN
(* solve real linear equations *)
(* *)
(* | w x | | H[i,en] | | -r | *)
(* | | | | = | | *)
(* | y z | | H[i+1,en] | | -s | *)
x:=H[i,i+1];
y:=H[i+1,i];
q:= sqr(Wr[i] - p) + Wi[i]*Wi[i];
t:=(x*s-zz*r)/q;
H[i,en]:=t;
IF (ABS(x) <= ABS(zz)) THEN
H[i+1,en]:=(-s-y*t)/zz;
ELSE
H[i+1,en]:=(-r-w*t)/x;
END; (* IF *)
ELSE
t:=w;
IF (t = 0.0) THEN
tst1:=norm;
t:=tst1;
REPEAT
t:=0.01*t;
tst2:=norm + t;
UNTIL (tst2 <= tst1);
END; (* IF *)
H[i,en]:= - r/t;
END; (* IF *)
(* overflow control *)
t:=ABS(H[i,en]);
IF (t <> 0.0) THEN
tst1:=t;
tst2:=tst1 + 1.0/tst1;
IF (tst2 <= tst1) THEN
FOR j:=i TO en DO
H[j,en]:=H[j,en] / t;
END; (* FOR *)
END; (* IF *)
END; (* IF *)
ELSE
zz:=w;
s:=r;
END; (* IF real vector *)
END; (* FOR ii *)
END; (* IF *)
END; (* IF *)
(* end complex vector *)
END; (* FOR nn (end back substitution) *)
(* vectors of isolated roots *)
FOR i:=1 TO N DO
IF ((i < Low) OR (i > High)) THEN
FOR j:=i TO N DO
Z[j,i] := H[i,j];
END; (* FOR *)
END; (* IF *)
END; (* FOR *)
(* multiply by transformation matrix to give *)
(* vectors of original full matrix *)
(* for j:=n step -1 until low do *)
FOR jj:=Low TO N DO
j:=N + Low - jj;
m:=MIN0(j,High);
FOR i:=Low TO High DO
zz:=0.0;
FOR k:=Low TO m DO
zz:=zz + Z[k,i]*H[k,j];
END; (* FOR *)
Z[j,i]:=zz;
END; (* FOR *)
END; (* FOR *)
END; (* IF *)
EXIT; (******)
ELSE
its:=0;
na := en - 1;
enm2 := na - 1;
END; (* IF en *)
done:=FALSE;
WHILE (NOT done) DO
done:=TRUE;
(* look for single small sub-diagonal element *)
(* for l:=en step -1 until low+1 do *)
(* if abs(h[l,l-1]) <= macheps*(abs(h[l-1,l-1]) + abs(h[l,l])) *)
(* then *)
(* goto cont1; *)
(* l:=low; *)
(* cont1: *)
l:=en; found:=FALSE;
WHILE (l > Low) AND NOT found DO
s := ABS(H[l-1,l-1]) + ABS(H[l,l]);
IF (s = 0.0) THEN
s:=norm;
END;
tst1:=s;
tst2:=tst1 + ABS(H[l,l-1]);
found := (tst2 = tst1);
IF NOT found THEN
l:=l - 1; (* DEC(l); *)
END;
END; (* WHILE *)
IF NOT found THEN l := Low; END;
(* form shift *)
x:=H[en,en];
IF (l = en) THEN
(* one root found *)
H[en,en] := x + t; (* t : shift *)
Wr[en] := H[en,en];
Wi[en] := 0.0;
en := na;
done := TRUE;
ELSE
y:=H[na,na];
w:=H[en,na]*H[na,en];
IF (l = na) THEN
(* two roots found *)
p := 0.5*(y - x);
q := p*p + w;
zz := sqrt(ABS(q));
H[en,en] := x + t; (* t : shift *)
x := H[en,en];
H[na,na] := y + t; (* t : shift *)
IF (q < 0.0) THEN
(* complex pair *)
Wr[na] := x + p;
Wr[en] := x + p;
Wi[na] := zz;
Wi[en] := - zz;
ELSE
(* real pair *)
IF (p < 0.0) THEN zz:=p - zz; ELSE zz:=p + zz; END;
Wr[na] := x + zz;
Wr[en] := Wr[na];
IF (zz <> 0.0) THEN
Wr[en] := x - w/zz;
END; (* IF *)
Wi[na]:=0.0;
Wi[en]:=0.0;
x := H[en,na];
s := ABS(x) + ABS(zz);
p := x/s;
q := zz/s;
r := Pythag(p,q); (* r := sqrt(p*p + q*q); *)
p := p/r;
q := q/r;
(* row modification *)
FOR j:=na TO N DO
zz := H[na,j];
H[na,j] := q*zz + p*H[en,j];
H[en,j] := q*H[en,j] - p*zz;
END; (* FOR *)
(* column modification *)
FOR i:=1 TO en DO
zz := H[i,na];
H[i,na] := q*zz + p*H[i,en];
H[i,en] := q*H[i,en] - p*zz;
END; (* FOR *)
(* accumulate transformations *)
FOR i:=Low TO High DO
zz := Z[na,i];
Z[na,i] := q*zz + p*Z[en,i];
Z[en,i] := q*Z[en,i] - p*zz;
END; (* FOR *)
END; (* IF q *)
en:=enm2;
ELSIF (itn = 0) THEN
(* set error -- all eigenvalues have not *)
(* converged after 30*n iterations *)
iErr:=en;
EXIT; (******)
ELSE
(* Wilkinsons original ad hoc shift *)
IF ((its MOD 10) = 0) AND (its <> 30) THEN
(* form exceptional shift *)
t:=t + x; (* t : shift *)
FOR i:=Low TO en DO
H[i,i]:=H[i,i] - x;
END; (* FOR *)
s := ABS(H[en,na]) + ABS(H[na,enm2]);
x := 0.75*s;
y := x;
w := - 0.4375*s*s;
END; (* IF *)
(* MATLABs new ad hoc shift *)
IF (its = 30) THEN
s := (y - x) / 2.0;
s := s*s + w;
IF (s > 0.0) THEN
s := sqrt(s);
IF (y < x) THEN
s := -s;
END;
s := x - w / ((y - x) / 2.0 + s);
FOR i:=Low TO N DO
H[i,i] := H[i,i] - s;
END;
t:=t + s; (* t : shift *)
x := 0.964; y := x; w := y;
END;
END;
its:=its + 1;
itn:=itn - 1;
(* look for two consecutive small *)
(* sub-diagonal elements *)
(* for m:=en-2 step -1 until l do *)
m:=en-1;
LOOP
m:=m - 1;
IF (m < l) THEN EXIT; END;
zz := H[m,m]; r := x - zz; s := y - zz;
p := (r*s - w) / H[m+1,m] + H[m,m+1];
q := H[m+1,m+1] - zz - r - s;
r := H[m+2,m+1];
s := ABS(p) + ABS(q) + ABS(r);
p := p/s; q := q/s; r := r/s;
IF (m = l) THEN EXIT; END;
tst1:=ABS(p)*(ABS(H[m-1,m-1]) + ABS(zz) + ABS(H[m+1,m+1]));
tst2:=tst1 + ABS(H[m,m-1])*(ABS(q) + ABS(r));
IF (tst2 = tst1) THEN EXIT; END;
(*
IF ABS(h[m,m-1])*(ABS(q) + ABS(r)) <=
MachEps*ABS(p)*(ABS(h[m-1,m-1]) + ABS(zz) + ABS(h[m+1,m+1]))
THEN
EXIT;
END;
*)
END; (* LOOP *)
mp2:=m + 2;
FOR i:=mp2 TO en DO
H[i,i-2]:=0.0;
IF (i <> mp2) THEN
H[i,i-3]:=0.0;
END; (* IF *)
END; (* FOR *)
(* double qr step involving rows l to en and *)
(* columns m to en *)
FOR k:=m TO na DO
notlas:=k <> na;
IF (k <> m) THEN
p := H[k,k-1];
q := H[k+1,k-1];
r := 0.0;
IF (notlas) THEN
r := H[k+2,k-1];
END; (* IF *)
x:=ABS(p) + ABS(q) + ABS(r);
IF (x <> 0.0) THEN
p:=p/x;
q:=q/x;
r:=r/x;
END; (* IF *)
END; (* IF *)
IF (x <> 0.0) THEN
s:=sqrt(p*p + q*q + r*r);
IF (p < 0.0) THEN s := - s; END;
IF (k = m) THEN
IF (l <> m) THEN
H[k,k-1]:= - H[k,k-1];
END;
ELSE
H[k,k-1]:= - s*x;
END; (* IF *)
p := p + s;
x := p/s;
y := q/s;
zz := r/s;
q := q/p;
r := r/p;
IF (notlas) THEN
(* row modification *)
FOR j:=k TO N DO
p := H[k,j] + q*H[k+1,j] + r*H[k+2,j];
H[k ,j]:=H[k,j ] - p*x;
H[k+1,j]:=H[k+1,j] - p*y;
H[k+2,j]:=H[k+2,j] - p*zz;
END; (* FOR *)
j:=MIN0(en,k+3);
(* column modification *)
FOR i:=1 TO j DO
p := x*H[i,k] + y*H[i,k+1] + zz*H[i,k+2];
H[i,k ]:=H[i,k ] - p;
H[i,k+1]:=H[i,k+1] - p*q;
H[i,k+2]:=H[i,k+2] - p*r;
END; (* FOR *)
(* accumulate transformations *)
FOR i:=Low TO High DO
p := x*Z[k,i] + y*Z[k+1,i] + zz*Z[k+2,i];
Z[k ,i]:=Z[k ,i] - p;
Z[k+1,i]:=Z[k+1,i] - p*q;
Z[k+2,i]:=Z[k+2,i] - p*r;
END; (* FOR *)
ELSE
(* row modification *)
FOR j:=k TO N DO
p := H[k,j] + q*H[k+1,j];
H[k ,j]:=H[k,j ] - p*x;
H[k+1,j]:=H[k+1,j] - p*y;
END; (* FOR *)
j:=MIN0(en,k+3);
(* column modification *)
FOR i:=1 TO j DO
p := x*H[i,k] + y*H[i,k+1];
H[i,k ]:=H[i,k ] - p;
H[i,k+1]:=H[i,k+1] - p*q;
END; (* FOR *)
(* accumulate transformations *)
FOR i:=Low TO High DO
p := x*Z[k,i] + y*Z[k+1,i];
Z[k ,i]:=Z[k ,i] - p;
Z[k+1,i]:=Z[k+1,i] - p*q;
END; (* FOR *)
END; (* IF *)
END; (* IF *)
END; (* FOR k *)
done:=FALSE;
END; (* IF (l = na) *)
END; (* IF (l = en) *)
END; (* WHILE NOT done *)
END; (* LOOP *)
END HQR2;
(*================= NUMAL Routinen =========================================*)
PROCEDURE TfmReaHes(VAR a : MATRIX;
n : INTEGER;
VAR norm : LONGREAL;
VAR int : INTVEKTOR);
PROCEDURE matvec( l,u : INTEGER;
i : INTEGER;
VAR a : MATRIX;
VAR b : ARRAY OF LONGREAL) : LONGREAL;
VAR k : INTEGER;
s : LONGREAL;
BEGIN
s:=0.0;
FOR k:=l TO u DO s:=s + a[i,k]*b[k]; END;
RETURN s;
END matvec;
VAR i,j,j1,k,l : INTEGER;
s,t,machtol : LONGREAL;
b : ARRAY [0..MaxDim-1] OF LONGREAL;
BEGIN
norm:=0.0;
FOR i:=1 TO n DO
s:= 0.0;
FOR j:=1 TO n DO s:=s + ABS(a[i,j]); END;
IF s > norm THEN norm:= s; END;
END;
machtol:= norm*MachEps; int[1]:= 0;
FOR j:=2 TO n DO
j1:=j-1; l:=0; s:= machtol;
FOR k:=j+1 TO n DO
t:= ABS(a[k,j1]);
IF t > s THEN l:= k; s:= t; END;
END;
IF l <> 0 THEN
IF ABS(a[j,j1]) < s THEN
ichrow(1,n,j,l,a);
ichcol(1,n,j,l,a)
ELSE
l := j;
END;
t:= a[j,j1];
FOR k:=j+1 TO n DO
a[k,j1]:= a[k,j1] / t;
END;
ELSE
FOR k:=j+1 TO n DO a[k,j1]:=0.0; END;
END;
FOR i:=1 TO n DO
(*
b[i - 1]:= a[i,j]:= a[i,j] +
(if l = 0 then 0 else MatMat(j+1,n,i,j1,a,a))-
matvec(1, if j1 < i-2 then j1 else i-2, i,a,b);
*)
IF (l # 0) THEN
a[i,j]:=a[i,j] + MatMat(j+1,n,i,j1,a,a);
END;
IF (j1 < i-2) THEN
a[i,j]:=a[i,j] - matvec(1,j1,i,a,b);
ELSE
a[i,j]:=a[i,j] - matvec(1,i-2,i,a,b);
END;
b[i-1] := a[i,j];
END;
int[j] := l;
END;
END TfmReaHes;
PROCEDURE BakReaHes1(VAR A : MATRIX;
N : CARDINAL;
VAR Int : CARDVEKTOR;
VAR V : VEKTOR);
VAR i,k,l : CARDINAL;
s,w : LONGREAL;
x : VEKTOR;
BEGIN
FOR i:=2 TO N DO x[i-1]:= V[i]; END;
FOR i:=N TO 2 BY -1 DO
(* v[i]:= v[i] + MatVec(1,i-2,i,a,x) *)
s:=0.0;
FOR k:=1 TO i-2 DO s:=s + A[i,k]*x[k]; END;
V[i]:= V[i] + s;
l := Int[i];
IF (l > i) THEN
w := V[i]; V[i] := V[l]; V[l] := w;
END;
END;
END BakReaHes1;
PROCEDURE RealVecHes(VAR a : MATRIX;
n : INTEGER;
lambda : LONGREAL;
VAR v : VEKTOR;
Norm : LONGREAL; (* em[1] *)
Toler : LONGREAL; (* em[6] = -10 *)
MaxIt : INTEGER; (* em[8] = 5 *)
VAR Iter : INTEGER; (* em[9] *)
VAR Resid : LONGREAL); (* em[7] *)
VAR i,i1,j,count,max : INTEGER;
m,r,machtol,tol : LONGREAL;
p : ARRAY [1..MaxDim] OF BOOLEAN;
BEGIN
machtol:= MachEps*Norm;
tol:= Toler*Norm;
max:= MaxIt;
a[1,1]:= a[1,1] - lambda;
FOR i:=1 TO n-1 DO (* gauss *)
i1 := i + 1; r := a[i,i]; m := a[i1,i];
IF (ABS(m) < machtol) THEN m:= machtol;END;
p[i] := ABS(m) <= ABS(r);
IF p[i] THEN
m:= m / r;
a[i1,i]:= m;
FOR j:=i1 TO n DO
(* a[i1,j]:= (if j > i1 then a[i1,j]
else a[i1,j] - lambda) - m * a[i,j] *)
IF (j > i1) THEN
a[i1,j]:= a[i1,j] - m*a[i,j];
ELSE
a[i1,j]:= a[i1,j] - m*a[i,j] - lambda;
END;
END;
ELSE
a[i,i]:= m;
m := r / m;
a[i1,i]:= m;
FOR j:= i1 TO n DO
(* r:= (if j > i1 then a[i1,j] else a[i1,j] - lambda); *)
IF (j > i1) THEN r:=a[i1,j] ELSE r:= a[i1,j] - lambda;END;
a[i1,j] := a[i,j] - m*r; a[i,j]:= r
;END;
END;
END; (* gauss *)
IF (ABS(a[n,n]) < machtol) THEN a[n,n]:= machtol; END;
FOR j:=1 TO n DO v[j]:=1.0; END;
count:=0;
REPEAT (* forward *)
FOR i:=1 TO n-1 DO
i1 := i + 1;
IF p[i] THEN
v[i1]:= v[i1] - a[i1,i]*v[i]
ELSE
r:= v[i1]; v[i1]:= v[i] - a[i1,i]*r;
v[i]:=r;
END;
END; (* forward *)
FOR i:=n TO 1 BY -1 DO (* backward *)
v[i]:= (v[i] - MatVec(i+1,n,i,a,v)) / (a[i,i] + MachEps);
END;
r:= 1.0 / sqrt(VecVec(1,n,0,v,v));
FOR j:=1 TO n DO v[j]:= v[j]*r; END;
count:=count + 1;
UNTIL (r <= tol) OR (count > max);
Resid := r;
Iter := count;
END RealVecHes;
PROCEDURE ComVecHes(VAR a : MATRIX;
n : INTEGER;
lambda : LONGREAL;
mu : LONGREAL;
VAR u,v : VEKTOR;
Norm : LONGREAL; (* em[1] *)
Toler : LONGREAL; (* em[6] = -10 *)
MaxIt : INTEGER; (* em[8] = 5 *)
VAR Iter : INTEGER; (* em[9] *)
VAR Resid : LONGREAL); (* em[7] *)
VAR i,j,i1,count,max : INTEGER;
aa,bb,d,m,r,s,w,x,y : LONGREAL;
rAii : LONGREAL;
machtol,tol : LONGREAL;
g,f : VEKTOR;
p : ARRAY [1..MaxDim] OF BOOLEAN;
BEGIN
machtol:= MachEps*Norm;
tol:= Toler*Norm;
max:= MaxIt;
FOR i:=2 TO n DO
f[i-1]:= a[i,i-1];
a[i,1]:=0.0;
END;
aa := a[1,1] - lambda; bb:= -mu;
FOR i:=1 TO n-1 DO
i1:= i + 1; m := f[i];
IF (ABS(m) < machtol) THEN m:= machtol; END;
a[i,i]:= m; d:= sqr(aa) + sqr(bb); p[i] := ABS(m) < sqrt(d);
IF p[i] THEN (* a[i,j] * factor and a[i1,j] - a[i,j] *)
r := m*aa / d; f[i] := r;
s := -m*bb / d; g[i] := s;
w:= a[i1,i]; x:= a[i,i1];
y := x*s + w*r; a[i1,i] := y;
x := x*r - w*s; a[i,i1] := x;
aa:= a[i1,i1] - lambda - x; bb:= -(mu + y);
FOR j:=i+2 TO n DO
w:= a[j,i]; x:= a[i,j];
y := x*s + w*r; a[j,i] := y;
x := x*r - w*s; a[i,j] := x;
a[j,i1]:= -y;
a[i1,j]:= a[i1,j] - x;
END;
ELSE
(* interchange a[i1,j] and a[i,j] - a[i1,j] * factor*)
r := aa / m; f[i]:= r;
s := bb / m; g[i]:= s;
w:= a[i1,i1] - lambda; aa:= a[i,i1] - r * w - s * mu;
a[i,i1]:= w; bb:= a[i1,i] - s*w + r*mu;
a[i1,i]:= -mu;
FOR j:=i+2 TO n DO
w := a[i1,j]; a[i1,j]:= a[i,j] - r*w;
a[i,j] := w;
a[j,i1]:= a[j,i] - s*w; a[j,i]:=0.0;
END;
END;
END;
p[n]:= TRUE;
d:= sqr(aa) + sqr(bb);
IF (d < sqr(machtol)) THEN
aa:= machtol; bb:=0.0; d:= sqr(machtol);
END;
a[n,n]:= d; f[n]:= aa; g[n]:= -bb;
FOR i:=1 TO n DO
u[i]:=1.0; v[i]:=0.0;
END;
count:= 0;
REPEAT (* forward *)
FOR i:=1 TO n DO
IF p[i] THEN
w:= v[i]; v[i]:= g[i] * u[i] + f[i] * w;
u[i]:= f[i]*u[i] - g[i]*w;
IF (i < n) THEN
v[i+1]:=v[i+1] - v[i];
u[i+1]:=u[i+1] - u[i];
END;
ELSE
aa:= u[i + 1]; bb:= v[i + 1];
u[i+1]:= u[i] - (f[i]*aa - g[i]*bb); u[i]:= aa;
v[i+1]:= v[i] - (g[i]*aa + f[i]*bb); v[i]:= bb;
END;
END; (* forward *)
FOR i:=n TO 1 BY -1 DO (* backward *)
IF (ABS(a[i,i]) <= Small*Small) THEN
rAii:=1.0 / Small;
ELSE
rAii:=1.0 / a[i,i];
END;
i1:= i + 1;
(* u[i]:= (u[i] - MatVec(i1,n,i,a,u) + (if p[i] then
TamVec(i1,n,i,a,v) else a[i1,i]*v[i1])) / a[i,i]; *)
IF p[i] THEN
u[i]:= (u[i] - MatVec(i1,n,i,a,u) + TamVec(i1,n,i,a,v))*rAii;
ELSE
u[i]:= (u[i] - MatVec(i1,n,i,a,u) + a[i1,i]*v[i1] )*rAii;
END;
(* v[i]:= (v[i] - MatVec(i1,n,i,a,v) - (if p[i] then
TamVec(i1,n,i,a,u) else a[i1,i]*u[i1])) / a[i,i] *)
IF p[i] THEN
v[i]:= (v[i] - MatVec(i1,n,i,a,v) - TamVec(i1,n,i,a,u))*rAii;
ELSE
v[i]:= (v[i] - MatVec(i1,n,i,a,v) - a[i1,i]*u[i1] )*rAii;
END;
END; (* backward *)
(* normalise *)
w:= 1.0 / sqrt(VecVec(1,n,0,u,u) + VecVec(1,n,0,v,v));
FOR j:=1 TO n DO
u[j]:=u[j]*w; v[j]:=v[j]*w;
END;
count:= count + 1;
UNTIL (w <= tol) OR (count > max);
Resid := w;
Iter := count;
END ComVecHes;
PROCEDURE EigenRGM(VAR A : MATRIX;
N : CARDINAL;
VAR Wr,Wi : VEKTOR;
VAR V : MATRIX;
ifBal : BOOLEAN;
norm : CARDINAL;
VAR iFehl : INTEGER);
VAR low,upp : CARDINAL;
iErr : INTEGER;
Index : POINTER TO CARDVEKTOR;
Scale : POINTER TO VEKTOR;
BEGIN
iFehl:=0;
low:=1; upp:=N;
IF TRUE THEN (* Berechnung mit ElmHess/ElmHessTrans/HQR2 *)
ALLOCATE(Index,N*TSIZE(CARDINAL));
IF (Index = NIL) THEN
iFehl := -1;
RETURN;
END;
IF ifBal THEN
ALLOCATE(Scale,N*TSIZE(LONGREAL));
IF (Scale = NIL) THEN
iFehl := -1;
RETURN;
ELSE
Balance(A,N,low,upp,Scale^);
END;
END;
ElmHess(A,N,low,upp,Index^);
ElmHessTrans(V,A,Index^,N,low,upp);
HQR2(N,low,upp,A,V,Wr,Wi,iFehl);
IF ifBal THEN
BalBak(N,low,upp,Scale^,N,V);
END;
END;
IF (norm = 1) THEN
NormOne(N,V,Wi,iErr);
ELSIF (norm = 2) THEN
NormTwo(N,V,Wi,iErr);
END;
IF ifBal THEN
DEALLOCATE(Scale,N*TSIZE(LONGREAL));
END;
DEALLOCATE(Index,N*TSIZE(CARDINAL));
END EigenRGM;
PROCEDURE EigenRGMsel(VAR A : MATRIX;
N : CARDINAL;
VAR EW : CVEKTOR;
VAR V : MATRIX;
ifBal : BOOLEAN;
norm : CARDINAL;
ord : INTEGER;
ifAll : BOOLEAN;
VAR Select : ARRAY OF BOOLEAN;
VAR iFehl : INTEGER);
CONST eps = 10.0*MachEps;
eps4 = MachEpsR4;
VAR i,k,n,m,ii : CARDINAL;
low,upp : CARDINAL;
Iter,iErr : INTEGER;
Atr : MATRIX;
Index,I : CARDVEKTOR;
cnt,IndexI : INTVEKTOR;
Wr,Wi,Scale : VEKTOR;
Norm,Resid,eps3 : LONGREAL;
Lambda,Zw : LONGCOMPLEX;
gefunden : BOOLEAN;
BEGIN
iFehl := 0;
n := N;
Resid := 0.0; (* Wg. Compilerwarnung *)
IF ifAll THEN
FOR i:=0 TO n-1 DO Select[i]:=TRUE; END;
END;
IF ifBal THEN
Balance(A,n,low,upp,Scale);
ELSE
low:=1; upp:=n;
END;
TfmReaHes(A,n,Norm,IndexI);
eps3:=Norm*MachEps;
FOR i:=1 TO n DO Index[i]:=VAL(CARDINAL,IndexI[i]); END;
CopyMat(Atr,A,n,n); (* Sichere A in Atr *)
HQR(A,EW,n,low,upp,cnt);
IF (ord = 1) THEN
(* Stelle sicher dass bei konjugiert komplexen Eigenwerten *)
(* derjenige mit positivem Imaginaerteil zuerst kommt *)
CSortVek(EW,I,n,ord);
FOR i:=2 TO n DO
IF (IM(EW[i]) # 0.0) THEN
IF (ABS(RE(EW[i]) - RE(EW[i-1])) < eps4) AND
(ABS(IM(EW[i]) + IM(EW[i-1])) < eps4)
THEN
IF (IM(EW[i]) > 0.0) THEN
Zw:=EW[i]; EW[i]:=EW[i-1]; EW[i-1]:=Zw;
END;
END;
END;
END;
END;
FOR i:=1 TO n DO
Wr[i]:=RE(EW[i]);
Wi[i]:=IM(EW[i]);
END;
i:=1; ii:=1;
WHILE (i <= n) DO
(* Diese Codefragment stammt aus der Routine "CHessInvIter" *)
Lambda := EW[i];
IF (i > 1) THEN
(* perturb eigenvalue if it is close to any previous eigenvalue *)
LOOP (* test *)
k:=i-1;
IF (k = 0) THEN EXIT; END;
REPEAT
gefunden:= (ABS(RE(EW[k]) - RE(Lambda)) < eps3) AND
(ABS(IM(EW[k]) - IM(Lambda)) < eps3);
DEC(k);
UNTIL gefunden OR (k = 0);
IF gefunden THEN
Lambda:=CMPLX(RE(Lambda) + eps3,IM(Lambda)); (* GOTO test: *)
ELSE
EXIT;
END;
END; (* LOOP *)
EW[i] := CMPLX(RE(Lambda),IM(EW[i]));;
Wr[i]:=RE(Lambda); Wi[i]:=IM(Lambda);
END; (* IF *)
IF Select[i-1] THEN
CopyMat(A,Atr,n,n);
IF (IM(EW[i]) = 0.0) THEN
RealVecHes(A,n,Wr[i],V[ii],Norm,eps,5,Iter,Resid);
BakReaHes1(Atr,n,Index,V[ii]);
INC(ii);
ELSIF (IM(EW[i]) > 0.0) THEN
(* Berechne nur einen der Eigenvektoren *)
(* einer konjugiert komplexen Wurzel *)
ComVecHes(A,n,Wr[i],Wi[i],V[ii],V[ii+1],Norm,eps,5,Iter,Resid);
BakReaHes1(Atr,n,Index,V[ii]); (* Realteil *)
BakReaHes1(Atr,n,Index,V[ii+1]); (* Imaginaerteil *)
INC(ii,2);
END;
IF (Resid > MachEpsR4) THEN iFehl:=1; END;
END; (* IF Select *)
INC(i);
END;
IF ifBal THEN
m := ii-1;
BalBak(N,low,upp,Scale,m,V);
END;
IF (norm = 1) THEN
NormOne(n,V,Wi,iErr);
ELSIF (norm = 2) THEN
NormTwo(n,V,Wi,iErr);
END;
END EigenRGMsel;
PROCEDURE ZTransQL(VAR EW : ARRAY OF LONGREAL;
VAR ND : ARRAY OF LONGREAL;
dim : CARDINAL;
u,o : CARDINAL;
VAR Z : ARRAY OF ARRAY OF LONGCOMPLEX;
genau : LONGREAL;
VAR iFehl : INTEGER);
PROCEDURE AbsMin(x,y : LONGREAL) : LONGREAL;
BEGIN
IF (ABS(x) <= ABS(y)) THEN x:=ABS(x) ELSE x:=ABS(y); END;
RETURN x;
END AbsMin;
CONST MaxIter = 30;
VAR Iter,m,l,i,k : CARDINAL;
m0 : INTEGER;
b,c,f,g,p,r,s : LONGREAL;
BEGIN
IF (o > dim) OR (u < 1) OR (u > o) THEN
Fehler := TRUE;
Fehlerflag:="Dimensionierungsfehler (ZTransQL)";
RETURN;
END;
Fehler := FALSE; iFehl:=0;
u:=u-1; o:=o-1;
FOR l:=u TO o DO
Iter:=0;
LOOP
IF (Iter <= MaxIter) THEN INC(Iter) ELSE
Fehler:=TRUE; Fehlerflag:='MaxIter Ueberschritten (TriQL).';
ErrOut(Fehlerflag);
iFehl := VAL(INTEGER,l);
RETURN;
END;
m0:=VAL(INTEGER,l) - 1;
REPEAT (* Suche kleine Subdiagonalelemente *)
INC(m0); m:=VAL(CARDINAL,m0);
(* Sollte das nich m = o heissen (war m = dim-1) ??? *)
UNTIL (m = o) OR (ABS(ND[m]) <= genau*AbsMin(EW[m],EW[m+1]));
p:=EW[l];
IF (m = l) THEN EXIT; END; (* Ausgang LOOP !!! *)
g:=(EW[l+1] - p) / (2.0*ND[l]);
r:=Pythag(g,1.0);
IF (g < 0.0) THEN
g := EW[m] - p + ND[l] / (g - r);
ELSE
g := EW[m] - p + ND[l] / (g + r);
END;
s:=1.0; c:=1.0; p:=0.0;
FOR i:=m-1 TO l BY -1 DO
f := s*ND[i];
b := c*ND[i];
IF (ABS(f) < ABS(g)) THEN
s := f / g;
r := sqrt(s*s + 1.0);
ND[i+1] := g*r;
c := 1.0 / r;
s := s*c;
ELSE
c := g/f;
r := sqrt(c*c+1.0);
ND[i+1] := f*r;
s := 1.0 / r;
c := c*s;
END; (* IF *)
g := EW[i+1] - p;
r := s*(EW[i] - g) + 2.0*c*b;
p := s*r;
EW[i+1] := g + p;
g := c*r - b;
FOR k:=0 TO dim-1 DO (* Imaginaerteil wird nicht beruehrt *)
f := RE(Z[i+1,k]);
Z[i+1,k] := CMPLX( s*RE(Z[i,k]) + c*f, 0.0);
Z[i ,k] := CMPLX( c*RE(Z[i,k]) - s*f, 0.0);
END;
END; (* FOR i *)
EW[l]:=EW[l]-p;
ND[l]:=g;
ND[m]:=0.0;
END; (* LOOP *)
END; (* FOR l *)
END ZTransQL;
PROCEDURE CHTriDi(VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
n : INTEGER;
VAR Tau : ARRAY OF LONGCOMPLEX;
VAR D,E : ARRAY OF LONGREAL;
VAR E2 : ARRAY OF LONGREAL);
VAR i,j,k,l,jp1 : INTEGER;
fr,gr,h,hh,si,scale : LONGREAL;
f,g,ctmp : LONGCOMPLEX;
BEGIN
Tau[n-1] := one;
si := 0.0;
FOR i:=0 TO n-1 DO
D[i] := RE(A[i,i]);
END;
FOR i:=n-1 TO 0 BY -1 DO
l := i - 1;
h := 0.0;
scale := 0.0;
IF (l < 0) THEN
E[i] := 0.0;
E2[i] := 0.0;
ELSE
FOR k:=0 TO l DO
(* scale row (algol tol then not needed) *)
scale:=scale + ABS(RE(A[i,k])) + ABS(IM(A[i,k]));
END;
END;
IF (scale = 0.0) THEN
IF (l >= 0) THEN
Tau[l] := one;
END;
ELSE
FOR k:=0 TO l DO
A[i,k] := scalarMult((1.0 / scale), A[i,k]);
h:=h + sqr(RE(A[i,k])) + sqr(IM(A[i,k]));
END;
E2[i] := scale*scale*h;
gr := sqrt(h);
E[i] := scale*gr;
fr := CABS(A[i,l]);
(* form next diagonal element of matrix t *)
IF (fr # 0.0) THEN
ctmp := A[i,l]*Tau[i];
Tau[l] := CMPLX(-RE(ctmp) / fr, IM(Tau[l]));
si := IM(ctmp) / fr;
h := h + fr*gr;
gr := 1.0 + gr / fr;
A[i,l] := scalarMult(gr,A[i,l]);
END;
IF (l # 0) OR (fr = 0.0) THEN
IF (fr = 0.0) THEN
Tau[l] := CMPLX(-RE(Tau[i]),IM(Tau[l]));
si := IM(Tau[i]);
A[i,l] := CMPLX(gr,IM(A[i,l]));
END;
fr := 0.0;
FOR j:=0 TO l DO
g := zero;
(* form element of a*u *)
FOR k:=0 TO j DO
g:=g + A[j,k]*conj(A[i,k]);
END;
jp1 := j + 1;
IF (l >= jp1) THEN
FOR k:=jp1 TO l DO
g:=g + conj(A[k,j]*A[i,k]);
END;
END;
(* form element of p *)
E[j] := RE(g) / h;
Tau[j] := CMPLX(RE(Tau[j]), IM(g) / h);
fr:=fr + E[j]*RE(A[i,j]) - IM(Tau[j])*IM(A[i,j]);
END;
hh := fr / (h + h);
(* form reduced a *)
FOR j:=0 TO l DO
f := conj(A[i,j]);
g := CMPLX(E[j] - hh*RE(f), IM(Tau[j]) - hh*IM(f));
E[j] := RE(g);
Tau[j] := CMPLX(RE(Tau[j]), -IM(g));
FOR k:=0 TO j DO
A[j,k]:=A[j,k] - f*CMPLX(E[k],IM(Tau[k])) - g*A[i,k];
END;
END;
END;
FOR k:=0 TO l DO
A[i,k] := scalarMult(scale,A[i,k]);
END;
Tau[l] := CMPLX(RE(Tau[l]), -si);
END;
hh := D[i];
D[i] := RE(A[i,i]);
A[i,i] := CMPLX(hh,scale*sqrt(h));
END; (* FOR i *)
END CHTriDi;
PROCEDURE CHTriBak(VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR Tau : ARRAY OF LONGCOMPLEX;
N : INTEGER;
M : INTEGER;
VAR Z : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR iFehl : INTEGER);
VAR i,j,k,l : INTEGER;
h : LONGREAL;
s : LONGCOMPLEX;
BEGIN
IF (M = 0) THEN iFehl:=1; RETURN; END;
iFehl := 0;
(* transform the eigenvectors of the real symmetric tridiagonal *)
(* matrix to those of the hermitian tridiagonal matrix. *)
FOR j:=0 TO M-1 DO
FOR k:=0 TO N-1 DO
Z[j,k] := conj(scalarMult(RE(Z[j,k]),Tau[k]));
END;
END;
IF (N > 1) THEN
(* recover and apply the householder matrices *)
FOR i:=1 TO N-1 DO
l := i - 1;
h := IM(A[i,i]);
IF (h # 0.0) THEN
FOR j:=0 TO M-1 DO
s:=zero;
FOR k:=0 TO l DO s:=s + A[i,k]*Z[j,k]; END;
(* double divisions avoid possible underflow *)
s := scalarMult((1.0 / h),s);
s := scalarMult((1.0 / h),s);
FOR k:=0 TO l DO Z[j,k]:=Z[j,k] - s*conj(A[i,k]); END;
END;
END; (* IF h *)
END; (* FOR i *)
END; (* IF *)
END CHTriBak;
PROCEDURE EigenCHM2(VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
dim : CARDINAL;
m : CARDINAL;
VAR EW : ARRAY OF LONGREAL;
VAR EV : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR iFehl : INTEGER);
VAR i,j : CARDINAL;
N,M,iErr : INTEGER; (* iErr fuer CHTriBak *)
D0,E,E2 : POINTER TO ARRAY [0..MAX(INTEGER)-1] OF LONGREAL;
Tau : POINTER TO ARRAY [0..MAX(INTEGER)-1] OF LONGCOMPLEX;
BEGIN
Fehler:=FALSE;
ALLOCATE(E ,dim*TSIZE(LONGREAL));
ALLOCATE(E2 ,dim*TSIZE(LONGREAL));
ALLOCATE(Tau,dim*TSIZE(LONGCOMPLEX));
IF (E = NIL) OR (E2 = NIL) (*OR (Tau = NIL)*) THEN
iFehl:=-2;
Fehler:=TRUE;
Fehlerflag:="Kein Freispeicher vorhanden (EigenLib2.EigenCHM)";
ErrOut(Fehlerflag);
RETURN;
END;
ALLOCATE(D0,dim*TSIZE(LONGREAL));
IF (D0 = NIL) THEN
iFehl:=-1;
ELSE
(* Sichere die Haupdiagonale *)
FOR i:=0 TO dim-1 DO D0^[i]:=RE(A[i,i]); END;
END;
N:=VAL(INTEGER,dim);
m:=dim;
M:=VAL(INTEGER,m);
Fehler := FALSE;
iFehl := 0;
(* Auf EW liegt temporaer die Diagonle der Trilinearmatrix *)
CHTriDi(A,N,Tau^,EW,E^,E2^);
(* Eventuell CHTriDi oder ZTransQL anpassen *)
(* um diesen Schritt eleminieren zu koennen *)
FOR i:=1 TO dim-1 DO E^[i-1] := E^[i]; END; E^[dim-1]:=0.0;
(* Eventuell noch Code mit Bisect einfuegen wenn nur wenige *)
(* Eigenvektoren benoetigt werden *)
FOR i:=0 TO dim-1 DO
FOR j:=0 TO dim-1 DO EV[i,j]:=zero END;
EV[i,i]:=one;
END;
ZTransQL(EW,E^,N,1,N,EV,MachEps,iErr);
IF Fehler OR (iErr # 0) THEN
ErrOut(Fehlerflag);
(* Noch nicht gut da iFehl = -1 so nicht signalisiert wird *)
(* Alledings steht Fehler in diesem Fall schon auf "wahr" *)
iFehl:=iErr;
END;
CHTriBak(A,Tau^,N,M,EV,iErr);
IF (D0 # NIL) THEN
(* Stelle die Matrix A wieder her *)
A[0,0]:=CMPLX(D0^[0],0.0);
FOR i:=1 TO dim-1 DO
FOR j:=0 TO i-1 DO
A[i,j]:=A[j,i];
END;
A[i,i]:=CMPLX(D0^[i],0.0); (* Diagonale reel !!! *)
END;
DEALLOCATE(D0,dim*TSIZE(LONGREAL));
END;
DEALLOCATE(Tau,dim*TSIZE(LONGCOMPLEX));
DEALLOCATE(E2 ,dim*TSIZE(LONGREAL));
DEALLOCATE(E ,dim*TSIZE(LONGREAL));
END EigenCHM2;
PROCEDURE CJacobi1D( N : INTEGER;
VAR A : ARRAY OF LONGCOMPLEX; (* SUPERVEKTOR *)
VAR EW : ARRAY OF LONGREAL;
VAR EV : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR iFehl : INTEGER);
VAR i,j,k,l,nrot,Iter,MaxIt : INTEGER;
ii,ij,kk,ll,kl,ik,ki,il,li : INTEGER;
norm,tresh : LONGREAL;
delta,s,s0,t0,t1,w0,Are : LONGREAL;
c0,c1,Aki,Aik,EVki : LONGCOMPLEX;
rot : BOOLEAN;
Index : POINTER TO ARRAY [0..MAX(INTEGER)-1]
OF INTEGER;
BEGIN
IF (N < 1) THEN iFehl:=1; RETURN; END;
IF (N = 1) THEN
EW[0]:=RE(A[0]); EV[0,0]:=one;
RETURN;
END;
ALLOCATE(Index,(N+1)*TSIZE(INTEGER));
ii:=0; (* Index haelt die Dreieckszahlen *)
FOR i:=0 TO N DO (* Wichtig BIS N !!! *)
Index^[i]:=ii; INC(ii,i+2);
END;
FOR i:=0 TO N-1 DO (* EV ist Einheitsmatrix zum Start *)
FOR j:=0 TO N-1 DO EV[i,j]:=zero; END;
EV[i,i]:=one;
END;
ij:=0;
s := ABS(RE(A[0])) + ABS(IM(A[0]));
FOR i:=1 TO N-1 DO (* calculate norm of matrix A *)
FOR j:=0 TO i-1 DO
INC(ij);
s:=s + 2.0*ABS(RE(A[ij])) + ABS(IM(A[ij]));
END;
INC(ij);
s:=s + ABS(RE(A[ij])) + ABS(IM(A[ij]));
END;
norm:=s;
MaxIt := N*(N DIV 2); IF (MaxIt < 60) THEN MaxIt:=60; END;
iFehl:=0; Iter:=0; nrot:=0;
LOOP
INC(Iter);
IF (Iter > MaxIt) THEN
iFehl:=-Iter;
EXIT;
END;
s:=0.0; ij:=0;
FOR i:=1 TO N-1 DO (* calculate norm of off-diagonal elements *)
FOR j:=0 TO i-1 DO
INC(ij);
s:=s + ABS(RE(A[ij])) + ABS(IM(A[ij]));
END;
INC(ij); (* Ueberspringe Diagonalelement *)
END;
s:=s / VAL(LONGREAL,(N*N));
IF (s < Small*norm) THEN EXIT; END;
(*********************************)
tresh:= 0.2*s;
rot:=FALSE;
kk:=0;
FOR k:=0 TO N-1 DO
l := k + 1;
ll := ((l*l + 3*l) DIV 2);
kl := Index^[k+1] - 1;
FOR l:=k+1 TO N-1 DO
IF (CABS(A[kl]) > tresh) THEN
rot:=TRUE;
INC(nrot);
delta:=sqr(RE(A[ll]) - RE(A[kk])) + 4.0*sqr(CABS(A[kl]));
t0 := RE(A[ll]) - RE(A[kk]) + sqrt(delta);
t1 := RE(A[ll]) - RE(A[kk]) - sqrt(delta);
IF (ABS(t0) >= ABS(t1)) THEN w0:=t0 ELSE w0:=t1;END;
s0 := ABS(w0) / sqrt(sqr(w0) + 4.0*sqr(CABS(A[kl])));
t0 := 2.0*s0 / w0;
c0 := scalarMult(t0,A[kl]);
c1 := conj(c0);
ik := Index^[k] - k;
il := Index^[l] - l;
FOR i:=0 TO k-1 DO
Aik := A[ik];
A[ik] := c0*Aik + scalarMult(s0,A[il]);
A[il] := c1*A[il] - scalarMult(s0,Aik);
INC(ik); INC(il);
END;
ki := Index^[k] + k + 1;
il := Index^[l-1] + k + 2;
FOR i:=k+1 TO l-1 DO
Aki := A[ki];
A[ki] := c1*Aki + scalarMult(s0,conj(A[il]));
A[il] := c1*A[il] - scalarMult(s0,conj(Aki));
INC(ki,i+1); INC(il);
END;
ki := Index^[l] + k + 1;
li := Index^[l+1] - 1;
FOR i:=l+1 TO N-1 DO
Aki := A[ki];
A[ki] := c1*Aki + scalarMult(s0,A[li]);
A[li] := c0*A[li] - scalarMult(s0,Aki);
INC(ki,i+1); INC(li,i+1);
END;
t0 := RE(A[kk]);
t1 := 4.0*sqr(s0*CABS(A[kl])) / w0;
Are := sqr(CABS(c0))*t0 + t1 + sqr(s0)*RE(A[ll]);
A[kk] := CMPLX(Are,IM(A[kk]));
Are := sqr(s0)*t0 - t1 + sqr(CABS(c0))*RE(A[ll]);
A[ll] := CMPLX(Are,IM(A[ll]));
A[kl] := zero;
FOR i:=0 TO N-1 DO (* Berechne Eigenvektoren *)
EVki := EV[k,i]; (* c1,c0 vertauscht gegen 2D routine *)
EV[k,i] := c1*EVki + scalarMult(s0,EV[l,i]);
EV[l,i] := c0*EV[l,i] - scalarMult(s0,EVki);
END;
END; (* IF *)
INC(ll,l+2);
INC(kl,l+1);
END; (* FOR l *)
INC(kk,k+2);
END; (* FOR k *)
IF (NOT rot) THEN EXIT; END;
END; (* LOOP *)
ii:=0;
FOR i:=0 TO N-1 DO EW[i]:=RE(A[ii]); INC(ii,i+2); END;
DEALLOCATE(Index,(N+1)*TSIZE(INTEGER));
IF debug THEN
TIO.WrStr(" MaxIt = "); TIO.WrInt(MaxIt,3); TIO.WrLn;
TIO.WrStr(" nrot = "); TIO.WrInt(nrot ,3); TIO.WrLn;
TIO.WrStr(" Iter = "); TIO.WrInt(Iter ,3); TIO.WrLn;
TIO.WrStr(" s(ij) = "); TIO.WrLngReal(s, 9,-3); TIO.WrLn;
TIO.WrLn;
END;
END CJacobi1D;
PROCEDURE CJacobi2D( N : INTEGER;
VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR EW : ARRAY OF LONGREAL;
VAR EV : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR iFehl : INTEGER);
VAR i,j,k,l,nrot,Iter,MaxIt : INTEGER;
tresh,norm : LONGREAL;
delta,s,s0,t0,t1,w0,Are : LONGREAL;
c0,c1,Aki,Aik,EVki : LONGCOMPLEX;
rot : BOOLEAN;
BEGIN
IF (N < 1) THEN iFehl:=1; RETURN; END;
IF (N = 1) THEN
EW[0]:=RE(A[0,0]); EV[0,0]:=one;
RETURN;
END;
FOR i:=0 TO N-1 DO
FOR j:=0 TO N-1 DO EV[i,j]:=zero; END;
EV[i,i]:=one;
END;
s:=0.0;
FOR i:=0 TO N-1 DO
FOR j:=0 TO N-1 DO
s:=s + ABS(RE(A[i,j])) + ABS(IM(A[i,j]));
END;
END;
norm:=s;
MaxIt := N*(N DIV 2); IF (MaxIt < 60) THEN MaxIt:=60; END;
iFehl:=0; Iter:=0; nrot:=0;
LOOP
INC(Iter);
IF (Iter > MaxIt) THEN
iFehl:=-MaxIt;
EXIT;
END;
s:=0.0;
FOR i:=0 TO N-2 DO (* calculate norm of off-diagonal elements *)
FOR j:=i+1 TO N-1 DO s:=s + ABS(RE(A[i,j])) + ABS(IM(A[i,j])); END;
END;
s:=s / VAL(LONGREAL,(N*N));
IF (s < Small*norm) THEN EXIT; END;
(*********************************)
tresh:= 0.2*s;
rot:=FALSE;
FOR k:=0 TO N-1 DO
FOR l:=k+1 TO N-1 DO
IF (CABS(A[k,l]) > tresh) THEN
rot:=TRUE;
INC(nrot);
delta:=sqr(RE(A[l,l]) - RE(A[k,k])) + 4.0*sqr(CABS(A[k,l]));
t0 := RE(A[l,l]) - RE(A[k,k]) + sqrt(delta);
t1 := RE(A[l,l]) - RE(A[k,k]) - sqrt(delta);
IF (ABS(t0) >= ABS(t1)) THEN w0:=t0 ELSE w0:=t1;END;
s0 := ABS(w0) / sqrt(sqr(w0) + 4.0*sqr(CABS(A[k,l])));
t0 := 2.0*s0 / w0;
c0 := scalarMult(t0,A[k,l]);
c1 := conj(c0);
FOR i:=0 TO k-1 DO
Aik := A[i,k];
A[i,k] := c0*Aik + scalarMult(s0,A[i,l]);
A[i,l] := c1*A[i,l] - scalarMult(s0,Aik);
END;
FOR i:=k+1 TO l-1 DO
Aki := A[k,i];
A[k,i] := c1*Aki + scalarMult(s0,conj(A[i,l]));
A[i,l] := c1*A[i,l] - scalarMult(s0,conj(Aki));
END;
FOR i:=l+1 TO N-1 DO
Aki := A[k,i];
A[k,i] := c1*Aki + scalarMult(s0,A[l,i]);
A[l,i] := c0*A[l,i] - scalarMult(s0,Aki);
END;
t0 := RE(A[k,k]);
t1 := 4.0*sqr(s0*CABS(A[k,l])) / w0;
Are := sqr(CABS(c0))*t0 + t1 + sqr(s0)*RE(A[l,l]);
A[k,k] := CMPLX(Are,IM(A[k,k]));
Are := sqr(s0)*t0 - t1 + sqr(CABS(c0))*RE(A[l,l]);
A[l,l] := CMPLX(Are,IM(A[l,l]));
A[k,l] := zero;
FOR i:=0 TO N-1 DO (* Berechne Eigenvektoren *)
EVki := EV[k,i];
EV[k,i] := c0*EVki + scalarMult(s0,EV[l,i]);
EV[l,i] := c1*EV[l,i] - scalarMult(s0,EVki);
END;
END; (* IF *)
END; (* FOR k *)
END; (* FOR l *)
IF (NOT rot) THEN EXIT; END;
END; (* LOOP *)
FOR i:=0 TO N-1 DO EW[i]:=RE(A[i,i]); END;
IF debug THEN
TIO.WrStr(" MaxIt = "); TIO.WrInt(MaxIt,3); TIO.WrLn;
TIO.WrStr(" nrot = "); TIO.WrInt(nrot ,3); TIO.WrLn;
TIO.WrStr(" Iter = "); TIO.WrInt(Iter ,3); TIO.WrLn;
TIO.WrStr(" s(ij) = "); TIO.WrLngReal(s, 9,-3); TIO.WrLn;
TIO.WrLn;
END;
END CJacobi2D;
END EigenLib2.