/*+T */
S=4D00;
/**********************************************************************/
/*                                                                    */
/*  Eigenwert- und Eigenvektorberechnung einer Matrix                 */
/*  =================================================                 */
/*                                                                    */
/*  Mit dem Unterprogramm EIGEN koennen die Eigenwerte und die Eigen- */
/*  vektoren einer (n,n)-Matrix bestimmt werden. Diese werden in      */
/*  vielfaeltiger Weise in der Berechnung von regelungstechnischen    */
/*  Problemen benoetigt.                                              */
/*  Es gilt:     A * s = V * s                                        */
/*       mit  A: (n,n)-Matrix                                         */
/*            V: (n)-Eigenvektor                                      */
/*            s: Eigenwert (Skalar)                                   */
/*                                                                    */
/*  Ueber die Eigenwertberechnung koennen z.B. auch die Nullstellen   */
/*  eines Polynoms bestimmt werden. Der Zusammenhang zwischen den     */
/*  Matrixkomponenten und den Polynomkoeffizienten besteht darin,     */
/*  dass wenn eine Matrix in der Regelungsnormalform vorliegt, die    */
/*  unterste Zeile mit den Polynomkoeffizienten des charakteris-      */
/*  tischen Polynoms belegt ist. Die Eigenvektoren der Matrix sind    */
/*  dann die Nullstellen des charakteristischen Polynoms.             */
/*  Die Nullstellen des Polynoms koennen mit dem Unterprogramm NULLST */
/*  berechnet werden.                                                 */
/*                                                                    */
/*  Die Programme wurden teilweise aus Fortran uebertragen. Leider    */
/*  war dies nicht bei allen benoetigten Unterprogrammen der Einfach- */
/*  heit halber moeglich, deshalb haben die Unterprogramme HQR2 und   */
/*  BALANC keine Blockstruktur, sondern arbeiten mit Labels und Gotos.*/
/*                                                                    */
/*  Die Programme wurden folgendem Buch entnommen:                    */
/*     Jordan-Engeln,G.,Reutter,F., Formelsammlung zur numerischen    */
/*     Mathematik mit Standard-Fortran-Programmen,                    */
/*     BI Hochschultaschenbuecher Band 106   (gibts nicht mehr!)      */
/*  HQR2 wurde in dieser Form der Regelungstechnischen Programmsamm-  */
/*  lung RASP des DFVLR entnommen, da mit ihr auch die Schur-Form     */
/*  einer oberen Hessenberg-Matrix berechenbar ist.                   */
/*  Die Programme stammen eigentlich aus :                            */
/*   - Peters, G. and Wilkinson, J.H., Eigenvectors of real and com-  */
/*     plex Matrices by LR and QR triangulations, Num. Math. 16,      */
/*     1970                                                           */
/*   - Peters, G. and Wilkinson, J.H., similarity reduction of a      */
/*     general matrix to Hessenberg form, Num. Math. 12, 1968         */
/*   - Parlett, B.N. and Reinsch, C., Balancing a matrix for calcu-   */
/*     lations of eigenvalues and eigenvectors, Num. Math. 13, 1969   */
/*                                                                    */
/*--------------------------------------------------------------------*/
/*                                                                    */
/*  Version 1.0   vom 02.11.86                                        */
/*                                                                    */
/*--------------------------------------------------------------------*/
/*                                                                    */
/*  Geschrieben von Matthias Schade, Korbacher Str. 96,               */
/*                  3501 Schauenburg                                  */
/*                                                                    */
/**********************************************************************/

MODULE EIGENW;

SYSTEM;
  DIALOG:A1<->;

PROBLEM;
  LENGTH FLOAT(55);
  SPC DIALOG DATION INOUT ALPHIC CONTROL(ALL);

/**********************************************************************/
/* FELDINIT: Feld loeschen.                                           */
/**********************************************************************/
FELDINIT:PROC (F(,) FLOAT IDENT);
FOR I FROM 1 TO (1 UPB F) REPEAT;
  FOR J FROM 1 TO (2 UPB F) REPEAT;
    F(I,J) := 0.0;
  END;
END;
END;

/**********************************************************************/
/* COMABS: Betragsbildung einer komplexen Zahl.                       */
/**********************************************************************/
COMABS:PROC ((X,Y) FLOAT) RETURNS (FLOAT);
DCL COMABS1 FLOAT;
IF (X/=0 OR Y/=0)
  THEN IF (ABS(X) < ABS(Y))
         THEN COMABS1 := ABS(Y) * SQRT(X/Y * X/Y + 1.0);
         ELSE COMABS1 := ABS(X) * SQRT(Y/X * Y/X + 1.0);
       FIN;
  ELSE COMABS1 := 0.0;
FIN;
RETURN (COMABS1);
END;

