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 / MATINV.LZB / MATINV.ÌIB
Text File  |  2000-06-30  |  3KB  |  101 lines
  1. { Gauss Jordan matrix inversion }
  2. procedure matinv (
  3.             var b: NBYN;        { square matrix of coefficients inverted}
  4.                 n: SUBS;        { order of matrix }
  5.           var det: real;        { determinant of matrix }
  6.         var error: boolean);    { true, i.e.=1, if matrix singular }
  7.  
  8. label 99;
  9.  
  10. var
  11.    index : array[SUBS,1..3] of integer;
  12.    i,j,k,l,irow,icol,l1 : integer;
  13.    pivot,hold,sum,t,ab,big : real;
  14.  
  15. procedure swap(var a,b: real);
  16. var     hold : real;
  17. begin
  18.   hold:=a; a:=b; b:=hold
  19. end;
  20.  
  21. procedure minv2;
  22. label   98;
  23. var     i,j,k,l,l1 : integer;
  24. begin   { procedure minv2 actual start of matinv }
  25.   error:=false; det:=1.0;
  26.   for i:=1 to n do index[i,3]:=0;
  27.   for i:=1 to n do
  28.     begin
  29.     { search for largest element }
  30.     big:=0.0;
  31.     for j:=1 to n do
  32.       begin
  33.         if index[j,3]<>1 then
  34.           begin
  35.             for k:=1 to n do
  36.               begin
  37.                 if index[k,3]>1 then
  38.                   begin
  39.                     writeln('ERROR: matrix is singular');
  40.                     error:=true; goto 98 { abort }
  41.                   end;
  42.                if index[k,3]<1 then
  43.                   if abs(b[j,k])>big then
  44.                     begin
  45.                       irow:=j; icol:=k; big:=abs(b[j,k])
  46.                     end
  47.              end        { k-loop }
  48.         end
  49.   end;          { j-loop }
  50.   index[icol,3]:=index[icol,3]+1; index[i,1]:=irow; index[i,2]:=icol;
  51.   { interchange rows to put pivot on diagonal }
  52.   if irow<>icol then
  53.     begin
  54.       det:=-det;
  55.       for l:=1 to n do
  56.         swap(b[irow,l],b[icol,l]);
  57.   end;
  58.   { divide pivot row by pivot column }
  59.   pivot:=b[icol,icol]; det:=det*pivot; b[icol,icol]:=1.0;
  60.   for l:=1 to n do b[icol,l]:=b[icol,l]/pivot;
  61.  
  62.   { reduce nonpivot rows }
  63.   for l1:=1 to n do
  64.     begin
  65.       if l1<>icol then
  66.         begin
  67.           t:=b[l1,icol]; b[l1,icol]:=0.0;
  68.           for l:=1 to n do b[l1,l]:=b[l1,l]-b[icol,l]*t;
  69.       end       { if l1<>icol }
  70.      end
  71.     end;        { i-loop }
  72. 98:
  73. end;            { minv2 }
  74.  
  75. begin           { gaus-jordan main program }
  76.   minv2;        { first half of matinv }
  77.   if error then goto 99;
  78.      { interchange columns }
  79.   for i:=1 to n do
  80.     begin
  81.       l:=n-i+1;
  82.       if index[l,1]<>index[l,2] then
  83.         begin
  84.           irow:=index[l,1]; icol:=index[l,2];
  85.           for k:=1 to n do swap(b[k,irow],b[k,icol]);
  86.       end       { if index }
  87.   end;          { i-loop }
  88. for k:=1 to n do
  89.   if index[k,3]<>1 then
  90.     begin
  91.       writeln(chr(7),'ERROR: matrix is singular');
  92.       error:=true; goto 99 { abort }
  93.     end;
  94. 99:
  95. end;    { procedure matinv }
  96.  
  97. ot column }
  98.   pivot:=b[icol,icol]; det:=det*pivot; b[icol,icol]:=1.0;
  99.   for l:=1 to n do b[icol,l]:=b[icol,l]/pivot;
  100.  
  101.   { r