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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
|
unit DGCMatrices;
{$mode delphi}
{$DEFINE WITH_GHTOPO}
{.$UNDEF WITH_GHTOPO}
interface
uses
Classes, SysUtils, math;
type TDGCMatrixRow = array of double;
type
{ TDGCMatrix }
TDGCMatrixSquared = record
strict private
function getNbRows(): integer;
function getNbCols(): integer;
private
FName : string;
FRows : array of TDGCMatrixRow;
procedure MultiplyRowByScalar(const ARow: integer; const Value: double);
procedure SwapRows(const ARow1, ARow2: integer);
procedure SwapColumns(const ACol1, ACol2: integer);
public
property Name: string read FName;
property NbRows: integer read getNbRows;
property NbCols: integer read getNbCols;
procedure Redim(const Name: string; const NbRows: integer);
procedure setValue(const ARow, ACol: integer; const V: double);
function getValue(const ARow, ACol: integer): double;
procedure FillRandom();
procedure LoadIdentity();
function Inverse(): boolean;
function ProduitBy(const M: TDGCMatrixSquared): TDGCMatrixSquared;
function DebugString(): string;
function DebugStringOfRow(const ARow: integer): string;
end;
implementation
uses
{$IFDEF WITH_GHTOPO}
Common, // pour AfficherMessageErreur,
{$ENDIF WITH_GHTOPO}
DGCDummyUnit; // unité vide sinon erreur 'Fin du code source non trouvée'
procedure Swap(var V1, V2: Double); overload;
var Tmp: Double;
begin
Tmp := V1;
V1 := V2;
V2 := Tmp;
end;
{$IFNDEF WITH_GHTOPO}
procedure AfficherMessageErreur(const Msg: string);
begin
WriteLn(Msg);
end;
{$ENDIF}
{ TDGCMatrixSquared }
function TDGCMatrixSquared.getNbRows(): integer;
begin
Result := Length(FRows);
end;
function TDGCMatrixSquared.getNbCols(): integer;
begin
Result := Length(FRows[0]);
end;
procedure TDGCMatrixSquared.Redim(const Name: string; const NbRows: integer);
var
i: Integer;
begin
FName := Name;
setLength(FRows, NbRows);
for i := 0 to NbRows-1 do SetLength(FRows[i], 2 * NbRows);
end;
(*
procedure TDGCMatrix.CopyTo(out MR: TDGCMatrix);
var
i, j: Integer;
begin
MR.Redim(self.name + '(Copy)', self.NbRows, self.NbCols);
for i := 1 to self.NbRows do
for j := 1 to self.NbCols do
MR.setValue(i, j, self.getValue(i, j));
end;
//*)
procedure TDGCMatrixSquared.setValue(const ARow, ACol: integer; const V: double);
var
R: TDGCMatrixRow;
begin
R := FRows[ARow-1];
R[ACol-1] := V;
end;
function TDGCMatrixSquared.getValue(const ARow, ACol: integer): double;
var
R: TDGCMatrixRow;
begin
R := FRows[ARow-1];
Result := R[ACol-1];
end;
procedure TDGCMatrixSquared.FillRandom();
var
i, j: Integer;
begin
AfficherMessageErreur('FillRandom');
for i := 1 to self.NbRows do
for j := 1 to self.NbRows do
self.setValue(i, j, Random);
// matrice augmentée mise à l'identité
for i := 1 to self.NbRows do
for j := 1 to self.NbRows do
begin
self.setValue(i, j + NbRows, Ifthen(i = j, 1.0, 0.0));
end;
end;
procedure TDGCMatrixSquared.LoadIdentity();
var
i, j: Integer;
begin
for i := 1 to self.NbRows do
for j := 1 to self.NbCols do
self.setValue(i, j, Ifthen(i=j, 1.0, 0.0));
end;
function TDGCMatrixSquared.Inverse(): boolean;
var
p, r, j, k, i: Integer;
Q1, Q2: Double;
function IdxMaxOfAij(const Col, p, q: integer): integer;
var
QMax: Extended;
ii: Integer;
ww: Double;
begin
AfficherMessageErreur(Format('-- Passe dans IdxMaxOfAij: Col= %d, p=%d, q=%d', [Col, p, q]));
Result := 0;
QMax := -Infinity;
for ii := p to q do
begin
ww := abs(self.getValue(ii, Col));
if (ww > QMax)then
begin
QMax := ww;
Result := ii;
end;
AfficherMessageErreur(Format('-- IdxMaxOfAij: Col= %d, ii: %d, ww = %.3f, Result: %d, p=%d, q=%d', [Col, ii, ww, Result, p, q]));
end;
end;
begin
result := false;
//self.CopyTo(MR);
self.DebugString();
// (r est l'indice de ligne du dernier pivot trouvé)
for j := 1 to self.getNbCols() do
begin
r := 0;
k := IdxMaxOfAij(j, r+1, self.getNbRows());
if (not (IsZero(self.getValue(k, j)))) then
begin
r += 1;
self.MultiplyRowByScalar(k, 1 / self.getValue(k, j)); // Diviser la ligne k par A[k,j] (On normalise la ligne de pivot de façon que le pivot prenne la valeur 1)
if (k <> r) then self.SwapRows(k, r);
for i := 1 to self.getNbRows() do
begin
if (i <> r) then
begin
for p := 1 to self.getNbCols() do // Soustraire à la ligne i la ligne r multipliée par A[i,j] (de façon à annuler A[i,j])
begin
Q1 := self.getValue(i, p);
Q2 := self.getValue(r, p) * self.getValue(i, j);
self.setValue(i, p, Q1 - Q2);
end;
end;
end;
end;
end;
(*
Gauss-Jordan
r = 0 (r est l'indice de ligne du dernier pivot trouvé)
Pour j de 1 jusqu'à m (j décrit tous les indices de colonnes)
| Rechercher max(|A[i,j]|, r+1 ≤ i ≤ n). Noter k l'indice de ligne du maximum
| (A[k,j] est le pivot)
| Si A[k,j]≠0 alors (A[k,j] désigne la valeur de la ligne k et de la colonne j)
| | r=r+1 (r désigne l'indice de la future ligne servant de pivot)
| | Diviser la ligne k par A[k,j] (On normalise la ligne de pivot de façon que le pivot prenne la valeur 1)
| | Si k≠r alors
| | | Échanger les lignes k et r (On place la ligne du pivot en position r)
| | Fin Si
| | Pour i de 1 jusqu'à n (On simplifie les autres lignes)
| | | Si i≠r alors
| | | | Soustraire à la ligne i la ligne r multipliée par A[i,j] (de façon à annuler A[i,j])
| | | Fin Si
| | Fin Pour
| Fin Si
Fin Pour
Fin Gauss-Jordan
//*)
self.DebugString();
end;
function TDGCMatrixSquared.ProduitBy(const M: TDGCMatrixSquared): TDGCMatrixSquared;
var
s: Extended;
i, j, k: Integer;
begin
Result.Redim(Format('%sx%s', [self.Name, M.Name]), self.NbRows);
for i := 1 to self.NbRows do
for j := 1 to self.NbRows do
begin
s := 0.00;
for k := 1 to self.NbCols do s += self.getValue(i, k) * M.getValue(k, j);
Result.setValue(i, j, s);
end;
end;
function TDGCMatrixSquared.DebugString(): string;
var
EWE: String;
i, j: Integer;
begin
result := '';
AfficherMessageErreur(format('%s: %d lignes, %d colonnes', [self.Name, self.getNbRows(), self.getNbCols()]));
for i := 1 to self.getNbRows() do
begin
EWE :='';
for j := 1 to self.getNbCols() do EWE += format('%.8f ', [self.getValue(i, j)]);
AfficherMessageErreur(EWE);
end;
end;
function TDGCMatrixSquared.DebugStringOfRow(const ARow: integer): string;
var
j: Integer;
begin
Result := format('Ligne %d: ', [ARow]);
for j := 1 to self.getNbCols() - 1 do Result += format('%.8f ', [self.getValue(ARow, j)]);
end;
procedure TDGCMatrixSquared.MultiplyRowByScalar(const ARow: integer; const Value: double);
var
j: Integer;
WU: Double;
begin
for j := 1 to self.NbCols do
begin
WU := self.getValue(ARow, j);
self.setValue(ARow, j, WU * Value);
end;
end;
procedure TDGCMatrixSquared.SwapRows(const ARow1, ARow2: integer);
var
c: Integer;
v1, v2: Double;
begin
for c := 1 to self.getNbCols() do
begin
v1 := self.getValue(ARow1, c);
v2 := self.getValue(ARow2, c);
self.setValue(ARow1, c, v2);
self.setValue(ARow2, c, v1);
end;
end;
procedure TDGCMatrixSquared.SwapColumns(const ACol1, ACol2: integer);
var
v1, v2: Double;
r: Integer;
begin
for r := 1 to self.getNbRows() do
begin
v1 := self.getValue(r, ACol1);
v2 := self.getValue(r, ACol2);
self.setValue(r, ACol1, v2);
self.setValue(r, ACol2, v1);
end;
end;
end. |
