home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / DISK_425.ZIP / GAUSSJ.LIB < prev    next >
Encoding:
Text File  |  1986-06-20  |  3.4 KB  |  167 lines
  1.    {87}
  2.  
  3. procedure gaussj
  4.        (var    b: ary2s;    { square matrix of coefficients }
  5.         y: arys;    { constant vector }
  6.     var  coef: arys;    { solution vector }
  7.          ncol: integer;    { order of matrix }
  8.     var error: boolean);    { true if matrix singular }
  9.  
  10. { Gauss Jordan matrix inversion and solution }
  11.  
  12. { B(n,n) coefficient matrix becomes inverse }
  13. { Y(n) original constant vector }
  14. { W(n,m) constant vector(s) become solution vector }
  15. { DETERM is the determinant }
  16. { ERROR=1 if singular }
  17. { INDEX(n,3) }
  18. { NV is number of constant vectors }
  19.  
  20. label 99;
  21.  
  22. var
  23.     w    : array[1..maxc,1..maxc] of real;
  24.     index    : array[1..maxc,1..3] of integer;
  25.     i,j,k,l,nv,
  26.     irow,icol,
  27.     n,l1    : integer;
  28.     determ,pivot,
  29.     hold,sum,t,
  30.     ab,big    : real;
  31.  
  32.  
  33.  
  34.  
  35. procedure swap(var a,b: real);
  36. var    hold    : real;
  37.  
  38. begin    { swap }
  39.   hold:=a;
  40.   a:=b;
  41.   b:=hold
  42. end;        { procedure swap }
  43.  
  44.  
  45. procedure gausj2;
  46. label    98;
  47. var    i,j,k,l,l1    : integer;
  48.  
  49.  
  50. procedure gausj3;
  51.  
  52. var    l    : integer;
  53.  
  54. begin    { procedure gausj3 }
  55.     { interchange rows to put pivot on diagonal }
  56.   if irow<>icol then
  57.     begin
  58.       determ:=-determ;
  59.       for l:=1 to n do
  60.     swap(b[irow,l],b[icol,l]);
  61.       if nv>0 then
  62.         for l:=1 to nv do
  63.       swap(w[irow,l],w[icol,l])
  64.   end    { if iroe<>icol }
  65. end;    { gausj3 }
  66.  
  67. begin    { procedure gausj2 }
  68.     { actual start of gaussj }
  69.   error:=false;
  70.   nv:=1;    { single constant vector }
  71.   n:=ncol;
  72.   for i:=1 to n do
  73.     begin
  74.       w[i,1]:=y[i];    { copy constant vector }
  75.       index[i,3]:=0
  76.     end;
  77.   determ:=1.0;
  78.   for i:=1 to n do
  79.     begin
  80.     { search for largest element }
  81.     big:=0.0;
  82.     for j:=1 to n do
  83.       begin
  84.     if index[j,3]<>1 then
  85.       begin
  86.         for k:=1 to n do
  87.           begin
  88.         if index[k,3]>1 then
  89.           begin
  90.             writeln('ERROR: matrix is singular');
  91.             error:=true;
  92.             goto 98        { abort }
  93.                   end;
  94.                if index[k,3]<1 then
  95.           if abs(b[j,k])>big then
  96.             begin
  97.               irow:=j;
  98.               icol:=k;
  99.               big:=abs(b[j,k])
  100.             end
  101.          end    { k-loop }
  102.     end
  103.   end;        { j-loop }
  104.  
  105.   index[icol,3]:=index[icol,3]+1;
  106.   index[i,1]:=irow;
  107.   index[i,2]:=icol;
  108.  
  109.   gausj3;    { further subdivision of gaussj }
  110.         { divide pivot row by pivot column }
  111.   pivot:=b[icol,icol];
  112.   determ:=determ*pivot;
  113.   b[icol,icol]:=1.0;
  114.  
  115.   for l:=1 to n do
  116.     b[icol,l]:=b[icol,l]/pivot;
  117.   if nv>0 then
  118.     for l:=1 to nv do
  119.       w[icol,l]:=w[icol,l]/pivot;
  120.  
  121.     { reduce nonpivot rows }
  122.  
  123.   for l1:=1 to n do
  124.     begin
  125.       if l1<>icol then
  126.     begin
  127.       t:=b[l1,icol];
  128.       b[l1,icol]:=0.0;
  129.       for l:=1 to n do
  130.         b[l1,l]:=b[l1,l]-b[icol,l]*t;
  131.       if nv>0 then
  132.         for l:=1 to nv do
  133.           w[l1,l]:=w[l1,l]-w[icol,l]*t;
  134.       end    { if l1<>icol }
  135.      end
  136.     end;    { i-loop }
  137. 98:
  138. end;        { gausj2 }
  139.  
  140. begin        { gaus-jordan main program }
  141.   gausj2;    { first half of gaussj }
  142.   if error then goto 99;
  143.      { interchange columns }
  144.   for i:=1 to n do
  145.     begin
  146.       l:=n-i+1;
  147.       if index[l,1]<>index[l,2] then
  148.     begin
  149.       irow:=index[l,1];
  150.       icol:=index[l,2];
  151.       for k:=1 to n do
  152.         swap(b[k,irow],b[k,icol])
  153.       end    { if index }
  154.   end;        { i-loop }
  155. for k:=1 to n do
  156.   if index[k,3]<>1 then
  157.     begin
  158.       writeln(chr(7),'ERROR: matrix is singular');
  159.       error:=true;
  160.       goto 99        { abort }
  161.     end;
  162.   for i:=1 to n do
  163.     coef[i]:=w[i,1];
  164. 99:
  165. end;    { procedure gaussj }
  166.  
  167.