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.
