DEFINITION MODULE EigenLib1;
(*------------------------------------------------------------------------*)
(* Stellt Routinen zur Eigenwert und Vektorenberechnung symmetrischer *)
(* reeller Matritzen zur Verf\"ugung. *)
(* Routines for calculating eigensystems of symmetric real matrices *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: EigenLib1.def,v 1.7 2017/09/27 07:51:14 mriedl Exp mriedl $ *)
FROM Deklera IMPORT MATRIX,VEKTOR,SUPERVEKTOR,PMATRIX;
TYPE ListenZeiger = POINTER TO Liste;
Liste = RECORD
Ungr : CARDINAL;
Obgr : CARDINAL;
next : ListenZeiger;
END;
(*--------------------------------------------------------------------*)
(* Der Type Liste wird benutzt um Blockungen eine Matrix anzuzeigen. *)
(*--------------------------------------------------------------------*)
PROCEDURE Gerschgorin(VAR A : ARRAY OF LONGREAL;
dim : CARDINAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Berechnet das Intervall, in dem sich alle Eigenwerte der *)
(* symetrischen, superlinear gespeicherten Matrix A befinden. *)
(* - Gerschgorin(A,dim) <= E_i <= Gerschgorin(A,dim), *)
(* E_i Eigenwert von A. *)
(*----------------------------------------------------------------*)
PROCEDURE SortEwEv(VAR EW : ARRAY OF LONGREAL;
VAR EV : ARRAY OF ARRAY OF LONGREAL;
dim,M : CARDINAL;
Ordnung : INTEGER);
(*----------------------------------------------------------------*)
(* Sortiert Eigenwerte und Eigenvektoren nach Gr\"o\3e der EW. *)
(* Wenn Ordnung = 1 in absteigender Reihenfolge *)
(* Wenn Ordnung <> 1 in ansteigender Reihenfolge *)
(* genau : Genauigkeit, mit der die EW bekannt sind *)
(*----------------------------------------------------------------*)
PROCEDURE PSortEwEv(VAR EW : ARRAY OF LONGREAL;
VAR EV : PMATRIX;
dim : CARDINAL;
Ordnung : INTEGER);
(*----------------------------------------------------------------*)
(* Sortiert Eigenwerte und Eigenvektoren. Parameter wie SortEwEV *)
(*----------------------------------------------------------------*)
PROCEDURE SortEigen( n : CARDINAL;
VAR d : ARRAY OF LONGREAL;
VAR r : ARRAY OF CARDINAL);
(*----------------------------------------------------------------*)
(* On exit d[r[k]] is the k.th eigenvalue in descending order *)
(* and v[i, r[k]] is the i.th component of the corresponding *)
(* eigenvector. This routine is an adopton of the "Handbook" *)
(* Algol60 procedure sorteig *)
(*----------------------------------------------------------------*)
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);
(*----------------------------------------------------------------*)
(* Diagonalisiert die quadratische Matrix A der Dimension n und *)
(* gibt die Eigenwerte in "d" zurueck. Wenn eivec = wahr werden *)
(* in V auch die Eigenvektoren berechnet. Die Anzahl der be- *)
(* noetigen Jacobi-Rotation wird in rot zurueckgegeben. *)
(* Die Routine ist eine 1:1 Uebersetzung der Algol-60 Routine *)
(* jacobi aus "dem Handbuch". *)
(* *)
(* Computes to eigenvalues d and - if "eivec" = true - the - *)
(* corresponding eigenvectors V of square symmetric matrix A. *)
(* Parameter rot gives on output the number of Jacobi rotations *)
(* needed. *)
(* *)
(* [1] Jacobi, C.G.J.; "Uber ein leichtes Verfahren, die in der *)
(* der Theorie der S"akularst"orungen vorkommenden Gleichun- *)
(* gen numerisch aufzul"osen. *)
(* Crelle's Journal 30, 51-94 (1846). *)
(* [2] Ruithauser, Heinz; "The Jacobi Method for Real Symmetric *)
(* Matrices", Numerische Mathematik 9, 1-10 (1966) *)
(*----------------------------------------------------------------*)
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 *)
(*----------------------------------------------------------------*)
(* Diagonalisiert die in A in Superlinearform "ubergebene Matrix *)
(* (A[ij] = A'[i][j] mit ij=i*(i-1) + j). Die Matrix A wird bei *)
(* der Diagonalisierung zerst"ort. *)
(* *)
(* EW : Eigenwerte der Matrix A. *)
(* EV : Eigenvektoren der Matrix A. *)
(* dim : Die aktuelle Dimension der Matrix A. *)
(* MaxIter : Maximalzahl der Jacobi-Iterationen (wenn 0, wird *)
(* MaxIter auf max(dim**3/2,15) gesetzt). *)
(* genau : Abbruchkriterium, das angibt, ab wann die Nichtdiag- *)
(* gonalelemente der Matrix A als Null anzusehen sind *)
(* (wenn 0.0, wird der Vorgabewert 1.0E-10 verwandt). *)
(* Ordnung : Wenn 1 weden die Eigenwerte in absteigender *)
(* wenn -1 aussteigender Reihenfolge sortiert. Ist der *)
(* Betrag von Ordnung ungleich 1, wird nicht sortiert. *)
(*----------------------------------------------------------------*)
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 *)
(*----------------------------------------------------------------*)
(* Erkl"arung wie Jacobi, nur die Form der "Ubergabe der Eigen- *)
(* vektormatrix ist verschieden, so da3 beleibig gro3e Matrizen *)
(* diagonalisiert werden k"onnen. *)
(*----------------------------------------------------------------*)
PROCEDURE Kaiser(VAR A : ARRAY OF ARRAY OF LONGREAL;
dim : CARDINAL;
VAR EW : ARRAY OF LONGREAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Eigenvalues and vectors of a symmetric positive (semi-) *)
(* definite matrix using Kaisers method. *)
(* Any symmetric matrix may be input, but if it is not positive *)
(* (semi-)definite, the absolute values of the eigenvalues will *)
(* be found. *)
(* *)
(* Arguments: *)
(* *)
(* A = on input, an array containing the matrix *)
(* on output, the columns of A contain the normalized *)
(* eigenvectors of A. *)
(* dim = The order of A *)
(* EW = on output the ordered eigenvalues *)
(* iFehl = error indicator *)
(* = 0 no error *)
(* = 1 n > nrows or n < 1 *)
(* = 2 failed to converge in 10 iterations *)
(* > 2 if Trace = SumEW, then all of the eigenvalues *)
(* are positive or zero. If SumEW > Trace, the *)
(* difference is twice the sum of the eigenvalues *)
(* which have been given the wrong signs ! *)
(* *)
(* This version is a translation of Alan Miller's Fortran version *)
(* to Modula-2. *)
(* *)
(* Kaiser, H. F.; 'The JK method: A procedure for finding the *)
(* eigenvalues of a real symmetric matrix', Comput. Journal 15, *)
(* 271-273 (1972) *)
(*----------------------------------------------------------------*)
PROCEDURE Wielandt( dim : CARDINAL;
VAR A : MATRIX;
VAR EVi : VEKTOR;
VAR EWi : LONGREAL;
genau : LONGREAL;
MaxIter : CARDINAL;
VAR nIter : CARDINAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Berechnet Eigenwert und Eigenvektor der Matrix A nach dem *)
(* Algorithmus von Wielandt. *)
(* *)
(* Wenn die zu iterirende Matrix singul\"ar oder sonst nicht *)
(* l\"o\3bar ist, wird Fehler auf TRUE gesetzt und iFehl ent- *)
(* sprechend gesetzt. *)
(* *)
(* iFehl = 1 : A mit Zeilenpivotierung nicht LU-zerlegbar *)
(* iFehl = 2 : A kann mit Vollpivotierung nicht LU-zerlegt werden *)
(* iFehl = 3 : Maximalzahl der Iterationen ueberschritten *)
(* die Loesung ist eventuell nicht brauchbar *)
(* iFehl = 4 : A kann nicht Householder-reduziert werden - *)
(* keine brauchbare Loesung wurde berechenet *)
(*----------------------------------------------------------------*)
PROCEDURE Tred2(VAR A : MATRIX; (* Zu Trilinearisierende Matrix *)
VAR Z : MATRIX; (* Reduzierte Matrix A *)
VAR HD : VEKTOR; (* Hauptdiagonale d. reduzierten Matrix A *)
VAR ND : VEKTOR; (* Nebendiagonale d. reduzierten Matrix A *)
dim : CARDINAL; (* Dimension von A *)
genau : LONGREAL); (* Genauigkeitswert *)
(*-----------------------------------------------------------------*)
(* Transformiert die Matrix A nach dem Householder-Algorithmus *)
(* auf Trilinearform. Die Informationen "uber die Transformationen *)
(* werden dabei in A gespeichert, die Trilineare Matrix in HD,ND *)
(* abgespalten. Die Eigenvektoren von HD,ND k"onnen mit d. Routine *)
(* EWEV dierekt berechnet werden, wenn Z unver"andert an EWEV als *)
(* Parameter EV "ubergeben wird. *)
(* Der Eingabeparameter Z kann mit A identisch sein, dann wird *)
(* die urspr"ungliche Matrix in EWEV mit ihren Eigenvektoren *)
(* "uberschrieben. *)
(* *)
(* Bearbeitet nach : *)
(* R.S.Matrin et.al. Numerische Mathematik 11, 181 (1968) *)
(* Handbook - Algolprozedur 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 *)
(*----------------------------------------------------------------*)
(* <== Fehler : Wird auf TRUE gesetzt, wenn ein Fehler *)
(* auftritt. *)
(* <== Fehlerflag : Fehlermeldung. *)
(* <==> HD : ==> Hauptdiagonale d. reduzierten Matrix A. *)
(* *)
(* Dient der Berechnung der Eigenwerte und Vektoren d. mit Tred2 *)
(* transformierten Matrix A. Es ist der Ausgabeparameter Z von *)
(* Tred2 als EV zu "ubergeben. *)
(* *)
(* Lit.: ??? *)
(* Handbook - Algolprozedur ??? *)
(*----------------------------------------------------------------*)
PROCEDURE ImTQL2( dim : CARDINAL;
VAR D,E : ARRAY OF LONGREAL;
VAR Z : ARRAY OF ARRAY OF LONGREAL;
VAR iErr : INTEGER);
(*----------------------------------------------------------------*)
(* This procedure finds the eigenvalues and eigenvectors of a *)
(* symmetric tridiagonal matrix by the implicit QL method. The *)
(* eigenvectors of a full symmetric matrix can also be found if *)
(* tred2 has been used to reduce this matrix to tridiagonal form. *)
(* *)
(* on entry *)
(* *)
(* dim : order N of the matrix. *)
(* D : the diagonal elements of the input matrix. *)
(* E : the subdiagonal elements of the input matrix *)
(* in its last N-1 positions. e(1) is arbitrary. *)
(* Z : contains the transformation matrix produced in the *)
(* reduction by tred2, if performed. If the *)
(* eigenvectors of the tridiagonal matrix are desired, *)
(* Z must contain the identity matrix. *)
(* *)
(* on exit *)
(* *)
(* D : contains the eigenvalues in ascending order. If an *)
(* error exit is made, the eigenvalues are correct but *)
(* unordered for indices 1,2,...,iErr-1. *)
(* E : has been destroyed. *)
(* Z : contains orthonormal eigenvectors of the symmetric *)
(* tridiagonal (or full) matrix. If an error exit is *)
(* made, Z contains the eigenvectors associated with the *)
(* stored eigenvalues. *)
(* iErr : zero - for normal return, *)
(* i - if the i-th eigenvalue has not been determined *)
(* after 30 iterations. *)
(* *)
(* This procedure is a translation of the Algol procedure imtql2 *)
(* *)
(* Ref: Martin,R.S.; Wilkinson,J.H.; Num. Math. 12, 377 (1968) *)
(* Dubrulle, Num. Math. 15, 450 (1970) *)
(* Handbook for Auto. Comp., Vol. II - Linear Algebra, *)
(* 241-248(1971). *)
(*----------------------------------------------------------------*)
PROCEDURE GivTri(VAR A : SUPERVEKTOR; (* ==> Marix *)
VAR HD : VEKTOR; (* <== Hauptdiagonale *)
VAR ND : VEKTOR; (* <== Nebendiagonale *)
dim : CARDINAL;
genau : LONGREAL);
(*----------------------------------------------------------------*)
(* Transformiert die symmetrische Matrix A nach dem Givens- *)
(* Algorithmus auf Tridiagonalgestalt. *)
(* Lit.: H.R.Schwarz, Numerische Mathematik *)
(* B.G.Teubner, Stuttgart 1988 *)
(*----------------------------------------------------------------*)
PROCEDURE GivTriBak(VAR EV : MATRIX; (* Eigenvektoren *)
VAR A : SUPERVEKTOR; (* Transformierte Matrix *)
maxEV : CARDINAL; (* Anzahl der EV *)
VAR genau : LONGREAL;
dim : CARDINAL); (* Dimension von A *)
(*----------------------------------------------------------------*)
(* A : In GivTri transformierte symmetrische Matrix A. *)
(* Berechnet aus den Eigenvektoren der trilinearisierten Matrix *)
(* A mit den Informationen "uber sin und cos der Givens- *)
(* transformation die Eigenvektoren der untransformierten Matrix *)
(* A zur"uck. *)
(* Dabei gilt : J = Qt'AQ' , EV(A) = Q*EV(J) *)
(* Q'=U(1)*U(2)*...*U(n) *)
(* Q=U(n)*U(n-1)*...*U(1) *)
(* mit U(i) : Elemantare Jacobi-Rotationsmatrix, um Element i *)
(* in A zu Null zu machen.' *)
(*----------------------------------------------------------------*)
PROCEDURE HhQL(VAR A : MATRIX; (* Zu Diagonalisierende Matrix *)
VAR EV : MATRIX; (* Eigenvektoren von A. *)
VAR EW : ARRAY OF LONGREAL; (* Eigenwerte von A. *)
dim : CARDINAL; (* Dimension von A *)
genau : LONGREAL; (* Genauigkeitsanforderung. *)
ord : INTEGER); (* Sortieroption, siehe SortEwEv *)
(*----------------------------------------------------------------*)
(* Transformiert die Matrix A nach dem Householder-Algorithmus *)
(* auf Trilinearform. Die Informationen \"uber die *)
(* Transformationen werden dabei in A gespeichert, die Trilineare *)
(* Matrix in EW,ND abgespalten. Die Eigenvektoren von EW,ND *)
(* werden dann nach dem QL-Algoritmus berechnet und die Eigen- *)
(* vektoren aus den Information \"uber die Transformationen *)
(* ermittelt. *)
(* Der Eingabeparameter EV kann mit A identisch sein, dann wird *)
(* die urspr\"ungliche Matrix in HhQL mit ihren Eigenvektoren *)
(* \"uberschrieben. *)
(* *)
(* Bearbeitet nach : *)
(* [1] Matrin, R.S. et.al; Numer. Math. 11, 181 (1968) *)
(* [2] Martin, R.S.; Wilkinson, J.H.; Numer. Math. 12, 377 (1968) *)
(* Zusammenfassung der *)
(* Handbook - Algolprozeduren tred2 [1] und imtql2 [2] *)
(*----------------------------------------------------------------*)
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 *)
(*----------------------------------------------------------------*)
(* Transformiert die Matrix A nach dem Householder-Algorithmus *)
(* auf Trilinearform. Die Informationen \"uber die *)
(* Transformationen werden dabei in A gespeichert, die trilineare *)
(* Matrix in HD,ND abgespalten. *)
(* abgespalten. *)
(* Die Eigenvektoren von HD,ND k\"onnen mit d. Routine HhTriBak *)
(* auf diejenigen von A zur\"ucktransformiert werden, wenn A *)
(* unver\"andert an HhTriBak \"ubergeben wird. *)
(* *)
(* This proeedure reduces the given lower triangle of a symmetrie *)
(* matrix, A, stored row by row in the array a[1:n(n+1)/2], to *)
(* tridiagonal form using Householders reduction. The diagonal *)
(* of the result is stored in the array HD[1:n] and the *)
(* sub-diagonal in the last n-1 stores of the array ND[1:n] (with *)
(* the additional element nd[1] = 0). SqrND[i] is set to equal *)
(* ND[i]^2. The array A is used to store sufficient information *)
(* for the details of the transformation to be recoverable in the *)
(* procedure HhTirBak *)
(* *)
(* Bearbeitet nach : *)
(* R.S.Matrin et.al. Numerische Mathematik 11, 181 (1968) *)
(* Handbook - Algolprozedur tred3 *)
(*----------------------------------------------------------------*)
PROCEDURE HhTriBak(VAR EV : MATRIX;
VAR A : SUPERVEKTOR;
dim : CARDINAL;
m1,m2 : CARDINAL);
(*----------------------------------------------------------------*)
(* Transformiert die Eigenvektoren der mit TRED3 trilinearisier- *)
(* ten Matrix A auf diejenigen der urspr\"unglichen Matrix A *)
(* *)
(* EV : Zu transformierende Eigenvekoren *)
(* A : Housholdertransformierte Matrix *)
(* dim : Dimension von A & EV *)
(* m1,m2 : EV[m1] bis EV[m2] werden transformiert. *)
(* *)
(* This procedure performs, on the matrix of eigenvectors, Z(EV), *)
(* stored in the array EV[1:n,m1:m2], a backtransformation to *)
(* form the eigenvectors of the original symmetrie matrix from *)
(* the eigenvectors of the tridiagonal matrix. The new vectors *)
(* are overwritten on the old ones. *)
(* The details of the Householder reduction must be stored in *)
(* the array a[1:n(n+1)/2], as left by the procedure HhTri. *)
(* If z denotes any column of the resultant matrix Z, then z *)
(* satisfies z^t x z = Z(input)^T * z(input); *)
(* *)
(* Bearbeitet nach : *)
(* *)
(* R.S.Matrin et.al. Numerische Mathematik 11, 181 (1968) *)
(* Handbook - Algolprozedur trbak3 *)
(*----------------------------------------------------------------*)
PROCEDURE Bisect(VAR HD : VEKTOR; (* Hauptdiagonale der Matrix. *)
ND : VEKTOR; (* Nebendiagonale der Matrix. *)
SqrND : VEKTOR; (* Quadrate von ND. *)
VAR EW : VEKTOR; (* Eigenwerde 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. *)
(*----------------------------------------------------------------*)
(* Berichnet die Eigenwerte einer in HD,ND gespeicherten tri- *)
(* digonalen Matrix durch Bisection und Sturmsequenzen. *)
(* *)
(* HD is the diagonal, ND the sub-diagonal and SqrND the squared *)
(* subdiagonal of a symmetric tridiagonal matrix of order dim. *)
(* The eigenvalues EW[m1] , ... , EW[m2], where m2 is not less *)
(* than m1 and EW[i+1] is not less than EW[i], are calculated by *)
(* the method of bisection and stored in the vector x. Bisection *)
(* is continued until the upper and lower bounds for an eigen- *)
(* value differ by less than eps1, unless at some earlier stage, *)
(* the upper and lower bounds differ only in the least *)
(* significant digits, eps2 gives an extreme upper bound for the *)
(* error in any eigenvalue, but for certain types of matrices the *)
(* small eigenvalues are determined to a very much higher *)
(* accuracy. In this case, ebs1 should be set equal to the error *)
(* to be tolerated in the smallest eigenvalue. It must not be set *)
(* equal to zero. *)
(* *)
(* Bearbeite nach : W.Barth et.al., Num. Math. 9, 386 (1967) *)
(* Handbook - Algolprozedur bisect. *)
(*----------------------------------------------------------------*)
PROCEDURE TriQR( HD,ND : VEKTOR; (* ==> Trilineare Matrix *)
VAR EW : VEKTOR; (* <== Eigenwerte von A *)
dim : CARDINAL; (* Dimension von A *)
VAR UngrObgr : ListenZeiger; (* <== Blockungen *)
genau : LONGREAL);
(*----------------------------------------------------------------*)
(* Berechnet die Eigenwerte einer tridiagonalen symmetrischen *)
(* Matrix A nach dem QR-Algoritmus. *)
(* A wird dabei in in HD (Hauptdiagonale) und ND (Nebendiagonale) *)
(* \"ubergeben. *)
(* Tritt ein Fehler in der Routine auf, so wird die globale *)
(* Variable Fehler auf TRUE gesetzt. Die Argumente HD,ND werden *)
(* aber nicht ver\"andert, so da\3 in diesem Fall z.B. die *)
(* Routine Bisect benutzt werden kann, da UngrObgr nicht wieder *)
(* gel\"oscht wird. *)
(* genau : Alle Werte > genau werden 'wegrotiert'. *)
(* Lit. : H.R.Schwarz, Numerische Mathematik, *)
(* B.G.Teubner, Stuttgart 1988 *)
(*----------------------------------------------------------------*)
PROCEDURE TriQL( HD,ND : VEKTOR; (* ==> Trilineare Matrix *)
VAR EW : VEKTOR; (* <== Eigenwerte von A *)
dim : CARDINAL; (* Dimension von A *)
VAR UngrObgr : ListenZeiger; (* <== Blockungen *)
genau : LONGREAL);
(*----------------------------------------------------------------*)
(* Berechnet die Eigenwerte einer tridiagonalen symmetrischen *)
(* Matrix nach dem QL-Algoritmus. *)
(* Diese wird dabei in in HD (Hauptdiagonale) und ND (Nebendia- *)
(* gonale) \"ubergeben. *)
(* genau : Alle Werte > genau werden 'wegrotiert'. *)
(* Fehlerbehandlung wie in PROCEDURE TriQR ! *)
(* *)
(* Lit.: Martin,R.S.; Wilkinson,J.H.; Numer. Math. 12, 377 (1968) *)
(* Handbook - Algolprozedur imtql1. *)
(*----------------------------------------------------------------*)
PROCEDURE EwEvQR(VAR HD : VEKTOR;
ND : VEKTOR;
VAR EW : VEKTOR;
VAR EV : MATRIX;
dim : CARDINAL;
genau : LONGREAL);
(*----------------------------------------------------------------*)
(* Berechnet die Eigenwerte einer tridiagonalen symmetrischen *)
(* Matrix nach dem QR-Algoritmus. *)
(* Diese wird dabei in in HD (Hauptdiagonale) und ND (Neben- *)
(* diagonale) \"ubergeben. *)
(* Die Eigenvektoren der tridiagonalen Matrix werden dabei in EV *)
(* aus den QR-Schritten aufgebaut. *)
(* genau : Alle Werte > genau werden 'wegrotiert'. *)
(* Lit. : H.R.Schwarz, Numerische Mathematik *)
(* B.G.Teubner, Stuttgart 1988 *)
(*----------------------------------------------------------------*)
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 *)
(*----------------------------------------------------------------*)
(* Berechnet den Eigenvektor EVi von A ( in HD,ND ) durch Vektor- *)
(* iteration nach Wielandt, wobei die Nachiteration der Eigen- *)
(* werte von A unterbleibt. *)
(* *)
(* A*X=C X : Vektor der m+1. Iteration *)
(* C : Vektor der m. Iteration *)
(*----------------------------------------------------------------*)
PROCEDURE RSymEwEv(VAR EV : MATRIX;
VAR EW : VEKTOR;
VAR A : SUPERVEKTOR;
dim : CARDINAL;
maxEW : CARDINAL;
maxEV : CARDINAL;
Ortho : CARDINAL;
genau : LONGREAL;
ord : INTEGER);
(*----------------------------------------------------------------*)
(* Berechent die Eigenvektoren und -werte einer symmetrischen *)
(* reellen Matrix in Superlinearform (Treiberroutine). *)
(* *)
(* <== EV : Eigenvektoren von A. *)
(* <== EW : Eigenwerte von A. *)
(* ==> A : quadratische Matrix *)
(* ==> dim : Dimension von A *)
(* ==> maxEW : Anzahl der zu bestimmenden Eigenwerte. *)
(* ==> maxEV : Anzahl der zu bestimmenden Eigenvektoren. *)
(* ==> Ortho : Wenn 1 und maxEW und maxEV gleich dim *)
(* dann werden die Eigenvektoren so berechent, *)
(* da\3 deren Orthogonalit\"at sicher ist. *)
(* ==> genau : Geforderte Genauigkeit. *)
(* ==> ord : wenn *)
(* -1 : EW[1] <= EW[2] ... <= EW[n] *)
(* 0 : keine Ordnung. *)
(* 1 : EW[1] >= EW[2] ... >= EW[n] *)
(*----------------------------------------------------------------*)
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);
(*----------------------------------------------------------------*)
(* L\"o\3t das generalisierte Eigenwertproblem Hx=\Lambda Sx nach *)
(* dem Algorithmus von L\"owdin. *)
(* *)
(* \lambda : Eigenwert des Problems (Aus EW). *)
(* ord : Wenn ord = *)
(* -1 : die EW,EV werden nach aufsteigendem EW *)
(* 1 : die EW,EV werden nach absteigendem Eigenwert *)
(* 0 : nicht *)
(* geordnet. *)
(* Neu : Wenn Neu = *)
(* 1 : wird die Tansformationsmatrix neu berechnet. *)
(* 0 : wird die Transformationsmatrix des vorherigen *)
(* Aufrufs benutzt. *)
(* Norm : Wenn Norm = *)
(* 0 : werden die Eigenvektoren nicht normiert. *)
(* 1 : werden die Eigenvektoren auf 1 normiert. *)
(* Transform : Wenn 1 werden die Eigenvektoren von *)
(* H'= V^{-1/2)' H V^{-1/2) auf diejenigen von H *)
(* transformiert. *)
(* Wenn 2 wird nur H' berechnet, ohne dieses zu *)
(* diagonalisieren. *)
(* *)
(* On the Non-Orthogonality Problem Connected with the use of *)
(* Atomic Wave Functions in the Theory of Molecules and Crystals, *)
(* J. Chem. Phys. 18, 367-370, 1950 *)
(*----------------------------------------------------------------*)
PROCEDURE Transform(VAR EV : MATRIX;
dim : CARDINAL;
Form : CARDINAL);
(*----------------------------------------------------------------*)
(* Wenn Form = *)
(* 0 : SV'= {S^{-1/2}}^+ SV S^{-1/2} *)
(* 1 : SV'= {S^{ 1/2}}^+ SV S^{ 1/2} *)
(* wobei S diejenige prositiv definite Matrix ist, die bei dem *)
(* letzten Aufruf der Prozedur Loewdin mit dem Parameter Neu = 1 *)
(* an diese \"ubergeben wurde. *)
(*----------------------------------------------------------------*)
PROCEDURE TransfSV(VAR SV : SUPERVEKTOR;
dim : CARDINAL;
Form : CARDINAL);
(*----------------------------------------------------------------*)
(* Wenn Form = *)
(* 0 : SV'= {S^{-1/2}}^+ SV S^{-1/2} *)
(* 1 : SV'= {S^{ 1/2}}^+ SV S^{ 1/2} *)
(* wobei S diejenige prositiv definite Matrix ist, die bei dem *)
(* letzten Aufruf der Prozedur Loewdin mit dem Parameter Neu = 1 *)
(* an diese \"ubergeben wurde. *)
(*----------------------------------------------------------------*)
PROCEDURE Reduce1(VAR A,B : ARRAY OF ARRAY OF LONGREAL;
VAR dL : ARRAY OF LONGREAL;
dim : CARDINAL; (* n in original *)
transB : BOOLEAN); (* FALSE instead of dim < 0 *)
(*-----------------------------------------------------------------*)
(* ==> getestet <== *)
(*-----------------------------------------------------------------*)
(* Reduction of the general symmetrie eigenvalue problem *)
(* A*x = lambda*B*x, *)
(* with symmetrie matrix A and symmetrie positive definite matrix *)
(* B, to the equivalent standard problem P*z = lambda*z. *)
(* The upper triangle, inc1uding diagonal elements, of A and B are *)
(* given in the arrays a[1:n,1:n] and b[1:n,1:n]. *)
(* L (B=L xLT) is formed in the remaining strietly lower triangle *)
(* of the array b with its diagonal elements in the array dl[i:n], *)
(* and the lower tri angle of the symmetrie matrix *)
(* P (P=inv(L) x A x inv(LT)) *)
(* is formed in the lower triangle of the array a, inc1uding the *)
(* diagonal elements. Hence the diagonal elements of A are lost. *)
(* If n < 0, it is assumed that L has already been formed. *)
(* The proeedure will fail if B, perhaps on aeeount of rounding *)
(* errors, is not positive definite *)
(* *)
(* Lit.: *)
(* Martin, R,S.; Wilkinson, J.H.; Numer. Math. 11(2) 99 (1968) *)
(* Handbook - Algolprozedur reduc 1. *)
(*-----------------------------------------------------------------*)
PROCEDURE ReBakA(VAR EV : ARRAY OF ARRAY OF LONGREAL; (* Z in original *)
dim : CARDINAL; (* n in original *)
m1,m2 : CARDINAL;
VAR B : ARRAY OF ARRAY OF LONGREAL;
VAR dL : ARRAY OF LONGREAL);
(*-----------------------------------------------------------------*)
(* ==> getestet <== *)
(*-----------------------------------------------------------------*)
(* This procedure performs, on the matrix of eigenvectors, Z, *)
(* stored in the array z[1:n,m1:m2], a backward substitution *)
(* LT xX =Z, overwriting X on Z. *)
(* The diagonal elements of L must be stored in the array dL[1:n], *)
(* and the remaining triangle in the strictly lower triangle of *)
(* the array B[1:n,1:n]. The procedures reduc1 and reduc2 leave L *)
(* in this desired form. *)
(* If x denotes any column of the resultant matrix X, then x *)
(* satisfies xT*B*x=zTXz, where B=L*LT *)
(* *)
(* Lit.: Martin, R,S.; Wilkinson, J.H.; Numer. Math. 11, 99 (1968) *)
(* Handbook - Algolprozedur rebaka. *)
(*-----------------------------------------------------------------*)
END EigenLib1.