GENERIC MODULE EigenN(VRep, MRep);

(* 
  Generic implementation of "EigenN(VRep, MRep)"
  assuming that 
    "VRep.T" is "ARRAY OF VRep.ElemT",
    "MRep.T" is "ARRAY OF ARRAY OF MRep.ElemT", AND 
    "VRep.ElemT" and "MRep.ElemT" are either "REAL" or "LONGREAL"; 
  so that "m[i,j]", "a[j]" and "FLOAT(x, MRep.ElemT)" are legal.
  
  Created by J. Stolfi with contributions by Marcellus R. Macêdo. *)

FROM Math IMPORT sqrt;

PROCEDURE SymEigen(
    VAR m: MRep.T;
    VAR d: ARRAY OF LONGREAL;
    VAR v: ARRAY OF VRep.T;
    VAR b: ARRAY OF LONGREAL;
    VAR z: ARRAY OF LONGREAL;
  ) =
  CONST MaxIter = 50;
  VAR sm: LONGREAL;
      tresh, theta, tau, tan, sin, h, g, cos: LONGREAL;
  BEGIN
    WITH
      N = NUMBER(m),
      nSqr = FLOAT(N*N, LONGREAL)
    DO
      <* ASSERT NUMBER(m[0]) = N *>
      (* Initialize "v" with the canonical basis: *)
      FOR i := 0 TO N-1 DO
        FOR j := 0 TO N-1 DO v[i][j] := FLOAT(0.0d0, VRep.ElemT) END;
        v[i][i] := FLOAT(1.0d0, VRep.ElemT)
      END;
 
      FOR i := 0 TO N-1 DO
        b[i] := FLOAT(m[i][i], LONGREAL);
        d[i] := b[i];
        z[i] := 0.0d0;
      END;

      FOR iter := 0 TO MaxIter-1 DO
        (* Add elements off the diagonal of "m": *)
        sm := 0.0d0;
        FOR i := 0 TO N-2 DO
          FOR j := i+1 TO N-1 DO 
            WITH fmij = FLOAT(m[i][j], LONGREAL) DO sm := sm + ABS(fmij) END
          END
        END;
        IF sm = 0.0d0 THEN 
          RETURN 
        ELSIF iter < 4 THEN
          tresh := 0.2d0 * sm/nSqr
        ELSE
          tresh := 0.0d0
        END;

        FOR i := 0 TO N-2 DO
          FOR k := i+1 TO N-1 DO
            WITH fmik = FLOAT(m[i][k], LONGREAL) DO
              g := 100.0d0 * ABS(fmik);
              IF iter > 4 
              AND ABS(d[i])+g = ABS(d[i]) 
              AND ABS(d[k])+g = ABS(d[k])
              THEN
                m[i][k] := FLOAT(0.0d0, MRep.ElemT);
              ELSIF ABS(fmik) <= tresh THEN
                (* Ok *)
              ELSE
                h := d[k] - d[i];
                IF ABS(h) + g = ABS(h) THEN
                  tan := fmik/h;
                ELSE
                  theta := 0.5d0 * h/fmik;
                  tan := 1.0d0 /(ABS(theta) + sqrt(1.0d0 + theta*theta));
                  IF theta < 0.0d0  THEN tan := -tan END;
                END;

                cos := 1.0d0/sqrt(1.0d0 + tan*tan);
                sin := tan*cos;
                tau := sin/(1.0d0 + cos);
                WITH
                  h = tan*fmik
                DO
                  z[i] := z[i] - h;
                  z[k] := z[k] + h;
                  d[i] := d[i] - h;
                  d[k] := d[k] + h;
                END;
                m[i][k] := FLOAT(0.0d0, MRep.ElemT);

                FOR j := 0 TO i-1 DO
                  WITH
                    g = FLOAT(m[j][i], LONGREAL), 
                    h = FLOAT(m[j][k], LONGREAL)
                  DO
                    m[j][i] := FLOAT(g - sin*(h + g*tau), MRep.ElemT);
                    m[j][k] := FLOAT(h + sin*(g - h*tau), MRep.ElemT);
                  END
                END;

                FOR j := i+1 TO k-1 DO
                  WITH
                    g = FLOAT(m[i][j], LONGREAL), 
                    h = FLOAT(m[j][k], LONGREAL)
                  DO
                    m[i][j] := FLOAT(g - sin*(h + g*tau), MRep.ElemT);
                    m[j][k] := FLOAT(h + sin*(g - h*tau), MRep.ElemT);
                  END
                END;

                FOR j := k+1 TO N-1 DO
                  WITH
                    g = FLOAT(m[i][j], LONGREAL), 
                    h = FLOAT(m[k][j], LONGREAL)
                  DO
                    m[i][j] := FLOAT(g - sin*(h + g*tau), MRep.ElemT);
                    m[k][j] := FLOAT(h + sin*(g - h*tau), MRep.ElemT);
                  END
                END;

                FOR j := 0 TO N-1 DO
                  WITH
                    g = FLOAT(v[i][j], LONGREAL), 
                    h = FLOAT(v[k][j], LONGREAL)
                  DO
                    v[i][j] := FLOAT(g - sin*(h + g*tau), MRep.ElemT);
                    v[k][j] := FLOAT(h + sin*(g - h*tau), MRep.ElemT);
                  END
                END;

              END; (* IF iter > 4 AND... *)
            END 
          END; (* FOR k *)
        END; (* FOR i *)
        FOR i := 0 TO N-1 DO
          b[i] := b[i] + z[i];
          d[i] := b[i];
          z[i] := 0.0d0;
        END;
      END; (* FOR iter *)
      (* Too many iterations: *)
      <* ASSERT FALSE *> 
    END (* WITH *)
  END SymEigen;

BEGIN
END EigenN.