home *** CD-ROM | disk | FTP | other *** search
/ POINT Software Programming / PPROG1.ISO / pascal / swag / math.swg / 0080_Matrix Algebra Unit.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1995-03-26  |  7.6 KB  |  284 lines
  1. {
  2. From: "Victor B. Wagner" <vitus@agropc.msk.su>
  3. ************************************************************
  4. *                   ALGEBRA.PAS                            *
  5. *           a simple matrix algebra unit                   *
  6. *           Turbo Pascal 4.0 or higher                     *
  7. *           Copyright (c) by Vitus Wagner,1990             *
  8. ************************************************************
  9. }
  10. unit algebra;
  11. interface
  12.    const MaxN=30;{You can increase it up to 100,not greater
  13.                  but each matrix variable would have size of
  14.                  sqr(MaxN)*sizeof(Real). It is possible to write
  15.                  unit for work with dinamically sized matrices,
  16.                  but i have no needs to do this.
  17.                  You can work with matrices with size less that MaxN,
  18.                  but while you work with this unit you must allocate
  19.                  memory for matrix MaxN x MaxN and leave rest of space
  20.                  unised}
  21.    type vector=array[1..MaxN]of real;
  22.         matrix=array[1..MaxN,1..MaxN]of real;
  23.         sett=set of 1..MaxN;
  24.  var algebrerr:boolean;
  25.  function scalar(a,b:vector;n:integer):real;
  26.  {Scalar multiplication of vectors a and b of n components}
  27.  procedure systeq(a:matrix;b:vector;var x:vector;n:integer);
  28.  { solving n line equation system A*X=B by Gauss method}
  29.  { sets algebrerr to True if system cannot be solved}
  30.  procedure matmult(a,b:matrix;var c:matrix;n:integer);
  31.  { multiplication of two NxN matrixs A and B.Result - matrix C AxB=C}
  32.  procedure matadd(a,b:matrix;var c :matrix;n:integer);
  33.  { addition of two NxN matrixs A+B=C }
  34.  procedure matconst(c:real;a:matrix;var b:matrix;n:integer);
  35.  { multiplication matrix A on constant c  cxA=B }
  36.  procedure matcadd(c1:real;a1:matrix;c2:real;a2:matrix;var b:matrix;n:integer);
  37.  { addition of two NxN matrixs with multiplication each of them on constant
  38.    c1xA1+c2xA2=B }
  39.  procedure matinv(a:matrix;var ainv:matrix;n:integer);
  40.  { inversion of NxN matrix A}
  41.  { sets algebrerr to True if matrix cannot be inverted}
  42.  procedure matvec(a:matrix;b:vector;var c:vector;n:integer);
  43.  { multiplication NxN matrix A to N-component vector B AxB=C}
  44.  function det(a:matrix;n:integer):real;
  45.  { determinant of NxN matrix }
  46.  procedure compress(a:matrix;var b:matrix;n:integer;s:sett);
  47.  { converse triangle matrix to simmetric,exclude rows and columns that is not
  48.    in set s (type sett=set of 0..maxN)}
  49.  function distance(a,b:vector;n:integer):real;
  50.  { Calculate Euclide distantion in N-dimensioned space between A & B }
  51.  Procedure Transpose(var A:Matrix;M,N:Integer);
  52.  { Transpose MxN Matrix. Put result in the same place}
  53.  Procedure EMatrix(var A:Matrix;N:Integer);
  54.  {Fills matrix by 0 and main diagonal by 1}
  55. implementation
  56.  function scalar(a,b:vector;n:integer):real;
  57.   var i:integer;
  58.       r:real;
  59.   begin
  60.    r:=0.0;
  61.    for i:=1 to n do
  62.     r:=r+a[i]*b[i];
  63.    scalar:=r;
  64.   end;
  65.  procedure systeq(a:matrix;b:vector;var x:vector;n:integer);
  66.   var i,j,k:integer;
  67.       max:real;
  68.   begin
  69.   algebrerr:=false;
  70.    { Conversion matrix to triangle }
  71.    for i:=1 to n do
  72.     begin
  73.      max:=abs(a[i,i]);k:=i;
  74.      for j:=succ(i) to n do
  75.       if abs(a[j,i])>max then
  76.        begin
  77.         max:=abs(a[j,i]);k:=j
  78.        end;
  79.       if max<1E-10 then begin algebrerr:=true;exit end;
  80.       if k<>i then
  81.        begin
  82.         for  j:=i to n do
  83.          begin
  84.           max:=a[k,j];
  85.           a[k,j]:=a[i,j];
  86.           a[i,j]:=max;
  87.          end;
  88.         max:=b[k];
  89.         b[k]:=b[i];
  90.         b[i]:=max;
  91.        end;
  92.       for j:=succ(i) to n do
  93.        a[i,j]:=a[i,j]/a[i,i];
  94.        b[i]:=b[i]/a[i,i];
  95.       for j:=succ(i) to n do
  96.        begin
  97.         for k:=succ(i) to n do
  98.          a[j,k]:=a[j,k]-a[i,k]*a[j,i];
  99.         b[j]:=b[j]-b[i]*a[j,i];
  100.        end;
  101.     end;
  102.      { X calculation}
  103.      x[n]:=b[n];
  104.      for i:=pred(n) downto 1 do
  105.       begin
  106.        max:=b[i];
  107.        for j:=succ(i) to n do
  108.         max:=max-a[i,j]*x[j];
  109.        x[i]:=max;
  110.       end;
  111.   end;
  112.  procedure matmult(a,b:matrix;var c:matrix;n:integer);
  113.   var i,j,k:integer;r:real;
  114.   begin
  115.    for i:=1 to n do
  116.     for j:=1 to n do
  117.      begin
  118.       r:=0.0;
  119.       for k:=1 to n do
  120.        r:=r+a[i,k]*b[k,j];
  121.       c[i,j]:=r;
  122.      end;
  123.   end;
  124.  procedure matadd(a,b:matrix;var c :matrix;n:integer);
  125.   var i,j:integer;
  126.   begin
  127.    for i:=1 to n do
  128.     for j:=1 to n do
  129.      c[i,j]:=a[i,j]+b[i,j];
  130.   end;
  131.  procedure matinv(a:matrix;var ainv:matrix;n:integer);
  132.   var i,j,k:integer;
  133.       e:matrix;
  134.       max:real;
  135.   begin
  136.    algebrerr:=false;
  137.    { creating single matrix }
  138.    for i:=1 to n do
  139.     for j:=1 to n do
  140.      e[i,j]:=0.0;
  141.    for i:=1 to n do
  142.     e[i,i]:=1.0;
  143.    { Conversion matrix to triangle }
  144.    for i:=1 to n do
  145. {1} begin
  146.      max:=abs(a[i,i]);k:=i;
  147.      for j:=succ(i) to n do
  148.       if abs(a[j,i])>max then
  149. {2}    begin
  150.         max:=abs(a[j,i]);k:=j
  151. {2}    end;
  152.       if max<1E-10 then begin algebrerr:=true;exit end;
  153.       if k<>i then
  154. {2}    begin
  155.         for  j:=i to n do
  156. {3}      begin
  157.           max:=a[k,j];
  158.           a[k,j]:=a[i,j];
  159.           a[i,j]:=max;
  160. {3}      end;
  161.       for j:=1 to n do
  162. {3}    begin
  163.         max:=e[k,j];
  164.         e[k,j]:=e[i,j];
  165.         e[i,j]:=max;
  166. {3}    end;
  167. {2}   end;
  168.       for j:=succ(i) to n do
  169.        a[i,j]:=a[i,j]/a[i,i];
  170.        for k:=1 to n do
  171.        e[i,k]:=e[i,k]/a[i,i];
  172.       for j:=succ(i) to n do
  173. {2}    begin
  174.         for k:=succ(i) to n do
  175.          a[j,k]:=a[j,k]-a[i,k]*a[j,i];
  176.         for k:=1 to n do
  177.         e[j,k]:=e[j,k]-e[i,k]*a[j,i];
  178. {2}    end;
  179. {1} end;
  180.      { ainv calculation}
  181.     for k:=1 to n do
  182. {1} begin
  183.      ainv[n,k]:=e[n,k];
  184.      for i:=pred(n) downto 1 do
  185. {2}   begin
  186.        max:=e[i,k];
  187.        for j:=succ(i) to n do
  188.         max:=max-a[i,j]*ainv[j,k];
  189.        ainv[i,k]:=max;
  190. {2}   end;
  191. {1}  end;
  192.   end;
  193.  procedure matvec(a:matrix;b:vector;var c:vector;n:integer);
  194.   var i,j:integer;r:real;
  195.   begin
  196.    for i:=1 to n do
  197.     begin
  198.      r:=0.0;
  199.      for j:=1 to n do
  200.       r:=r+a[i,j]*b[j];
  201.      c[i]:=r;
  202.     end;
  203.   end;
  204.  function det(a:matrix;n:integer):real;
  205.   var i,j,k:integer;d:real;
  206.   begin
  207.    for i:=1 to pred(n) do
  208.     begin
  209.      if abs(a[i,i])<1E-10 then begin det:=0.0;exit end;
  210.      for j:=succ(i) to n do
  211.       begin
  212.        d:=a[j,i]/a[i,i];
  213.        for k:=i to n do
  214.         a[j,k]:=a[j,k]-d*a[i,k];
  215.       end;
  216.     end;
  217.    d:=1.0;
  218.    for i:=1 to n do
  219.     d:=d*a[i,i];
  220.    det:=d;
  221.   end;
  222.  procedure matconst(c:real;a:matrix;var b:matrix;n:integer);
  223.  var i,j:integer;
  224.  begin
  225.   for i:=1 to n do
  226.    for j:=1 to n do
  227.     b[i,j]:=c*a[i,j];
  228.  end;
  229.  procedure matcadd(c1:real;a1:matrix;c2:real;a2:matrix;var b:matrix;n:integer);
  230.  var i,j:integer;
  231.  begin
  232.   for i:=1 to n do
  233.    for j:=1 to n do
  234.     b[i,j]:=c1*a1[i,j]+c2*a2[i,j];
  235.  end;
  236.  procedure compress(a:matrix;var b:matrix;n:integer;s:sett);
  237.  var i,j,k,l:integer;
  238.   begin
  239.    k:=1;
  240.    for i:=1 to pred(n) do
  241.     if i in s then
  242.      begin
  243.       l:=1;
  244.       b[k,k]:=a[i,i];
  245.       for j:=succ(i) to n do
  246.        if j in s then
  247.         begin
  248.          b[k,l]:=a[i,j];
  249.          b[l,k]:=a[i,j];
  250.          inc(l);
  251.         end;
  252.       inc(k);
  253.      end;
  254.   end;
  255.  function distance(a,b:vector;n:integer):real;
  256.  var i:integer;r:real;
  257.  begin
  258.   r:=0;
  259.   for i:=1 to n do
  260.    r:=r+sqr(a[i]-b[i]);
  261.   distance:=sqrt(r);
  262.  end;
  263. Procedure Transpose(var A:Matrix;M,N:Integer);
  264. var i,j:Integer;Tmp:Real;
  265. begin
  266.  For i:=1 to n do
  267.   for j:=i+1 to m do
  268.    begin
  269.     Tmp:=A[i,j];
  270.     A[i,j]:=A[j,i];
  271.     A[J,i]:=Tmp;
  272.    end;
  273. end;
  274. Procedure EMatrix(var A:Matrix;N:Integer);
  275. var I:Integer;
  276. begin
  277.   FillChar(A,SizeOf(A),0);
  278.   For i:=1 to n do
  279.    A[i,i]:=1.0;
  280. end;
  281.  
  282. end.
  283.  
  284.