MODE R = REAL;
PROC eigen = (REF[,] R m, REF[]R l, REF[,] R vc) VOID: (
  INT n = 1 UPB m; [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);

([3, 3]R g; FOR i TO 3 DO FOR j TO 3 DO g[i, j] := (i=j|1|1/2) OD OD;
[3, 3]R vc; [3]R e; eigen(g, e, vc);
FOR i TO 3 DO FOR j TO 3 DO print(g[i, j]) OD; print(newline) OD;
FOR i TO 3 DO print(e[i]); FOR j TO 3 DO print(vc[i, j]) OD; print(newline) OD)

#
+1.00000000000000e  +0+5.00000000000000e  -1+5.00000000000000e  -1
+5.00000000000000e  -1+1.00000000000000e  +0+5.00000000000000e  -1
+5.00000000000000e  -1+5.00000000000000e  -1+1.00000000000000e  +0
+2.00000000000000e  +0+5.77350269189626e  -1+5.77350269189626e  -1+5.77350269189626e  -1
+5.00000000000000e  -1-7.07106781186547e  -1+7.07106781186548e  -1+0.00000000000000e  +0
+5.00000000000000e  -1-4.08248290463863e  -1-4.08248290463863e  -1+8.16496580927726e  -1
#