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); INT s = 38; [s,s]R m; FOR j TO s DO FOR k FROM j TO s DO m[j, k] := next random; m[ k, j] := m[j, k] OD OD; [s]R evl; [s, s]R evc; eigen(m, evl, evc); R me := 0; # below: i is which eigenvector we consider; j is which component of that vector we consider; k is which summand of that component. # FOR i TO s DO FOR j TO s DO R t := 0; FOR k TO s DO t +:= m[j, k] * evc[i, k] OD; R er = ABS(evl[i]*evc[i, j] - t); IF er > me THEN me := er FI OD OD; print((me, newline, evl))