--- a
+++ b/EigenLib3.mod
@@ -0,0 +1,869 @@
+IMPLEMENTATION MODULE EigenLib3;
+
+  (*------------------------------------------------------------------------*)
+  (* Routinen zu Berechnung der Eigensysteme komplexer genereller Matrizen  *)
+  (* Routines to calculate the eigensystems of complex general matrices     *)
+  (*------------------------------------------------------------------------*)
+  (* Letzt Bearbeitung:                                                     *)
+  (*                                                                        *)
+  (* 07.08.93, MRi: Durchsicht                                              *)
+  (* 06.02.15, MRi: Ausgetestete Routinen aus EisPack herausgenommen und    *)
+  (*                in EigenLib3.mod integriert (CSortEwEv,ComEig,ComHess   *)
+  (*                CHessLR,CHessInvIter)                                   *)
+  (* 15.07.15, MRi: CSortEwEV,RHessQR,CHessLR auf ISO COMPLEX umgestellt    *)
+  (* 18.07.15, MRi: Umstellung vom ComEig auf ISO-COMPLEX                   *)
+  (* 01.04.16, MRi: Aufraeumen wg. Compilerwarnungen, kleine Code-          *)
+  (*                veraenderung in CHessInvIter bei der Bestimmung von     *)
+  (*                Norm & VNorm                                            *)
+  (* 20.09.17, MRi: Umstellung vom ComEig auf offene Felder                 *)
+  (* 21.09.17, MRi: Umstellen von ComHess und CHessLR auf ISO-COMPLX        *)
+  (* 23.09.17, MRi: Umstellen von CHessInvIter auf ISO-COMPLX (neu)         *)
+  (* 28.09.17, MRi: CBalance und CBalBak aufgenommen. In CBalance den       *)
+  (*                Parameter iFehl ergaenzt                                *)
+  (* 19.10.17, MRi: In ComEig die Transformationsschleifen von Komponenten- *)
+  (*                arithmetik zur expliziter komplexer Arithmetik umge-    *)
+  (*                stellt, ebenfalls die Berechnung von H                  *)
+  (* 22.10.17, MRi: In ComEig MaxIter fix auf 70 begrenzt                   *)
+  (*------------------------------------------------------------------------*)
+  (* Offene Punkte                                                          *)
+  (*                                                                        *)
+  (* - Umstellen auf offene Felder                                          *)
+  (* - ComEig auf komplexe Arithmethikg durchsehen und vereinfachen         *)
+  (* - Paramer "select" f��r CHessInvIter und ComEigEvEw einfuehren ?        *)
+  (* - CHessInvIter laeuft fuer einen der beiden Eigenwert (1.0,1.0) aus    *)
+  (*   Num. Math. 16, 181-204 (1970) - BSP 2 auf einen Fehler - warum ???   *)
+  (* - CBalance liefert fuer "Test unbekannter Herkunft" unsinnige Werte    *)
+  (*------------------------------------------------------------------------*)
+  (* Testprogramm ist TstCRM.mod                                            *)
+  (*------------------------------------------------------------------------*)
+  (* Implementation : Michael Riedl                                         *)
+  (* Licence        : GNU Lesser General Public License (LGPL)              *)
+  (*------------------------------------------------------------------------*)
+
+  (* $Id: EigenLib3.mod,v 1.4 2017/10/22 17:52:02 mriedl Exp mriedl $ *)
+
+FROM SYSTEM          IMPORT TSIZE,ADR;
+FROM Storage         IMPORT ALLOCATE,DEALLOCATE;
+FROM Errors          IMPORT Fehler,Fehlerflag,ErrOut,FatalError;
+                     IMPORT Errors;
+FROM Deklera         IMPORT MaxDim,VEKTOR,CARDVEKTOR,CVEKTOR,CMATRIX;
+FROM LongMath        IMPORT sqrt;
+                     IMPORT LongComplexMath;
+FROM LongComplexMath IMPORT zero,one,scalarMult;
+FROM LMathLib        IMPORT MachEps,MachEpsR4;
+FROM EigenLibAux     IMPORT CNormMat;
+                     IMPORT TIO;
+
+CONST debug      = FALSE;
+
+(*================= Hilfsroutinen ==========================================*)
+
+PROCEDURE CEucNorm(VAR X    : ARRAY OF LONGCOMPLEX; (* unchanged *)
+                       n    : INTEGER;
+                   VAR Norm : LONGREAL);
+
+          (*--------------------------------------------------------*)
+          (* Berechnet die euclidische Norm des komplexen Vektors X *)
+          (* Modula-2 Transkription der LAPACK- Routine dznrm2      *)
+          (*--------------------------------------------------------*)
+
+          VAR scale,ssq,temp,t1 : LONGREAL;
+              i                 : INTEGER;
+BEGIN
+      IF (n < 1) THEN
+        Norm := 0.0;
+      ELSE
+        scale := 0.0;
+        ssq   := 1.0;
+        FOR i:=0 TO n-1 DO
+          IF (RE(X[i]) # 0.0)THEN
+            temp := ABS(RE(X[i]));
+            IF (scale < temp) THEN
+              t1    := scale / temp;
+              ssq   := 1.0 + ssq*t1*t1;
+              scale := temp;
+            ELSE
+              t1  := temp / scale;
+              ssq :=ssq + t1*t1;
+            END;
+          END;
+          IF (IM(X[i]) # 0.0) THEN
+            temp := ABS(IM(X[i]));
+            IF (scale < temp) THEN
+              t1    := scale / temp;
+              ssq   := 1.0 + ssq*t1*t1;
+              scale := temp;
+            ELSE
+              t1 := temp/scale;
+              ssq := ssq+t1*t1;
+            END;
+          END;
+        END; (* FOR i *)
+        Norm  := scale*sqrt(ssq);
+      END; (* IF *)
+END CEucNorm;
+
+PROCEDURE CNormMat2(   dim  : CARDINAL;
+                   VAR CMat : CMATRIX);
+
+          VAR Norm    : LONGREAL;
+              i,j     : CARDINAL;
+BEGIN
+      FOR i:=1 TO dim DO
+        CEucNorm(CMat[i],dim,Norm);
+        IF (Norm # 0.0) THEN
+          Norm := 1.0 / Norm;
+        END;
+        FOR j:=1 TO dim DO
+          CMat[i,j] := scalarMult(Norm,CMat[i,j]);
+        END;
+      END;
+END CNormMat2;
+
+(*================= Ende Hilfsroutinen =====================================*)
+
+PROCEDURE CSortEwEv(VAR EW  : CVEKTOR;  (* Eigenwerte    *)
+                    VAR EV  : CMATRIX;  (* Eigenvektoren *)
+                        dim : CARDINAL;
+                        ord : INTEGER);
+
+          VAR  i,l,iminmax : CARDINAL;
+               MinMax      : LONGREAL;
+               Zw          : LONGCOMPLEX;
+               AbsEW       : VEKTOR;
+BEGIN
+      FOR i:=1 TO dim DO AbsEW[i]:=LongComplexMath.abs(EW[i]); END;
+      FOR l:=1 TO dim-1 DO
+        MinMax:=AbsEW[l]; iminmax:=l;
+        FOR i:=l+1 TO dim DO
+          IF (ord = 1) THEN  (* Berechne Maximum von AbsEW[dim] *)
+            IF (AbsEW[i] > MinMax) THEN  MinMax:=AbsEW[i]; iminmax:=i; END;
+          ELSE               (* Berechne Minimum von AbsEW[dim] *)
+            IF (AbsEW[i] < MinMax) THEN  MinMax:=AbsEW[i]; iminmax:=i; END;
+          END; (* IF Ordnung *)
+        END; (* FOR i *)
+        IF (iminmax <> l) THEN (* Vertausche EWi mit EWimax *)
+          Zw:=EW[l]; EW[l]:=EW[iminmax]; EW[iminmax]:=Zw;
+          AbsEW[iminmax]:=AbsEW[l];
+          FOR i:=1 TO dim DO (* Vertausche Eigenvektoren entsprechend *)
+            Zw:=EV[l,i]; EV[l,i]:=EV[iminmax,i]; EV[iminmax,i]:=Zw;
+          END; (* i *)
+        END; (* IF *)
+      END; (* l *)
+END CSortEwEv;
+
+PROCEDURE CBalance(VAR A        : ARRAY OF ARRAY OF LONGCOMPLEX;
+                       dim      : CARDINAL;
+                   VAR low,high : CARDINAL;
+                   VAR Scale    : ARRAY OF LONGREAL;
+                   VAR iFehl    : INTEGER);
+
+          VAR j,k,l,N : INTEGER;
+
+          PROCEDURE exc(m : INTEGER); (* lokal zu CBalance *)
+
+                    VAR t : LONGCOMPLEX;
+                        i : INTEGER;
+          BEGIN
+                Scale[m]:=LFLOAT(j);
+                IF (j # m) THEN
+                  FOR i:=0 TO k   DO t:=A[i,j]; A[i,j]:=A[i,m]; A[i,m]:=t; END;
+                  FOR i:=l TO N-1 DO t:=A[j,i]; A[j,i]:=A[m,i]; A[m,i]:=t; END;
+                END;
+          END exc;
+
+          CONST radix = 2.0; (* Basis der Flie\3kommadarstellung *)
+                sqrdx = radix*radix;
+
+          VAR   s,r,g,f,c       : LONGREAL;
+                i               : INTEGER;
+                Konvergenz,cont,flag,found : BOOLEAN;
+BEGIN
+      N:=VAL(INTEGER,dim);
+      iFehl:=0;
+
+      l := 0;
+      k := N-1;
+      REPEAT
+        (* Search for rows isolating an eigenvalue and push them down *)
+        j := k;
+        REPEAT
+          i := 0;
+          REPEAT
+            flag := (i # j) AND (A[j,i] # zero);
+            INC(i);
+          UNTIL flag OR (i > k);
+          found := NOT flag;
+          IF found THEN
+            exc(k);
+            DEC(k)
+          END;
+          DEC(j);
+        UNTIL found OR (j < 0);
+      UNTIL (NOT found) OR (k < 0);
+
+      IF (k < 0) THEN k := 0; END;
+      IF (k = 0) THEN
+        Errors.ErrOut("high = 0 (CBalance)");
+        low :=1;
+        high:=VAL(CARDINAL,k+1);
+        iFehl:=-1;
+        RETURN;
+      END;
+
+      REPEAT
+        (* Search for columns isolating an eigenvalue and push them left *)
+        j := l;
+        REPEAT
+          i := l;
+          REPEAT
+            flag := (i # j) AND (A[i,j] # zero);
+            INC(i);
+          UNTIL flag OR (i > k);
+          found := NOT flag;
+          IF found THEN
+            exc(l);
+            INC(l);
+          END;
+          INC(j);
+        UNTIL found OR (j > k);
+      UNTIL (NOT found);
+
+      low :=VAL(CARDINAL,l+1);
+      high:=VAL(CARDINAL,k+1);
+
+      (* Balanciere die Submatrix A[i,l..k] *)
+      FOR i:=l TO k DO Scale[i]:=1.0; END;
+      REPEAT
+        Konvergenz:=TRUE; (* Wird getestet *)
+        FOR i:=l TO k DO
+          c:=0.0; r:=0.0;
+          FOR j:=l TO k DO
+            IF (j # i) THEN
+              c:=c + ABS(RE(A[j,i])) + ABS(IM(A[j,i]));
+              r:=r + ABS(RE(A[i,j])) + ABS(IM(A[i,j])); (* ??? *)
+            END;
+          END;
+          IF (c # 0.0) AND (r # 0.0) THEN
+            g := r / radix;
+            f := 1.0;
+            s := c + r;
+            WHILE (c < g)  DO f:=f * radix; c:=c * sqrdx; END;
+            g := r*radix;
+            WHILE (c >= g) DO f:=f / radix; c:=c / sqrdx; END;
+            IF (((c+r) / f) < 0.95*s) THEN
+              Konvergenz:=FALSE;
+              g := 1.0 / f;
+              Scale[i]:=Scale[i] * f;
+              FOR j:=l TO N-1 DO A[i,j]:=scalarMult(g,A[i,j]); END;
+              FOR j:=0 TO k   DO A[j,i]:=scalarMult(f,A[j,i]); END;
+            END;
+          END; (* IF *)
+        END; (* FOR i *)
+      UNTIL Konvergenz;
+END CBalance;
+
+PROCEDURE CBalBak(    dim      : CARDINAL; (* Ordnung der Matrix Z *)
+                      low,high : CARDINAL;
+                  VAR Scale    : ARRAY OF LONGREAL;
+                      m        : CARDINAL;
+                  VAR Z        : ARRAY OF ARRAY OF LONGCOMPLEX);
+
+          VAR  i,j,k : CARDINAL;
+               S     : LONGREAL;
+               Zw    : LONGCOMPLEX;
+
+BEGIN
+      IF (m = 0) THEN RETURN END;
+
+      DEC(low); DEC(high);
+
+      IF (low # high) THEN
+        FOR i:=low TO high DO
+          S:=Scale[i]; (* Linksvektoren durch S:=1.0 / Scale[i] ! *)
+          FOR j:=0 TO m-1 DO
+            Z[j,i]:=scalarMult(S,Z[j,i]);
+          END;
+        END;
+      END;
+
+      (* for i=low-1 step -1 until 1,high+1 step 1 until dim do *)
+
+      IF (low > 0) THEN
+        FOR i:=low-1 TO 0 BY -1 DO
+          k := TRUNC(Scale[i]);
+          IF (k # i) THEN
+            FOR j:=0 TO m-1 DO Zw:=Z[j,i]; Z[j,i]:=Z[j,k]; Z[j,k]:=Zw; END;
+          END;
+        END;
+      END;
+      FOR i:=high+1 TO dim-1 DO
+        k := TRUNC(Scale[i]);
+        IF (k # i) THEN
+          FOR j:=0 TO m-1 DO Zw:=Z[j,i]; Z[j,i]:=Z[j,k]; Z[j,k]:=Zw; END;
+        END;
+      END;
+END CBalBak;
+
+PROCEDURE ComHess(VAR A     : CMATRIX;
+                      dim   : CARDINAL;
+                      k,l   : CARDINAL;
+                  VAR Index : CARDVEKTOR);
+
+          VAR i,j,m : CARDINAL;
+              X,Y   : LONGCOMPLEX;
+BEGIN
+      Index[k]:=k; Index[l]:=l;
+      FOR m:=k+1 TO l-1 DO
+        i:=m; X:=CMPLX(0.0,0.0);
+        FOR j:=m TO l DO
+          IF ABS(RE(A[j,m-1]))+ABS(IM(A[j,m-1])) > ABS(RE(X))+ABS(IM(X)) THEN
+            X:=A[j,m-1]; i:=j;
+          END;
+        END;
+        Index[m]:=i;
+        IF (i # m) THEN (* interchange rows and columns of arrays a *)
+          FOR j:=m-1 TO dim DO Y:=A[i,j]; A[i,j]:=A[m,j]; A[m,j]:=Y; END;
+          FOR j:=1   TO l   DO Y:=A[j,i]; A[j,i]:=A[j,m]; A[j,m]:=Y; END;
+        END;
+        IF (RE(X) # 0.0) OR (IM(X) # 0.0) THEN
+          FOR i:=m+1 TO l DO
+            Y:=A[i,m-1];
+            IF (RE(Y) # 0.0) OR (IM(Y) # 0.0) THEN
+              Y:=Y / X; A[i,m-1]:=Y;
+              FOR j:=m TO dim DO A[i,j]:=A[i,j] - Y*A[m,j]; END;
+              FOR j:=1 TO l   DO A[j,m]:=A[j,m] + Y*A[j,i]; END;
+            END; (* IF Y *)
+          END; (* FOR i *)
+        END; (* IF X *)
+      END; (* FOR m *)
+END ComHess;
+
+PROCEDURE ComHessBak(VAR EV    : CMATRIX;
+                     VAR A     : CMATRIX;
+                         Index : CARDVEKTOR;
+                         dim   : CARDINAL;
+                         EvAnz : CARDINAL;
+                         k,l   : CARDINAL);
+
+          VAR  i,j,m : CARDINAL;
+               X     : LONGCOMPLEX;
+BEGIN
+      IF (l > dim) THEN
+        Errors.FatalError("upp > dim (ComHessBak)");
+      END;
+      INC(k);
+      FOR m:=l-1 TO k BY -1 DO
+        FOR i:=m+1 TO l DO
+          X:=A[i,m-1];
+          IF (RE(X) # 0.0) OR (IM(X) # 0.0) THEN
+            FOR j:=1 TO EvAnz DO
+              EV[j,i]:=EV[j,i] + X*EV[j,m];
+            END;
+          END;
+        END; (* FOR i *)
+        i:=Index[m];
+        IF (i # m) THEN
+          FOR j:=1 TO EvAnz DO X:=EV[j,i]; EV[j,i]:=EV[j,m]; EV[j,m]:=X; END;
+        END;
+      END; (* FOR m *)
+END ComHessBak;
+
+PROCEDURE CHessLR(VAR H     : CMATRIX;   (* komplexe Hessenbergmatrix     *)
+                  VAR EW    : CVEKTOR;   (* komplexe Eigenwerte von H     *)
+                      dim   : CARDINAL); (* Dimension von H               *)
+
+          CONST MaxIter = 30;
+
+          VAR   n,i,j,l,m,ll,mm,Iter : CARDINAL;
+                S,T,X,Y,Z            : LONGCOMPLEX;
+                x,y,tmp              : LONGREAL;
+BEGIN
+      X:=zero; Y:=zero; Z:=zero;
+      T:=zero; n:=dim;
+      WHILE (n > 0) DO
+        Iter:=0;
+        LOOP
+          LOOP (* Suche kleine Subdiagonalelemente *)
+            FOR l:=n TO 2 BY -1 DO
+              x:=ABS(RE(H[l  ,l-1])) + ABS(IM(H[l  ,l-1]));
+              y:=ABS(RE(H[l-1,l-1])) + ABS(IM(H[l-1,l-1])) +
+                 ABS(RE(H[l  ,l  ])) + ABS(IM(H[l  ,l  ]));
+              IF (x <= MachEps*y) THEN ll:=l; EXIT; END;
+            END;
+            ll:=1; EXIT;
+          END;
+          l:=ll;
+          IF (l = n) THEN EXIT END; (* Wurzel gefunden ! *)
+
+          IF (Iter = MaxIter) THEN
+            Fehler:=TRUE; Fehlerflag:='MaxIter ueberschritten (CHessLR)';
+            RETURN;
+          END;
+          IF ((Iter MOD 10) = 0) AND (n > 2) THEN
+            S:=CMPLX(ABS(RE(H[n,n-1])) + ABS(RE(H[n-1,n-2])),
+                     ABS(IM(H[n,n-1])) + ABS(IM(H[n-1,n-2])));
+          ELSE
+            S := H[n,n];
+            X := H[n-1,n]*H[n,n-1];
+            IF (RE(X) # 0.0) OR (IM(X) # 0.0) THEN
+              Y:=scalarMult(0.5,(H[n-1,n-1] - S));
+              Z:=CMPLX(
+                   RE(Y)*RE(Y) - IM(Y)*IM(Y) + RE(X),
+                             2.0*RE(Y)*IM(Y) + IM(X)
+                 );
+              Z:=LongComplexMath.sqrt(Z);
+              IF ((RE(Y)*RE(Z) + IM(Y)*IM(Z)) < 0.0) THEN
+                Z:= - Z;
+              END;
+              X:=X / (Y + Z);
+              S:=S - X;
+            END;
+          END; (* IF Iter *)
+
+          FOR i:=1 TO n DO H[i,i]:=H[i,i] - S; END;
+
+          (* Suche zwei aufeinander folgende kleine Subdiagonalelemente *)
+
+          X := CMPLX(ABS(RE(H[n-1,n-1])) + ABS(IM(H[n-1,n-1])),IM(X));
+          Y := CMPLX(ABS(RE(H[n  ,n-1])) + ABS(IM(H[n  ,n-1])),IM(Y));
+          Z := CMPLX(ABS(RE(H[n  ,n  ])) + ABS(IM(H[n  ,n  ])),IM(Z));
+
+          T:=T + S;
+          j:=l+1;
+          INC(Iter); (* !!! *)
+          LOOP
+            FOR m:=n-1 TO j BY -1 DO
+              Y := CMPLX(ABS(RE(H[m  ,m-1])) + ABS(IM(H[m  ,m-1])),RE(Y));
+              tmp := RE(X);
+              X := CMPLX(ABS(RE(H[m-1,m-1])) + ABS(IM(H[m-1,m-1])),RE(Z));
+              Z := CMPLX(tmp,IM(Z));
+              IF (RE(Y) < MachEps*(RE(Z) + RE(X) + IM(X))*(RE(Z)/IM(Y))) THEN
+                mm:=m; EXIT;
+              END;
+            END;
+            mm:=l; EXIT;
+          END;
+          m:=mm;
+          FOR i:=m+1 TO n DO
+            X:=H[i-1,i-1]; Y:=H[i,i-1];
+            IF ABS(RE(X)) + ABS(IM(X)) < ABS(RE(Y)) + ABS(IM(Y)) THEN
+              FOR j:=i-1 TO n DO Z:=H[i-1,j]; H[i-1,j]:=H[i,j]; H[i,j]:=Z; END;
+              Z := X / Y; EW[i]:=CMPLX(+1.0,IM(EW[i]));
+            ELSE
+              Z := Y / X; EW[i]:=CMPLX(-1.0,IM(EW[i]));
+            END;
+            H[i,i-1]:=Z;
+            FOR j:=i TO n DO H[i,j]:=H[i,j] - Z*H[i-1,j]; END;
+          END; (* FOR i *)
+          (* H = R*L *)
+          FOR j:=m+1 TO n DO
+            X:=H[j,j-1]; H[j,j-1]:=zero;
+            IF (RE(EW[j]) > 0.0) THEN
+              FOR i:=l TO j DO Z:=H[i,j-1]; H[i,j-1]:=H[i,j]; H[i,j]:=Z; END;
+            END;
+            FOR i:=l TO j DO H[i,j-1]:=H[i,j-1] + X*H[i,j]; END;
+          END; (* FOR j *)
+        END; (* LOOP *)
+        EW[n]:=H[n,n] + T;
+        DEC(n);
+      END; (* WHILE n > 0 *)
+END CHessLR;
+
+PROCEDURE CHessInvIter(VAR EV    : CMATRIX;   (* Eigenvektoren    *)
+                       VAR EW    : CVEKTOR;   (* Eigenwerte von A *)
+                       VAR H     : CMATRIX;   (* Hessenbergmatrix *)
+                           dim   : CARDINAL;  (* Dimension von H  *)
+                       VAR IErr  : ARRAY OF CARDINAL;
+                       VAR iFehl : INTEGER);
+
+          CONST MinEps = 100.0*MachEps;
+
+          VAR   X,Y,eps3,Norm,NormV,GrowTo,SqrtUk : LONGREAL;
+                i,j,k,m,uk,Iter,ukk               : CARDINAL;
+                Zw,Lambda                         : LONGCOMPLEX;
+                B                                 : POINTER TO CMATRIX;
+                EVk                               : POINTER TO CVEKTOR;
+                gefunden                          : BOOLEAN;
+BEGIN
+      NEW(B);
+      IF (B = NIL) THEN
+        FatalError("kein Freispeicher vorhanden (EigenLib3.CHessInvIter) !");
+      END;
+      eps3:=0.0; GrowTo:=0.0; (* Wg. Compilerwarnungen *)
+
+      iFehl:=0; uk:=0; Fehler:=FALSE;
+      FOR k:=1 TO dim DO
+        EVk := ADR(EV[k]); (* Zeiger auf k-ten Eigenvektor *)
+        IF (uk < k) THEN
+          LOOP (* split *)
+            FOR uk:=k TO dim-1 DO
+              IF (RE(H[uk+1,uk]) = 0.0) AND (IM(H[uk+1,uk]) = 0.0) THEN
+                ukk:=uk;
+                EXIT;
+              END;
+            END;
+            ukk:=dim;
+            EXIT;
+          END;
+          uk:=ukk;
+          Norm:=0.0; m:=1;
+          FOR i:=1 TO uk DO (* norm of leading uk X uk matrix *)
+            X:=0.0;
+            FOR j:=m TO uk DO X:=X + LongComplexMath.abs(H[i,j]); END;
+            IF (X > Norm) THEN Norm:=X; END;
+            m:=i;
+          END; (* FOR *)
+
+          IF (Norm = 0.0) THEN Norm:=1.0; END;
+          eps3:=MachEps*Norm;
+          SqrtUk := sqrt(VAL(LONGREAL,uk));
+          GrowTo := 0.1 / SqrtUk
+        END; (* IF uk *)
+
+        Lambda := EW[k];
+        (* perturb eigenvalue if it is close to any previous eigenvalue *)
+        IF (k > 1) THEN
+          LOOP (* test *)
+            i:=k-1;
+            IF (i = 0) THEN EXIT; END;
+            REPEAT
+              gefunden:= (ABS(RE(EW[i]) - RE(Lambda)) < eps3) AND
+                         (ABS(IM(EW[i]) - IM(Lambda)) < eps3);
+              DEC(i);
+            UNTIL gefunden OR (i = 0);
+            IF gefunden THEN
+              Lambda:=CMPLX(RE(Lambda) + eps3,IM(Lambda)); (* GOTO test: *)
+            ELSE
+              EXIT;
+            END;
+          END; (* LOOP *)
+          EW[k] := CMPLX(RE(Lambda),IM(EW[k]));;
+        END; (* IF *)
+
+        (* generate B = A -lambda x I, the matrix to be reduced, *)
+        (* and initial vector u + iv (EVk)                       *)
+        m:=1;
+        FOR i:=1 TO uk DO
+          FOR j:=m TO uk DO
+            B^[i,j] := H[i,j];
+          END;
+          B^[i,i]:=B^[i,i] - Lambda;
+          m:=i;
+          EVk^[i]:=CMPLX(eps3,0.0);
+        END; (* FOR *)
+
+        IF (uk # 1) THEN
+          (* reduce B to LU form *)
+          FOR i:=2 TO uk DO
+            m := i - 1;
+            IF (LongComplexMath.abs(B^[i,m])   >
+                LongComplexMath.abs(B^[m,m]))
+            THEN
+              FOR j:=m TO uk DO
+                Zw:=B^[i,j]; B^[i,j]:=B^[m,j]; B^[m,j]:=Zw;
+              END;
+            END; (* IF *)
+
+            IF (ABS(RE(B^[m,m])) <= MinEps) AND
+               (ABS(IM(B^[m,m])) <= MinEps) THEN
+              B^[m,m]:=CMPLX(eps3,IM(B^[m,m]));
+            END;
+            Zw := B^[i,m] / B^[m,m];
+            IF ((RE(Zw) # 0.0) OR (IM(Zw) # 0.0)) THEN
+              FOR j:=i TO uk DO
+                B^[i,j]:=B^[i,j] - Zw*B^[m,j];
+              END;
+            END; (* IF Zw = 0.0 *)
+          END; (* FOR i *)
+        END; (* IF uk # 1 *)
+        IF (RE(B^[uk,uk]) = 0.0) AND (IM(B^[uk,uk]) = 0.0) THEN
+          B^[uk,uk]:=CMPLX(eps3,0.0);
+        END;
+        Iter:=0;
+
+        LOOP (* back substitution *)
+          FOR i:=uk TO 1 BY -1 DO
+            Zw:=CMPLX(RE(EVk^[i]),0.0);
+            IF (i # uk) THEN
+              FOR j:=i+1 TO uk DO
+                Zw:=Zw - B^[i,j]*EVk^[j];
+              END; (* FOR *)
+            END; (* IF i # uk *)
+            EVk^[i] := Zw / B^[i,i];
+          END; (* FOR i *)
+
+          INC(Iter);
+
+          Norm:=0.0; NormV:=Norm; j:=1;
+          FOR i:=1 TO uk DO
+            X:= LongComplexMath.abs(EVk^[i]);
+            IF (NormV < X) THEN NormV:=X; j:=i; END;
+            Norm:=Norm + X;
+          END;
+
+          IF (Iter > 1) AND (Norm > GrowTo) THEN (* Akzeptabler Eigenvektor *)
+            Zw:=EVk^[j];
+            FOR i:=1    TO uk  DO EVk^[i]:=EVk^[i] / Zw;   END;
+            FOR i:=uk+1 TO dim DO EVk^[i]:=CMPLX(0.0,0.0); END;
+            j := uk + 1;
+            EXIT; (* !!! *)
+          END; (* IF *)
+
+          IF (Iter >= uk) THEN
+            IErr[iFehl]:=k; INC(iFehl); Fehler:=TRUE;
+            Fehlerflag:=
+            "nicht akzeptabler Eigenvektor (EigenLib3.CHessInvIter) !";
+            Errors.ErrOut(Fehlerflag);
+            Errors.WriteString(" uk = ");
+            Errors.WriteCard(uk);
+            Errors.WriteLn;
+            FOR i:=1 TO dim DO EVk^[i]:=CMPLX(0.0,0.0); END;
+            EXIT;
+          END;
+          (* Neuer Testvektor *)
+          X := sqrt(VAL(LONGREAL,dim));
+          Y := eps3 / (X + 1.0);
+          EVk^[1] := CMPLX(eps3,IM(EVk^[1]));
+          FOR i:=2 TO dim DO EVk^[i]:=CMPLX(Y,IM(EVk^[i])); END;
+          EVk^[j]:=EVk^[j] - CMPLX(eps3,0.0);
+        END; (* LOOP *)
+      END; (* FOR k *)
+      DISPOSE(B);
+END CHessInvIter;
+
+PROCEDURE ComEigEvEw(VAR A      : CMATRIX;
+                         dim    : CARDINAL;
+                     VAR EW     : CVEKTOR;
+                     VAR EV     : CMATRIX;
+                         Norm   : CARDINAL;
+                         ord    : INTEGER;
+                     VAR iFehl  : INTEGER);
+
+          VAR i,j,u,low,upp : CARDINAL;
+              Index,IErr    : CARDVEKTOR;
+              H             : CMATRIX;
+
+BEGIN
+          low := 1; upp := dim;
+
+          FOR i:=1 TO dim DO Index[i]:=i; END;
+
+          ComHess(A,dim,low,upp,Index);
+
+          u:=1;
+          FOR i:=1 TO dim DO
+            FOR j:=1 TO u-1 DO H[i,j]:=CMPLX(0.0,0.0); END;
+            FOR j:=u TO dim DO H[i,j]:=A[i,j]; END;
+            IF (i > 2) THEN INC(u); END;
+          END;
+          CHessLR(H,EW,dim);
+
+          CHessInvIter(EV,EW,A,dim,IErr,iFehl);
+
+          ComHessBak(EV,A,Index,dim,dim,low,upp);
+
+          IF (Norm = 1) THEN
+            CNormMat(EV,dim,TRUE);
+          ELSE
+            CNormMat2(dim,EV);
+          END;
+
+          IF (ord # 0) THEN
+            IF (ord > 0) THEN ord:=1; ELSE ord:=-1; END;
+            CSortEwEv(EW,EV,dim,ord);
+          END;
+END ComEigEvEw;
+
+PROCEDURE ComEig(VAR A     : ARRAY OF ARRAY OF LONGCOMPLEX;
+                 VAR REV   : ARRAY OF ARRAY OF LONGCOMPLEX; (* Rechts-EV  *)
+                 VAR LEV   : ARRAY OF ARRAY OF LONGCOMPLEX; (* Links-EV   *)
+                 VAR EW    : ARRAY OF LONGCOMPLEX;          (* Eigenwerte *)
+                     dim   : CARDINAL;
+                     genau : LONGREAL;
+                     ord   : INTEGER);
+
+          VAR    i,j,k,m,Iter,MaxIter                        : CARDINAL;
+                 eta,tse,tem,tep,tau,Max,Root1,Root2,Root    : LONGREAL;
+                 Hi,Hr,Hj,g,Te,Tee,Br,Bi,Er,Ei,Dr,Di,sig,isw : LONGREAL;
+                 b,c,Ca,Cb,Ch,Cx,nC,e,s,Sa,Sb,Sh,Sx,d,De,nD  : LONGREAL;
+                 Cot2X,CotX,Cos2a,Sin2a,TanH                 : LONGREAL;
+                 Aik,Aki,Aim,Ami,EVim,EVik,Zw,s1,s2,c1,c2,H  : LONGCOMPLEX;
+                 s1re,s1im,s2re,s2im                         : LONGREAL;
+                 re,im                                       : LONGREAL;
+                 Mark                                        : BOOLEAN;
+                 En  : POINTER TO ARRAY [0..MAX(INTEGER)-1] OF LONGREAL;
+BEGIN
+      ALLOCATE(En,dim*TSIZE(LONGREAL));
+      IF (genau > MachEpsR4) THEN genau:=MachEpsR4; END;
+      IF (genau < MachEps  ) THEN genau:=MachEps;   END;
+
+      MaxIter := 70; (* Sollte reichen, wie im Originalatikel *)
+
+      FOR i:=0 TO dim-1 DO            (* Vorbelege Eigenvektormatritzen *)
+        FOR j:=0 TO dim-1 DO REV[i,j]:=zero; LEV[i,j]:=zero; END;
+        REV[i,i]:=one; LEV[i,i]:=one;
+      END;
+
+      Mark:=FALSE; Fehler:=FALSE; Iter:=0;
+      LOOP
+        INC(Iter);
+        IF Mark OR (Iter > MaxIter) THEN EXIT END; (* !!! *)
+        tau:=0.0;
+        FOR k:=0 TO dim-1 DO
+          tem:=0.0;
+          FOR i:=0 TO dim-1 DO
+            IF (i # k) THEN tem:=tem + ABS(RE(A[i,k])) + ABS(IM(A[i,k])); END;
+          END;
+          tau:=tau + tem;
+          En^[k]:=tem + ABS(RE(A[k,k])) + ABS(IM(A[k,k]));
+        END; (* FOR k *)
+
+        FOR k:=0 TO dim-2 DO (* Vertausche Zeilen und Spalten *)
+          Max:=En^[k]; i:=k;
+          FOR j:=k+1 TO dim-1 DO
+            IF (En^[j] > Max) THEN Max:=En^[j]; i:=j; END;
+          END;
+          IF (i # k) THEN
+            En^[i]:=En^[k];
+            FOR j:=0 TO dim-1 DO Zw:=A[k,j]; A[k,j]:=A[i,j]; A[i,j]:=Zw; END;
+            FOR j:=0 TO dim-1 DO
+              Zw:=A  [j,k]; A  [j,k]:=A  [j,i]; A  [j,i]:=Zw;
+              Zw:=REV[k,j]; REV[k,j]:=REV[i,j]; REV[i,j]:=Zw;
+              Zw:=LEV[k,j]; LEV[k,j]:=LEV[i,j]; LEV[i,j]:=Zw;
+            END;
+          END; (* IF i # k *)
+        END; (* FOR k *)
+
+        IF debug THEN
+          TIO.WrStr(" ComEig MaxIter,Iter, tau = ");
+          TIO.WrCard(Iter,4); 
+          TIO.WrLngReal(tau,12,-5);
+          TIO.WrLn;
+        END;
+
+        IF (tau < genau) THEN EXIT END; (* !!! *)
+
+        Mark:=TRUE;
+        FOR k:=0 TO dim-1-1 DO
+          FOR m:=k+1 TO dim-1 DO
+            Hj:=0.0; g:=0.0; H:=zero;
+            FOR i:=0 TO dim-1 DO
+              IF (i # k) AND (i # m) THEN
+                H:=H + A[k,i]*A[m,i] - A[i,k]*A[i,m];
+                Te :=    RE(A[i,k])*RE(A[i,k]) + IM(A[i,k])*IM(A[i,k])
+                       + RE(A[m,i])*RE(A[m,i]) + IM(A[m,i])*IM(A[m,i]);
+                Tee:=    RE(A[i,m])*RE(A[i,m]) + IM(A[i,m])*IM(A[i,m])
+                       + RE(A[k,i])*RE(A[k,i]) + IM(A[k,i])*IM(A[k,i]);
+                g  :=    g  + Te + Tee;
+                Hj :=    Hj - Te + Tee;
+              END; (* IF *)
+            END; (* FOR i *)
+            Hr:=RE(H); Hi:=IM(H);
+            Br:=RE(A[k,m]) + RE(A[m,k]); Bi:=IM(A[k,m]) + IM(A[m,k]);
+            Er:=RE(A[k,m]) - RE(A[m,k]); Ei:=IM(A[k,m]) - IM(A[m,k]);
+            Dr:=RE(A[k,k]) - RE(A[m,m]); Di:=IM(A[k,k]) - IM(A[m,m]);
+            Te  := Br*Br + Ei*Ei + Dr*Dr;
+            Tee := Bi*Bi + Er*Er + Di*Di;
+            IF (Te >= Tee) THEN
+              isw:= 1.0; c:=Br; s:= Ei; d:=Dr; De:=Di; Root2:=sqrt(Te);
+            ELSE
+              isw:=-1.0; c:=Bi; s:=-Er; d:=Di; De:=Dr; Root2:=sqrt(Tee);
+            END;
+            Root1:= LongComplexMath.abs(CMPLX(s,c)); (* M.R. *)
+            Sa:=0.0;
+            IF (c >= 0.0) THEN Ca :=1.0; ELSE Ca :=-1.0; END;
+            IF (d >= 0.0) THEN sig:=1.0; ELSE sig:=-1.0; END;
+            IF (Root1 < MachEps) THEN
+              Sx:=0.0; Sa:=0.0; Cx:=1.0; Ca:=1.0;
+              IF (isw > 0.0) THEN e:=Er; b:=Bi; ELSE e:=Ei; b:=-Br; END;
+              nD:=d*d + De*De;
+              (* GOTO enter1; *)
+            ELSE
+              IF (ABS(s) > MachEps) THEN Ca := c / Root1; Sa := s / Root1; END;
+              Cot2X:= d / Root1;
+              CotX := Cot2X + (sig*sqrt(1.0 + Cot2X*Cot2X));
+              Sx   := sig / sqrt(1.0 + CotX*CotX);
+              Cx   := Sx*CotX;
+              eta  := (Er*Br + Bi*Ei) / Root1;
+              tse  := (Br*Bi - Er*Ei) / Root1;
+              Te   := sig*( - Root1*De + tse*d) / Root2;
+              Tee  := (d*De + Root1*tse) / Root2;
+              nD   := Root2*Root2 + Tee*Tee;
+              Tee  := Hj*Cx*Sx;
+              Cos2a:= Ca*Ca - Sa*Sa;
+              Sin2a:= 2.0*Ca*Sa;
+              tem  := Hr*Cos2a + Hi*Sin2a;
+              tep  := Hi*Cos2a - Hr*Sin2a;
+              Hr   := Cx*Cx*Hr - Sx*Sx*tem - Ca*Tee;
+              Hi   := Cx*Cx*Hi + Sx*Sx*tep - Sa*Tee;
+              b    := isw*Te*Ca + eta*Sa;
+              e    := eta*Ca    - isw*Te*Sa;
+            END; (* IF (Root1 < MachEps) *)
+(*enter1:*) s := Hr - sig*Root2*e;
+            c := Hi - sig*Root2*b;
+            Root := LongComplexMath.abs(CMPLX(c,s)); (* M.R.*)
+            IF (Root < MachEps) THEN
+              Cb:=1.0; Ch:=1.0; Sb:=0.0; Sh:=0.0;
+              (* GOTO trans; *)
+            ELSE
+              Cb  := - c / Root;
+              Sb  := s / Root;
+              Tee := Cb*b - e*Sb;
+              nC  := Tee*Tee;
+              TanH:= Root / (g + 2.0*(nC + nD));
+              Ch  := 1.0 / sqrt(1.0 - TanH*TanH);
+              Sh  := Ch*TanH;
+            END; (* IF *)
+(*trans:*)  tem   := Sx*Sh*(Sa*Cb - Sb*Ca);
+            im := - Sx*Sh*(Ca*Cb + Sa*Sb);
+            c1 := CMPLX(Cx*Ch - tem,im);
+            c2 := CMPLX(Cx*Ch + tem,im);
+            tep   := Sx*Ch*Ca;
+            tem   := Cx*Sh*Sb;
+            s1re :=   tep - tem;
+            s2re := - tep - tem;
+            tep   := Sx*Ch*Sa;
+            tem   := Cx*Sh*Cb;
+            s1im := tep + tem;
+            s2im := tep - tem;
+            s1    := CMPLX(s1re,s1im);
+            s2    := CMPLX(s2re,s2im);
+            tem   := LongComplexMath.abs(s1); (* M.R.*)
+            tep   := LongComplexMath.abs(s2);
+            IF (tem > MachEps) OR (tep > MachEps) THEN
+              Mark:=FALSE;
+              FOR i:=0 TO dim-1 DO (* Linksseitige Transformation *)
+                Aki:=A[k,i]; Ami:=A[m,i];
+                A[k,i] := c1*Aki + s1*Ami;
+                A[m,i] := s2*Aki + c2*Ami;
+              END;
+              FOR i:=0 TO dim-1 DO (* Linksvektoren *)
+                EVik:=LEV[k,i]; EVim:=LEV[m,i];
+                LEV[k,i] := c1*EVik + s1*EVim;
+                LEV[m,i] := s2*EVik + c2*EVim;
+              END;
+              FOR i:=0 TO dim-1 DO (* Rechtsseitige Transformation *)
+                Aik:=A[i,k]; Aim:=A[i,m];
+                A[i,k] :=   c2*Aik - s2*Aim;
+                A[i,m] := - s1*Aik + c1*Aim;
+              END;
+              FOR i:=0 TO dim-1 DO (* Rechtsvektoren *)
+                EVik:=REV[k,i]; EVim:=REV[m,i];
+                REV[k,i] :=   c2*EVik - s2*EVim;
+                REV[m,i] := - s1*EVik + c1*EVim;
+              END;
+            END; (* IF Transformation *)
+          END; (* FOR m *)
+        END; (* FOR k *)
+        (* tau2:=tau1; tau1:=tau; *)
+      END; (* LOOP *)
+      FOR i:=0 TO dim-1 DO EW[i]:=A[i,i]; END;
+      IF (Iter > MaxIter) THEN
+        Fehler:=TRUE; Fehlerflag:='MaxIter ueberschritten (ComEig).';
+        ErrOut(Fehlerflag);
+      END;
+      DEALLOCATE(En,dim*TSIZE(LONGREAL));
+END ComEig;
+
+END EigenLib3.