DEFINITION MODULE EigenLib2;
(*------------------------------------------------------------------------*)
(* Stellt Routinen fuer die Eigenwertberechung von unsymmetrischen *)
(* realwertigen Matrizen und symmetrischer hermitische Matrizen zur *)
(* Verfuegung. *)
(* Einige der Prozeduren zur Berechnung der Eigenwerte und -vektoren *)
(* reelen generellen Matrix stammen aus der NumAl Algol Bibiliothek die *)
(* des Stichting Centrum Wiskunde & Informatica ( *)
(* *)
(* Routines for the calculation of eigensystems of unsymmetric real *)
(* general and hermitian matrices *)
(* *)
(* This Module contains some Routines for calculatiing Eigenvalues and *)
(* Eigenvectors of a real general matrices where parts of the routines *)
(* stem from the NumAl Algol library which had been provided by *)
(* Stichting Centrum Wiskunde & Informatica. *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Lizenz : GNU public licence *)
(*------------------------------------------------------------------------*)
(* $Id: EigenLib2.def,v 1.5 2017/10/22 17:51:52 mriedl Exp mriedl $ *)
FROM Deklera IMPORT VEKTOR,CVEKTOR,MATRIX,INTVEKTOR,CARDVEKTOR;
PROCEDURE Balance(VAR A : ARRAY OF ARRAY OF LONGREAL;
dim : CARDINAL;
VAR low,high : CARDINAL;
VAR Scale : ARRAY OF LONGREAL);
(*----------------------------------------------------------------*)
(* Balanciert eine quadratische Matrix (Als Vorbereitung zur *)
(* Eigensystemberechnung) *)
(* Werden die Eigenvektoren von A ben"otigt, ist nach deren *)
(* Berechnung die Prozedur BalBak anzuwenden ! *)
(* *)
(* Reduce the norm of a [1:n,1:n] by exact diagonal similarity *)
(* transformations stored in Scale[1:n] *)
(* *)
(* Lit.: Parlett,B.N.; Reinsch,C.; Numer. Math 13, 293 (1969) *)
(* Handbook - Algolprozedur balance. *)
(*----------------------------------------------------------------*)
PROCEDURE BalBak( dim : CARDINAL; (* Ordnung der Matrix Z *)
low,high : CARDINAL; (* Von Balance ermittelt *)
VAR Scale : ARRAY OF LONGREAL; (* Von Balance ermittelt *)
m : CARDINAL;
VAR Z : ARRAY OF ARRAY OF LONGREAL); (* Eigenvekt. *)
(*----------------------------------------------------------------*)
(* Formt die Eigenvektoren einer reelen, unsymmetrischen Matrix *)
(* durch R"ucktransformation der Eigenvektoren einer durch die *)
(* Prozedur Balance 'balancierten' Matrix. *)
(* *)
(* ==> Scale : enth"alt Informationen "uber die Permutationen *)
(* und Skalierungsfaktoren, die von der Prozedur *)
(* Balance benutzt wurden. *)
(* m : Anzahl der Vektoren in Z, die R"ucktrans- *)
(* formiert werden sollen. *)
(* <==> Z : Eigenvektoren, die R"ucktransformiert werden *)
(* sollen (in den ersten m Zeilen). *)
(* *)
(* Backward transformation of a set of right-hand eigenvectors *)
(* of a balanced matrix into the eigenvectors of the original *)
(* matrix from which the balanced matrix was derived by a call of *)
(* procedure balance *)
(* *)
(* Lit.: Parlett; Reinsch; Num. Math. 13, 293-304 (1969) *)
(* Handbook - Algolprozedur balbak. *)
(*----------------------------------------------------------------*)
PROCEDURE GivHess(VAR A : MATRIX;
dim : CARDINAL;
genau : LONGREAL);
(*----------------------------------------------------------------*)
(* Transformiert die unsymmetrische Matrix A nach dem Givens- *)
(* Algorithmus auf obere Hessenbergform. *)
(* Lit.: H.R.Schwarz, Numerische Mathematik *)
(* B.G.Teubner, Stuttgart 1988 *)
(*----------------------------------------------------------------*)
PROCEDURE ElmHess(VAR A : ARRAY OF ARRAY OF LONGREAL;
dim : CARDINAL;
k,l : CARDINAL; (* Low,High *)
VAR Index : ARRAY OF CARDINAL);
(*----------------------------------------------------------------*)
(* Transformiert die unsymmetrische reelle Matrix A auf obere *)
(* Hessenberg. *)
(* *)
(* Given the unsymmetric matrix, A, stored in the array *)
(* a[1:n,1:n], this procedure reduces the sub-matrix of order *)
(* l-k+1, which starts at the element A[k,k] and finishes at the *)
(* element a[l,l], to Hessenberg form, H, by non-orthogonal *)
(* elementary transformations. The matrix H is overwritten on A *)
(* with details of the transformations stored in the remaining *)
(* triangle under H and in the array int[k:l] *)
(* *)
(* Lit.: Martin, R.S. et. al. Numer. Math. 12, 349 (1968) *)
(* Handbook - Algolprozedur elmhes. *)
(*----------------------------------------------------------------*)
PROCEDURE ElmHessTrans(VAR V : ARRAY OF ARRAY OF LONGREAL;
VAR H : ARRAY OF ARRAY OF LONGREAL;
VAR Index : ARRAY OF CARDINAL;
dim : CARDINAL;
low,upp : CARDINAL);
(*-----------------------------------------------------------------*)
(* Akkumuliert die in der Prozedur ElmHess verwandten Transformat- *)
(* ionsschritte (in H und Index von ElmHess unver"andert "uber- *)
(* geben). low und upp werden von der Prozedur Balance berechnet, *)
(* wird diese nicht verwandt, setzt man low=1 und upp=dim. *)
(* Die Matrix V kann dann als Eingabe f"ur die Matrix EV (Eigen- *)
(* vektormatrix) der Prozedur HessQREwEv verwandt werden. *)
(* *)
(* Form the matrix of accumulated transformations in the array *)
(* V[1:n,1:n] from the information left by procedure elmhes *)
(* below the upper Hessenberg matrix, H, in the array H[1:n,1:n] *)
(* and in the integer array int[1:n] *)
(* *)
(* Lit.: Peters, G. et. al. Numer. Math. 16, 181 (1970) *)
(* Handbook - Algolprozedur elmtrans. *)
(*-----------------------------------------------------------------*)
PROCEDURE RHessQR(VAR H : MATRIX; (* Hessenbergmatrix *)
VAR EW : ARRAY OF LONGCOMPLEX; (* Eigenwerte *)
dim : CARDINAL); (* Getestet *)
(*----------------------------------------------------------------*)
(* Berechnet die Eigenwerte der oberen Hessenbergmatrix H. *)
(* Die Matrix H wird dabei ver"andert ! *)
(* Lit.: Matrin, R.S. et. al; Numer. Math. 14, 219 (1970) *)
(* Handbook - Algolprozedure hqr. *)
(*----------------------------------------------------------------*)
PROCEDURE HQR(VAR H : MATRIX; (* Wird zerst"ort *)
VAR EW : ARRAY OF LONGCOMPLEX; (* Eigenwerte *)
n : CARDINAL; (* Dimension von H *)
low,upp : CARDINAL; (* Ausgabe von Balence *)
VAR cnt : ARRAY OF INTEGER);
(*----------------------------------------------------------------*)
(* Berechnet die Eigenwerte einer oberen Hessenbergmatrix. *)
(* *)
(* Die Routine ist eine Uebersetzung der Algol routine hqr2 *)
(* wobei alle Eigenvektoranteile herausgenommen wurden. *)
(* *)
(* Lit. : Peters, G. et. al; Numer. Math. 16, 181 (1970) *)
(* Handbook - Algolprozedur hqr2. *)
(*----------------------------------------------------------------*)
PROCEDURE HQR2( N : INTEGER;
Low : INTEGER;
High : INTEGER;
VAR H : MATRIX;
VAR Z : MATRIX;
VAR Wr : VEKTOR;
VAR Wi : VEKTOR;
VAR iErr : INTEGER);
(*----------------------------------------------------------------*)
(* hqr2 computes eigenvalues and eigenvectors of a *)
(* real upper hessenberg matrix *)
(* *)
(* This procedure is a translation of the Fortran subroutine hqr2 *)
(* which was a translation of the algol procedure hqr2. *)
(* Num. Math. 16, 181-204 (1970) by Peters and Wilkinson. *)
(* handbook f. auto. comp., vol.ii-linear algebra, 372-395 (1971) *)
(* *)
(* This procedure finds the eigenvalues and eigenvectors *)
(* of a real upper hessenberg matrix by the qr method. The *)
(* eigenvectors of a real general matrix can also be found *)
(* if elmhes and eltran or orthes and ortran have *)
(* been used to reduce this general matrix to hessenberg *)
(* form and to accumulate the similarity transformations. *)
(* *)
(* on input *)
(* *)
(* n is the order of the matrix. *)
(* *)
(* Low and High are integers determined by the balancing *)
(* routine balanc. if balanc has not been used, *)
(* set Low:=1, High:=n. *)
(* *)
(* h contains the upper hessenberg matrix. *)
(* *)
(* z contains the transformation matrix produced by *)
(* eltran after the reduction by elmhes, or by ortran *)
(* after the reduction by orthes, if performed. If the *)
(* eigenvectors of the hessenberg matrix are desired, *)
(* z must contain the identity matrix. *)
(* *)
(* on output *)
(* *)
(* h has been destroyed. *)
(* *)
(* wr and wi contain the real and imaginary parts, *)
(* respectively, of the eigenvalues. the eigenvalues *)
(* are unordered except that complex conjugate pairs *)
(* of values appear consecutively with the eigenvalue *)
(* having the positive imaginary part first. if an *)
(* error exit is made, the eigenvalues should be correct *)
(* for indices ierr+1, ,n. *)
(* *)
(* z contains the real and imaginary parts of the eigen- *)
(* vectors. If the i-th eigenvalue is real, the i-th *)
(* column of z contains its eigenvector. If the i-th *)
(* eigenvalue is complex with positive imaginary part, *)
(* the i-th and (i+1)-th columns of z contain the real *)
(* and imaginary parts of its eigenvector. The eigen- *)
(* vectors are unnormalized. If an error exit is made, *)
(* none of the eigenvectors has been found. *)
(* *)
(* ierr is set to *)
(* *)
(* zero for normal return, *)
(* j if the limit of 30*n iterations is *)
(* exhausted while the j-th eigenvalue is *)
(* being sought. *)
(* *)
(* calls cdiv for complex division. *)
(* *)
(* This version dated August 1983. *)
(*----------------------------------------------------------------*)
PROCEDURE TfmReaHes(VAR a : MATRIX;
n : INTEGER;
VAR norm : LONGREAL;
VAR int : INTVEKTOR);
(*----------------------------------------------------------------*)
(* transforms a matrix into a similar upper-hessenberg matrix by *)
(* means of wilkinsons transformation. *)
(* *)
(* n : the order of the given matrix *)
(* a : a[1:n,1:n] *)
(* entry : the matrix to be transformed *)
(* exit : the upper-hessenberg matrix is delivered in *)
(* the upper triangle and the first subdiagonal *)
(* of a, the (nontrivial elements of the) *)
(* transforming matrix, l, in the remaining part *)
(* of a, i.e. a[i,j] = l[i,j + 1], for i=3,...,n *)
(* and j = 1,...,i - 2 *)
(* norm : on exit the infinity norm of the original matrix *)
(* int : int[1:n] *)
(* exit : the pivotal indices defining the stabilizing *)
(* row and column interchanges *)
(* *)
(* Source: NUMAL Numerical Algol library (Stichting CWI) *)
(*----------------------------------------------------------------*)
PROCEDURE BakReaHes1(VAR A : MATRIX;
N : CARDINAL;
VAR Int : CARDVEKTOR;
VAR V : VEKTOR);
(*----------------------------------------------------------------*)
(* tfmreahes transforms a matrix into a similar upper-hessenberg *)
(* matrix by means of wilkinsons transformation. *)
(* bakreahes1 performs the corresponding back transformation *)
(* on a vector and should be called after tfmreahes. *)
(* *)
(* N : the length of the vector to be transformed *)
(* A : the (nontrivial elements of the) transforming matrix, *)
(* l, as produced by tfmreahes must be given in the part *)
(* below the first subdiagonal of A. *)
(* i.e. a[i,j] = l[i,j+1], for i=3,...,n and *)
(* j = 1,...,i - 2 *)
(* Int : pivotal indices defining the stabilizing row and *)
(* column interchanges as produced by tfmreahes *)
(* V : on entry the vector to be transformed *)
(* on exit the transformed vector. *)
(* *)
(* Source: NUMAL Numerical Algol library (Stichting CWI) *)
(*----------------------------------------------------------------*)
PROCEDURE RealVecHes(VAR a : MATRIX;
n : INTEGER;
lambda : LONGREAL;
VAR v : VEKTOR;
Norm : LONGREAL; (* em[1] *)
Toler : LONGREAL; (* em[6] = -10 *)
MaxIt : INTEGER; (* em[8] = 5 *)
VAR Iter : INTEGER; (* em[9] *)
VAR Resid : LONGREAL); (* em[7] *)
(*----------------------------------------------------------------*)
(* RealVecHes calculates an eigenvector corresponding to a given *)
(* approximate real eigenvalue of a real upper-hessenberg matrix, *)
(* by means of inverse iteration *)
(* *)
(* formal parameters *)
(* *)
(* a : on entry: the elements of the real upper *)
(* hessenberg matrix must be given in the *)
(* upper triangle and the first subdiagonal *)
(* of array a[1:n,1:n] *)
(* on exit : the hessenberg part of array a is altered *)
(* n the order of the given matrix *)
(* lambda : the given real eigenvalue of the upper hessenberg *)
(* matrix *)
(* v : the calculated eigenvector is delivered in v[1:n] *)
(* Norm : a norm of the given matrix *)
(* Toler : the tolerance used for the eigenvector. The inverse *)
(* iteration ends if the euclidian norm of the residue *)
(* vector is smaller than Norm*Tol *)
(* MaxIt : on entry the maximum allowed number of iterations *)
(* Iter : the number of inverse iterations performed. *)
(* If Resid remains larger than Norm*Tol during *)
(* MaxIt iterations, the value MaxIt+1 is delivered *)
(* Resid : on exit the euclidian norm of the residue vector of *)
(* the calculated eigenvector *)
(* *)
(* Source: NUMAL Numerical Algol library (Stichting CWI) *)
(*----------------------------------------------------------------*)
PROCEDURE ComVecHes(VAR a : MATRIX;
n : INTEGER;
lambda : LONGREAL;
mu : LONGREAL;
VAR u,v : VEKTOR;
Norm : LONGREAL; (* em[1] *)
Toler : LONGREAL; (* em[6] = -10 *)
MaxIt : INTEGER; (* em[8] = 5 *)
VAR Iter : INTEGER; (* em[9] *)
VAR Resid : LONGREAL); (* em[7] *)
(*----------------------------------------------------------------*)
(* ComVecHes calculates an eigenvector corresponding to a given *)
(* approximate real eigenvalue of a real upper-hessenberg matrix, *)
(* by means of inverse iteration *)
(* *)
(* formal parameters *)
(* *)
(* a : on entry: the elements of the real upper *)
(* hessenberg matrix must be given in the *)
(* upper triangle and the first subdiagonal *)
(* of array a[1:n,1:n] *)
(* on exit : the hessenberg part of array a is altered *)
(* n the order of the given matrix *)
(* lambda, *)
(* mu : the real and imaginary part of the given eigenvalue *)
(* of the upper hessenberg matrix *)
(* u,v : the real and imaginary parts of the calculated *)
(* eigenvector are delivered in the arrays u,v[1:n] *)
(* Norm : a norm of the given matrix *)
(* Toler : the tolerance used for the eigenvector. The inverse *)
(* iteration ends if the euclidian norm of the residue *)
(* vector is smaller than Norm*Tol *)
(* MaxIt : on entry the maximum allowed number of iterations *)
(* Iter : the number of inverse iterations performed. *)
(* If Resid remains larger than Norm*Tol during *)
(* MaxIt iterations, the value MaxIt+1 is delivered *)
(* Resid : on exit the euclidian norm of the residue vector of *)
(* the calculated eigenvector *)
(* *)
(* Source: NUMAL Numerical Algol library (Stichting CWI) *)
(*----------------------------------------------------------------*)
PROCEDURE EigenRGM(VAR A : MATRIX;
N : CARDINAL;
VAR Wr,Wi : VEKTOR;
VAR V : MATRIX;
ifBal : BOOLEAN;
norm : CARDINAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Compute all eigenvalues and -vectors of a real general matrix *)
(* *)
(* A : input matrix *)
(* N : size of matrix *)
(* Wr,Wi : real and maginary parts of eigenvalues *)
(* V : Eigenvectors *)
(* ifBal : use balancing of the input matrix A *)
(* norm : normalize Eigenvectors if norm = 1 or 2 *)
(* iFehl : return code *)
(*----------------------------------------------------------------*)
PROCEDURE EigenRGMsel(VAR A : MATRIX;
N : CARDINAL;
VAR EW : CVEKTOR;
VAR V : MATRIX;
ifBal : BOOLEAN;
norm : CARDINAL;
ord : INTEGER;
ifAll : BOOLEAN;
VAR Select : ARRAY OF BOOLEAN;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Berechne die Eigenwerte und -vektoren der reelen generellen *)
(* Matrix A der Dimension N. *)
(* *)
(* <==> A : Eingabematrix (wird zerstoert) *)
(* ==> N : Dimension der Matrix *)
(* <== EW : Berechnete Eigenwerte [1..N] *)
(* ==> ifBal : Verwende Balance/BalBak J/N *)
(* ==> norm : Wenn 1 werden die Eigenvektoren 1-Normiert *)
(* Wenn 2 werden die Eigenvektoren 2-Normiert *)
(* ==> ord : Wenn 1 sortiere die Eigenwerte nach Absolutwert *)
(* absteigend, ansonsten keine Sortierung *)
(* ==> ifAll : Wenn wahr werden alle Eigenvektoren berechnet *)
(* und Select wird ignoriert. In diesem Fall werden *)
(* alle Elemente von Select auf "wahr" gesetzt *)
(* <==> Select : Wenn Select[i] = wahr wird der zu EW[i] *)
(* gehoerende Eigenvektor berechnet, sonst nicht *)
(* <== iFehl : Rueckgabecode *)
(* 0 : Alles in Ordnung *)
(* 1 : Residium mindestens eines Eigenvektors ist *)
(* grosser als MachEps(REAL) *)
(* *)
(* Die Realteile nahezu entarteter Eigenwerte werden in der *)
(* Routine solange um den Betrag Norm(A)*MachEps verschoben bis *)
(* eine getrennte Berechnung der Eigenvektoren sichergestellt ist *)
(* *)
(* Benutzt die Routinen *)
(* Balance,TrmReaHes,HQR,RealVecHes, ComVecHes BakReaHes1,BalBak *)
(* zur Berechnung der Eigenwerte und -vektoren *)
(*----------------------------------------------------------------*)
PROCEDURE CHTriDi(VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
n : INTEGER;
VAR Tau : ARRAY OF LONGCOMPLEX;
VAR D,E : ARRAY OF LONGREAL;
VAR E2 : ARRAY OF LONGREAL);
(*---------------------------------------------------------------*)
(* Diese Routine ist eine Uebersetzung der Eispack Fortran *)
(* Subroutine htridi *)
(* *)
(* This routine is a Modula-2 translation of the eispack *)
(* subroutine htridi *)
(* *)
(* A : Complex hermitian matrix to be bi-diagonalized *)
(* N : Dimension of A *)
(* Tau : Complex vector of diagonal elements *)
(* D : diagonal of the calculated trilinear real matrix *)
(* E : subdiagonal of the calculated trilinear real matrix *)
(* E2 : squares of subdiagonal elements *)
(*---------------------------------------------------------------*)
PROCEDURE CHTriBak(VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR Tau : ARRAY OF LONGCOMPLEX;
N : INTEGER;
M : INTEGER;
VAR Z : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR iFehl : INTEGER);
(*---------------------------------------------------------------*)
(* Diese Routine ist eine Uebersetzung der Eispack Fortran *)
(* Subroutine htribk *)
(* *)
(* This routine is a Modula-2 translation of the eispack *)
(* subroutine htribk completly in direct complex arithmetic *)
(* *)
(* A : Matrix as transformed by routine CHTriDi *)
(* Tau : Vector as as calculeated by by routine CHTriDi *)
(* N : Dimension of A *)
(* M : Number of eigenvectors to be back-transformed *)
(* Z : eigenvectors of the trilinear matrix calculated by *)
(* ny CHTriDi (imaginary part is to be set to zero) *)
(* iFehl : Error indicator *)
(*---------------------------------------------------------------*)
PROCEDURE ZTransQL(VAR EW : ARRAY OF LONGREAL;
VAR ND : ARRAY OF LONGREAL;
dim : CARDINAL;
u,o : CARDINAL;
VAR Z : ARRAY OF ARRAY OF LONGCOMPLEX;
genau : LONGREAL;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Die Procedure TransQL erledigt die QL-Transformation der *)
(* Ttilinearmatrix (HD,ND) in den Grenzen u:o *)
(* *)
(* QL algorithm with implicit shifts for a real tridiaginal *)
(* symmetric matrix. *)
(* *)
(* input *)
(* *)
(* EW : is the vector of diagonal elements *)
(* ND : subdiagonal elements, e[1] is arbitrary *)
(* dim : dimension of original matrix *)
(* u,o : starting- and endpoint of matrix *)
(* Z : is identity matrix for initially tridiagonal *)
(* input matrix or z is output from tred2 or *)
(* similar routine for real symmetric matrix *)
(* genau : all elements > genau are rotates to "zero" *)
(* *)
(* output *)
(* *)
(* EW : contains eigenvalues *)
(* Z : The k-th column Z[k,i] returns the normalized *)
(* eigenvectors, corresponding to EW[k] *)
(* iFehl : Fehlercode *)
(* -2 : Speicherfehler *)
(* -1 : Matrix A wurde zerstoert *)
(* 0 : alles in Ordnung *)
(* >= 1 : Maximalzahl der Iterationen bei der *)
(* Berechnung des Eigenwerter mit Index *)
(* iFehl wurde ueberschritten - Eigenwerte *)
(* und Vektoren bis Index iFehl sollten in *)
(* Ordnung sein *)
(*----------------------------------------------------------------*)
PROCEDURE EigenCHM2(VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
dim : CARDINAL;
m : CARDINAL;
VAR EW : ARRAY OF LONGREAL;
VAR EV : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Berechnet die Eigenwerte- und Vektoren der hermitischen *)
(* Matrix A die als quadratisches Feld vorgegeben wird *)
(* *)
(* ==> A : Hermitische Matrix (2-dimensionales complexes *)
(* Feld. *)
(* ==> dim : Dimension von A *)
(* ==> m : Anzahl der zu berechnenden Eigenvektoren *)
(* (ungenutzt) *)
(* <== EW : Eigenwerte von A (reelwertig) *)
(* <== EV : Eigenvektoren von A (complxwertig) *)
(* <== iFehl : Fehlercode *)
(* -2 : Speicherfehler *)
(* -1 : Matrix A wurde zerstoert *)
(* 0 : alles in Ordnung *)
(* >= 1 : Maximalzahl der Iterationen bei der *)
(* Berechnung des Eigenwerter mit Index *)
(* iFehl wurde ueberschritten - Eigenwerte *)
(* und Vektoren bis Index iFehl sollten in *)
(* Ordnung sein *)
(*----------------------------------------------------------------*)
PROCEDURE CJacobi1D( N : INTEGER;
VAR A : ARRAY OF LONGCOMPLEX; (* SUPERVEKTOR *)
VAR EW : ARRAY OF LONGREAL;
VAR EV : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Calculate the eigensystem of a complex hermitian matrix *)
(* *)
(* This routine is based on an idea found on the web page of *)
(* Jean-Pierre Moreau (www.jpmoreau.fr) and the literture he *)
(* cites there for his Pascal routine EPHJ. *)
(* *)
(* Adopted to Modula-2 and heavily modified by M.Riedl *)
(* *)
(* ==> N : Dimension of matrix A *)
(* <==> A : Matrix to be diagonalized in packed form *)
(* <== EW : Eigenvalues of A *)
(* <== EV : Eigenvectors of A *)
(* <== iFehl : 0 if all fine *)
(* : 1 if N < 1 *)
(* : < 0 exeeded maximal number of iterations. The *)
(* negative number of iterations is returned in that *)
(* case *)
(* *)
(* Lit: Jardrin, Jean-Louis; "ALGEBRE: Algorithmes et programmes *)
(* en Pascal", Dunod, Paris (1988) *)
(* Sec. 6.7, pp. 209-218, "Calcul des elements propres d une *)
(* matrice hermitienne par la methode de Jacobi *)
(*----------------------------------------------------------------*)
PROCEDURE CJacobi2D( N : INTEGER;
VAR A : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR EW : ARRAY OF LONGREAL;
VAR EV : ARRAY OF ARRAY OF LONGCOMPLEX;
VAR iFehl : INTEGER);
(*----------------------------------------------------------------*)
(* Calculate the eigensystem of a complex hermitian matrix *)
(* *)
(* This routine is essential the same as CJacob1D but for *)
(* unpacked matrices. *)
(*----------------------------------------------------------------*)
END EigenLib2.