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 / REGGAUSS.LZB / REGGAUSS.ÌIB
Text File  |  2000-06-30  |  3KB  |  103 lines
  1. { Solution of Simultaneous Linear Equations  Ax=b Using Gaussian Elimination
  2.   Method with No Pivoting - Suitable for Symmetric matrices but less efficient
  3.   than Choelski decomp for positive semi definite matrices - outputs LU
  4.   factorization in A suitable for InverseDiagonal subroutine below - slicing
  5.   used for speed.
  6.  
  7.   The types are assumed defined in the calling program:
  8.  
  9.   Const
  10.      N=...; (maximum dimensions of coefficient matrix A)
  11.   Type
  12.      SUBS= 1..N;
  13.      RVEC = Array [SUBS] Of Real;
  14.      NBYN = Array [SUBS] Of RVEC;}
  15.  
  16. { function used by both routines }
  17. Function Dotprod(u, v: RVEC; k, n: SUBS): Real;
  18. Var
  19.   i: SUBS;
  20.   sum: Real;
  21. Begin
  22.   sum:=0; For i:=k To n Do sum:=sum+u[i]*v[i]; Dotprod:=sum;
  23. End; { of Dotprod }
  24.  
  25. Procedure Gausslice(
  26.      n     : SUBS;    { Dimension of Square matrix A }
  27. Var  a     : NBYN;    { Coefficient Matrix A returns LU decomposition }
  28.      b     : RVEC;    { Initial b's }
  29. Var  x     : RVEC;    { Solved X Values }
  30. Var  error : Boolean);{ error - not solved }
  31. label 999;
  32. const
  33.   assumezero=1E-20;
  34. Var
  35.   i, j: SUBS;
  36.   mult, pivot: Real;
  37.   state: (incrementingj, jisnminus1, zeropivot);
  38.  
  39. Procedure SubtractRow(Var u: RVEC; v: RVEC; m: Real; k, n:SUBS);
  40. Var
  41.   i: SUBS;
  42. Begin
  43.   For i:=k To n Do u[i]:=u[i]-m*v[i];
  44. End; { of SubtractRow }
  45.  
  46. Begin { of GAUSSLICE main }
  47.   j:=1; state:=incrementingj;
  48.   Repeat
  49.     pivot:=a[j,j];
  50.     If Abs(pivot)<=assumezero Then state:=zeropivot Else
  51.       Begin
  52.       For i:=j+1 To n Do
  53.         Begin
  54.         mult:=a[i,j]/pivot;
  55.         If Abs(mult)>assumezero Then
  56.           Begin
  57.           a[i,j]:=mult; SubtractRow(a[i],a[j],mult,j+1,n);
  58.           b[i]:=b[i]-mult*b[j];
  59.           End;
  60.         End;
  61.       If j=n-1 Then state:=jisnminus1 Else j:=j+1;
  62.     End;
  63.   Until state<>incrementingj;
  64.   Case state Of
  65.     zeropivot: error:=true;
  66.     jisnminus1:
  67.       Begin
  68.       error:=Abs(a[n,n])<=assumezero;
  69.       If not error Then
  70.         Begin
  71.         x[n]:=b[n]/a[n,n];
  72.         For i:=n-1 Downto 1 Do x[i]:=(b[i]-Dotprod(a[i],x,i+1,n))/a[i,i];
  73.         End;
  74.       End;
  75.     End; { of case };
  76. End; { of GAUSSLICE }
  77.  
  78. { Companion program to GAUSSLICE uses repeated calls To LUSolving with Ax=b
  79.   coming from columns of AinvA=Identity to solve for inverse matrix columns in
  80.   this version only diagonals are computed for statistical tests on x's }
  81.  
  82. Procedure InverseDiagonals(
  83.      n     : SUBS;    { Dimension of Square matrix A }
  84.      a     : NBYN;    { LU decomposition }
  85. Var  x     : RVEC);   { Inverse diagonals }
  86. Var
  87.   i, j: SUBS;
  88.   b: RVEC;
  89. Begin
  90.   For i:=1 To n Do
  91.     Begin
  92.     b[i]:=1;
  93.     For j:=i+1 To n Do b[j]:=b[j]-Dotprod(a[j],b,i,j-1);
  94.     x[n]:=b[n]/a[n,n];
  95.     For j:=n-1 Downto i Do x[j]:=(b[j]-Dotprod(a[j],x,j+1,n))/a[j,j];
  96.     End;
  97. End; { of InverseDiagonals }
  98. ngj;
  99.   Case state Of
  100.     zeropivot: error:=true;
  101.     jisnminus1:
  102.       Begin
  103.       erro