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.