/**********************************************************************/
/* COMDIV: Division zweier komplexer Zahlen.                          */
/*         1.Zahl: A = Re, B = Im                                     */
/*         2.Zahl: C = Re, D = Im                                     */
/*         Ergebnis: X = Re, Y = Im                                   */
/**********************************************************************/
COMDIV:PROC ((A,B,C,D) FLOAT,(X,Y) FLOAT IDENT);
DCL (U,AM,AN,V,P,Q,F) FLOAT;
IF (C/=0 OR D/=0)
  THEN IF (A==0 AND B==0)
         THEN X := 0.0;
              Y := 0.0;
         ELSE IF (ABS(A) <= ABS(B))
                THEN U := B;
                     AM := A/B;
                     AN := 1.0;
                ELSE U := A;
                     AM := 1.0;
                     AN := B/A;
              FIN;
              IF (ABS(C) <= ABS(D))
                THEN V := D;
                     P := C/D;
                     Q := 1.0;
                ELSE V := C;
                     P := 1.0;
                     Q := D/C;
              FIN;
              F := U/V;
              V := P * P + Q * Q;
              U := (AM * P + AN * Q) / V;
              X := U * F;
              U := (-AM * Q + AN * P) / V;
              Y := U * F;
       FIN;
  ELSE PUT 'COMDIV: Division durch Null!' TO DIALOG BY A,SKIP;
FIN;
END;

/**********************************************************************/
/* NORMAL:                                                            */
/**********************************************************************/
NORMAL:PROC (N FIXED,V(,) FLOAT IDENT,(WR,WI)() FLOAT IDENT) GLOBAL;
DCL (X,Y,RETEIL,IMTEIL,GROSS) FLOAT, KOMPL BIT(1);
KOMPL := '0'B1;
FOR I FROM 1 TO N REPEAT;
  IF NOT KOMPL
    THEN IF WI(I) == 0.0
           THEN GROSS := V(1,I);
                FOR J FROM 2 TO N REPEAT;
                  IF (ABS(V(J,I)) > ABS(GROSS))
                    THEN GROSS := V(J,I);
                  FIN;
                END;
                FOR J FROM 1 TO N REPEAT;
                  V(J,I) := V(J,I) / GROSS;
                END;
           ELSE KOMPL := '1'B1;
                RETEIL := V(1,I);
                IMTEIL := V(1,I+1);
                FOR J FROM 2 TO N REPEAT;
                  IF COMABS(V(J,I),V(J,I+1)) > COMABS(RETEIL,IMTEIL)
                    THEN RETEIL := V(J,I);
                         IMTEIL := V(J,I+1);
                  FIN;
                END;
                FOR J FROM 1 TO N REPEAT;
                  CALL COMDIV(V(J,I),V(J,I+1),RETEIL,IMTEIL,X,Y);
                  V(J,I) := X;
                  V(J,I+1) := Y;
                END;
         FIN;
    ELSE KOMPL := '0'B1;
  FIN;
END;
END;

/**********************************************************************/
/* BALBAK: Bestimmung der Eigenvektoren einer Matrix, die durch das   */
/*         Unterprogramm BALANC transformiert wurde.                  */
/* Eingangsargumente:                                                 */
/*  Z     : Matrix, enthaelt die Eigenvektoren die zuruecktransfor-   */
/*          miert werden sollen.                                      */
/*  N     : Zeilen und Spaltenzahl der aktuellen Matrix Z.            */
/*  LOW,IHI: Unterer bzw. oberer Spaltenindex fuer die Hessenberg-    */
/*          transformation.                                           */
/*  D     : enthaelt Information ueber die Balancierung.              */
/*  M     : Anzahl der Vektoren, die zuruecktransformiert werden      */
/*          sollen.                                                   */
/* Ausgangsargumente:                                                 */
/*  Z     : Matrix mit den transformierten Eigenvektoren.             */
/**********************************************************************/
BALBAK:PROC (N FIXED,Z(,) FLOAT IDENT,D() FLOAT IDENT,
             (LOW,IHI,M) FIXED) GLOBAL;
DCL S FLOAT, (K,I) FIXED;
IF IHI /= LOW
  THEN FOR G FROM LOW TO IHI REPEAT;
         S := D(G);
         FOR J FROM 1 TO M REPEAT;
           Z(G,J) := Z(G,J) * S;
         END;
       END;
FIN;
FOR L FROM 1 TO N REPEAT;
  IF (L<=IHI AND L>= LOW)
    THEN ;
    ELSE I := LOW - L;
         IF L>IHI
           THEN I := L;
         FIN;
         K := ENTIER(D(I)+0.0001);
         IF K/=I
           THEN FOR J FROM 1 TO M REPEAT;
                  S := Z(I,J);
                  Z(I,J) := Z(K,J);
                  Z(K,J) := S;
                END;
         FIN;
    FIN;
END;
END;

