INT n = 2;
# 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);

PROC ec = (REF[,] 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 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);

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

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 assert = (INT j, BOOL b) VOID: (b |~| print(("Exiting: ", j)); quit);

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

MODE PX = STRUCT(REF[] R lambda, REF[,] R vecr);
PROC le = ([,]R c, []R d) PX: (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;
assert(1, e*e = md(d));
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);
# 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: #
[n,n]R ex; FOR j TO n DO ex[j,] := e*x[j,] OD;
[n,n]R dc; FOR j TO n DO dc[j,] := d[j] * c[j,] OD;
assert(6, dc = md(d)*c);
ec(dc, lambda, ex, 9);
(lambda, HEAP [n,n]R := ex)));

CASE 1 IN (PX px = (n-1|(R p=2, q=1; le(((-p, p), (p, -p-q)), (1/10, 1/3))),
   le(((3, 4, 6), (4, -1, 5), (6, 5, 2)), (5, 6, 2)));
[,]R x = vecr OF px;
print((lambda OF px, newline, x, newline));
IF n=2 THEN [,]R r = inv(x); []R y = (8, 14);
print((y, newline, r*y, newline, x*r*y, newline)) FI),

(PROC ldesolv = ([,]R c, []R cap, []R init) VOID: (
  [n]R d; FOR j TO n DO d[j] := 1/cap[j] OD;
  PX r = le(c, d); [,]R x = vecr OF r; []R lam = lambda OF r;
  print(("lam ", lam, newline));
  []R it = inv(x)*init;
  FOR m FROM 0 TO 20 DO R t = m/10;
    [n]R dv; FOR j TO n DO dv[j] := exp(lam[j]*t)*it[j] OD;
    []R yv = x*dv;
    print((t, "  ", yv[1], "  ", yv[2], newline)) OD);

(R p=2, q=1; 
ldesolv(((-p, p), (p, -p-q)), (10, 3), (8, 14)))) ESAC;

 quit: ~

# n=2
-5.83974396909360e-2  -1.14160256030906e+0

+2.94834610067746e-1  +2.08746678266587e-1
-1.14335264490883e-1  +5.38291888891147e-1
#

# n=3
+4.46483308008168e+1  -2.14240821614804e+1  -1.02242486393363e+1

+1.55740572903600e+0  +1.18264956574805e+0  +7.50716297686185e-1
-7.89042907035766e-1  +2.14233171396218e+0  -4.70215686814650e-1
-1.39710367762110e+0  +1.08419702950204e-1  +1.10241654934139e+0
#