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