/**********************************************************************/
/* HQR2: Berechnung der Eigenwerte, Eigenvektoren oder oberen Schur-  */
/*       Form einer oberen Hessenberg-Matrix mit Hilfe des QR-Ver-    */
/*       fahrens.                                                     */
/* Eingangsargumente:                                                 */
/*  H      : Matrix in der oberen Hessenbergform (wird ueberschrieben)*/
/*  N      : Aktuelle Ordnung der Matrix H.                           */
/*  LOW,IHI: Oberer bzw. unterer Spaltenindex fuer die Hessenberg-    */
/*           Transformation. Wenn H nicht balanciert ist: LOW=1 IHI=N */
/*  IOP = 0: Nur Eigenwertberechnung.                                 */
/*      = 1: Berechnung der Schur-Form von H.                         */
/*      = 2: Eigenwert und Eigenvektorberechnung.                     */
/*  VECS   : IOP=0, Dummyarray, wird nicht veraendert.                */
/*           IOP=1,2 Transformationsmatrix der Hessenberg-Transforma- */
/*           tion.                                                    */
/* Ausgangsargumente:                                                 */
/*  H      : Bei IOP=1 ist H in Schur-Form.                           */
/*  WR,WI  : Vektoren der Laenge N, enthalten die Real- und Imaginaer-*/
/*           teile der Eigenwerte.                                    */
/*  Z      : Bei IOP=1 ist Z die Transformationsmatrix                */
/*           H(ausgang) = Z^-1 * A * Z                                */
/*           A = urspruengliche Matrix vor Hessenbergtransformation   */
/*           und Balancieren.                                         */
/*           War Z(eingang) orthonormal, so ist es auch Z(ausgang).   */
/*           Bei IOP=2 enthaelt Z die Eigenvektoren.                  */
/*           Ist der I-te Eigenvektor komplex, so ist in der Matrix Z */
/*           die I-te Spalte der Realteil und die I+1-te Spalte der   */
/*           Imaginaerteil des zugehoerenden Eigenvektors.            */
/*  ITERI  : Anzahl der Iterationen bei den Eigenvektoren.            */
/*  FAIL   : Falls Wahr, ist keine Konvergenz aufgetreten.            */
/**********************************************************************/
HQR2:PROC (N FIXED,H(,) FLOAT IDENT,ITERI() FIXED IDENT,
     (LOW,IHI) FIXED,VECS(,) FLOAT IDENT,(WR,WI)() FLOAT IDENT,
     IOP FIXED,FAIL BIT IDENT) GLOBAL;
DCL EPSI FLOAT INIT (1.0E-30); /* relative Maschinengenauigkeit */
DCL (I,IP1,ITS,J,K1,L,M,MEND,MP2,MP3,NA,NE,NAM1) FIXED;
DCL (F,NORM,P,Q,R,S,T,VI,VR,W,X,Y,Z,RA,SA) FLOAT;
DCL NOLAST BIT(1);
FAIL := '0'B1;
/* Store roots isolated by BALANC and compute matrix norm */
     NORM := 0.0;
     K1 := 1;
     FOR I1 FROM 1 TO N REPEAT;
       FOR J1 FROM K1 TO N REPEAT;
         NORM := NORM + ABS(H(I1,J1));
       END;
       K1 := I1;
       ITERI(I1) :=0;
       IF (I1>=LOW AND I1<=IHI)
         THEN GOTO L100;
       FIN;
       WR(I1) := H(I1,I1);
       WI(I1) := 0.0;
L100:END;
     NE := IHI;
     T := 0.0;
/* Search for next eigenvalues */
L150:IF (NE < LOW)
       THEN GOTO L700;
     FIN;
     ITS := 0;
     NA := NE - 1;
/* Look for single small sub-diagonal element */
L200:FOR LL FROM LOW TO NE REPEAT;
       L := NE + LOW - LL;
       IF L==LOW
         THEN GOTO L250;
       FIN;
       S := ABS(H(L-1,L-1)) + ABS(H(L,L));
       IF S==0.0 THEN S := NORM; FIN;
       IF (ABS(H(L,L-1))<=EPSI*S)
         THEN GOTO L250;
       FIN;
     END;
/* Form shift */
L250:X := H(NE,NE);
     IF L==NE
       THEN GOTO L550;
     FIN;
     Y := H(NA,NA);
     W := H(NE,NA) * H(NA,NE);
     IF L==NA THEN GOTO L600; FIN;
     IF ITS <=30 THEN GOTO L270; FIN;
     ITERI(NE) := 31;
     FAIL := '1'B1;
     GOTO L990;
L270:IF (ITS/=10 AND ITS/=20) THEN GOTO L300; FIN;
/* Form exceptional shift */
     T := T + X;
     FOR I1 FROM LOW TO NE REPEAT;
       H(I1,I1) := H(I1,I1) - X;
     END;
     S := ABS(H(NE,NA)) + ABS(H(NA,NE-2));
     X := 0.75 * S;
     Y := X;
     W := -0.4375 * S * S;
