home *** CD-ROM | disk | FTP | other *** search
/ ftp.barnyard.co.uk / 2015.02.ftp.barnyard.co.uk.tar / ftp.barnyard.co.uk / cpm / walnut-creek-CDROM / CPM / TURBOPAS / MAPSTATF.LBR / EIGEN.LZB / EIGEN.ÌIB
Text File  |  2000-06-30  |  2KB  |  60 lines
  1. Procedure eigen(Var a, v:NBYN; Var e:RVEC; t: Real;
  2.                 Var dv, nf, ot:Integer);
  3. (* eigenvalues and vectors by iteration and exhaustion - a is input
  4.   matrix of dimension dv, nf the number of vectors to extract, e is
  5.   returned eigenvalues and v is eigenvectors - if ip=1 below then matrix
  6.   is squared to refine iteration, eigenvalues only are then returned *)
  7. Var
  8.   x,y:RVEC;
  9.   j,k,l,it,ip:Integer;
  10.   s1,s2:Real;
  11.  
  12. Procedure MatSQ(Var b, a:NBYN; n:SUBS);
  13. var
  14.   i,j,k: SUBS;
  15.   sigma: Real;
  16. begin
  17.   for i:=1 to n do
  18.       for j:=1 to n do begin
  19.       b[i,j]:=0.0; for k:=1 to n do b[i,j]:=b[i,j]+a[i,k]*a[j,k]; end;
  20.   a:=b;
  21. end; { of MatSQ }
  22.  
  23. Begin
  24.   ip:=1; { set to zero if eigenvectors are reqd or to avoid powering }
  25.   l:=1;
  26.   If(ip=1) Then MatSQ(v,a,dv);
  27.   While (dv>=l) and (nf>=l) Do
  28.     Begin
  29.     it:=0; s1:=1.0;
  30.     For j:=1 To dv Do y[j]:=1.0;
  31.     While (s1>t) Do Begin it:=it+1;
  32.       For j:=1 To dv Do Begin
  33.         x[j]:=0.0; For k:=1 To dv Do x[j]:=x[j]+(a[j,k]*y[k]); End;
  34.       e[l]:=x[1]; s1:=0.0;
  35.       For j:=1 To dv Do Begin
  36.         v[j,l]:=x[j]/x[1]; s1:=s1+(Abs(y[j]-v[j,l])); y[j]:=v[j,l]; End;
  37.       If (it=25) and (s1>=s2) Then Begin Writeln;
  38.         Writeln('*** Roots Not Converging in Eigen - Error Abort ***');
  39.         Bdos(0); End;
  40.       s2:=s1;
  41.       End;
  42.     s1:=0.0;
  43.     For j:=1 To dv Do s1:=s1+Sqr(v[j,l]);
  44.     s1:=Sqrt(s1);
  45.     For j:=1 To dv Do v[j,l]:=v[j,l]/s1;
  46.     For j:=1 To dv Do
  47.       For k:=1 To dv Do a[j,k]:=a[j,k]-(v[j,l]*v[k,l]*e[l]);
  48.     If(l>1) And (e[l-1]-e[l]<=t) And (l<nf) Then
  49.       Begin Writeln; nf:=l;
  50.       Writeln('** Number of Eigenvalues Extracted Limited to ',nf:2,' **');
  51.       Writeln('** Variance Explained by Factor too small For Accuracy **');
  52.       wait(ot);
  53.       End;
  54.     If(ip=1) Then e[l]:=Sqrt(e[l]);
  55.     Writeln('Eigenvalue ',l,' =',e[l],' obtained in ',it,' iterations');
  56.     l:=l+1;
  57.     End;
  58. End; (* of eigen *)
  59.   For j:=1 To dv Do Begin
  60.         x[j]:=0.0; For k:=1 To dv Do x[j]:=x[j]+(a[j,k]*y[k]);