IMPLEMENTATION MODULE EigenLibAux;
(*------------------------------------------------------------------------*)
(* Stellt Hilfsroutinen fuer die Eigenwertberechung zur Verfuegung. *)
(* Auxillary routines for the calculation of eigen systems. *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung: *)
(* *)
(* 20.09.17, MRi: Zusammenfassen aller Routinen in neuer EigenLibAux *)
(* 02.10.17, MRi: CmplxComRes, IstHermitisch eingefuegt *)
(* 01.06.18, MRi: CheckOrtho eingefuegt *)
(*------------------------------------------------------------------------*)
(* Offene Punkte *)
(* *)
(* - Alles auf offene Felder umstellen *)
(* - Komplexe Arithmetik konsequent einfueren *)
(*------------------------------------------------------------------------*)
(* Hint to sources of some routines *)
(* *)
(* - Some auxilliar routine dealing with matrices as provided *)
(* by HQR2 (e.g. NormOne) originate form the PMATH Pascal Library *)
(* PMATH is provided by Jean Debord, l'Universite de Limoges, France *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: EigenLibAux.mod,v 1.1 2017/09/27 08:13:54 mriedl Exp mriedl $ *)
FROM Deklera IMPORT VEKTOR,MATRIX;
FROM LongMath IMPORT sqrt;
FROM LMathLib IMPORT MachEps;
FROM LongComplexMath IMPORT zero,abs,conj,scalarMult;
FROM CmplxMath IMPORT CABS,CMUL,CDIV;
FROM LibDBlas IMPORT ddot;
CONST MinGenau = 10.0*MachEps;
CONST debug = FALSE;
PROCEDURE CSortEwEv(VAR Wr,Wi : VEKTOR;
VAR Vr,Vi : MATRIX;
eps : LONGREAL;
dim : INTEGER;
ord : INTEGER);
(*----------------------------------------------------------------*)
(* Sortiert Eigenwerte und Eigenvektoren nach Groesse der Eigen- *)
(* werte. Die komplexen Eigenwerde werden dabei in Wr,Wi *)
(* uebergeben, die komplexen Eigenvektoren in Vr,Vi. *)
(* *)
(* Wenn ord > 0 in absteigender Reihenfolge *)
(* Wenn ord < 0 in ansteigender Reihenfolge *)
(*----------------------------------------------------------------*)
VAR i,l,iminmax : INTEGER;
minmax,tmp : LONGREAL;
swap : BOOLEAN;
BEGIN
FOR l:=1 TO dim DO
minmax:=Wr[l];
iminmax:=l;
FOR i:=l+1 TO dim DO
IF (ord > 0) THEN (* Berechne Maximum von W[n] *)
IF (Wr[i] > minmax) THEN minmax:=Wr[i]; iminmax:=i; END;
ELSIF (ord < 0) THEN (* Berechne Minimum von W[n] *)
IF (Wr[i] < minmax) THEN minmax:=Wr[i]; iminmax:=i; END;
END; (* IF *)
END; (* FOR i *)
IF (iminmax <> l) THEN (* Vertausche W[i] mit W[imax] *)
tmp:=Wr[l]; Wr[l]:=Wr[iminmax]; Wr[iminmax]:=tmp;
tmp:=Wi[l]; Wi[l]:=Wi[iminmax]; Wi[iminmax]:=tmp;
FOR i:=1 TO dim DO (* Vertausche Eigenvektoren entspr. *)
tmp:=Vr[l,i]; Vr[l,i]:=Vr[iminmax,i]; Vr[iminmax,i]:=tmp;
tmp:=Vi[l,i]; Vi[l,i]:=Vi[iminmax,i]; Vi[iminmax,i]:=tmp;
END; (* FOR i *)
END; (* IF *)
END; (* FOR l *)
IF (eps > 0.0) THEN
(* Versuche bei konjugiert komplexen Eigenweten *)
(* den imaginaeren Part zu beruecksichtigen *)
swap:=FALSE;
l:=1;
WHILE (l <= dim) DO
IF (ABS(Wi[l]) > MachEps*ABS(Wr[l])) AND
(ABS(Wr[l] - Wr[l+1]) < eps*ABS(Wr[l]))
THEN
IF (ord > 0) THEN (* Berechne Maximum von W[n] *)
swap := (Wi[l] < Wi[l+1]);
ELSIF (ord < 0) THEN (* Berechne Minimum von W[n] *)
swap := (Wi[l] > Wi[l+1]);
END; (* IF *)
IF swap THEN
tmp:=Wr[l]; Wr[l]:=Wr[l+1]; Wr[l+1]:=tmp;
tmp:=Wi[l]; Wi[l]:=Wi[l+1]; Wi[l+1]:=tmp;
FOR i:=1 TO dim DO (* Vertausche Eigenvektoren entspr. *)
tmp:=Vr[l,i]; Vr[l,i]:=Vr[l+1,i]; Vr[l+1,i]:=tmp;
tmp:=Vi[l,i]; Vi[l,i]:=Vi[l+1,i]; Vi[l+1,i]:=tmp;
END; (* FOR i *)
END; (* IF *)
l:=l + 2;
ELSE
l:=l + 1;
END; (* IF ABS(Wi[l]) *)
END; (* WHILE l *)
END; (* IF eps *)
END CSortEwEv;
PROCEDURE NormOne( N : INTEGER;
VAR V : ARRAY OF ARRAY OF LONGREAL;
Wi : ARRAY OF LONGREAL;
VAR iErr : INTEGER);
(*----------------------------------------------------------------*)
(* Normalize eigenvectors to have norm 1 *)
(* *)
(* N : Dimension of matrix *)
(* <==> V : Matrix with eigenvektors *)
(* ==> Wi : Imaginary parts of eigenvalues *)
(* <== iErr : return code *)
(*----------------------------------------------------------------*)
VAR i,j : INTEGER;
maxi,tr,ti : LONGREAL;
BEGIN
IF (N < 1) THEN
iErr := 1;
ELSE
iErr := 0;
j:=0;
WHILE (j < N) DO
IF (Wi[j] = 0.0) THEN
maxi := V[j,0];
FOR i:=1 TO N-1 DO
IF (ABS(V[j,i]) > ABS(maxi)) THEN maxi := V[j,i]; END;
END;
IF (maxi <> 0.0) THEN
maxi := 1.0 / maxi;
FOR i:=0 TO N-1 DO V[j,i]:=V[j,i]*maxi; END;
END;
INC(j);
ELSE
tr := V[j ,0];
ti := V[j+1,0];
FOR i:=1 TO N-1 DO
IF (CABS(V[j,i],V[j+1,i]) > CABS(tr,ti)) THEN
tr := V[j ,i];
ti := V[j+1,i];
END;
END;
IF (tr <> 0.0) OR (ti <> 0.0) THEN
FOR i:=0 TO N-1 DO
CDIV(V[j,i],V[j+1,i], tr,ti, V[j,i],V[j+1,i]);
END;
END;
INC(j,2);
END;
END; (* j *)
END; (* IF N *)
END NormOne;
PROCEDURE NormTwo( N : INTEGER;
VAR V : MATRIX;
VAR Wi : VEKTOR;
VAR iErr : INTEGER);
(*----------------------------------------------------------------*)
(* Normalize eigenvectors to have one norm 2 *)
(* *)
(* N : Dimension of matrix *)
(* <==> V : Matrix with eigenvektors *)
(* ==> Wi : Imaginary parts of eigenvalues *)
(* <== iErr : return code *)
(*----------------------------------------------------------------*)
VAR i,j : INTEGER;
sr,normr : LONGREAL;
ctmpre,ctmpim : LONGREAL;
BEGIN
iErr := 0;
IF (N < 1) THEN
iErr := -1;
ELSIF (N < 2) THEN
iErr := 1;
ELSE
j:=1;
WHILE (j <= N) DO (* for j:=1 to N do *)
IF (Wi[j] = 0.0) THEN
sr:=0.0;
FOR i:=1 TO N DO
sr:=sr + V[j,i]*V[j,i];
END;
IF (sr <> 0.0) THEN
normr := 1.0 / sqrt(sr);
FOR i:=1 TO N DO V[j,i]:=V[j,i]*normr; END;
END;
j:=j + 1;
ELSE
sr:=0.0;
FOR i:=1 TO N DO
CMUL(V[j,i],V[j+1,i], V[j,i],-V[j+1,i], ctmpre,ctmpim);
sr:=sr + CABS(ctmpre,ctmpim);
END;
IF (sr <> 0.0) THEN
normr := 1.0 / sqrt(sr);
FOR i:=1 TO N DO
V[j ,i] := V[j ,i]*normr;
V[j+1,i] := V[j+1,i]*normr;
END;
END;
j:=j + 2;
END;
END; (* j *)
END; (* IF N *)
END NormTwo;
PROCEDURE ComRes( N,ii : INTEGER;
VAR Wr,Wi : VEKTOR;
VAR Vr,Vi : MATRIX;
VAR Ar,Ai : MATRIX;
VAR norm : LONGREAL);
(*---------------------------------------------------------------*)
(* Calculates the residual for eigenvector/value ii *)
(* *)
(* N : Dimension of the problem *)
(* ii : Index of Eigenvalue/vektor pair to be checked *)
(* (Wr,Wi) : Complex eigenvalues of the system *)
(* (Wr,Wi) : Complex eigenvalues of the system *)
(* (Vr,Vi) : Complex eigenvector of the system *)
(* (Ar,Ai) : Complex matrix under investigation *)
(*---------------------------------------------------------------*)
VAR j,k : INTEGER;
g,s,ss : LONGREAL;
BEGIN
ss:=0.0;
FOR j:=1 TO N DO (* (s,g) = W[i]*(T,Vi) - (A,Ai)*(T,Vi) *)
s := - Wr[ii]*Vr[ii,j] + Wi[ii]*Vi[ii,j];
g := - Wi[ii]*Vr[ii,j] - Wr[ii]*Vi[ii,j];
FOR k:=1 TO N DO
s:=s + Ar[j,k]*Vr[ii,k] - Ai[j,k]*Vi[ii,k];
g:=g + Ar[j,k]*Vi[ii,k] + Ai[j,k]*Vr[ii,k];
END;
ss:=ss + s*s + g*g;
END;
norm:=sqrt(ss) / LFLOAT(N);
END ComRes;
PROCEDURE CheckRGM( N : INTEGER;
VAR V : MATRIX;
VAR A : MATRIX;
VAR Wr,Wi : VEKTOR;
VAR Norm : LONGREAL;
VAR NormM : LONGREAL);
(*----------------------------------------------------------------*)
(* Calulates the norm of the residial of the eigensysem of A *)
(* *)
(* N : Dimension of matrix *)
(* <==> V : Matrix with eigenvektors as provided by HQR2 *)
(* ==> A : Original Matrix for which the eigensystem had *)
(* been calculated *)
(* ==> Wr : Real parts of eigenvalues *)
(* ==> Wi : Imaginary parts of eigenvalues *)
(* <== Norm : Caluclated weighted norm of the matrix *)
(* Z = (A x V) - (W x V) *)
(*----------------------------------------------------------------*)
VAR i,j,k : INTEGER;
sr,si,ss,max,tmp : LONGREAL;
BEGIN
ss:=0.0; max:=0.0; k:=1;
WHILE (k <= N) DO (* FOR k:=1 TO N DO *)
IF (Wi[k] = 0.0) THEN
FOR i:=1 TO N DO
sr := - Wr[k]*V[k,i];
FOR j:=1 TO N DO
sr:=sr + A[i,j]*V[k,j];
END;
tmp := sr*sr;
IF (tmp > max) THEN max:=tmp; END;
ss:=ss + tmp;
END;
INC(k);
ELSE
FOR i:=1 TO N DO
sr := - Wr[k]*V[k,i] + Wi[k]*V[k+1,i];
si := - Wi[k]*V[k,i] - Wr[k]*V[k+1,i];
FOR j:=1 TO N DO
sr:=sr + A[i,j]*V[k ,j];
si:=si + A[i,j]*V[k+1,j];
END;
tmp := 2.0*(sr*sr + si*si);
IF (tmp > max) THEN max:=tmp; END;
ss:=ss + tmp;
END;
INC(k,2);
END;
END; (* WHILE *)
Norm := sqrt(ss) / LFLOAT(N);
NormM := sqrt(max);
END CheckRGM;
PROCEDURE KonvCmplxC1( N : CARDINAL;
VAR Zreal : ARRAY OF ARRAY OF LONGREAL;
VAR Zcmplx : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR Wi : ARRAY OF LONGREAL);
VAR j,k : CARDINAL;
BEGIN
FOR k:=0 TO N-1 DO
IF (Wi[k] = 0.0) THEN
FOR j:=0 TO N-1 DO
Zcmplx[k,j]:=CMPLX(Zreal[k,j],0.0);
END;
ELSIF (Wi[k] < 0.0) THEN (* konjugiert komplexe Wurzel, k > 1 *)
FOR j:=0 TO N-1 DO
Zcmplx[k,j]:=conj(Zcmplx[k-1,j]);
END;
ELSE
FOR j:=0 TO N-1 DO
Zcmplx[k,j]:=CMPLX(Zreal[k,j],Zreal[k+1,j]);
END;
END; (* IF *)
END;
END KonvCmplxC1;
PROCEDURE KonvCmplxC2( N : CARDINAL;
VAR Zreal : ARRAY OF ARRAY OF LONGREAL;
VAR Zr,Zi : ARRAY OF ARRAY OF LONGREAL;
VAR Wi : ARRAY OF LONGREAL);
VAR j,k : CARDINAL;
BEGIN
FOR k:=0 TO N-1 DO
IF (Wi[k] = 0.0) THEN
FOR j:=0 TO N-1 DO
Zr[k,j] := Zreal[k,j];
Zi[k,j] := 0.0;
END;
ELSIF (Wi[k] < 0.0) THEN (* konjugiert komplexe Wurzel, k > 1 *)
FOR j:=0 TO N-1 DO
Zr[k,j] := Zr[k-1,j];
Zi[k,j] := -Zi[k-1,j];
END;
ELSE
FOR j:=0 TO N-1 DO
Zr[k,j] := Zreal[k ,j];
Zi[k,j] := Zreal[k+1,j];
END;
END; (* IF *)
END;
END KonvCmplxC2;
PROCEDURE InfMatNorm( N : CARDINAL;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR imax : CARDINAL;
VAR Norm : LONGREAL);
(*----------------------------------------------------------------*)
(* Berechnet die Zeilensummennorm der Matrix A *)
(* Calculates the infinity norm of the Matrix A *)
(*----------------------------------------------------------------*)
VAR i,j : CARDINAL;
r : LONGREAL;
BEGIN
Norm:=0.0;
FOR i:=0 TO N-1 DO
r:=0.0;
FOR j:=0 TO N-1 DO r:=r + ABS(A[i,j]); END;
IF (r > Norm) THEN Norm:=r; imax:=i; END;
END;
END InfMatNorm;
PROCEDURE CNormMat(VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
dim : CARDINAL;
MaxNorm : BOOLEAN);
(*----------------------------------------------------------------*)
(* Normiert die in A "ubergebene komplex Matrix. *)
(* Wenn MaxNorm = TRUE, wird die Maximumnorm berechnet, *)
(* ansonsten die euklidische Norm. *)
(*----------------------------------------------------------------*)
VAR i,j : CARDINAL;
Amax,Norm : LONGREAL;
Sum : LONGCOMPLEX;
BEGIN
FOR i:=0 TO dim-1 DO
IF MaxNorm THEN
Amax:=0.0;
FOR j:=0 TO dim-1 DO
IF (ABS(RE(A[i,j])) > Amax) THEN
Amax := ABS(RE(A[i,j]));
END;
IF (ABS(IM(A[i,j])) > Amax) THEN
Amax := ABS(IM(A[i,j]));
END;
END;
Norm := 1.0 / Amax;
ELSE
Sum:=zero;
FOR j:=0 TO dim-1 DO
Sum := Sum + A[i,j]*conj(A[i,j]);
END;
Norm := 1.0 / sqrt(abs(Sum));
END;
FOR j:=0 TO dim-1 DO
A[i,j]:=scalarMult(Norm,A[i,j]);
END;
END; (* FOR i *)
END CNormMat;
PROCEDURE SplitHess( N : CARDINAL;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR H : ARRAY OF ARRAY OF LONGREAL;
oben : BOOLEAN);
(*----------------------------------------------------------------*)
(* Trennt die Hessenbergmatrix H aus dem Ergebniss A' der Trans- *)
(* formations mit ElmHess oder einer vergleichbaren Routine ab. *)
(* Wenn oben = wahr wird die obere, ansonsten die untere Hessen- *)
(* bergmatrix abgetrennt. *)
(*----------------------------------------------------------------*)
VAR i,j,js : CARDINAL;
BEGIN
IF (oben) THEN
FOR i:=0 TO N-1 DO
IF (i > 1) THEN js:=i-1; ELSE js:=0; END;
IF (js > 0) THEN
FOR j:=0 TO js-1 DO H[i,j]:=0.0; END;
END;
FOR j:=js TO N-1 DO H[i,j]:=A[i,j]; END;
END;
ELSE
FOR i:=0 TO N-1 DO
IF (i <= N-3) THEN js:=i+1; ELSE js:=N-1; END;
FOR j:=0 TO js DO H[i,j]:=A[i,j]; END;
FOR j:=js+1 TO N-1 DO H[i,j]:=-0.0; END;
END;
END;
END SplitHess;
PROCEDURE CmplxComRes( N,ii : INTEGER;
VAR W : ARRAY OF LONGCOMPLEX;
VAR V : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR norm : LONGREAL);
VAR j,k : INTEGER;
s,ss : LONGCOMPLEX;
BEGIN
DEC(ii);
ss:=CMPLX(0.0,0.0);
FOR j:=0 TO N-1 DO (* s) = W[i]*V - A*V *)
s := - W[ii]*V[ii,j];
FOR k:=0 TO N-1 DO
s:=s + A[j,k]*V[ii,k];
END;
ss:=ss + s*s;
END;
norm:=sqrt(abs(ss)) / LFLOAT(N-1);
END CmplxComRes;
PROCEDURE IstHermitisch(VAR H : ARRAY OF ARRAY OF LONGCOMPLEX;
dim : CARDINAL;
VAR I,J : CARDINAL) : BOOLEAN;
VAR i,j : CARDINAL;
BEGIN
FOR i:=0 TO dim-1 DO
IF (i > 0) THEN
FOR j:=0 TO i-1 DO
IF (RE(H[i,j]) # RE(H[j,i])) THEN
I:=i+1; J:=j+1;
RETURN FALSE;
END;
IF ((IM(H[i,j]) + IM(H[j,i])) # 0.0) THEN
I:=i+1; J:=j+1;
RETURN FALSE;
END;
END;
END;
IF (IM(H[i,i]) # 0.0) THEN
I:=i+1; J:=i+1;
RETURN FALSE;
END;
END;
I:=0; J:=0;
RETURN TRUE;
END IstHermitisch;
PROCEDURE CheckOrtho(VAR A : ARRAY OF ARRAY OF LONGREAL;
N : CARDINAL;
VAR S : ARRAY OF LONGREAL;
VAR Smax : LONGREAL;
VAR Ssqr : LONGREAL);
VAR i,j,ij : CARDINAL;
BEGIN
Smax := 0.0; Ssqr := 0.0; ij:=1;
S[0]:= ddot(N,A[0,0],1,A[0,0],1);
FOR i:=1 TO N-1 DO
FOR j:=0 TO i-1 DO
S[ij]:= ddot(N,A[i,0],1,A[j,0],1);
IF (ABS(S[ij]) > Smax) THEN Smax := ABS(S[ij]); END;
Ssqr:=Ssqr + S[ij]*S[ij];
INC(ij);
END;
S[ij]:= ddot(N,A[i,0],1,A[i,0],1); INC(ij); (* Diagonalelement *)
END;
Ssqr := Ssqr / LFLOAT(N*(N-1) DIV 2);
END CheckOrtho;
END EigenLibAux.