L300:ITS := ITS + 1;
/* Look for two consecutive small sub-diagonal elements */
     MEND := NE - 2;
     FOR MM FROM L TO MEND REPEAT;
       M := MEND + L - MM;
       Z := H(M,M);
       R := X - Z;
       S := Y - Z;
       P := (R * S - W) / H(M+1,M) + H(M,M+1);
       Q := H(M+1,M+1) - Z - R - S;
       R := H(M+2,M+1);
       S := ABS(P) + ABS(Q) + ABS(R);
       P := P / S;
       Q := Q / S;
       R := R / S;
       IF M==L THEN GOTO L400; FIN;
       IF (ABS(H(M,M-1))*(ABS(Q)+ABS(R))<=EPSI*ABS(P)*(ABS(H(M-1,M-1))
          +ABS(Z)+ABS(H(M+1,M+1)))) THEN GOTO L400; FIN;
     END;
L400:MP2 := M + 2;
     FOR I1 FROM MP2 TO NE REPEAT;
       H(I1,I1-2) := 0.0;
       IF I1==MP2 THEN GOTO L410; FIN;
       H(I1,I1-3) := 0.0;
L410:END;
/* Double QR step involving rows L to NE and columns M to NE */
     FOR K FROM M TO NA REPEAT;
       NOLAST := K /= NA;
       IF K==M THEN GOTO L430; FIN;
       P := H(K,K-1);
       Q := H(K+1,K-1);
       R := 0.0;
       IF NOLAST THEN R := H(K+2,K-1); FIN;
       X := ABS(P) + ABS(Q) + ABS(R);
       IF X==0.0 THEN GOTO L500; FIN;
       P := P / X;
       Q := Q / X;
       R := R / X;
L430:  S := SQRT(P * P + Q * Q + R * R);
       IF P<0.0 THEN S := -S; FIN;
       IF K==M THEN GOTO L435; FIN;
       H(K,K-1) := -S * X;
       GOTO L440;
L435:  IF L/=M THEN H(K,K-1) := -H(K,K-1); FIN;
L440:  P := P + S;
       X := P / S;
       Y := Q / S;
       Z := R / S;
       Q := Q / P;
       R := R / P;
/* Row modification */
       FOR J1 FROM K TO N REPEAT;
         P := H(K,J1) + Q * H(K+1,J1);
         IF NOT NOLAST THEN GOTO L455; FIN;
         P := P + R * H(K+2,J1);
         H(K+2,J1) := H(K+2,J1) - P * Z;
L455:    H(K+1,J1) := H(K+1,J1) - P * Y;
         H(K  ,J1) := H(K  ,J1) - P * X;
       END;
       J := NE;
       IF K+3<NE THEN J := K + 3; FIN;
/* column modification */
       FOR I1 FROM 1 TO J REPEAT;
         P := X * H(I1,K) + Y * H(I1,K+1);
         IF NOT NOLAST THEN GOTO L465; FIN;
         P := P + Z * H(I1,K+2);
         H(I1,K+2) := H(I1,K+2) - P * R;
L465:    H(I1,K+1) := H(I1,K+1) - P * Q;
         H(I1,K  ) := H(I1,K  ) - P;
       END;
       IF IOP<1 THEN GOTO L500; FIN;
       FOR I1 FROM LOW TO IHI REPEAT;
         P := X * VECS(I1,K) + Y * VECS(I1,K+1);
         IF NOT NOLAST THEN GOTO L475; FIN;
         P := P + Z * VECS(I1,K+2);
         VECS(I1,K+2) := VECS(I1,K+2) - P * R;
L475:    VECS(I1,K+1) := VECS(I1,K+1) - P * Q;
L480:    VECS(I1,K  ) := VECS(I1,K  ) - P;
       END;
L500:END;
     GOTO L200;
/* One root found */
L550:H(NE,NE) := X + T;
     WR(NE) := H(NE,NE);
     WI(NE) := 0.0;
     ITERI(NE) := ITS;
     NE := NA;
     GOTO L150;
/* Two roots found */
L600:P := (Y - X) / 2;
     Q := P * P + W;
     Z := SQRT(ABS(Q));
     H(NE,NE) := X + T;
     X := H(NE,NE);
     H(NA,NA) := Y + T;
     ITERI(NA) := ITS;
     ITERI(NE) := - ITS;
     IF Q<=0.0 THEN GOTO L635; FIN;
