[Back to MATH SWAG index] [Back to Main SWAG index] [Original]

    {
    From: "Victor B. Wagner" <vitus@agropc.msk.su>
    ************************************************************
    *                   ALGEBRA.PAS                            *
    *           a simple matrix algebra unit                   *
    *           Turbo Pascal 4.0 or higher                     *
    *           Copyright (c) by Vitus Wagner,1990             *
    ************************************************************
    }
    unit algebra;
    interface
       const MaxN=30;{You can increase it up to 100,not greater
                     but each matrix variable would have size of
                     sqr(MaxN)*sizeof(Real). It is possible to write
                     unit for work with dinamically sized matrices,
                     but i have no needs to do this.
                     You can work with matrices with size less that MaxN,
                     but while you work with this unit you must allocate
                     memory for matrix MaxN x MaxN and leave rest of space
                     unised}
       type vector=array[1..MaxN]of real;
            matrix=array[1..MaxN,1..MaxN]of real;
            sett=set of 1..MaxN;
     var algebrerr:boolean;
     function scalar(a,b:vector;n:integer):real;
     {Scalar multiplication of vectors a and b of n components}
     procedure systeq(a:matrix;b:vector;var x:vector;n:integer);
     { solving n line equation system A*X=B by Gauss method}
     { sets algebrerr to True if system cannot be solved}
     procedure matmult(a,b:matrix;var c:matrix;n:integer);
     { multiplication of two NxN matrixs A and B.Result - matrix C AxB=C}
     procedure matadd(a,b:matrix;var c :matrix;n:integer);
     { addition of two NxN matrixs A+B=C }
     procedure matconst(c:real;a:matrix;var b:matrix;n:integer);
     { multiplication matrix A on constant c  cxA=B }
     procedure matcadd(c1:real;a1:matrix;c2:real;a2:matrix;var b:matrix;n:integer);
     { addition of two NxN matrixs with multiplication each of them on constant
       c1xA1+c2xA2=B }
     procedure matinv(a:matrix;var ainv:matrix;n:integer);
     { inversion of NxN matrix A}
     { sets algebrerr to True if matrix cannot be inverted}
     procedure matvec(a:matrix;b:vector;var c:vector;n:integer);
     { multiplication NxN matrix A to N-component vector B AxB=C}
     function det(a:matrix;n:integer):real;
     { determinant of NxN matrix }
     procedure compress(a:matrix;var b:matrix;n:integer;s:sett);
     { converse triangle matrix to simmetric,exclude rows and columns that is not
       in set s (type sett=set of 0..maxN)}
     function distance(a,b:vector;n:integer):real;
     { Calculate Euclide distantion in N-dimensioned space between A & B }
     Procedure Transpose(var A:Matrix;M,N:Integer);
     { Transpose MxN Matrix. Put result in the same place}
     Procedure EMatrix(var A:Matrix;N:Integer);
     {Fills matrix by 0 and main diagonal by 1}
    implementation
     function scalar(a,b:vector;n:integer):real;
      var i:integer;
          r:real;
      begin
       r:=0.0;
       for i:=1 to n do
        r:=r+a[i]*b[i];
       scalar:=r;
      end;
     procedure systeq(a:matrix;b:vector;var x:vector;n:integer);
      var i,j,k:integer;
          max:real;
      begin
      algebrerr:=false;
       { Conversion matrix to triangle }
       for i:=1 to n do
        begin
         max:=abs(a[i,i]);k:=i;
         for j:=succ(i) to n do
          if abs(a[j,i])>max then
           begin
            max:=abs(a[j,i]);k:=j
           end;
          if max<1E-10 then begin algebrerr:=true;exit end;
          if k<>i then
           begin
            for  j:=i to n do
             begin
              max:=a[k,j];
              a[k,j]:=a[i,j];
              a[i,j]:=max;
             end;
            max:=b[k];
            b[k]:=b[i];
            b[i]:=max;
           end;
          for j:=succ(i) to n do
           a[i,j]:=a[i,j]/a[i,i];
           b[i]:=b[i]/a[i,i];
          for j:=succ(i) to n do
           begin
            for k:=succ(i) to n do
             a[j,k]:=a[j,k]-a[i,k]*a[j,i];
            b[j]:=b[j]-b[i]*a[j,i];
           end;
        end;
         { X calculation}
         x[n]:=b[n];
         for i:=pred(n) downto 1 do
          begin
           max:=b[i];
           for j:=succ(i) to n do
            max:=max-a[i,j]*x[j];
           x[i]:=max;
          end;
      end;
     procedure matmult(a,b:matrix;var c:matrix;n:integer);
      var i,j,k:integer;r:real;
      begin
       for i:=1 to n do
        for j:=1 to n do
         begin
          r:=0.0;
          for k:=1 to n do
           r:=r+a[i,k]*b[k,j];
          c[i,j]:=r;
         end;
      end;
     procedure matadd(a,b:matrix;var c :matrix;n:integer);
      var i,j:integer;
      begin
       for i:=1 to n do
        for j:=1 to n do
         c[i,j]:=a[i,j]+b[i,j];
      end;
     procedure matinv(a:matrix;var ainv:matrix;n:integer);
      var i,j,k:integer;
          e:matrix;
          max:real;
      begin
       algebrerr:=false;
       { creating single matrix }
       for i:=1 to n do
        for j:=1 to n do
         e[i,j]:=0.0;
       for i:=1 to n do
        e[i,i]:=1.0;
       { Conversion matrix to triangle }
       for i:=1 to n do
    {1} begin
         max:=abs(a[i,i]);k:=i;
         for j:=succ(i) to n do
          if abs(a[j,i])>max then
    {2}    begin
            max:=abs(a[j,i]);k:=j
    {2}    end;
          if max<1E-10 then begin algebrerr:=true;exit end;
          if k<>i then
    {2}    begin
            for  j:=i to n do
    {3}      begin
              max:=a[k,j];
              a[k,j]:=a[i,j];
              a[i,j]:=max;
    {3}      end;
          for j:=1 to n do
    {3}    begin
            max:=e[k,j];
            e[k,j]:=e[i,j];
            e[i,j]:=max;
    {3}    end;
    {2}   end;
          for j:=succ(i) to n do
           a[i,j]:=a[i,j]/a[i,i];
           for k:=1 to n do
           e[i,k]:=e[i,k]/a[i,i];
          for j:=succ(i) to n do
    {2}    begin
            for k:=succ(i) to n do
             a[j,k]:=a[j,k]-a[i,k]*a[j,i];
            for k:=1 to n do
            e[j,k]:=e[j,k]-e[i,k]*a[j,i];
    {2}    end;
    {1} end;
         { ainv calculation}
        for k:=1 to n do
    {1} begin
         ainv[n,k]:=e[n,k];
         for i:=pred(n) downto 1 do
    {2}   begin
           max:=e[i,k];
           for j:=succ(i) to n do
            max:=max-a[i,j]*ainv[j,k];
           ainv[i,k]:=max;
    {2}   end;
    {1}  end;
      end;
     procedure matvec(a:matrix;b:vector;var c:vector;n:integer);
      var i,j:integer;r:real;
      begin
       for i:=1 to n do
        begin
         r:=0.0;
         for j:=1 to n do
          r:=r+a[i,j]*b[j];
         c[i]:=r;
        end;
      end;
     function det(a:matrix;n:integer):real;
      var i,j,k:integer;d:real;
      begin
       for i:=1 to pred(n) do
        begin
         if abs(a[i,i])<1E-10 then begin det:=0.0;exit end;
         for j:=succ(i) to n do
          begin
           d:=a[j,i]/a[i,i];
           for k:=i to n do
            a[j,k]:=a[j,k]-d*a[i,k];
          end;
        end;
       d:=1.0;
       for i:=1 to n do
        d:=d*a[i,i];
       det:=d;
      end;
     procedure matconst(c:real;a:matrix;var b:matrix;n:integer);
     var i,j:integer;
     begin
      for i:=1 to n do
       for j:=1 to n do
        b[i,j]:=c*a[i,j];
     end;
     procedure matcadd(c1:real;a1:matrix;c2:real;a2:matrix;var b:matrix;n:integer);
     var i,j:integer;
     begin
      for i:=1 to n do
       for j:=1 to n do
        b[i,j]:=c1*a1[i,j]+c2*a2[i,j];
     end;
     procedure compress(a:matrix;var b:matrix;n:integer;s:sett);
     var i,j,k,l:integer;
      begin
       k:=1;
       for i:=1 to pred(n) do
        if i in s then
         begin
          l:=1;
          b[k,k]:=a[i,i];
          for j:=succ(i) to n do
           if j in s then
            begin
             b[k,l]:=a[i,j];
             b[l,k]:=a[i,j];
             inc(l);
            end;
          inc(k);
         end;
      end;
     function distance(a,b:vector;n:integer):real;
     var i:integer;r:real;
     begin
      r:=0;
      for i:=1 to n do
       r:=r+sqr(a[i]-b[i]);
      distance:=sqrt(r);
     end;
    Procedure Transpose(var A:Matrix;M,N:Integer);
    var i,j:Integer;Tmp:Real;
    begin
     For i:=1 to n do
      for j:=i+1 to m do
       begin
        Tmp:=A[i,j];
        A[i,j]:=A[j,i];
        A[J,i]:=Tmp;
       end;
    end;
    Procedure EMatrix(var A:Matrix;N:Integer);
    var I:Integer;
    begin
      FillChar(A,SizeOf(A),0);
      For i:=1 to n do
       A[i,i]:=1.0;
    end;

    end.

[Back to MATH SWAG index] [Back to Main SWAG index] [Original]