EigenLib1.mod.m2pp
2399 lines (2239 with data), 86.2 kB
IMPLEMENTATION MODULE EigenLib1;
(*========================================================================*)
(* HINWEIS : Bitte nur die Datei EigenLib1.mod.m2pp editieren *)
(*========================================================================*)
(* Es sind 2 Versionen enthalten die mit *)
(* *)
(* m2pp -D __{Parameter}__ < EigenLib1.mod.m2pp > EigenLib1.mod *)
(* *)
(* mit Parameter = {inline|BLAS} erzeugt werden koennen. *)
(* *)
(* inline : Alle Schleifen werden expliziet ausgefuehrt *)
(* BLAS : Schleifen werden durch Routinen aus LibDBlas ausgefuehrt *)
(* *)
(* Mit dem zusatzlichen Parameter {DEBUG} werden in einigen Routinen *)
(* zusaetzliche Ausgaben auf dem Standardausgabekanal erzeugt. *)
(*========================================================================*)
(* Stellt Routinen zur Eigenwert- und Vektorenberechnung symmetrischer *)
(* reeller Matritzen zur Verf"ugung. *)
(* *)
(* Routines for calculating eigensystems of symmetric real matrices *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung: *)
(* *)
(* 03.08.93, MRi: Erstellend der ersten Version von Reduce1 und ReBakA *)
(* 27.11.94, MRi: Durchsicht *)
(* 03.08.15, MRi: Ersetzen von MinGenau durch MachEps, von MachEps *)
(* durch Small um konsistent mit anderen Definitionen zu *)
(* werden. Anpassen der Unterlaufdedektion in *)
(* TriQR.TransQR, *)
(* benutzung Pythag statt sqrt(x*x + y*y) *)
(* 13.10.15, MRi: Ersetzen von LFLOAT(#) durch VAL(LONGREAL,#) *)
(* 01.04.16, MRi: Korrektur bei der Berechung des Indexs des letzten *)
(* Elements von HD in GivTri, in GivTriBak.SinCos die Ab- *)
(* frage ELSIF ABS(rho) >= srh THEN in ELSE ... geaendert, *)
(* 14.04.16, MRi: Umstellen von SortEwEv auf offenen Felder *)
(* 17.12.16, MRi: In Jacobi.Rotiere die Berechnung der Drehwinkel umge- *)
(* stellt. *)
(* 18.12.16, MRi: Routine Jacobi2D und SortEigen eingefuegt. *)
(* 01.01.17, MRi: Leicht ueberarbeitete Routine Wielandt eingefuegt. *)
(* 02.01.17, MRi: Jacobi (1D) vollstaendig auf offene Felder umgestellt. *)
(* PJacobi leicht angepasst damit Prozedur Rotier weiter *)
(* fuer beide Routinen verwendet werden kann. *)
(* 10.01.17, MRi: Routine Kaiser aufgenommen *)
(* 11.01.17, MRi: In Jacobi & PJacobi genau auf MachEps/MachEpsR4 bezogen *)
(* 17.01.17, MRi: In Rotier bei der Berechnung der Drehwinkel Division *)
(* durch Null abgefangen *)
(* 26.02.17, MRi: Erste Version von ImTQL2 - tuts noch nicht *)
(* 21.05.17, MRi: Korrekturen an Reduce1, Parameter transB eingefuegt, *)
(* Reduce1 und ReBakA auf offene Felder umgestellt *)
(* 22.05.17, MRi: In Reduce1 Aufrufe von BLAS level 1 ddot eingefuegt *)
(* und zusammen mit ReBakA in Modul EigenLib1 eingefuegt *)
(* 04.08.17, MRi: Korrekturen an ImTQL2 und Umstellung auf offene Felder *)
(* - sieht gut aus *)
(* 27.09.17, MRi: Tred2,EWEV,GivTri,GivTriBak aus EigenLib2 uebernommen *)
(* 24.10.17, MRi: MatSvProd durch MatSvProdNN ersetzt *)
(*------------------------------------------------------------------------*)
(* Offene Punkte: *)
(* *)
(* - Prozedure {P}SortEwEV auf schnelleren Sortieralgorithmus umstellen *)
(* - Prozeduren auf offene Felder umstellen *)
(* - Prozedure TransfSV auf neue Matrixmultiplikationsroutinen umstellen *)
(*------------------------------------------------------------------------*)
(* Testroutinen fuer Jacobi(1D,dyn,2D), HhQL, RSymEwEv in TestEigenRS *)
(* Testroutinen fuer Loewdin und Reduce1 in AxEQlBx *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: EigenLib1.mod.m2pp,v 1.12 2018/01/16 09:19:51 mriedl Exp mriedl $ *)
FROM SYSTEM IMPORT ADR,TSIZE;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
FROM LongMath IMPORT sqrt;
FROM LMathLib IMPORT MachEps,MachEpsR4,Small,MaxCard,Pythag;
FROM Errors IMPORT Fehler,Fehlerflag;
IMPORT Errors; (* F"ur Prozeduren. *)
FROM Deklera IMPORT VEKTOR,MATRIX,SUPERVEKTOR,CARDVEKTOR,PMATRIX,
MaxPMat,MaxDim;
FROM MatLib IMPORT MatVekProd,NormMat,QuickSort,SkalarProd,
MatSvProdNN,MatMatProd,MatToSv,Stuerz;
FROM LinLib IMPORT Gauss,Householder,Invert;
<* IF (__BLAS__) THEN *>
FROM LibDBlas IMPORT ddot;
<* END *>
FROM FileSystem IMPORT StdErr;
IMPORT FIO2;
<* IF (__DEBUG__) THEN *>
IMPORT TIO;
<* END *>
CONST MinGenau = 10.0*MachEps;
PROCEDURE Gerschgorin(VAR A : ARRAY OF LONGREAL;
dim : CARDINAL) : LONGREAL;
VAR i,j,ii,ij : CARDINAL;
Max,Sum : LONGREAL;
BEGIN
Max:=0.0; ii:=0;
FOR i:=1 TO dim DO
Sum:=0.0; ij:=ii;
FOR j:=1 TO i DO Sum:=Sum + ABS(A[ij]); INC(ij); END; INC(ij,i-1);
FOR j:=i TO dim DO Sum:=Sum + ABS(A[ij]); INC(ij,j); END;
IF (Sum > Max) THEN Max:=Sum; END;
INC(ii,i);
END;
RETURN Max;
END Gerschgorin;
PROCEDURE SortEwEv(VAR EW : ARRAY OF LONGREAL;
VAR EV : ARRAY OF ARRAY OF LONGREAL;
dim,M : CARDINAL;
Ordnung : INTEGER);
VAR i,l,iMinMax : CARDINAL;
MinMax,Zw : LONGREAL;
BEGIN
FOR l:=0 TO M-1 DO
MinMax:=EW[l]; iMinMax:=l;
FOR i:=l+1 TO dim-1 DO
IF (Ordnung = 1) THEN (* Berechne Maximum von EW[dim] *)
IF (EW[i] > MinMax) THEN MinMax:=EW[i]; iMinMax:=i; END;
ELSE (* Berechne Minimum von EW[dim] *)
IF (EW[i] < MinMax) THEN MinMax:=EW[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;
FOR i:=0 TO dim-1 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 SortEwEv;
PROCEDURE PSortEwEv(VAR EW : ARRAY OF LONGREAL;
VAR EV : PMATRIX;
dim : CARDINAL;
Ordnung : INTEGER);
VAR i,l,iMinMax : CARDINAL;
MinMax,Zw : LONGREAL;
BEGIN
FOR l:=1 TO dim DO
MinMax:=EW[l-1]; iMinMax:=l;
FOR i:=l+1 TO dim DO
IF (Ordnung = 1) THEN (* Berechne Maximum von EW[dim] *)
IF (EW[i-1] > MinMax) THEN MinMax:=EW[i-1]; iMinMax:=i; END;
ELSE (* Berechne Minimum von EW[dim] *)
IF (EW[i-1] < MinMax) THEN MinMax:=EW[i-1]; iMinMax:=i; END;
END; (* IF Ordnung *)
END; (* FOR i *)
IF (iMinMax # l) THEN (* Vertausche EWi mit EWimax *)
Zw:=EW[l-1]; EW[l-1]:=EW[iMinMax-1]; EW[iMinMax-1]:=Zw;
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 PSortEwEv;
PROCEDURE SortEigen( n : CARDINAL;
VAR d : ARRAY OF LONGREAL;
VAR r : ARRAY OF CARDINAL);
VAR k,l,aux : CARDINAL;
BEGIN
FOR k:=0 TO n-1 DO r[k]:=k+1; END;
FOR k:=0 TO n-2 DO
FOR l:=k+1 TO n-1 DO
IF (d[r[k]-1] > d[r[l]-1]) THEN
aux:=r[k]; r[k]:=r[l]; r[l]:=aux;
END;
END;
END;
END SortEigen;
PROCEDURE Jacobi2D( n : CARDINAL;
eivec : BOOLEAN;
VAR A : ARRAY OF ARRAY OF LONGREAL;
VAR d : ARRAY OF LONGREAL;
VAR V : ARRAY OF ARRAY OF LONGREAL;
VAR rot : CARDINAL);
VAR sm,c,s,t,h,g : LONGREAL;
tau,theta,tresh : LONGREAL;
p,q,j,iter : CARDINAL;
b,z : POINTER TO
ARRAY [0..MAX(INTEGER)-1] OF LONGREAL;
MaxIter : CARDINAL;
BEGIN
ALLOCATE(b,n*TSIZE(LONGREAL));
ALLOCATE(z,n*TSIZE(LONGREAL));
IF (b = NIL) OR (z = NIL) THEN
IF (b # NIL) THEN DEALLOCATE(b,n*TSIZE(LONGREAL)); END;
Fehler:=TRUE;
Fehlerflag:="Kein Freispeicher vorhanden (Jacobi2D)";
RETURN;
END;
IF eivec THEN
FOR p:=0 TO n-1 DO
FOR q:=0 TO n-1 DO V[p,q]:=0.0; END;
V[p,p]:=1.0;
END;
END;
FOR p:=0 TO n-1 DO
b^[p] := A[p,p];
d [p] := b^[p];
z^[p] := 0.0;
END;
Fehler:=FALSE;
MaxIter:=n*n*(n DIV 2); IF (MaxIter < 30) THEN MaxIter:=30; END;
rot := 0;
iter := 0;
LOOP
INC(iter);
sm:=0.0;
FOR p:=0 TO n-2 DO
FOR q:=p+1 TO n-1 DO sm:=sm + ABS(A[p,q]); END;
END;
IF (sm = 0.0) THEN EXIT; END;
(***************************)
tresh:= 0.0;
IF (iter < 4) THEN tresh:= 0.2*sm / VAL(LONGREAL,(n*n)); END;
FOR p:=0 TO n-2 DO
FOR q:=p+1 TO n-1 DO
g:= 100.0*ABS(A[p,q]);
IF (iter > 4) AND (ABS(d[p]) + g = ABS(d[p]))
AND (ABS(d[q]) + g = ABS(d[q]))
THEN
A[p,q] := 0.0;
ELSIF (ABS(A[p,q]) > tresh) THEN
h := d[q] - d[p];
IF (ABS(h) + g = ABS(h)) THEN
t := A[p,q] / h;
ELSE
theta := 0.5*h / A[p,q];
t := 1.0 / (ABS(theta) + sqrt(1.0 + theta*theta));
IF (theta < 0.0) THEN t:= -t; END;
END; (* computing tan of rotation angle *)
c := 1.0 / sqrt(1.0 + t*t);
s := t*c;
tau:= s / (1.0 + c);
h:= t*A[p,q];
z^[p]:=z^[p] - h;
z^[q]:=z^[q] + h;
d [p]:=d [p] - h;
d [q]:=d [q] + h;
A[p,q]:=0.0;
IF (p > 0) THEN
FOR j:=0 TO p-1 DO
g := A[j,p];
h := A[j,q];
A[j,p] := g - s*(h + g*tau);
A[j,q] := h + s*(g - h*tau)
END; (* of case 1<=j<p *)
END;
IF (q > 0) THEN
FOR j:=p+1 TO q-1 DO
g := A[p,j];
h := A[j,q];
A[p,j] := g - s*(h + g*tau);
A[j,q] := h + s*(g - h*tau)
END; (* of case p<j<q *)
END;
FOR j:=q+1 TO n-1 DO
g := A[p,j];
h := A[q,j];
A[p,j] := g - s*(h + g*tau);
A[q,j] := h + s*(g - h*tau)
END; (* of case q<j<=n *)
IF eivec THEN
FOR j:=0 TO n-1 DO
g := V[p,j];
h := V[q,j];
V[p,j] := g - s*(h + g*tau);
V[q,j] := h + s*(g - h*tau)
END; (* of case v *)
END;
INC(rot);
END; (* rotate *)
END; (* q *)
END; (* p *)
FOR p:=0 TO n-1 DO
b^[p] := b^[p] + z^[p];
d [p] := b^[p];
z^[p] := 0.0;
END; (* p *)
IF (iter > MaxIter) THEN
Fehler:=TRUE;
Fehlerflag:="Maximalzahl der Iterationen ueberschritten (Jacobi2D)";
EXIT;
END;
END; (* LOOP *)
DEALLOCATE(z,n*TSIZE(LONGREAL));
DEALLOCATE(b,n*TSIZE(LONGREAL));
END Jacobi2D;
PROCEDURE Rotier(VAR A : ARRAY OF LONGREAL;
VAR EW : ARRAY OF LONGREAL;
VAR EVp,EVq : ARRAY OF LONGREAL; (* Eigenvektoren p,q *)
pq : CARDINAL;
p,q : CARDINAL;
dim : CARDINAL;
VAR Schranke : LONGREAL; (* Schranke f"ur Pivot *)
pk,qk : INTEGER);
(*----------------------------------------------------------------*)
(* Berechnet den Drehwinkel sin/cos phi der entsprechenden *)
(* Jacobiiteration und transformiert A. *)
(* Baut dabei auch die Eigenvektormatrix der entsprechenden Iter- *)
(* ation. *)
(* Cos \in [\sqrt 1/2,1], Sin \in [-\sqrt 1/2,\sqrt 1/2] *)
(* Unterroutine zu Jacobi und PJacobi. *)
(*----------------------------------------------------------------*)
CONST sqrt05 = 0.7071067811865475; (* 1.0 / sqrt(0.5) *)
VAR Apq,Apk,EVpk,TanApq : LONGREAL;
Cos,Sin,Tan,Omega,Delta,R : LONGREAL;
k : CARDINAL;
BEGIN
DEC(pk); DEC(qk); (* Wegen "Ubergabe von A als offenes Feld *)
Apq:=A[pq]; A[pq]:=0.0;
Delta := EW[q] - EW[p];
(*
(* Warum geht das in sehr seltenen Faellen schief ??? *)
IF (ABS(Delta) < ABS(Apq)*Schranke) THEN
IF (Apq < 0.0) THEN Sin:=-sqrt05; ELSE Sin:=sqrt05; END;
Cos:=sqrt05;
Tan:=1.0;
ELSE (* Berechne Sinus und Cosinus des Matrixdrehwinkels *)
Omega := Delta / (2.0*Apq);
IF (Omega >= 0.0) THEN
Tan := 1.0 / (Omega + sqrt(1.0 + Omega*Omega));
ELSE
Tan := 1.0 / (Omega - sqrt(1.0 + Omega*Omega));
END;
Cos := 1.0 / sqrt(1.0 + Tan*Tan);
Sin := Cos*Tan;
END;
*)
(*
IF (ABS(Delta) > Small) THEN
IF (ABS(Delta) < ABS(Apq)*Schranke) THEN
Tan := Apq / Delta;
ELSE
Omega := Delta / (2.0*Apq);
Tan := 1.0 / (ABS(Omega) + sqrt(1.0 + Omega*Omega));
IF (Omega < 0.0) THEN Tan := - Tan; END;
END;
Cos := 1.0 / sqrt(1.0 + Tan*Tan);
Sin := Tan*Cos;
ELSE (* Vermeide Division durch Null *)
Tan:=1.0; Sin:=sqrt05; Cos:=sqrt05;
END;
*)
IF (ABS(Delta) + 100.0*ABS(Apq) = ABS(Delta)) THEN
Tan := Apq / Delta;
ELSE
Omega := 0.5*Delta / Apq;
Tan := 1.0 / (ABS(Omega) + sqrt(1.0 + Omega*Omega));
IF (Omega < 0.0) THEN Tan:= -Tan; END;
END; (* computing tan of rotation angle *)
Cos := 1.0 / sqrt(1.0 + Tan*Tan);
Sin := Tan*Cos;
TanApq:=Tan*Apq;
EW[p]:=EW[p] - TanApq; (* Transformation der *)
EW[q]:=EW[q] + TanApq; (* Schnittpunktelemente. *)
IF (ABS(Sin) > 1.0E-05) THEN (* Sin,Cos vergleichbare Gr"o3enordnung *)
FOR k:=1 TO q DO
INC(pk); INC(qk); Apk:=A[pk];
A[pk]:=Apk*Cos - A[qk]*Sin;
A[qk]:=Apk*Sin + A[qk]*Cos;
END;
INC(pk); INC(qk);
FOR k:=q+1 TO p-1 DO
INC(pk); INC(qk,k); Apk:=A[pk];
A[pk]:=Apk*Cos - A[qk]*Sin;
A[qk]:=Apk*Sin + A[qk]*Cos;
END;
INC(pk); INC(qk,p);
FOR k:=p+1 TO dim-1 DO
INC(pk,k); INC(qk,k); Apk:=A[pk];
A[pk]:=Apk*Cos - A[qk]*Sin;
A[qk]:=Apk*Sin + A[qk]*Cos;
END;
FOR k:=0 TO dim-1 DO (* Berechne neue Eigenvektormatrix *)
EVpk := EVp[k];
EVp[k] := EVpk*Cos - EVq[k]*Sin;
EVq[k] := EVpk*Sin + EVq[k]*Cos;
END;
ELSE (* Wenn Sin << Cos, dann numerisch stabiler, aber langsamer *)
R:= Sin / (1.0 + Cos);
FOR k:=0 TO dim-1 DO
EVpk := EVp[k]; (* Berechne neue Eigenvektormatrix *)
EVp[k] := EVpk - Sin*(EVq[k] + R*EVpk );
EVq[k] := EVq[k] + Sin*(EVpk - R*EVq[k]);
END;
FOR k:=1 TO q DO
INC(pk); INC(qk); Apk:=A[pk];
A[pk] := Apk - Sin*(A[qk] + R*Apk );
A[qk] := A[qk] + Sin*(Apk - R*A[qk]);
END;
INC(pk); INC(qk);
FOR k:=q+1 TO p-1 DO
INC(pk); INC(qk,k); Apk:=A[pk];
A[pk] := Apk - Sin*(A[qk] + R*Apk );
A[qk] := A[qk] + Sin*(Apk - R*A[qk]);
END;
INC(pk); INC(qk,p);
FOR k:=p+1 TO dim-1 DO
INC(pk,k); INC(qk,k); Apk:=A[pk];
A[pk] := Apk - Sin*(A[qk] + R*Apk );
A[qk] := A[qk] + Sin*(Apk - R*A[qk]);
END;
END; (* IF Sin *)
END Rotier;
PROCEDURE Jacobi(VAR EV : ARRAY OF ARRAY OF LONGREAL; (* Eigenvektoren *)
VAR EW : ARRAY OF LONGREAL; (* Eigenwert-Vektor *)
VAR A : ARRAY OF LONGREAL; (* Zu diag. sym. Matrix *)
dim : CARDINAL; (* Gr"o3e der Matrix *)
MaxIter : CARDINAL; (* Maximale Iterationszahl *)
genau : LONGREAL; (* Geforderte Genauigkeit *)
Ordnung : INTEGER); (* Sortieroption *)
VAR (* Schwa,Schwi : Schwllwerte absolut / der i.-Iteration *)
SumQ,Anzahl,Schwi,Schwa : LONGREAL;
p,q,pq,Iter : CARDINAL;
Rotation : BOOLEAN;
Tab : ARRAY [0..MaxDim-1] OF INTEGER;
BEGIN
FOR p:=0 TO dim-1 DO FOR q:=0 TO dim-1 DO EV[p][q]:=0.0; END; END;
SumQ:=0.0; pq:=0; (* Ber. die Quadratsumme der nicht-Diag.-elem.*)
FOR p:=0 TO dim-1 DO
EV[p,p]:=1.0; (* E-Matrix als 1. N"aherung der Eigenvektormatrix *)
Tab[p]:=VAL(INTEGER,pq); (* Tab[p]:=(p*(p-1) DIV 2) *)
INC(pq);
IF (p > 0) THEN
FOR q:=0 TO p-1 DO SumQ:=SumQ + A[pq-1]*A[pq-1]; INC(pq); END;
END;
EW[p]:=A[pq-1]; (* Diagonalelemente von A auf EW legen *)
END;
(*
* IF (genau <= 0.0) THEN genau:=1.0E-10; END;
* IF (genau >= 1.0E-04) THEN genau:=1.0E-04; END;
*)
IF (genau < MachEps* 100.0) THEN genau:=MachEps* 100.0; END;
IF (genau > MachEpsR4*10.0) THEN genau:=MachEpsR4*10.0; END;
IF (MaxIter = 0) THEN MaxIter:=(dim DIV 2)*(dim*dim); END;
IF (MaxIter < 15) THEN MaxIter:= 15; END;
Anzahl:=2.0 / VAL(LONGREAL,dim*(dim-1));
Schwa:=genau*sqrt(Anzahl*SumQ);
Errors.Fehler:=FALSE; Iter:=0;
LOOP
REPEAT
Rotation:=FALSE;
IF (SumQ <= 0.0) THEN
Schwi:=Schwa;
ELSE
Schwi:=sqrt(Anzahl*SumQ);
IF (Schwi < Schwa) THEN Schwi:=Schwa; END;
END;
pq:=0;
FOR p:=1 TO dim-1 DO
INC(pq);
FOR q:=0 TO p-1 DO
IF (ABS(A[pq]) > Schwi) THEN (* Finde Pivot-Element *)
Rotation:=TRUE; INC(Iter);
SumQ:=SumQ - A[pq]*A[pq];
Rotier(A,EW,EV[p],EV[q],pq,p,q,dim,Schwa,Tab[p],Tab[q]);
IF (Iter >= MaxIter) THEN Errors.Fehler:=TRUE; EXIT; END;
END; (**************************)
INC(pq);
END; (* FOR q *)
END; (* FOR p *)
UNTIL NOT Rotation;
SumQ:=0.0; (* Ein zus"atzlicher Durchlauf ! *)
IF (Schwi <= Schwa) THEN EXIT END;
END; (* LOOP *)
IF Errors.Fehler THEN
Errors.ErrOut('Iter > MaxIter (EigenLib1.Jacobi) !');
END;
IF (ABS(Ordnung) = 1) THEN SortEwEv(EW,EV,dim,dim,Ordnung); END;
END Jacobi;
PROCEDURE PJacobi(VAR EV : PMATRIX; (* Eigenvektor-Matrix *)
VAR EW : ARRAY OF LONGREAL; (* Eigenwert-Vektor *)
VAR A : ARRAY OF LONGREAL; (* Zu diag. sym. Matrix *)
dim : CARDINAL; (* Gr"o3e der Matrix *)
MaxIter : CARDINAL; (* Max. Iterationszahl *)
genau : LONGREAL; (* Geford. Genauigkeit *)
Ordnung : INTEGER); (* Sortieroption *)
VAR (* Schwa,Schwi : Schwllwerte absolut / der i.-Iteration *)
SumQ,Anzahl,Schwi,Schwa : LONGREAL;
p,q,pq,Iter : CARDINAL;
Rotation : BOOLEAN;
Tab : POINTER TO
ARRAY [1..MaxPMat] OF INTEGER;
BEGIN
ALLOCATE(Tab,dim*TSIZE(INTEGER));
FOR p:=1 TO dim DO FOR q:=1 TO dim DO EV^[p]^[q]:=0.0; END; END;
SumQ:=0.0; pq:=0; (* Ber. die Quadratsumme der nicht-Diag.-elem.*)
FOR p:=1 TO dim DO
EV^[p]^[p]:=1.0; (* E-Matrix als 1. N"aherung der Eigenvektormatrix *)
Tab^[p]:=VAL(INTEGER,pq); (* Tab[p]:=(p*(p-1) DIV 2) *)
INC(pq);
FOR q:=1 TO p-1 DO SumQ:=SumQ + A[pq-1]*A[pq-1]; INC(pq); END;
EW[p-1]:=A[pq-1]; (* Diagonalelemente von A auf EW legen *)
END;
(*
* IF (genau <= 0.0) THEN genau:=1.0E-10; END;
* IF (genau >= 1.0E-04) THEN genau:=1.0E-04; END;
*)
IF (genau < MachEps* 100.0) THEN genau:=MachEps* 100.0; END;
IF (genau > MachEpsR4*10.0) THEN genau:=MachEpsR4*10.0; END;
IF (MaxIter = 0) THEN MaxIter:=(dim DIV 2)*(dim*dim); END;
IF (MaxIter < 15) THEN MaxIter:= 15; END;
Anzahl:=2.0 / VAL(LONGREAL,dim*(dim-1));
Schwa:=genau*sqrt(Anzahl*SumQ);
Errors.Fehler:=FALSE; Iter:=0;
LOOP
REPEAT
Rotation:=FALSE;
IF (SumQ <= 0.0) THEN
Schwi:=Schwa;
ELSE
Schwi:=sqrt(Anzahl*SumQ);
IF (Schwi < Schwa) THEN Schwi:=Schwa; END;
END;
pq:=0;
FOR p:=2 TO dim DO
INC(pq);
FOR q:=1 TO p-1 DO
IF (ABS(A[pq]) > Schwi) THEN (* Finde Pivot-Element *)
Rotation:=TRUE; INC(Iter);
SumQ:=SumQ - A[pq]*A[pq];
Rotier(A,EW,EV^[p]^,EV^[q]^,pq,p-1,q-1,dim,Schwa,Tab^[p],Tab^[q]);
IF (Iter >= MaxIter) THEN Errors.Fehler:=TRUE; EXIT END;
END; (*************************)
INC(pq);
END; (* FOR q *)
END; (* FOR p *)
UNTIL NOT Rotation;
SumQ:=0.0; (* Ein zus"atzlicher Durchlauf ! *)
IF (Schwi <= Schwa) THEN EXIT END;
END; (* LOOP *)
IF Errors.Fehler THEN
Errors.ErrOut('Iter > MaxIter (EigenLib1.Jacobi) !');
END;
IF (ABS(Ordnung) = 1) THEN PSortEwEv(EW,EV,dim,Ordnung); END;
DEALLOCATE(Tab,dim*TSIZE(INTEGER));
END PJacobi;
PROCEDURE Kaiser(VAR A : ARRAY OF ARRAY OF LONGREAL;
dim : CARDINAL;
VAR EW : ARRAY OF LONGREAL;
VAR iFehl : INTEGER);
PROCEDURE sqr(x : LONGREAL) : LONGREAL; BEGIN RETURN x*x; END sqr;
CONST small = 10.0*MachEps;
VAR i,j,k,nn,ncount,iter,MaxIter : CARDINAL;
absp,absq,halfp,p,q,ss,tmp,eps : LONGREAL;
ctn,CosA,SinA,TanA,Aik,Ajk : LONGREAL;
SumEW,Trace : LONGREAL;
BEGIN
IF (dim < 1) THEN
iFehl:=1;
RETURN;
END;
(* calculate convergence tolerance, eps. *)
(* calculate trace. initial settings. *)
Trace := 0.0; (* The trace of Matirx A *)
ss:=0.0;
FOR i:=0 TO dim-1 DO
Trace:=Trace + A[i,i];
FOR j:=0 TO dim-1 DO
ss:=ss + sqr(A[i,j]);
END;
END;
eps := small*ss / LFLOAT(dim);
nn := dim*(dim-1) DIV 2;
ncount := nn;
IF (dim < 20) THEN MaxIter:=10; ELSE MaxIter := dim DIV 2; END;
iFehl := 0;
iter := 0;
LOOP (* orthogonalize pairs of columns i & j, j > i. *)
INC(iter);
IF (iter > MaxIter) THEN iFehl:=2; EXIT; END;
FOR i:=0 TO dim-2 DO
FOR j:=i+1 TO dim-1 DO (* calculate planar rotation required *)
halfp:=0.0;
q:=0.0;
FOR k:=0 TO dim-1 DO
Aik := A[i,k];
Ajk := A[j,k];
halfp:=halfp + Aik*Ajk;
q:=q + (Aik + Ajk)*(Aik - Ajk);
END;
p := 2.0*halfp;
absp:=ABS(p);
IF (absp < eps) AND (q >= 0.0) THEN
(* if p is very small, the vectors are almost orthogonal. *)
(* skip the rotation if q >= 0 (correct ordering). *)
DEC(ncount);
IF (ncount = 0) THEN EXIT; END; (* Alle Elemente geprueft *)
ELSE (* rotation needed *)
absq:=ABS(q);
IF (absp <= absq) THEN
TanA := absp / absq;
CosA := 1.0 / sqrt(1.0 + TanA*TanA);
SinA := TanA*CosA;
ELSE
ctn := absq / absp;
SinA := 1.0 / sqrt(1.0 +ctn*ctn);
CosA := ctn*SinA;
END;
CosA := sqrt(0.5*(1.0 + CosA));
SinA := SinA / (2.0*CosA);
IF (q < 0.0) THEN
tmp := CosA;
CosA := SinA;
SinA := tmp;
END;
IF (p < 0.0) THEN
SinA:= - SinA;
END;
FOR k:=0 TO dim-1 DO (* perform rotation *)
Aik :=A[i,k];
Ajk :=A[j,k];
A[i,k] := Ajk*SinA + Aik*CosA;
A[j,k] := Ajk*CosA - Aik*SinA;
END;
END; (* IF *)
END; (* FOR k *)
END; (* FOR i *)
ncount:=nn;
END; (* LOOP *)
(* converged, or gave up after MaxIter iterations *)
SumEW := 0.0; (* SumEW = sum of the eigenvalues computed *)
FOR i:=0 TO dim-1 DO
tmp:=0.0;
FOR j:=0 TO dim-1 DO tmp:=tmp + sqr(A[i,j]); END;
EW[i] := sqrt(tmp);
SumEW := SumEW + EW[i];
END;
(* scale columns to have unit length *)
FOR i:=0 TO dim-1 DO
IF (EW[i] > 0.0) THEN tmp := 1.0 / EW[i]; ELSE tmp := 0.0; END;
FOR j:=0 TO dim-1 DO A[i,j]:=A[i,j]*tmp; END;
END;
IF ((SumEW - ABS(Trace)) > 0.1*SumEW*MachEpsR4) THEN
iFehl:=iFehl + 3;
Fehler:=TRUE;
Fehlerflag:="Matrix nicht positiv (semi-)definit (EigenLib1.Kaiser)";
END;
<* IF (__DEBUG__) THEN *>
TIO.WrLn;
TIO.WrStr("eps = "); TIO.WrLngReal(eps ,24,-9); TIO.WrLn;
TIO.WrStr("SumEW = "); TIO.WrLngReal(SumEW,24,-9); TIO.WrLn;
TIO.WrStr("Trace = "); TIO.WrLngReal(Trace,24,-9); TIO.WrLn;
TIO.WrStr("diff = "); TIO.WrLngReal(SumEW-ABS(Trace),24,-9); TIO.WrLn;
TIO.WrStr("iter = "); TIO.WrCard (iter ,10 ); TIO.WrLn;
TIO.WrStr("iFehl = "); TIO.WrInt (iFehl,10 ); TIO.WrLn;
TIO.WrLn;
<* END *>
END Kaiser;
PROCEDURE Wielandt( dim : CARDINAL;
VAR A : MATRIX;
VAR EVi : VEKTOR;
VAR EWi : LONGREAL;
genau : LONGREAL;
MaxIter : CARDINAL;
VAR nIter : CARDINAL;
VAR iFehl : INTEGER);
VAR Norm,TestNorm : LONGREAL;
i,Iter : CARDINAL;
iErr : INTEGER;
Regleigh : LONGREAL;
SkalarProdukt : LONGREAL;
Determ : LONGREAL;
Interim,Aii : VEKTOR;
PivStrat : CARDINAL;
BEGIN
iFehl:=0;
Fehler:=FALSE;
IF (MaxIter < 1) THEN MaxIter:=dim; END;
IF (genau <= 1.0E-10) THEN genau:=1.0E-10; END;
Norm:=0.0;
FOR i:=1 TO dim DO
Aii[i]:=A[i,i]; (* Sichere die Hauptdiagonale *)
A[i,i]:=A[i,i]-EWi;
Norm:=Norm+EVi[i]*EVi[i];
END;
Norm:=1.0/sqrt(Norm);
PivStrat:=0; (* Zeilenpivot *)
Iter:=0;
LOOP
INC(Iter);
FOR i:=1 TO dim DO
EVi[i]:=EVi[i]*Norm; (* Normiere Vektor *)
Interim[i]:=EVi[i]; (* Sicher Vektor *)
END;
(* Hier waere es gut wenn die LU-Zerlegung mit den neuen *)
(* Diagonalelementen nur aktualisiert werden muesste ... *)
Gauss(A,EVi,EVi,dim,Determ,1,PivStrat,iErr);
IF Fehler THEN
IF (PivStrat = 0) THEN iFehl:=1; ELSE iFehl:=2; END;
PivStrat:=1; (* Maxpivot in naechster Iteration *)
FOR i:=1 TO dim DO EVi[i]:=Interim[i]; END;
EWi:=EWi*(1.0+2.5E-16); (* Neuer Versuch ? *)
Fehler:=FALSE;
Householder(A,EVi,EVi,dim,dim);
IF Fehler THEN
FOR i:=1 TO dim DO EVi[i]:=Interim[i]; END;
Fehlerflag:='Eigenvektor nicht zu bestimmen (Wielandt)';
iFehl:=4;
EXIT;
END;
ELSE
TestNorm:=Norm;
Norm:=0.0;
FOR i:=1 TO dim DO Norm:=Norm+EVi[i]*EVi[i]; END;
SkalarProdukt:=Norm;
Norm:=1.0/sqrt(Norm);
MatVekProd(Interim,A,EVi,dim,dim,FALSE);
Regleigh:=SkalarProd(Interim,EVi,dim)/SkalarProdukt;
EWi:=EWi+Regleigh;
FOR i:=1 TO dim DO A[i,i]:=A[i,i]-Regleigh; END;
IF (ABS(Regleigh) < genau) THEN EXIT; END;
IF (ABS(TestNorm-Norm) < genau) THEN EXIT; END;
END; (* IF Fehler *)
IF (Iter >= MaxIter) THEN iFehl:=3; EXIT; END;
END; (* LOOP *)
nIter := Iter;
Norm:=0.0;
FOR i:=1 TO dim DO Norm:=Norm+EVi[i]*EVi[i]; END;
Norm:=1.0/sqrt(Norm);
FOR i:=1 TO dim DO
EVi[i]:=Norm*EVi[i]; (* Normiere den Vektor *)
A[i,i]:=Aii[i]; (* Schreibe die Haupdiagonale zur"uck *)
END;
END Wielandt;
PROCEDURE Tred2(VAR A : MATRIX; (* Zu Trilinearisierende Matrix *)
VAR Z : MATRIX; (* Reduzierte Matrix A *)
VAR HD : VEKTOR; (* Hauptdiag. der reduzierten Matrix A *)
VAR ND : VEKTOR; (* Nebendiag. der reduzierten Matrix A *)
dim : CARDINAL; (* Dimension von A *)
genau : LONGREAL); (* Genauigkeitswert *)
VAR i,j,k,l : CARDINAL;
f,g,h,hh,hr : LONGREAL;
BEGIN
IF (ADR(A) # ADR(Z)) THEN (* Kein unn"otiges Kopieren der Matrix *)
FOR i:=1 TO dim DO
FOR j:=1 TO dim DO
Z[i,j]:=A[i,j];
END;
END;
END;
IF (genau < MinGenau) THEN genau:=MinGenau; END;
FOR i:=dim TO 2 BY -1 DO
l:=i-2;
f:=Z[i-1,i];
g:=0.0;
FOR k:=1 TO l DO
g:=g+Z[k,i]*Z[k,i];
END;
h:=g+f*f;
IF (g <= genau) THEN
ND[i]:=f;
h:=0.0;
ELSE
INC(l);
g:=sqrt(h);
IF (f >= 0.0) THEN g:=-g; END;
ND[i]:=g;
h:=h-f*g;
Z[i-1,i]:=f-g;
f:=0.0;
hr:=1.0/h;
FOR j:=1 TO l DO
Z[i,j]:=Z[j,i]*hr;
g:=0.0;
FOR k:=1 TO j DO
g:=g+Z[k,j]*Z[k,i];
END;
FOR k:=j+1 TO l DO
g:=g+Z[j,k]*Z[k,i];
END;
ND[j]:=g*hr;
f:=f+g*Z[i,j];
END; (* FOR j *)
hh:=f/(h+h);
FOR j:=1 TO l DO
f:=Z[j,i];
g:=ND[j]-hh*f; ND[j]:=g;
FOR k:=1 TO j DO
Z[k,j]:=Z[k,j]-f*ND[k]-g*Z[k,i];
END;
END; (* FOR j *)
END; (* IF h *)
HD[i]:=h;
END; (* FOR i *)
HD[1]:=0.0; ND[1]:=0.0;
FOR i:=1 TO dim DO
l:=i-1;(**)
IF (HD[i] # 0.0) THEN
FOR j:=1 TO l DO
g:=0.0;
FOR k:=1 TO l DO
g:=g+Z[k,i]*Z[j,k];
END;
FOR k:=1 TO l DO
Z[j,k]:=Z[j,k]-g*Z[i,k];
END;
END; (* FOR j *)
END; (* IF D[j] *)
HD[i]:=Z[i,i];
Z[i,i]:=1.0;
FOR j:=1 TO l DO
Z[i,j]:=0.0; Z[j,i]:=0.0;
END;
END; (* FOR i *)
FOR i:=2 TO dim DO
ND[i-1]:=ND[i];
END;
ND[dim]:=0.0;
END Tred2;
PROCEDURE EWEV(VAR EV : MATRIX; (* Ausgabe : Eigenvektoren von A *)
VAR HD : VEKTOR; (* Ausgabe : Eigenwerte von A *)
VAR ND : VEKTOR; (* Nebendiagonale d. reduzierten Matrix A *)
dim : CARDINAL); (* Dimension von EV *)
PROCEDURE dSign(x,y : LONGREAL) : LONGREAL; (* FORTRAN DSIGN-Funktion *)
BEGIN
IF (y < 0.0) THEN x := -ABS(x); ELSE x := ABS(x); END; RETURN x;
END dSign;
CONST maxiter = 30;
epsi = MinGenau;
VAR ii,i,k,l,m,iter : CARDINAL;
b,c,f,g,p,r,s : LONGREAL;
BEGIN
Fehler:=FALSE;
IF (dim = 1) THEN RETURN END;
FOR l:=1 TO dim DO
iter:=0;
LOOP
INC(iter);
IF (iter > maxiter) THEN
Fehler:=TRUE;
Fehlerflag:='Maxiter ueberschritten (EigenLib2.EWEV) !';
Errors.ErrOut(Fehlerflag);
RETURN;
END; (* IF *)
m:=l-1;
REPEAT (* Suche kleine Subdiagonalelemente *)
INC(m);
UNTIL (m = dim) OR
(ABS(ND[m]) <= epsi*(ABS(HD[m])+ABS(HD[m+1])));
p:=HD[l];
IF (m = l) THEN EXIT END; (* Ausgang LOOP !!! *)
g:=(HD[l+1]-p)/(2.0*ND[l]);
r:=sqrt(g*g+1.0);
g:=HD[m]-p+ND[l]/(g+dSign(r,g));
s:=1.0;
c:=1.0;
p:=0.0;
FOR ii:=1 TO m-l DO
i:=m-ii;
(*ip1:=i+1;*)
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:=HD[i+1]-p;
r:=(HD[i]-g)*s+2.0*c*b;
p:=s*r;
HD[i+1]:=g+p;
g:=c*r-b;
FOR k:=1 TO dim DO (* Berechne Eigenvektoren *)
f:=EV[i+1,k];
EV[i+1,k]:=s*EV[i,k]+c*f;
EV[i,k]:=c*EV[i,k]-s*f;
END; (* FOR *)
END; (* FOR ii *)
HD[l]:=HD[l]-p;
ND[l]:=g;
ND[m]:=0.0;
END; (* LOOP *)
END; (* FOR l *)
END EWEV;
PROCEDURE ImTQL2( dim : CARDINAL;
VAR D,E : ARRAY OF LONGREAL;
VAR Z : ARRAY OF ARRAY OF LONGREAL;
VAR iErr : INTEGER);
CONST eps = 10.0*MachEps;
VAR N,i,l,k,m,its : INTEGER;
h,c,f,q,s,r,p : LONGREAL;
t1,t2 : LONGREAL;
underflow : BOOLEAN;
BEGIN
N:=VAL(INTEGER,dim);
iErr:=0;
FOR i:=1 TO N-1 DO E[i-1] := E[i]; END;
E[N-1]:= 0.0;
FOR l:=0 TO N-1 DO
its:= 0;
REPEAT
m:=l;
LOOP (* look for single small sub-diagonal element *)
IF (m = N-1) THEN EXIT; END;
t1:= ABS(D[m]) + ABS(D[m+1]);
t2:= t1 + ABS(E[m]);
IF (ABS(t1 - t2) < eps*t1) THEN EXIT; END;
INC(m);
END;
p:=D[l];
IF (m # l) THEN
IF (its = 30) THEN
iErr:=l; Fehler:=TRUE;
Fehlerflag:="MaxIter uerberschritten (EigenLib1.ImTQL2)";
Errors.ErrOut(Fehlerflag);
RETURN;
END;
INC(its);
(* form shift *)
q:= (D[l+1] - p) / (2.0*E[l]);
r:=Pythag(q,1.0);
IF (q < 0.0) THEN
q:= D[m] - p + E[l] / (q - r);
ELSE
q:= D[m] - p + E[l] / (q + r);
END;
s:=1.0;
c:=1.0;
p:=0.0;
underflow:=FALSE;
i := m - 1;
WHILE (i >= l) AND NOT underflow DO
f := s*E[i];
h := c*E[i];
r := Pythag(f,q);
E[i+1] := r;
IF (r = 0.0) THEN
underflow:=TRUE;
ELSE
s := f / r;
c := q / r;
q := D[i+1] - p;
r := (D[i] - q)*s + 2.0*c*h;
p := s*r;
D[i+1] := q + p;
q := c*r - h;
FOR k:=0 TO N-1 DO (* form vector *)
f:= Z[i+1,k];
Z[i+1,k]:= s*Z[i,k] + c*f;
Z[i ,k]:= c*Z[i,k] - s*f;
END;
DEC(i);
END; (* IF underflow *)
END; (* WHILE *)
IF NOT underflow THEN
D[l]:=D[l] - p;
E[l]:=q;
ELSE (* recover from underflow *)
D[i+1]:=D[i+1] - p;
END;
E[m]:=0.0;
END;
UNTIL (m = l);
END; (* l *)
FOR i:=0 TO N-1 DO
(* order eigenvalues and eigenvectors *)
k:=i; f:=D[i];
FOR l:=i+1 TO N-1 DO
IF (D[l] < f) THEN
k:=l; f:=D[l];
END;
END;
IF (k # i) THEN
D[k]:= D[i]; D[i]:= f;
FOR l:=0 TO N-1 DO
f:=Z[i,l]; Z[i,l]:=Z[k,l]; Z[k,l]:=f;
END;
END;
END; (* i *)
END ImTQL2;
PROCEDURE GivTri(VAR A : SUPERVEKTOR; (* ==> Marix *)
VAR HD : VEKTOR; (* <== Hauptdiagonale *)
VAR ND : VEKTOR; (* <== Nebendiagonale *)
dim : CARDINAL;
genau : LONGREAL);
CONST srh = 0.7071067811865475; (* û« *)
VAR i,j,k,p,ii,pp,kk,ik,kp : CARDINAL;
Omega,Delta,Interim,
sin,cos,sincos,sqrsin,sqrcos,
Aij,Aki,Aik,Apj,Aip,rho : LONGREAL;
Index : CARDVEKTOR;
BEGIN
ii:=0; FOR i:=1 TO dim DO Index[i]:=ii; INC(ii,i); END;
pp:=1; p:=1;
FOR j:=1 TO dim-2 DO
INC(p); ii:=Index[j+2];
FOR i:=j+2 TO dim DO
IF (A[ii+j] # 0.0) THEN
Aij:=A[ii+j]; Apj:=A[pp+j];
IF (ABS(Apj) < genau*ABS(Aij)) THEN
cos :=0.0; sin :=1.0;
sqrcos:=0.0; sqrsin:=1.0;
sincos:=0.0;
A[pp+j]:=-Aij;
ELSE
IF (Apj >= 0.0) THEN
Omega:= sqrt(Apj*Apj+Aij*Aij)
ELSE
Omega:=-sqrt(Apj*Apj+Aij*Aij);
END;
cos:=Apj/Omega; sin:=-Aij/Omega;
sincos:=sin*cos;
sqrsin:=sin*sin; sqrcos:=cos*cos;
A[pp+j]:=Omega;
END; (* IF *)
IF (sin = 1.0) THEN (* Berechne R"ucktransformations- *)
rho:=1.0 (* information. EV(J) ==> EV(A) *)
ELSIF (ABS(sin) < srh) THEN (* ³sin³ < cos *)
rho:=sin
ELSE
IF (sin >= 0.0) THEN rho:=1.0/cos ELSE rho:=-1.0/cos; END;
END; (* IF *)
A[ii+j]:=rho; (* A[ii+j]:=0.0; Abspeichern von rho *)
Aip:=A[ii+p];
Delta:=A[pp+p]-A[ii+i];
Interim:=Delta*sqrsin + 2.0*sincos*Aip;
A[pp+p]:=A[pp+p]-Interim;
A[ii+i]:=A[ii+i]+Interim;
A[ii+p]:=Delta*sincos+Aip*(sqrcos-sqrsin);
kk:=Index[j+2];
FOR k:=j+2 TO i-1 DO
kp:=kk+p; ik:=ii+k;
Aik:=A[ik];
A[ik]:=sin*A[kp]+cos*Aik;
A[kp]:=cos*A[kp]-sin*Aik;
INC(kk,k);
END; (* FOR k *)
kk:=ii+i;
FOR k:=i+1 TO dim DO
kp:=kk+p; ik:=kk+i;
Aki:=A[ik];
A[ik]:=sin*A[kp]+cos*Aki;
A[kp]:=cos*A[kp]-sin*Aki;
INC(kk,k);
END; (* FOR k *)
END; (* IF Aij *)
INC(ii,i);
END; (* FOR i *)
INC(pp,p);
END; (* FOR j *)
pp:=0;
FOR p:=1 TO dim-1 DO (* Spalte Haupt- und Nebendiagonale *)
INC(pp,p); (* von A ab. *)
HD[p]:=A[pp];
ND[p]:=A[pp+p];
END;
HD[dim]:=A[pp+dim]; (* A[dim*(dim+1) DIV 2] *)
ND[dim]:=0.0;
END GivTri;
PROCEDURE GivTriBak(VAR EV : MATRIX; (* Eigenvektoren *)
VAR A : SUPERVEKTOR; (* Transformierte Matrix *)
maxEV : CARDINAL; (* Anzahl der EV *)
VAR genau : LONGREAL;
dim : CARDINAL); (* Dimension von A *)
CONST srh = 0.707106781186547524; (* 1/sqrt(2) *)
PROCEDURE SinCos(VAR sin,cos : LONGREAL;
VAR rho : LONGREAL);
(*---------------------------------------------------------*)
(* Die Procedur SinCos berechnet aus dem Wert rho den Wert *)
(* von sin und cos zur"uck. *)
(*---------------------------------------------------------*)
BEGIN
IF (ABS(rho) = 1.0 ) THEN
sin:=1.0; cos:=0.0;
ELSIF (ABS(rho) < srh) THEN (* rho in (0,shr] *)
sin:=rho; cos:=sqrt(1.0-rho*rho);
ELSE (* ABS(rho) >= srh) , rho in (shr,2] *)
IF (rho >= 0.0) THEN
cos:=1.0/ABS(rho); sin:=sqrt(1.0-cos*cos);
ELSE
cos:=1.0/ABS(rho); sin:=-sqrt(1.0-cos*cos);
END;
END; (* IF *)
END SinCos;
VAR Q : MATRIX;
i,j,ii,k,p : CARDINAL;
Qik,Qpk,rho,sin,cos : LONGREAL;
Index : CARDVEKTOR;
BEGIN
ii:=0;
FOR i:=1 TO dim DO
Index[i]:=ii; INC(ii,i); Q[i,i]:=1.0;
FOR j:=1 TO i-1 DO Q[i,j]:=0.0; Q[j,i]:=0.0; END;
END;
FOR j:=dim-2 TO 1 BY -1 DO
p:=j+1;
FOR i:=dim TO j+2 BY -1 DO
rho:=A[Index[i]+j];
IF (ABS(rho) > genau) THEN
SinCos(sin,cos,rho);
FOR k:=p TO dim DO
Qpk:=Q[p,k]; Qik:=Q[i,k];
Q[p,k]:=Qpk*cos + Qik*sin;
Q[i,k]:=Qik*cos - Qpk*sin;
END;
END; (* IF rho *)
END; (* FOR i *)
END; (* FOR j *)
FOR i:=1 TO maxEV DO MatVekProd(EV[i],Q,EV[i],dim,dim,FALSE); END;
END GivTriBak;
PROCEDURE HhQL(VAR A : MATRIX; (* Zu Diagonalisierende Matrix *)
VAR EV : MATRIX; (* Eigenvektoren *)
VAR EW : ARRAY OF LONGREAL; (* Eigenwerte von A *)
dim : CARDINAL; (* Dimension von A *)
genau : LONGREAL; (* Genauigkeitswert *)
ord : INTEGER); (* Sortieroption *)
(* Berechnung der Eigenwerte/Vektoren der Trilinearmatrix nach dem *)
(* Algorithmus der Handbuch-Routine imtql2 *)
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;
PROCEDURE dSign(x,y : LONGREAL) : LONGREAL; (* FORTRAN DSIGN-Funktion *)
BEGIN
IF (y < 0.0) THEN x:=-ABS(x) ELSE x:=ABS(x) END; RETURN x;
END dSign;
CONST maxiter = 30;
epsi = MachEps;
VAR ii,i,j,k,l,m,iter : CARDINAL;
b,c,f,g,p,r,s,h,hh,hr : LONGREAL;
ND : VEKTOR;
BEGIN
IF (ADR(A) # ADR(EV)) THEN (* Kein unn"otiges Kopieren der Matrix *)
FOR i:=1 TO dim DO FOR j:=1 TO dim DO EV[i,j]:=A[i,j]; END; END;
END;
IF (genau < MachEps) THEN genau:=MachEps; END;
FOR i:=dim TO 2 BY -1 DO
l:=i-2;
f:=EV[i-1,i];
g:=0.0;
FOR k:=1 TO l DO g:=g+EV[k,i]*EV[k,i]; END;
h:=g+f*f;
IF (g <= genau) THEN
ND[i]:=f;
h:=0.0;
ELSE
INC(l);
g:=sqrt(h);
IF (f >= 0.0) THEN g:=-g; END;
ND[i]:=g;
h:=h-f*g;
EV[i-1,i]:=f-g;
f:=0.0;
hr:=1.0/h;
FOR j:=1 TO l DO
EV[i,j]:=EV[j,i]*hr;
g:=0.0;
FOR k:=1 TO j DO g:=g+EV[k,j]*EV[k,i]; END;
FOR k:=j+1 TO l DO g:=g+EV[j,k]*EV[k,i]; END;
ND[j]:=g*hr;
f:=f+g*EV[i,j];
END; (* FOR j *)
hh:=f/(h+h);
FOR j:=1 TO l DO
f:=EV[j,i];
g:=ND[j]-hh*f; ND[j]:=g;
FOR k:=1 TO j DO EV[k,j]:=EV[k,j]-f*ND[k]-g*EV[k,i]; END;
END; (* FOR j *)
END; (* IF h *)
EW[i-1]:=h;
END; (* FOR i *)
EW[0]:=0.0; ND[1]:=0.0;
l:=0;
FOR i:=1 TO dim DO (* Bereche die Transformationsmatrix *)
IF (EW[i-1] # 0.0) THEN
FOR j:=1 TO l DO
g:=0.0;
FOR k:=1 TO l DO g:=g+EV[k,i]*EV[j,k]; END;
FOR k:=1 TO l DO EV[j,k]:=EV[j,k]-g*EV[i,k]; END;
END; (* FOR j *)
END; (* IF EW[j] *)
EW[i-1]:=EV[i,i]; EV[i,i]:=1.0;
FOR j:=1 TO l DO EV[i,j]:=0.0; EV[j,i]:=0.0; END;
INC(l);
END; (* FOR i *)
FOR i:=2 TO dim DO ND[i-1]:=ND[i]; END;
ND[dim]:=0.0;
FOR l:=1 TO dim DO (* Berechnung der Eigenwerte und Vektoren *)
iter:=0;
LOOP
INC(iter);
IF (iter > maxiter) THEN
Fehler:=TRUE; Fehlerflag:='Maxiter Ueberschritten (HhQL).';
Errors.ErrOut(Fehlerflag);
RETURN;
END; (* IF *)
m:=l-1;
REPEAT (* Suche kleine Subdiagonalelemente *)
INC(m);
UNTIL (m = dim) OR (ABS(ND[m]) <= epsi*AbsMin(EW[m-1],EW[m]));
p:=EW[l-1];
IF (m = l) THEN EXIT END; (* Ausgang LOOP !!! *)
g := (EW[l] - p) / (2.0*ND[l]);
r := sqrt(g*g + 1.0);
g := EW[m-1] - p+ND[l] / (g + dSign(r,g));
s:=1.0; c:=1.0; p:=0.0;
FOR ii:=1 TO m-l DO
i:=m-ii;
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] - p;
r := (EW[i-1]-g)*s + 2.0*c*b;
p := s*r;
EW[i]:=g+p;
g:=c*r-b;
FOR k:=1 TO dim DO (* Berechne Eigenvektoren *)
f:=EV[i+1,k];
EV[i+1,k] := s*EV[i,k] + c*f;
EV[i,k] := c*EV[i,k] - s*f;
END; (* FOR *)
END; (* FOR ii *)
EW[l-1]:=EW[l-1] - p;
ND[l]:=g;
ND[m]:=0.0;
END; (* LOOP *)
END; (* FOR l *)
IF (ABS(ord) = 1) THEN SortEwEv(EW,EV,dim,dim,ord); END;
END HhQL;
PROCEDURE HhTri(VAR A : SUPERVEKTOR; (* <==> Matrix *)
VAR HD : VEKTOR; (* ==> Hauptdiagonale von A *)
VAR ND : VEKTOR; (* ==> Nebendiagonale von A *)
VAR SqrND : VEKTOR; (* ==> Quadrate von ND *)
dim : CARDINAL); (* Dimension von A *)
VAR i,j,k,l,ii,jj,ik,jk : CARDINAL;
f,g,h,hh : LONGREAL;
BEGIN
<* IF (__DEBUG__) THEN *>
TIO.WrStr("*** Aufruf von EigenLib1.HhTri (SV-Version) ***"); TIO.WrLn;
<* END *>
ii:=((dim*(dim+1)) DIV 2);
f:=0.0;
FOR i:=dim TO 1 BY -1 DO
DEC(ii,i);
l:=i-1;
h:=0.0; ik:=ii;
FOR k:=1 TO l DO INC(ik); f:=A[ik]; HD[k]:=f; h:=h+f*f; END;
IF (h <= Small) THEN
ND[i]:=0.0; SqrND[i]:=0.0;
h:=0.0;
ELSE
SqrND[i]:=h;
g:=sqrt(h);
IF (f >= 0.0) THEN g:=-g; END;
ND[i]:=g;
h:=h-f*g;
HD[l]:=f-g;
A[ii+l]:=HD[l];
f:=0.0;
jj:=1;
FOR j:=1 TO l DO
g:=0.0;
jk:=jj;
FOR k:=1 TO l DO
g:=g+A[jk]*HD[k];
IF (k < j) THEN INC(jk); ELSE INC(jk,k); END;
END;
g:=g/h; ND[j]:=g;
f:=f+g*HD[j];
INC(jj,j);
END; (* FOR j *)
hh:=f/(h+h);
jk:=0;
FOR j:=1 TO l DO
f:=HD[j];
g:=ND[j]-hh*f; ND[j]:=g;
FOR k:=1 TO j DO INC(jk); A[jk]:=A[jk]-f*ND[k]-g*HD[k]; END;
END; (* FOR j *)
END; (* IF h *)
HD[i]:=A[ii+i];
A[ii+i]:=h;
END; (* FOR i *)
FOR i:=2 TO dim DO
ND[i-1]:=ND[i];
SqrND[i-1]:=SqrND[i];
END;
ND[dim]:=0.0; SqrND[dim]:=0.0;
END HhTri;
PROCEDURE HhTriBak(VAR EV : MATRIX; (* Zu transformierende EV *)
VAR A : SUPERVEKTOR; (* Housholdertransf. Matrix *)
dim : CARDINAL;
m1,m2 : CARDINAL); (* EV[m1] bis EV[m2] werden trf. *)
VAR i,j,k,l,ii,ik : CARDINAL;
h,s : LONGREAL;
BEGIN
ii:=1;
FOR i:=2 TO dim DO
l:=i-1;
h:=A[ii+i];
IF (h # 0.0) THEN
FOR j:=m1 TO m2 DO
s:=0.0;
ik:=ii;
FOR k:=1 TO l DO INC(ik); s:=s+A[ik]*EV[j,k]; END;
s:=s/h;
ik:=ii;
FOR k:=1 TO l DO INC(ik); EV[j,k]:=EV[j,k]-s*A[ik]; END;
END; (* FOR j *)
END; (* IF h *)
INC(ii,i);
END; (* FOR i *)
END HhTriBak;
PROCEDURE Bisect(VAR HD : VEKTOR; (* Hauptdiagonale der Matrix. *)
ND : VEKTOR; (* Nebendiagonale der Matrix. *)
SqrND : VEKTOR; (* Quadrate von ND. *)
VAR EW : VEKTOR; (* Eigenwerte der Matrix. *)
dim : CARDINAL; (* Dimension der Matrix. *)
m1 : CARDINAL; (* EW[m1] bis EW[m2] werden *)
m2 : CARDINAL; (* berechnet. *)
genau : LONGREAL; (* Genauigkeitsanforderung. *)
AbsFeh : LONGREAL; (* Maschinengenauigkeit. *)
VAR RelFeh : LONGREAL; (* Fehlerabsch"atzung. *)
VAR z : CARDINAL; (* Anzahl Bisectionsschritte. *)
VAR Liste : ListenZeiger); (* Nicht Implementiert. *)
VAR i,j,k,MaxIter,Iter : CARDINAL;
q,x1,xu,xo : LONGREAL;
h,xmin,xmax : LONGREAL;
wu : VEKTOR;
BEGIN
FOR i:=1 TO dim DO EW[i]:=0.0; END;
IF (m2 < m1) THEN i:=m2; m2:=m1; m1:=i; END;
IF (m1 = 0) OR (m2 = 0) THEN RETURN END;
MaxIter:=dim*dim*MaxCard((m2-m1),4);
IF (genau < MachEps) THEN genau:=MachEps; END;
Liste:=NIL;
Fehler:=FALSE;
FOR i:=dim TO 2 BY -1 DO (* Speicher ND so um, da\3 ND[1]:=0.0 *)
SqrND[i]:=SqrND[i-1];
ND[i]:=ND[i-1];
END;
SqrND[1]:=0.0;
ND[1]:=0.0;
xmin:=HD[dim]-ABS(ND[dim]);
xmax:=HD[dim]+ABS(ND[dim]);
FOR i:=dim-1 TO 1 BY -1 DO
h:=ABS(ND[i])+ABS(ND[i+1]);
IF (HD[i]+h > xmax) THEN xmax:=HD[i]+h; END;
IF (HD[i]-h < xmin) THEN xmin:=HD[i]-h; END;
END; (* FOR i *)
IF ((xmin+xmax) > 0.0) THEN
RelFeh:=AbsFeh*xmax;
ELSE
RelFeh:=-AbsFeh*xmin;
END;
IF (genau < 0.0) THEN genau:=RelFeh; END;
RelFeh:=0.5*genau+7.0*RelFeh;
xo:=xmax;
FOR i:=m1 TO m2 DO
EW[i]:=xmax;
wu[i]:=xmin;
END;
z:=0;
FOR k:=m2 TO m1 BY -1 DO
xu:=xmin;
i:=k;
LOOP
IF (xu < wu[i]) THEN xu:=wu[i]; EXIT END;
IF (i = m1) THEN EXIT END;
DEC(i);
END;
IF (xo > EW[k]) THEN xo:=EW[k]; END;
Iter:=0;
LOOP
INC(Iter);
IF (xo-xu < 2.0*AbsFeh*(ABS(xu)+ABS(xo))) THEN EXIT END;
IF (Iter > MaxIter) THEN
Fehler:=TRUE;
Fehlerflag:='MaxIter "uberschritten (Bisect)';
Errors.ErrOut(Fehlerflag);
EXIT;
END;
x1:=0.5*(xu+xo);
INC(z);
j:=0;
q:=1.0;
FOR i:=1 TO dim DO
IF (q # 0.0) THEN
q:=HD[i]-x1-SqrND[i]/q;
ELSE
q:=HD[i]-x1-ABS(ND[i]/AbsFeh);
END;
IF (q < 0.0) THEN INC(j); END;
END; (* FOR i *)
IF (j < k) THEN
IF (j < m1) THEN
xu:=x1;
wu[m1]:=x1;
ELSE
xu:=x1;
wu[j+1]:=x1;
IF (EW[j] > x1) THEN EW[j]:=x1; END;
END;
ELSE
xo:=x1;
END;
END; (* LOOP *)
EW[k]:=0.5*(xo+xu);
END; (* FOR k *)
END Bisect;
PROCEDURE ListenEintrag(VAR UngrObgr : ListenZeiger;
u,o : CARDINAL);
(*------------------------------------------------------------*)
(* Tr"agt eventuelle Blockungen der Matrix in eine verkettete *)
(* Liste zum Gebrauch bei der Eigenvektorberechnung ein. *)
(*------------------------------------------------------------*)
VAR InterimZeiger : ListenZeiger;
BEGIN
NEW(InterimZeiger);
WITH InterimZeiger^ DO
Ungr:=u;
Obgr:=o;
next:=UngrObgr;
END;
UngrObgr:=InterimZeiger;
END ListenEintrag;
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;
PROCEDURE TriQR( HD,ND : VEKTOR; (* <== Trilineare Matrix *)
VAR EW : VEKTOR; (* ==> Eigenwerte von A *)
dim : CARDINAL; (* Dimension von A *)
VAR UngrObgr : ListenZeiger;
genau : LONGREAL);
CONST Eps = 10.0*MachEps;
Min = Small;
PROCEDURE TransQR(VAR EW : VEKTOR;
VAR ND : VEKTOR;
u,o : CARDINAL;
VAR genau : LONGREAL);
(*----------------------------------------------------------------*)
(* Die Procedure TransQR erledigt die eigendliche QR - *)
(* Transformation der im Hauptteil zerlegten Matrix (HD,ND). *)
(*----------------------------------------------------------------*)
CONST MaxIter = 60;
VAR sin,cos,sincos,sqrsin,sqrcos,omega,
interim,Delta,x,y,shift,EWo : LONGREAL;
Iter,p : CARDINAL;
BEGIN
Fehler:=FALSE;
WHILE (o > u) DO
Iter:=0;
LOOP
IF (Iter <= MaxIter) THEN INC(Iter) ELSE
Fehler:=TRUE; Fehlerflag:='MaxIter "uberschritten (TriQR)';
Errors.ErrOut(Fehlerflag); RETURN;
END;
Delta:=0.5*(EW[o-1]-EW[o]);
IF (ABS(Delta) < Eps*ABS(0.5*(EW[o-1]+EW[o]))) THEN
shift:=EW[o];
IF (ABS(shift) < Eps) THEN shift:=EW[o-1]; END; (* ?? *)
IF (shift = EW[u]) THEN shift:=(1.0 + 2.0*Eps)*shift; END; (*??*)
ELSE
shift:=EW[o]+Delta;
IF (Delta >= 0.0) THEN
shift:=shift - sqrt(Delta*Delta + ND[o-1]*ND[o-1]);
ELSE
shift:=shift + sqrt(Delta*Delta + ND[o-1]*ND[o-1]);
END;
END;
x:=EW[u]-shift; y:=ND[u]; EWo:=EW[o];
FOR p:=u TO o-1 DO
IF (ABS(x) <= genau*ABS(y)) THEN
omega:=-y;
sqrcos:=0.0; cos:=0.0;
sqrsin:=1.0; sin:=1.0;
sincos:=0.0;
ELSE (* Bereche Drehparameter *)
(* omega:=sqrt(x*x + y*y); Nicht allzu genau ! *)
omega:=Pythag(x,y);
IF (omega < Small*Small) THEN
Errors.Fehlerflag:=" Unterlauf in EigenLib1.TriQR !";
Fehler:=TRUE; RETURN;
END;
interim := 1.0 / omega;
cos:= x * interim; sqrcos:=cos*cos;
sin:=-y * interim; sqrsin:=sin*sin;
sincos:=sin*cos;
END;
Delta:=EW[p]-EW[p+1];
interim:=2.0*sincos*ND[p] + Delta*sqrsin;
EW[p] :=EW[p] - interim;
EW[p+1]:=EW[p+1] + interim;
ND[p]:=Delta*sincos + (sqrcos-sqrsin)*ND[p]; x:=ND[p];
IF (p > u) THEN ND[p-1]:=omega; END;
IF (p < o-1) THEN y:=-sin*ND[p+1]; ND[p+1]:=cos*ND[p+1]; END;
END;(* FOR p *)
(* Eigenwert isoliert ? *)
IF (ABS(ND[o-1]) < genau*AbsMin(EW[o],EW[o-1])) THEN EXIT END;
(* Konvergenz ? *)
IF (ABS(EWo - EW[o]) < Eps*ABS(EW[o])) THEN EXIT END;
(* Eigenwert numerisch Null oder fast Null ? *)
IF (ABS(EWo) + ABS(EW[o]) < 20.0*Eps) AND (Iter >= 2) THEN EXIT END;
END; (* LOOP *)
DEC(o);
END; (* WHILE ((o > u) *)
END TransQR;
VAR p,u : CARDINAL;
BEGIN (* TriQR *)
<* IF (__DEBUG__) THEN *>
TIO.WrLn();
TIO.WrStr("Aufruf von TriQR ... ");
TIO.WrLn();
TIO.WrLn();
<* END *>
FOR p:=1 TO dim DO (* Abfrage verhindert Unterlauf. *)
EW[p]:=HD[p]; IF (ABS(EW[p]) <= Min) THEN EW[p]:=Min; END;
END;
IF (genau < Eps) THEN genau:=Eps; END;
UngrObgr:=NIL; p:=0; u:=1;
REPEAT (* Zerlege (HD,ND) entsprechend der Blockung *)
INC(p);
IF (ABS(ND[p]) <= MachEps*AbsMin(HD[p],HD[p+1])) THEN
ListenEintrag(UngrObgr,u,p);
IF (u # p) THEN TransQR(EW,ND,u,p,genau); END;
IF Fehler THEN RETURN END;
u:=p+1;
END;
UNTIL (p = dim-1);
ListenEintrag(UngrObgr,u,dim);
IF (u < dim) THEN TransQR(EW,ND,u,dim,genau); END;
END TriQR;
PROCEDURE TriQL( HD,ND : VEKTOR; (* <== Trilineare Matrix *)
VAR EW : VEKTOR; (* ==> Eigenwerte von A *)
dim : CARDINAL; (* Dimension von A *)
VAR UngrObgr : ListenZeiger;
genau : LONGREAL);
CONST Eps = MachEps;
PROCEDURE TransQL(VAR EW : VEKTOR;
VAR ND : VEKTOR;
dim : CARDINAL;
u,o : CARDINAL;
VAR genau : LONGREAL);
(*-----------------------------------------------------------*)
(* Die Procedure TransQL erledigt die eigendliche QL - *)
(* Transformation der im Haupteile zerlegten Matrix (HD,ND). *)
(*-----------------------------------------------------------*)
CONST MaxIter = 30;
VAR Iter,m,l,i : CARDINAL;
b,c,f,g,p,r,s : LONGREAL;
BEGIN
FOR l:=u TO o DO
Iter:=0;
LOOP
IF (Iter <= MaxIter) THEN INC(Iter) ELSE
Fehler:=TRUE; Fehlerflag:='MaxIter Ueberschritten (TriQL).';
Errors.ErrOut(Fehlerflag);
RETURN;
END;
m:=l-1;
REPEAT (* Suche kleine Subdiagonalelemente *)
INC(m);
UNTIL (m = dim) 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:=sqrt(g*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:=(EW[i]-g)*s+2.0*c*b;
p:=s*r;
EW[i+1]:=g+p;
g:=c*r-b;
END; (* FOR i *)
EW[l]:=EW[l]-p;
ND[l]:=g;
ND[m]:=0.0;
END; (* LOOP *)
END; (* FOR l *)
END TransQL;
VAR p,u : CARDINAL;
BEGIN
FOR p:=1 TO dim DO EW[p]:=HD[p]; END;
IF (genau < Eps) THEN genau:=Eps; END;
UngrObgr:=NIL;
p:=0; u:=1;
REPEAT (* Zerlege (HD,ND) entsprechend der Blockung *)
INC(p);
IF (ABS(ND[p]) <= MachEps*AbsMin(HD[p],HD[p+1])) THEN
ListenEintrag(UngrObgr,u,p);
IF (u # p) THEN TransQL(EW,ND,dim,u,p,genau); END;
IF Fehler THEN RETURN END;
u:=p+1;
END;
UNTIL (p = dim-1);
ListenEintrag(UngrObgr,u,dim);
IF (u < dim) THEN TransQL(EW,ND,dim,u,dim,genau); END;
END TriQL;
PROCEDURE EwEvQR(VAR HD : VEKTOR;
ND : VEKTOR;
VAR EW : VEKTOR;
VAR EV : MATRIX;
dim : CARDINAL;
genau : LONGREAL);
CONST Eps = MachEps;
Min = Small*Small;
PROCEDURE TransQR(VAR EW : VEKTOR;
VAR EV : MATRIX;
VAR ND : VEKTOR;
dim : CARDINAL;
u,o : CARDINAL;
VAR genau : LONGREAL);
(*----------------------------------------------------------*)
(* Die Procedure Transform erledigt die eigendliche *)
(* QR - Transformation der im Hauptteil zerlegten Matrix A. *)
(*----------------------------------------------------------*)
CONST MaxIter = 30;
VAR sin,cos,sincos,sqrsin,sqrcos,omega,
interim,Delta,x,y,shift,EVqk,EVpk,EWo : LONGREAL;
Iter,p,q,k : CARDINAL;
BEGIN
Fehler:=FALSE;
WHILE (o > u) DO
Iter:=0;
LOOP
IF (Iter < MaxIter) THEN INC(Iter) ELSE
Fehler:=TRUE;
Fehlerflag:='MaxIter "uberschritten (EwEvQR)';
Errors.ErrOut(Fehlerflag);
RETURN;
END;
IF (ABS(ND[o-1]) < genau*AbsMin(EW[o],EW[o-1])) THEN EXIT END;
Delta:=0.5*(EW[o-1]-EW[o]);
IF ((ABS(Delta) > Eps*ABS(0.5*(EW[o-1]+EW[o]))) AND
(ABS(EW[o]) > Eps)) AND (* ?? *)
(EW[o] # EW[u])
THEN (* Berechen Shift *)
shift:=EW[o];
ELSE
shift:=EW[o]+Delta;
IF (Delta >= 0.0) THEN
shift:=shift - sqrt(Delta*Delta+ND[o-1]*ND[o-1]);
ELSE
shift:=shift + sqrt(Delta*Delta+ND[o-1]*ND[o-1]);
END;
END;
x:=EW[u]-shift;
y:=ND[u];
EWo:=EW[o];
FOR p:=u TO o-1 DO
IF (ABS(x) <= genau*ABS(y)) THEN (* ?? *)
omega:=-y;
sqrcos:=0.0; cos:=0.0;
sqrsin:=1.0; sin:=1.0;
sincos:=0.0;
ELSE
omega:=sqrt(x*x+y*y);
cos:= x/omega;
sin:=-y/omega;
sincos:=sin*cos;
sqrsin:=sin*sin;
sqrcos:=cos*cos;
END;
Delta:=EW[p]-EW[p+1];
interim:=2.0*sincos*ND[p] + Delta*sqrsin;
EW[p] :=EW[p] - interim;
EW[p+1]:=EW[p+1] + interim;
ND[p]:=Delta*sincos + (sqrcos-sqrsin)*ND[p];
x:=ND[p];
IF (p > u) THEN ND[p-1]:=omega; END;
IF (p < o-1) THEN
y:=-sin*ND[p+1];
ND[p+1]:=cos*ND[p+1];
END;
q:=p+1;
FOR k:=1 TO dim DO (* Berechne Eigenvektoren *)
EVqk:=EV[q,k]; EVpk:=EV[p,k];
EV[q,k]:=sin*EVpk + cos*EVqk;
EV[p,k]:=cos*EVpk - sin*EVqk;
END; (* FOR *)
END;(* FOR p *)
IF (EWo = EW[o]) THEN EXIT END;
END; (* LOOP *)
DEC(o);
END; (* WHILE ((o > u) *)
END TransQR;
VAR p,q,u : CARDINAL;
BEGIN
FOR p:=1 TO dim DO (* Abfrage verhindert Unterlauf. *)
EW[p]:=HD[p]; IF (ABS(EW[p]) <= Min) THEN EW[p]:=Min; END;
FOR q:=1 TO p-1 DO EV[p,q]:=0.0; EV[q,p]:=0.0; END;
EV[p,p]:=1.0;
END;
FOR p:=1 TO dim DO EW[p]:=HD[p]; END;
IF (genau < Eps) THEN genau:=Eps; END;
p:=0; u:=1;
REPEAT (* Zerlege A entsprechend der Blockung *)
INC(p);
IF (ABS(ND[p]) <= MachEps*AbsMin(HD[p],HD[p+1])) THEN
TransQR(EW,EV,ND,dim,u,p,genau);
u:=p+1;
END;
UNTIL (p = dim-1);
IF (u < dim) THEN
TransQR(EW,EV,ND,dim,u,dim,genau);
END;
END EwEvQR;
PROCEDURE BerechneEV(HD,ND : VEKTOR; (* Haupt- u. Nebendiagonale v. A *)
VAR EWi : LONGREAL; (* Sch"atzwert d. i. Eigenwerts *)
VAR EVi : VEKTOR; (* Eigenvektor *)
i : CARDINAL; (* Index von EWi,EVi *)
genau : LONGREAL; (* selbsterkl"arend *)
MaxIter : CARDINAL; (* MaximalZahl der Iterationen *)
UngrObgr : ListenZeiger;
dim : CARDINAL); (* Dimension der Matrix A *)
PROCEDURE BestimmeUngrObgr(VAR Untergr : CARDINAL;
VAR Obergr : CARDINAL;
UngrObgr : ListenZeiger;
i : CARDINAL;
dim : CARDINAL);
(*--------------------------------------------------------------*)
(* Berechnet aus den in QR in die Liste UngrObger eingetragenen *)
(* Blockungen der trilinearisierten zu diagonalisierenden *)
(* Matrix A die Werte fuer Untergrenzen und Obergrenzen zur"uck *)
(*--------------------------------------------------------------*)
VAR p : ListenZeiger; (* Hilfszeiger *)
BEGIN
p:=UngrObgr; (* p zeigt auf erstes Element *)
WHILE (p # NIL)
AND NOT ((p^.Ungr <= i) AND (p^.Obgr >= i)) DO
p:=p^.next;
END;
IF (p = NIL) THEN
Untergr:=1;
Obergr:=dim;
ELSE
WITH p^ DO
Untergr:=Ungr;
Obergr:=Obgr;
END;
END;
END BestimmeUngrObgr;
PROCEDURE GaussTriLinear(VAR L,HD,ND,T : VEKTOR;
VAR X : VEKTOR;
VAR C : VEKTOR;
VAR II : CARDVEKTOR;
u,o : CARDINAL);
(*--------------------------------------------------------------*)
(* L"o\3t das Gleichungssyste A*X = C f"ur eine in HD,ND *)
(* gespeicherte trilineare symmetrische Matrix. *)
(* HD : Hauptdiagonale von A *)
(* ND : Nebendiagonale von A (oberhalb von HD) *)
(* T : Nebendiagonale von ND (oberhalb von ND) *)
(* L : Nebendiagonale von A (unterhalb von HD) *)
(*--------------------------------------------------------------*)
VAR alpha,beta,z : LONGREAL;
i,ip : CARDINAL; (* ip = i+1 *)
BEGIN
Fehler:=FALSE;
FOR i:=u TO o DO T[i]:=0.0; II[i]:=0; END;
ip:=u;
FOR i:=u TO o-1 DO
INC(ip);
alpha := ABS(HD[i]) + ABS(ND[i]);
beta := ABS(L[i]) + ABS(HD[ip]) + ABS(ND[ip]);
IF (alpha < Small) THEN alpha:=Small; END;
IF (beta < Small) THEN beta :=Small; END;
IF (ABS(HD[i])/alpha < ABS(L[i])/beta) THEN (*??*)
II[i]:=i; (* Vertausche Elemente *)
z := HD[i]; HD[i] := L[i]; L[i] := z;
z := HD[ip]; HD[ip] := ND[i]; ND[i] := z;
z := ND[ip]; ND[ip] := T[i]; T[i] := z;
z := C[ip]; C[ip] := C[i]; C[i] := z;
END; (* IF *)
IF (HD[i] = 0.0) THEN
Fehler:=TRUE;
Fehlerflag:='Matrix singulaer (GaussTriLinear)';
RETURN;
END;
z := L[i] / HD[i];
HD[ip] := HD[ip] - z*ND[i];
ND[ip] := ND[ip] - z*T[i];
C[ip] := C[ip] - z*C[i];
L[i] := z;
END; (* FOR i *);
FOR i:=u TO o DO (* Test ! *)
IF (HD[i] = 0.0) THEN
Fehler:=TRUE;
Fehlerflag:='Diagonalelement Null (GaussTriDiagonal)';
RETURN;
END;
END;
X[o] := C[o] / HD[o];
X[o-1] := (C[o-1]-ND[o-1]*X[o]) / HD[o-1];
FOR i:=o-2 TO u BY -1 DO
X[i]:=(C[i]-T[i]*X[i+2]-ND[i]*X[i+1]) / HD[i];
END;
END GaussTriLinear;
PROCEDURE BerechneNeu( L,HD,ND,T : VEKTOR;
VAR X : VEKTOR;
VAR C : VEKTOR;
II : CARDVEKTOR;
u,o : CARDINAL);
(*--------------------------------------------------------------*)
(* Berechnet durch Vorwaerts- und Rueckwaertseinsetzen in die *)
(* LR-Zerlegte von A neuen Vektor der i-ten Wielandt-Iteration *)
(* (i > 1) *)
(*--------------------------------------------------------------*)
VAR i,ip : CARDINAL;
Z : LONGREAL;
BEGIN
ip:=u;
FOR i:=u TO o-1 DO
INC(ip);
IF (II[i] = i) THEN Z:=C[i]; C[i]:=C[ip]; C[ip]:=Z; END;
C[ip]:=C[ip]-L[i]*C[i];
END;
X[o]:=C[o] / HD[o];
X[o-1]:=(C[o-1]-ND[o-1]*X[o]) / HD[o-1];
FOR i:=o-2 TO u BY -1 DO
X[i]:=(C[i]-T[i]*X[i+2]-ND[i]*X[i+1]) / HD[i];
END;
END BerechneNeu;
VAR Norm,TestNorm : LONGREAL;
p,Iter,u,o : CARDINAL;
II : CARDVEKTOR;
L,T : VEKTOR;
BEGIN (* BerechneEV *)
BestimmeUngrObgr(u,o,UngrObgr,i,dim);
IF (u = o) THEN
FOR p:=1 TO dim DO EVi[p]:=0.0; END;
EVi[o]:=1.0;
ELSE
IF (MaxIter < 1) THEN MaxIter:=dim; END;
IF (MaxIter < 6) THEN MaxIter:=6; END;
IF (genau <= 0.0) THEN genau:=1.0E-10; END;
FOR p:=1 TO u-1 DO EVi[p]:=0.0; END;
FOR p:=o+1 TO dim DO EVi[p]:=0.0; END;
Norm:=1.0 / sqrt((VAL(LONGREAL,o-u))+1.0);
FOR p:=u TO o DO
HD[p]:=HD[p]-EWi;
IF (ABS(HD[p]) <= Small) THEN HD[p]:=Small; END;
EVi[p]:=Norm;
END;
L:=ND;
GaussTriLinear(L,HD,ND,T,EVi,EVi,II,u,o);
IF Fehler THEN RETURN; END; (* !!! *)
Norm:=0.0;
FOR p:=u TO o DO (* Einfache Normierung ! *)
IF (ABS(EVi[p]) > Norm) THEN Norm:=ABS(EVi[p]); END;
END;
Norm := 1.0 / Norm;
Iter:=1; (* Berechnung des Startvektors in GaussTriLinear *)
REPEAT
INC(Iter);
FOR p:=u TO o DO EVi[p]:=EVi[p]*Norm; END; (* Norm. Vektor *)
BerechneNeu(L,HD,ND,T,EVi,EVi,II,u,o);
TestNorm:=Norm;
Norm:=0.0;
FOR p:=u TO o DO (* Einfache Normierung ! *)
IF (ABS(EVi[p]) > Norm) THEN Norm:=ABS(EVi[p]); END;
END;
Norm := 1.0 / Norm;
UNTIL (ABS(TestNorm-Norm) < genau*Norm) OR (Iter >= MaxIter);
Norm:=0.0;
FOR p:=u TO o DO Norm:=Norm + EVi[p]*EVi[p]; END;
Norm := 1.0 / sqrt(Norm);
FOR p:=u TO o DO EVi[p]:=EVi[p]*Norm; END; (* Normiere Vektor *)
END; (* IF u = o *)
END BerechneEV;
PROCEDURE RSymEwEv(VAR EV : MATRIX;
VAR EW : VEKTOR;
VAR A : SUPERVEKTOR;
dim : CARDINAL;
maxEW : CARDINAL;
maxEV : CARDINAL;
Ortho : CARDINAL;
genau : LONGREAL;
ord : INTEGER);
VAR HD,ND,sqrND : VEKTOR;
i : CARDINAL;
EWi,x : LONGREAL;
UngrObgr : ListenZeiger;
p : ListenZeiger;
II : CARDVEKTOR;
BEGIN
UngrObgr:=NIL;
IF (maxEV > dim) THEN maxEV:=dim; END;
IF (maxEW > dim) THEN maxEW:=dim; END;
IF (maxEV > maxEW) THEN maxEW:=maxEV; END;
IF (maxEW = 0) THEN RETURN; END;
IF (Ortho > 1) THEN Ortho:=0; END;
HhTri(A,HD,ND,sqrND,dim);
IF (maxEW = dim) AND (maxEV = dim) AND (Ortho = 1) THEN
EwEvQR(HD,ND,EW,EV,dim,genau);
HhTriBak(EV,A,dim,1,dim);
SortEwEv(EW,EV,dim,dim,ord);
RETURN; (* !!! *)
END;
IF ((dim DIV maxEW) >= 4) AND (ord = -1) THEN
Bisect(HD,ND,sqrND,EW,dim,1,maxEW,genau,MachEps,x,i,UngrObgr);
FOR i:=1 TO maxEV DO II[i]:=i; END;
IF Fehler THEN
TriQR(HD,ND,EW,dim,UngrObgr,genau);
QuickSort(EW,II,dim,ord);
END;
ELSE
TriQR(HD,ND,EW,dim,UngrObgr,genau);
IF Fehler THEN
Bisect(HD,ND,sqrND,EW,dim,1,maxEW,genau,MachEps,x,i,UngrObgr);
FOR i:=1 TO maxEV DO II[i]:=i; END;
ELSE
QuickSort(EW,II,dim,ord);
END;
END;
FOR i:=1 TO maxEV DO
BerechneEV(HD,ND,EW[i],EV[i],II[i],genau,dim,UngrObgr,dim);
IF Fehler THEN (* Neuer Versuch mit schlechterem Eigenwert *)
Errors.ErrOut('Eigenvektoren unsicher (BerechneEV)');
IF (EW[i] = 0.0) THEN
EWi:=1.0E-14
ELSE
EWi:=(1.0+1.0E-14)*EW[i];
END;
BerechneEV(HD,ND,EWi,EV[i],II[i],genau,dim,UngrObgr,dim);
IF Fehler THEN Errors.FatalError(Fehlerflag); END;
END;
END;
IF (maxEV < dim) THEN
(* Initialisiere nicht berechnete EV *)
FOR i:=1 TO dim DO EV[maxEV+1,i]:=0.0; END;
FOR i:=maxEV+2 TO dim DO EV[i]:=EV[maxEV+1]; END;
END;
WHILE (UngrObgr # NIL) DO (* Gebe Speicher wieder frei *)
p:=UngrObgr;
UngrObgr:=p^.next;
DISPOSE(p);
END;
IF (maxEV > 0) THEN HhTriBak(EV,A,dim,1,maxEV); END;
END RSymEwEv;
MODULE LOEWDIN;
IMPORT ADR;
IMPORT ALLOCATE,DEALLOCATE;
IMPORT Errors;
IMPORT StdErr,FIO2;
IMPORT MATRIX,VEKTOR,SUPERVEKTOR;
IMPORT MachEps;
IMPORT NormMat,MatSvProdNN,MatMatProd,MatToSv,Stuerz;
IMPORT Invert;
IMPORT RSymEwEv,Jacobi;
IMPORT sqrt;
<* IF (__XDS__) THEN *>
IMPORT Loewdin,Transform,TransfSV; (* EXPORT *)
<* ELSE *>
EXPORT Loewdin,Transform,TransfSV;
<* END *>
(*----------------------------------------------------------------*)
(* Die Matrix V dient der Zwischenspeicherung der *)
(* Transformationsmatrix {S^{-1/2}}^t von Aufruf zu Aufruf. *)
(*----------------------------------------------------------------*)
VAR V : MATRIX;
PROCEDURE Loewdin(VAR H : SUPERVEKTOR;
VAR S : SUPERVEKTOR; (* positiv define Matrix *)
VAR EW : VEKTOR; (* Eigenwerte *)
VAR EV : MATRIX; (* Eigenvektoren *)
dim : CARDINAL; (* Dimension des Problems *)
maxEV : CARDINAL; (* Anzahl d. Eigenvektoren *)
genau : LONGREAL; (* geforderte Genauigkeit *)
ord : INTEGER;
Neu : CARDINAL;
Norm : CARDINAL;
Transform : CARDINAL);
VAR i,j,k,ii,ij,ik : CARDINAL;
Sum,EWi,x : LONGREAL;
Vek : VEKTOR; (* Hilfsvektor *)
BEGIN
IF (Neu > 1) OR (Norm > 1) OR (Transform > 2) THEN
Errors.ErrOut('Steuerparameter falsch gesetzt (Loewdin)');
IF (Neu > 1) THEN Neu:=1; END;
IF (Norm > 1) THEN Norm:=0; END;
IF (Transform > 1) THEN Transform:=1; END;
END;
IF (Neu > 0) THEN
IF (ADR(H) = ADR(S)) THEN
Errors.FatalError('H und S identisch und (Neu = TRUE) (Loewdin)');
END;
Jacobi(V,EW,S,dim,0,MachEps,0); (* Jacobi numerisch am stabilsten. *)
(* Auf V liegen die Eigenvektoren von S *)
FOR i:=1 TO dim DO
IF (EW[i] <= 0.0) THEN
Errors.WriteString("Nichtpositiver Eigenwert EW(S)[");
Errors.WriteCard(i); Errors.WriteString("] = ");
FIO2.WrLngReal(StdErr,EW[i],18,-9);
Errors.FatalError('S-Matrix nicht positiv definit (Loewdin)');
END;
EWi := 1.0 / sqrt(EW[i]);
FOR j:=1 TO dim DO V[i,j]:=EWi*V[i,j]; END; (* V = {S^{-1/2}}^t *)
END; (* FOR i *)
END; (* IF Neu *)
ii:=0;
FOR i:=1 TO dim DO
FOR j:=1 TO dim DO
Sum:=0.0; ik:=ii;
FOR k:=1 TO i DO INC(ik); Sum:=Sum + H[ik]*V[j,k]; END;
FOR k:=i+1 TO dim DO INC(ik,k-1); Sum:=Sum + H[ik]*V[j,k]; END;
EV[i,j]:=Sum; (* EV = S^{-1/2} H *)
END; (* FOR j *)
INC(ii,i);
END; (* FOR i *)
ij:=0;
FOR i:=1 TO dim DO
FOR j:=1 TO i DO
Sum:=0.0;
FOR k:=1 TO dim DO Sum:=Sum + V[i,k]*EV[k,j]; END;
INC(ij);
H[ij]:=Sum; (* H' = S^{-1/2} H {S^{-1/2}}^t *)
END;
END;
IF (Transform < 2) THEN (* Nur H' bilden ? *)
(* Jacobi(EV,EW,H,dim,0,genau,ord); *)
RSymEwEv(EV,EW,H,dim,maxEV,maxEV,0,genau,ord);
IF (Transform = 1) THEN
FOR i:=1 TO maxEV DO (* Berechne EV = EV*V *)
FOR j:=1 TO dim DO Vek[j]:=0.0; END;
FOR j:=1 TO dim DO
x:=EV[i,j];
IF (x # 0.0) THEN
FOR k:=1 TO dim DO Vek[k]:=Vek[k] + x*V[j,k]; END;
END;
END;
FOR j:=1 TO dim DO EV[i,j]:=Vek[j]; END;
END; (* FOR i *)
END; (* IF *)
IF (Norm = 1) THEN NormMat(EV,dim,maxEV,FALSE); END;
END; (* IF Transform < 2 *)
END Loewdin;
PROCEDURE Transform(VAR EV : MATRIX;
dim : CARDINAL;
Form : CARDINAL);
CONST FehlStr = 'kein Freispeicher vorhanden (EigenLib1.Transform) !';
VAR i,j,k : CARDINAL;
Sum,x : LONGREAL;
Vinv : POINTER TO MATRIX; (* Inverse zu V *)
Vek : VEKTOR;
iFehl : INTEGER;
BEGIN
IF (Form > 2) THEN Form:=0; END;
IF (Form = 0) THEN
FOR i:=1 TO dim DO
FOR j:=1 TO dim DO Vek[j]:=0.0; END;
FOR j:=1 TO dim DO
x:=EV[i,j];
IF (x # 0.0) THEN
FOR k:=1 TO dim DO Vek[k]:=Vek[k] + x*V[j,k]; END;
END;
END;
FOR j:=1 TO dim DO EV[i,j]:=Vek[j]; END;
END;
ELSIF (Form = 1) THEN (* Erweiter von rechts mit V^{-1}*)
NEW(Vinv); IF (Vinv = NIL) THEN Errors.FatalError(FehlStr); END;
Invert(Vinv^,V,x,dim,iFehl);
FOR i:=1 TO dim DO (* Berechne EV = EV * V^{-1} *)
FOR j:=1 TO dim DO
Sum:=0.0;
FOR k:=1 TO dim DO Sum:=Sum + EV[i,k]*Vinv^[k,j]; END;
Vek[j]:=Sum;
END; (* FOR j *)
FOR j:=1 TO dim DO EV[i,j]:=Vek[j]; END;
END; (* FOR i *)
NormMat(EV,dim,dim,FALSE);
DISPOSE(Vinv);
ELSE (* Form = 2 , Berechne EV = {S^{-1/2}}^+ EV S^{-1/2} *)
(*
* MatMatProd(EV,EV,V,dim);
* Stuerz(V,dim);
* MatMatProd(EV,V,EV,dim);
* Stuerz(V,dim);
*)
FOR i:=1 TO dim DO
FOR j:=1 TO dim DO
Sum:=0.0;
FOR k:=1 TO dim DO Sum:=Sum + V[i,k]*EV[k,j]; END;
Vek[j]:=Sum;
END;
FOR j:=1 TO dim DO EV[i,j]:=Vek[j]; END;
END;
FOR i:=1 TO dim DO
FOR j:=1 TO dim DO
Sum:=0.0;
FOR k:=1 TO dim DO Sum:=Sum + EV[i,k]*V[j,k]; END;
Vek[j]:=Sum;
END;
FOR j:=1 TO dim DO EV[i,j]:=Vek[j]; END;
END;
END; (* IF Form *)
END Transform;
PROCEDURE TransfSV(VAR SV : SUPERVEKTOR;
dim : CARDINAL;
Form : CARDINAL);
CONST FehlStr = 'kein Freispeicher vorhanden (EigenLib1.TransfSV) !';
VAR Det : LONGREAL;
Z,Vinv : POINTER TO MATRIX;
iFehl : INTEGER;
BEGIN
NEW(Z); IF (Z = NIL) THEN Errors.FatalError(FehlStr); END;
IF (Form = 0) THEN
MatSvProdNN(dim,Z^,V,SV);
Stuerz(V,dim);
MatMatProd(Z^,Z^,V,dim);
Stuerz(V,dim);
MatToSv(Z^,SV,dim,0);
ELSE
NEW(Vinv); IF (Vinv = NIL) THEN Errors.FatalError(FehlStr); END;
Invert(Vinv^,V,Det,dim,iFehl); (* Vinv = S^{1/2} = V^{-1} *)
MatSvProdNN(dim,Z^,Vinv^,SV);
Stuerz(Vinv^,dim);
MatMatProd(Z^,Z^,Vinv^,dim);
Stuerz(V,dim);
MatToSv(Z^,SV,dim,0);
DISPOSE(Vinv);
END;
DISPOSE(Z);
END TransfSV;
END LOEWDIN;
PROCEDURE Reduce1(VAR A,B : ARRAY OF ARRAY OF LONGREAL;
VAR dL : ARRAY OF LONGREAL;
dim : CARDINAL;
transB : BOOLEAN);
VAR i,j,k : CARDINAL;
x,y : LONGREAL;
BEGIN
IF transB THEN
y:=1.0; (* Wg. Compilerwarnung *)
FOR i:=0 TO dim-1 DO
FOR j:=i TO dim-1 DO
x:=B[i,j];
IF (i > 0) THEN
<* IF (__BLAS__) THEN *>
x:=x - ddot(i,B[i][0],1,B[j][0],1);
<* ELSE *>
FOR k:=i-1 TO 0 BY -1 DO x:=x - B[i,k]*B[j,k]; END;
<* END *>
END;
IF (i = j) THEN
IF (x <= 0.0) THEN
Fehler:=TRUE;
Fehlerflag:="B nicht positiv definit (EisPack.Reduce1) !";
Errors.ErrOut(Fehlerflag);
RETURN;
END;
y:=sqrt(x); dL[i]:=y;
ELSE
B[j,i]:= x / y;
END;
END; (* FOR j *)
END; (* FOR i *)
END; (* IF *)
FOR i:=0 TO dim-1 DO (* L has been formed in the array b *)
y := dL[i];
FOR j:=i TO dim-1 DO
x := A[i,j];
IF (i > 0) THEN
<* IF (__BLAS__) THEN *>
x:=x - ddot(i,B[i][0],1,A[j][0],1);
<* ELSE *>
FOR k:=i-1 TO 0 BY -1 DO x:=x - B[i,k]*A[j,k]; END;
<* END *>
END;
A[j,i] := x / y;
END;
END;
(* The transpose of the upper triangel of inv(L)xA has *)
(* been formed in the lower triangel of the array A *)
FOR j:=0 TO dim-1 DO
FOR i:=j TO dim-1 DO
x := A[i,j];
IF (i > 0) THEN
FOR k:=i-1 TO j BY -1 DO x:=x - A[k,j]*B[i,k]; END;
END;
IF (j > 0) THEN
<* IF (__BLAS__) THEN *>
x:=x - ddot(j,A[j][0],1,B[i][0],1);
<* ELSE *>
FOR k:=j-1 TO 0 BY -1 DO x:=x - A[j,k]*B[i,k]; END;
<* END *>
END;
A[i,j] := x / dL[i];
END;
END;
END Reduce1;
PROCEDURE ReBakA(VAR EV : ARRAY OF ARRAY OF LONGREAL;
dim : CARDINAL;
m1,m2 : CARDINAL;
VAR B : ARRAY OF ARRAY OF LONGREAL;
VAR dL : ARRAY OF LONGREAL);
VAR i,j,k : CARDINAL;
x : LONGREAL;
BEGIN
FOR j:=m1-1 TO m2-1 DO
FOR i:=dim-1 TO 0 BY -1 DO
x := EV[j,i];
FOR k:=i+1 TO dim-1 DO x:=x - B[k,i]*EV[j,k]; END;
EV[j,i]:= x / dL[i];
END;
END;
END ReBakA;
END EigenLib1.