/* Real pair */
     Z := ABS(Z) + ABS(P);
     IF P<0.0 THEN Z := -Z; FIN;
     WR(NA) := X + Z;
     WR(NE) := WR(NA);
     IF Z/=0.0 THEN WR(NE) := X - W / Z; FIN;
     WI(NA) := 0.0;
     WI(NE) := 0.0;
     X := H(NE,NA);
     S := ABS(X) + ABS(Z);
     P := X / S;
     Q := Z / S;
     R := SQRT(P * P + Q * Q);
     P := P / R;
     Q := Q / R;
     FOR J1 FROM NA TO N REPEAT;
       Z := H(NA,J1);
       H(NA,J1) := Q * Z + P * H(NE,J1);
       H(NE,J1) := Q * H(NE,J1) - P * Z;
     END;
     FOR I1 FROM 1 TO NE REPEAT;
       Z := H(I1,NA);
       H(I1,NA) := Q * Z + P * H(I1,NE);
       H(I1,NE) := Q * H(I1,NE) - P * Z;
     END;
     IF IOP<1 THEN GOTO L630; FIN;
     FOR I1 FROM LOW TO IHI REPEAT;
       Z := VECS(I1,NA);
       VECS(I1,NA) := Q * Z + P * VECS(I1,NE);
       VECS(I1,NE) := Q * VECS(I1,NE) - P * Z;
     END;
L630:GOTO L640;
/* complex pair */
L635:WR(NA) := X + P;
     WR(NE) := X + P;
     WI(NA) := Z;
     WI(NE) := -Z;
L640:NE := NE - 2;
     GOTO L150;
/* ?????? */
L700:IF IOP<2 THEN GOTO L990; FIN;
     IF NORM==0.0 THEN GOTO L990; FIN;
     FOR NX FROM 1 TO N REPEAT;
       NE := N + 1 - NX;
       P := WR(NE);
       Q := WI(NE);
       NA := NE - 1;
       IF Q/=0.0 THEN GOTO L810; FIN;
       M := NE;
       H(NE,NE) := 1.0;
       IF NA==0.0 THEN GOTO L900; FIN;
       FOR IX FROM 1 TO NA REPEAT;
         I := NA + 1 - IX;
         W := H(I,I) - P;
         R := H(I,NE);
         IF M>NA THEN GOTO L760; FIN;
         FOR J1 FROM M TO NA REPEAT;
           R := R + H(I,J1) * H(J1,NE);
         END;
L760:    IF WI(I)>=0.0 THEN GOTO L770; FIN;
         Z := W;
         S := R;
         GOTO L800;
L770:    M := I;
         IF WI(I)/=0.0 THEN GOTO L780; FIN;
         F := W;
         IF W==0.0 THEN F := EPSI * NORM; FIN;
         H(I,NE) := -R / F;
         GOTO L800;
L780:    X := H(I,I+1);
         Y := H(I+1,I);
         Q := (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I);
         T := (X * S - Z * R) / Q;
         H(I,NE) := T;
         IF ABS(X)<=ABS(Z) THEN GOTO L790; FIN;
         H(I+1,NE) := (-R - W * T) / X;
         GOTO L800;
L790:    H(I+1,NE) := (-S - Y * T) / Z;
L800:  END;
       GOTO L900;
L810:  IF Q>0.0 THEN GOTO L900; FIN;
       M := NA;
       IF ABS(H(NE,NA))<=ABS(H(NA,NE)) THEN GOTO L820; FIN;
       H(NA,NA) := Q / H(NE,NA);
       H(NA,NE) := -(H(NE,NE) - P) / H(NE,NA);
       GOTO L830;
L820:  CALL COMDIV(0.0,-H(NA,NE),H(NA,NA)-P,Q,H(NA,NA),H(NA,NE));
L830:  H(NE,NA) := 0.0;
       H(NE,NE) := 1.0;
       NAM1 := NA - 1;
       IF NAM1==0 THEN GOTO L900; FIN;
       FOR IX FROM 1 TO NAM1 REPEAT;
         I := NA - IX;
         W := H(I,I) - P;
         RA := 0.0;
         SA := H(I,NE);
         FOR J1 FROM M TO NA REPEAT;
           RA := RA + H(I,J1) * H(J1,NA);
           SA := SA + H(I,J1) * H(J1,NE);
         END;
         IF WI(I)>=0.0 THEN GOTO L870; FIN;
         Z := W;
         R := RA;
         S := SA;
         GOTO L890;
L870:    M := I;
         IF WI(I)/=0.0 THEN GOTO L880; FIN;
         CALL COMDIV(-RA,-SA,W,Q,H(I,NA),H(I,NE));
         GOTO L890;
L880:    X := H(I,I+1);
         Y := H(I+1,I);
         VR := (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I) - Q * Q;
         VI := (WR(I) - P) * 2.0 * Q;
         IF (VR==0.0 AND VI==0.0)
           THEN VR := EPSI * NORM * (ABS(W) + ABS(Q) + ABS(X) +
                      ABS(Y) + ABS(Z));
         FIN;
         CALL COMDIV(X*R - Z*RA + Q*SA,X*S - Z*SA - Q*RA,VR,VI,
                  H(I,NA),H(I,NE));
         IF ABS(X)<=ABS(Z)+ABS(Q) THEN GOTO L885; FIN;
         H(I+1,NA) := (-RA - W * H(I,NA) + Q * H(I,NE)) / X;
         H(I+1,NE) := (-SA - W * H(I,NE) - Q * H(I,NA)) / X;
         GOTO L890;
