1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
   |  
 
unit mat17_inv;
 
interface
 
const   C_17  = 17;
const   C_32  = 32;
 
type Matrice  = array[1..C_32,1..C_32] of Extended;
 
function inverse_M( rg : integer; M: matrice; var M_1 : Matrice) : boolean;
 
implementation
 
procedure Somme(rg : integer; MatA,MatB:Matrice;  var MatC:Matrice);
var i,j       : byte;
   begin
   for i:=1 to rg do  for j:=1 to rg do
      MatC[i,j]               :=  MatA[i,j] + MatB[i,j]
   end;
procedure Difference(rg : integer; MatA,MatB:Matrice;var MatC:Matrice);
var i,j       : byte;
   begin
   for i:=1 to rg do for j:=1 to rg do
      MatC[i,j]               :=  MatA[i,j] - MatB[i,j]
   end;
procedure Produit(rg : integer; MatA,MatB:Matrice;var MatC:Matrice);
var i,j,k     : byte;
   begin
   fillchar(MatC,sizeof(MatC),0);
   for i:=1 to rg do for j:=1 to rg do for k:=1 to rg do
      MatC[i,j]               :=  MatC[i,j]+MatA[i,k] * MatB[k,j]
   end;
function Norme_Ligne(rg : integer; B:Matrice):Extended;
var Tmp,Tmp2  : Extended;
    i,j       : byte;
   begin
   Tmp                        :=  0;
   for i:=1 to rg do
      begin
      Tmp2                    :=  0;
      for j:=1 to rg do
         Tmp2                 :=  Tmp2 + abs(B[i,j]);
      if Tmp2>Tmp then
         Tmp                  :=  Tmp2
      end;
   Norme_Ligne                :=  Tmp
   end;
function Norme_Colonne(rg : integer; B:Matrice):Extended;
var Tmp,Tmp2  : Extended;
    i,j       : byte;
   begin
   Tmp                        :=  0;
   for i:=1 to rg do
      begin
      Tmp2                    :=  0;
      for j:=1 to rg do
         Tmp2                 :=  Tmp2 + abs(B[j,i]);
      if Tmp2>Tmp then
         Tmp                  :=  Tmp2
      end;
   Norme_Colonne              :=  Tmp
   end;
procedure Calcule_E(rg : integer; var A, B, E, Id : matrice);
var C         : Matrice;
   begin
   Produit(rg, B, A, C);
   Difference(rg, Id, C, E)
   end;
procedure Init(rg : integer; A:Matrice;var B,Id :Matrice;var OK:boolean);
var i,j       : byte;
    T         : Extended;
   begin
   fillchar(Id,SizeOf(Id),0);
   T                          :=  Norme_Ligne(rg, A) * Norme_Colonne(rg, A);
   if abs(T)>1E-100 then
      begin
      OK                      :=  True;
      for i:=1 to rg do
         begin
         Id[i,i]              :=  1;
         B[i,i]               :=  A[i,i]/T;
         for j:=i+1 to rg do
            begin
            B[i,j]            :=  A[j,i]/T;
            B[j,i]            :=  A[i,j]/T
            end
         end
      end
   else
      OK                      :=  false
   end;
procedure Calcule_B( rg : integer; var B,E,Id : Matrice);
var C,D       : Matrice;
   begin
   Somme(rg,Id,E,D);
   Produit(rg,D,B,C);
   B                          :=  C
   end;
const Seuil_  = 1e-17;
const max_CPT = 225;
function Inverse_M(rg : integer; M : Matrice; var  M_1 : Matrice) : boolean;
var Id,E,B    : Matrice;
    Old,New   : Extended;
    Cpt       : byte;
    Ok        : boolean;
   begin
   Init(rg, M, B, Id, Ok);
   if  OK then
       begin
       Calcule_E(rg, M, B, E, Id );
       Old                    :=  Norme_Ligne(rg, E);
       New                    :=  Old;
       Cpt                    :=  0;
       while not ( (New < Seuil_) or ((New > Old) and (Cpt>max_CPT))) do
          begin
          Old                 :=  New;
          Calcule_B(rg, B, E, Id);
          Calcule_E(rg, M, B, E, Id );
          New                 :=  Norme_Ligne(rg, E);
          inc(Cpt)
          end;
       Ok                     :=  New < Seuil_;
       M_1                    :=  B
       end;
   Inverse_M                  :=  Ok
   end;
end. | 
Partager