| 12
 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