L885:    CALL COMDIV(-R-Y*H(I,NA),-S-Y*H(I,NE),Z,Q,H(I+1,NA),H(I+1,NE));
L890:  END;
L900:END;
     FOR I1 FROM 1 TO N REPEAT;
       IF (I1>=LOW AND I1<=IHI) THEN GOTO L940; FIN;
       FOR J1 FROM I1 TO N REPEAT;
         VECS(I1,J1) := H(I1,J1);
       END;
L940:END;
     FOR JX FROM LOW TO N REPEAT;
       J := N + LOW - JX;
       M := IHI;
       IF J<IHI THEN M := J; FIN;
       FOR I1 FROM LOW TO IHI REPEAT;
         Z := 0.0;
         FOR K FROM LOW TO M REPEAT;
           Z := Z + VECS(I1,K) * H(K,J);
         END;
         VECS(I1,J) := Z;
       END;
       GOTO L980;
L980:END;
L990:RETURN;
END;

/**********************************************************************/
/* ELMTRA: Berechnung der Transformationsmatrix T zu der in UP ELMHES */
/*         ausgefuehrten Hessenberg-Transformation: Ahess = Z^-1*A*Z  */
/* Eingangsargumente:                                                 */
/*  H     : Matrix (obere Hessenbergform) im oberen Dreieck und der   */
/*          unteren Nebendiagonalen. Auf den uebrigen Plaetzen und in */
/*          INF steht Information fuer die Ruecktransformation        */
/*          (s. UP ELMHES).                                           */
/*  N     : Aktuelle Ordnung der Matrix H.                            */
/*  LOW,IHI: Unterer bzw. oberer Spaltenindex fuer die Hessenberg-    */
/*          Transformation. Wenn A nicht balanciert: LOW=1 IHI=N      */
/*  INF   : Vektor Laenge IHI, enthaelt Information fuer die Rueck-   */
/*          transformation.                                           */
/* Ausgangsargumente:                                                 */
/*  V     : Transformationmatrix (n,n).                               */
/**********************************************************************/
ELMTRA:PROC (N FIXED,H(,) FLOAT IDENT,INF() FIXED IDENT,
             (LOW,IHI) FIXED,V(,) FLOAT IDENT) GLOBAL;
DCL (I,IEND,IP1,J) FIXED;
FOR I1 FROM 1 TO N REPEAT;
  FOR J1 FROM 1 TO N REPEAT;
    V(I1,J1) := 0.0;
  END;
  V(I1,I1) := 1.0;
END;
IEND := IHI - LOW - 1;
IF IEND >= 1
  THEN FOR L FROM 1 TO IEND REPEAT;
         I := IHI - L;
         J := INF(I);
         IP1 := I + 1;
         FOR K FROM IP1 TO IHI REPEAT;
           V(K,I) := H(K,I-1);
         END;
         IF I/=J
           THEN FOR K FROM I TO IHI REPEAT;
                  V(I,K) := V(J,K);
                  V(J,K) := 0.0;
                END;
                V(J,I) := 1.0;
         FIN;
       END;
FIN;
END;

/**********************************************************************/
/* ELMHES: Berechnung der oberen Hessenberg-Form einer Matrix.        */
/* Eingangsargumente:                                                 */
/*  A     : Matrix die transformiert werden soll. Wird ueberschrieben.*/
/*  N     : Aktuelle Ordnung der Matrix A.                            */
/*  LOW,IHI: Unterer bzw. oberer Spaltenindex fuer den Bereich, in dem*/
/*          die Transformation ausgefuehrt werden soll. Wenn A nicht  */
/*          balanciert ist: LOW=1 IHI=N.                              */
/* Ausgangsargumente:                                                 */
/*  A     : Transformierte Matrix im oberen Dreieck und der unteren   */
/*          Nebendiagonalen. Auf den sonstigen Plaetzen und in        */
/*  INF   : Vektor Laenge IHI, steht Information fuer die Ruecktrafo. */
/**********************************************************************/
ELMHES:PROC (N FIXED,A(,) FLOAT IDENT,INF() FIXED IDENT,
             (LOW,IHI) FIXED) GLOBAL;
DCL (X,Y) FLOAT,(I,IANF,JANF,MANF,MEND) FIXED;
MANF := LOW + 1;
MEND := IHI - 1;
IF MANF<=MEND
  THEN FOR M FROM MANF TO MEND REPEAT;
         I := M;
         X := 0.0;
         FOR J FROM M TO IHI REPEAT;
           IF (ABS(A(J,M-1))>ABS(X))
             THEN X := A(J,M-1);
                  I := J;
           FIN;
         END;
         INF(M) := I;
         IF I/=M
           THEN JANF := M - 1;
                FOR J FROM JANF TO N REPEAT;
                  Y := A(I,J);
                  A(I,J) := A(M,J);
                  A(M,J) := Y;
                END;
                FOR J FROM 1 TO IHI REPEAT;
                  Y := A(J,I);
                  A(J,I) := A(J,M);
                  A(J,M) := Y;
                END;
         FIN;
         IF X/=0.0
           THEN IANF := M + 1;
                FOR I1 FROM IANF TO IHI REPEAT;
                  Y := A(I1,M-1);
                  IF Y/=0.0
                    THEN Y := Y / X;
                         A(I1,M-1) := Y;
                         FOR J FROM M TO N REPEAT;
                           A(I1,J) := A(I1,J) - Y * A(M,J);
                         END;
                         FOR J FROM 1 TO IHI REPEAT;
                           A(J,M) := A(J,M) + Y * A(J,I1);
                         END;
                  FIN;
                END;
         FIN;
       END;
