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

    {
                 ÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛ
                 ÛÛÛÝÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÞÛÛÛ±±
                 ÛÛÛÝÛÛ                                      ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ           Matrix inversion           ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ                                      ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ           Aleksandar Dlabac          ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ    (C) 1995. Dlabac Bros. Company    ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ    ------------------------------    ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ      adlabac@urcpg.urc.cg.ac.yu      ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ      adlabac@urcpg.pmf.cg.ac.yu      ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ                                      ÛÛÞÛÛÛ±±
                 ÛÛÛÝßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßÞÛÛÛ±±
                 ÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛ±±
                   ±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±
    }

       {   Gausian algorithm for matrix inversion  }

       Const MaxN = 10;

       Type Row    = array [1..MaxN] of real;
            Matrix = array [1..MaxN] of Row;

       Var I, J        : integer;
           InversionOK : Boolean;
           A           : Matrix;

       { A is matrix to be inverted, and N is dimension of matrix. }
       { Result is return in A.                                    }

       Procedure MatrixInversion (Var A:Matrix; N:integer);
         Var I, J, K : integer;
             Factor  : real;
             Temp    : Row;
             B       : Matrix;
           Begin
             InversionOK:=False;
             For I:=1 to N do
               For J:=1 to N do
                 If I=J then
                   B [I,J]:=1
                        else
                   B [I,J]:=0;
             For I:=1 to N do
               Begin
                 For J:=I+1 to N do
                   If Abs (A [I,I])<Abs (A [J,I]) then
                     Begin
                       Temp:=A [I];
                       A [I]:=A [J];
                       A [J]:=Temp;
                       Temp:=B [I];
                       B [I]:=B [J];
                       B [J]:=Temp
                     End;
                 If A [I,I]=0 then Exit;
                 Factor:=A [I,I];
                 For J:=N downto 1 do
                   Begin
                     B [I,J]:=B [I,J]/Factor;
                     A [I,J]:=A [I,J]/Factor
                   End;
                 For J:=I+1 to N do
                   Begin
                     Factor:=-A [J,I];
                     For K:=1 to N do
                       Begin
                         A [J,K]:=A [J,K]+A [I,K]*Factor;
                         B [J,K]:=B [J,K]+B [I,K]*Factor
                       End
                   End
               End;
             For I:=N downto 2 do
               Begin
                 For J:=I-1 downto 1 do
                   Begin
                     Factor:=-A [J,I];
                     For K:=1 to N do
                       Begin
                         A [J,K]:=A [J,K]+A [I,K]*Factor;
                         B [J,K]:=B [J,K]+B [I,K]*Factor
                       End
                   End
               End;
             A:=B;
             InversionOK:=True
           End;

       Begin
         A [1,1]:=3; A [1,2]:=-2; A [1,3]:=0;
         A [2,1]:=1; A [2,2]:=4;  A [2,3]:=-1;
         A [3,1]:=7; A [3,2]:=5;  A [3,3]:=-3;
         Writeln ('Matrix A is: ');
         For I:=1 to 3 do
           Begin
             For J:=1 to 3 do
               Write (A [I,J]:6:2);
             Writeln
           End;
         MatrixInversion (A,3);
         If not InversionOK then
           Writeln ('Matrix cannot be inverted.')
                            else
           Begin
             Writeln ('Inverse matrix of A (A^(-1)) is: ');
             For I:=1 to 3 do
               Begin
                 For J:=1 to 3 do
                   Write (A [I,J]:6:2);
                 Writeln
               End
           End
       End.

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