INT case = 1; INT n = (case|2, 2, 3, 1);
# n is a compile time parameter to choose a test case. #
MODE R = REAL;

# Algol 68 Matrix Inversion #
PROC inv = ([,]R m) [,]R: (INT n = 1 UPB m;
   IF n /= 2 UPB m THEN print("not square") FI;
   [n, 2*n]R am; FOR i TO n DO FOR j TO n DO
      am[i, j] := m[i, j]; am[i, j+n] := (i=j|1|0) OD OD;
   FOR i TO n DO INT im := 0; R mv := ABS am[i, i];
     FOR j FROM i+1 TO n DO R ae = ABS am[j, i];
       IF mv < ae THEN mv := ae; im := j FI OD;
     IF mv = 0.0 THEN print("Singular") FI;
     IF im > i THEN FOR n FROM i TO 2*n DO
        R s = am[i, n]; am[i, n] := am[im, n]; am[im, n] := s OD FI;
   (R r = 1.0 / am[i, i]; FOR j FROM i TO 2*n DO am[i, j] *:= r OD);
   FOR k FROM i+1 TO n DO R f = am[k, i]; 
     FOR j FROM i+1 TO 2*n DO am[k, j] -:= f * am[i, j] OD OD OD;
   FOR i TO n DO FOR j FROM i+1 TO n DO R p = am[i, j];
     FOR k FROM i+1 TO 2*n DO am[i, k] -:= am[j, k] * p OD OD OD;
   am[,n+1:2*n]);

# Algol 68 eigenvalues and vectors #
PROC eigen = (REF[,] R m, REF[]R l, REF[,] R vc) VOID: (
  [n,n] R m2 := m;
  FOR i TO n DO FOR j TO n DO vc[i, j] := (i=j|1|0) OD OD;
  WHILE R mod := 0; INT i, j; FOR k TO n DO FOR m FROM k+1 TO n DO R q = ABS m2[k,m];
    IF q > mod THEN mod:=q; i:=k; j:=m FI OD OD; mod > 1.e-14 DO 
    R th = 0.5*(R d = m2[i,i] - m2[j,j]; (d=0|pi/2|atan(2.0*m2[i,j]/d)));
    R c = cos(th), s = sin(th);
    PROC twst = (REF[,]R m)VOID: FOR k TO n DO R t = c*m[i,k]+s*m[j,k];
      m[j,k] := -s*m[i,k]+c*m[j,k]; m[i,k]:=t OD;
    FOR k TO n DO R t = c*m2[k,i] + s*m2[k,j];
      m2[k,j] := -s*m2[k,i] + c*m2[k,j]; m2[k,i] := t OD;
    twst(m2); twst(vc) OD;
  FOR j TO n DO l[j] := m2[j,j] OD);

OP = = ([]R a, []R b)BOOL:
  (FOR j TO n DO (ABS(a[j] - b[j]) > 1.e-13 | ex) OD; 0=0 EXIT ex: 0=1);

OP * = (R s, []R v)[]R: ([n]R a; FOR j TO n DO a[j] := s*v[j] OD; a);

OP * = ([,]R m, []R v)[]R: ([n]R a;
  FOR j TO n DO R s := 0; FOR k TO n DO s +:= m[j,k]*v[k] OD; a[j] := s OD; a);

PROC assert = (INT j, BOOL b) VOID: (b |print((j, newline))
  | print(("Exiting: ", j)); quit);

PROC ec = ([,]R m, REF[]R l, [,] R vc, INT t) VOID:
 FOR i TO n DO assert(t, m*vc[i,] = l[i]*vc[i,]) OD;

OP * = ([,]R a, [,]R b) [,]R: ([n,n]R ans;
  FOR i TO n DO FOR j TO n DO R s := 0;
    FOR k TO n DO s +:= a[i,k]*b[k,j] OD; ans[i, j] := s OD OD;
  ans);

PROC md = ([]R d)[,]R: ([n,n]R m;
  FOR i TO n DO FOR j TO n DO m[i,j] := (i=j|d[j]|0) OD OD; m);

OP = = ([,]R x, [,]R y) BOOL:
  (FOR i TO n DO (x[i,] = y[i,] |~| ex) OD; 0=0 EXIT ex: 0=1);

OP + = ([,]R x)[,]R: ([n,n]R a; FOR j TO n DO a[j,] := x[,j] OD; a);

PROC sym = ([,]R c) VOID: assert(20, c = +c);

PROC pr = (STRING s, [,]R c) VOID: (print((s, newline));
  FOR j TO n DO print((c[j,], newline)) OD);

PROC unif = ([]INT a)VOID: FOR j TO UPB a DO assert(-1, a[j] = a[1]) OD;

MODE PY = STRUCT(REF[] R lambda, REF[,] R vec);
PROC le = ([,]R c, []R d) PY: (assert(0, n = 1 UPB c);
assert(0, n = 2 UPB c ANDF n = UPB d); sym(c);
[n,n]R e; FOR i TO n DO FOR j TO n DO e[i,j] := (i=j|sqrt(d[j])|0) OD OD;
REF[,]R ece = (HEAP [n,n]R := e*c*e); sym(ece);
HEAP [n]R lambda; [n,n]R x;
eigen(ece, lambda, x); ec(ece, lambda, x, 5);
FOR i TO n DO assert(4, ece*x[i,] = lambda[i]*x[i,]) OD;
# Now verify DCEx = λEx for each of the n vectors x. #
FOR j TO n DO assert(5, md(d)*c*e*x[j,] = lambda[j]*(e*x[j,])) OD;
( # eigenvectors ex of lopsided matrix: #
HEAP [n,n]R ex; FOR j TO n DO ex[j,] := e*x[j,] OD;
[,]R dc = md(d)*c; 
FOR j TO n DO assert(7, dc*ex[j,] = lambda[j]*ex[j,]) OD;
print((lambda, " lambda", newline));
pr("dc", dc); pr("e*x", e*x); pr("ex", ex);
FOR j TO n DO print(("nx", newline, dc*ex[j,], newline,
  lambda[j]*ex[j,], newline)) OD;
ec(dc, lambda, ex, 9);
(lambda, ex)));

MODE PX = STRUCT([n,n]R c, [n]R d, [n]R init);
(PX px = (case|
    (R p = 2, q=1; (((-p, p), (p, -p-q)), (1/10, 1/3), (8, 14))),
    (((1, 2), (2, 1)), (1, 1), (0, 0)),
    (((-3, 4, 6), (4, -1, -5), (6, -5, -2)), (5, 3.3, 2), (1, 1, 2)),
    (((1)), (1/2), (2)));
  [,]R c = c OF px; []R d = d OF px;
  PY py = le(c, d); []R lam = lambda OF py;
  unif((1 UPB c, 2 UPB c, UPB d, UPB init OF px));
  print((lam, "evals", newline));
  pr("evecs:", vec OF py);
  [,]R dc = md(d)*c, g = +vec OF py, gi = inv(g);
  pr("G", g); pr("DC(=M)", dc); pr("G^-1MG", gi*dc*g);
  ([]R a = gi*init OF px;
  FOR j FROM 0 TO 20 DO R time = j*0.1;
    [n]R xp; FOR k TO n DO xp[k] := a[k]*exp(lam[k]*time) OD;
    print((time, "  ", g*xp, newline)) OD))
EXIT quit: ~