FIN;
END;

/**********************************************************************/
/* BALANC: Balancieren einer Matrix vor Eigenwert-, Eigenvektorbe-    */
/*         rechnungen. (Normreduktion von A durch Aehnlichkeits-      */
/*         transformation mit Diagonalmatrizen - Balancieren).        */
/* Eingangsargumente:                                                 */
/*  A     : Matrix.                                                   */
/*  N     : Aktuelle Ordnung der Matrix A.                            */
/* Ausgangsargumente:                                                 */
/*  A     : Balancierte Matrix.                                       */
/*  LOW,IHI: Unterer bzw. oberer Spaltenindex fuer die Hessenberg-    */
/*          Transformation.                                           */
/*  D     : Vektor Laenge N, enthaelt Information fuer die Rueck-     */
/*          transformation durch BALBAK.                              */
/**********************************************************************/
BALANC:PROC (N FIXED,A(,) FLOAT IDENT,D() FLOAT IDENT,
             (LOW,IHI) FIXED IDENT) GLOBAL;
DCL BASIS FLOAT INIT (10);
DCL (NOCONV,WEG) BIT(1),(B2,C,F,G,R,S) FLOAT,(J,M) FIXED;
     B2 := BASIS * BASIS;
     LOW := 1;
     IHI := N;
     WEG := '1'B1;
B100:FOR JX FROM 1 TO IHI REPEAT;
       J := IHI + 1 - JX;
       FOR I FROM 1 TO IHI REPEAT;
         IF I==J THEN GOTO B180; FIN;
         IF A(I,J)/=0.0 THEN GOTO B200; FIN;
B180:  END;
       M := IHI;
       GOTO B450;
B200:END;
     WEG := '0'B1;
B300:FOR J1 FROM LOW TO IHI REPEAT;
       FOR I FROM LOW TO IHI REPEAT;
         IF I==J1 THEN GOTO B380; FIN;
         IF A(I,J1)/=0.0 THEN GOTO B400; FIN;
B380:  END;
       M := LOW;
       GOTO B450;
B400:END;
     GOTO B500;
B450:D(M) := J;
     IF J==M THEN GOTO B480; FIN;
     FOR I FROM 1 TO IHI REPEAT;
       F := A(I,J);
       A(I,J) := A(I,M);
       A(I,M) := F;
     END;
     FOR I FROM LOW TO N REPEAT;
       F := A(J,I);
       A(J,I) := A(M,I);
       A(M,I) := F;
     END;
B480:IF WEG THEN GOTO B490; FIN;
     LOW := LOW + 1;
     GOTO B300;
B490:IF IHI==1 THEN GOTO B950; FIN;
     IHI := IHI - 1;
     GOTO B100;
B500:FOR I FROM LOW TO IHI REPEAT;
       D(I) :=1.0;
     END;
B600:NOCONV := '0'B1;
     FOR I FROM LOW TO IHI REPEAT;
       C := 0.0;
       R := 0.0;
       FOR J1 FROM LOW TO IHI REPEAT;
         IF J1==I THEN GOTO B650; FIN;
         C := C + ABS(A(J1,I));
         R := R + ABS(A(I,J1));
B650:  END;
       G := R / BASIS;
       F := 1.0;
       S := C + R;
B660:  IF C>=G THEN GOTO B700; FIN;
       F := F * BASIS;
       C := C * B2;
       GOTO B660;
B700:  G := R * BASIS;
B760:  IF C<G THEN GOTO B800; FIN;
       F := F / BASIS;
       C := C / B2;
       GOTO B760;
B800:  IF ((C+R)/F >= 0.95*S) THEN GOTO B900; FIN;
       G := 1.0 / F;
       D(I) := D(I) * F;
       NOCONV := '1'B1;
       FOR J1 FROM LOW TO N REPEAT;
         A(I,J1) := A(I,J1) * G;
       END;
       FOR J1 FROM 1 TO IHI REPEAT;
         A(J1,I) := A(J1,I) * F;
       END;
B900:END;
     IF NOCONV THEN GOTO B600; FIN;
B950:RETURN;
END;

/**********************************************************************/
/* EIGEN: Berechnung von Eigenwerten und Eigenvektoren von A.         */
/* Eingangsargumente:                                                 */
/*  N    : Aktuelle Ordnung der Matrix A.                             */
/*  A    : Matrix, von der die Eigenwerte,-vektoren berechnet werden. */
/* Ausgangsargumente:                                                 */
/*  Z    : Matrix mit den Eigenvektoren                               */
/*  WR,WI: Eigenwerte (Real-,Imaginaerteil).                          */
/*  IFEHL: ==0 Berechnung erfolgreich, /= nicht erfolgreich beendet.  */
/* Die restlichen Felder werden nur deshalb uebergeben, damit die     */
/* Feldgroesse mit der vom aufrufenden Programm uebereinstimmt.       */
/**********************************************************************/
EIGEN:PROC (N FIXED,A(,) FLOAT IDENT,IOP FIXED,D() FLOAT IDENT,
            Z(,) FLOAT IDENT,(WR,WI)() FLOAT IDENT,
            (INF,ITERI)() FIXED IDENT,IFEHL FIXED IDENT) GLOBAL;
DCL FAIL BIT(1),(LOW,IHI) FIXED;
IFEHL := 1;
CALL BALANC(N,A,D,LOW,IHI);
CALL ELMHES(N,A,INF,LOW,IHI);
CALL ELMTRA(N,A,INF,LOW,IHI,Z);
IF N/=2
  THEN FOR I FROM 3 TO N REPEAT;
         FOR J FROM 3 TO I REPEAT;
           A(I,J-2) := 0.0;
         END;
       END;
FIN;
CALL HQR2(N,A,ITERI,LOW,IHI,Z,WR,WI,IOP,FAIL);
IF (NOT FAIL AND IOP/=0)
  THEN CALL BALBAK(N,Z,D,LOW,IHI,N);
       CALL NORMAL(N,Z,WR,WI);
FIN;
IF NOT FAIL THEN IFEHL := 0; FIN;
END;

/**********************************************************************/
/* RNFORM: Ein Polynom (Polynomkoeffizienten) in eine                 */
/*                      Matrix in Regelungsnormalform umsetzen.       */
/* Eingangsargumente:                                                 */
/*  N    : Grad des Polynoms.                                         */
/*  V    : Polynomkoeffizienten (Vektor Laenge N+1)                   */
/* Ausgangsargumente:                                                 */
/*  F    : Matrix in die die Koeffizienten eingesetzt werden sollen.  */
/*  KREISVERSTAERKUNG: Falls Polynom mit Koeffizient/=1 multipliziert.*/
/**********************************************************************/
RNFORM:PROC(N FIXED,V() FLOAT IDENT,F(,) FLOAT IDENT,
                         KREISVERSTAERKUNG FLOAT IDENT) GLOBAL;
IF V(N+1) /= 1
  THEN KREISVERSTAERKUNG := V(N+1);
       FOR I FROM 1 TO N+1 REPEAT;
         V(I) := V(I) / V(N+1);
       END;
  ELSE KREISVERSTAERKUNG := 1.0;
FIN;
CALL FELDINIT(F);
FOR I FROM 1 TO N-1 REPEAT;
  F(I,I+1) := 1.0;
END;
FOR I FROM 1 TO N REPEAT;
  F(N,I) := -V(I);
END;
END;

/**********************************************************************/
/* NULLST: Berechnung der Nullstellen eines Polynoms.                 */
/* Eingangsargumente:                                                 */
/*  N     : Grad des Polynoms.                                        */
/*  V     : Polynomkoeffizienten (Vektor Lange N+1).                  */
/* Ausgangsargumente:                                                 */
/*  VN    : Nullstellen des Polynoms (Matrix (n+1,2), in (i,1) Real-  */
/*          teil, in (i,2) Imaginaerteil).                            */
/*  KREISVERSTAERKUNG: Falls Polynom mit Koeffizient/=1 multipliziert.*/
/* Die restlichen Felder werden nur uebergeben, damit die Feldgroesse */
/* mit dem aufrufenden Programm uebereinstimmt.                       */
/**********************************************************************/
NULLST:PROC (N FIXED,V() FLOAT IDENT,VN(,) FLOAT IDENT,
             KREISVERSTAERKUNG FLOAT IDENT,(A,Z)(,) FLOAT IDENT,
             (D,WR,WI)() FLOAT IDENT,(INF,ITERI)() FIXED IDENT) GLOBAL;
DCL (IFEHL,IOP) FIXED INIT (1,0);  /* IOP=0: Nur Eigenwertberechnung */
CALL RNFORM(N,V,A,KREISVERSTAERKUNG);
CALL EIGEN(N,A,IOP,D,Z,WR,WI,INF,ITERI,IFEHL);
CALL FELDINIT(VN);
IF IFEHL == 0
  THEN FOR I FROM 1 TO N REPEAT;
         VN(I,1) := WR(I);
         VN(I,2) := WI(I);
       END;
FIN;
END;


MODEND.