Algo Pivot de Gauss : résultats aberrants

**JP CASSOU** · 25/04/2022, 16h08

Bjr à vous,

Je suis une fois de plus confronté à une transcription d'algo qui ne fonctionne pas: pivot de Gauss dont voici l'algo.
Je me prends un 'Violation d'accès' incompréhensible alors que mes indices ne débordent pas

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
 
 Gauss-Jordan
     r = 0                                       (r est l'indice de ligne du dernier pivot trouvé)
     Pour j de 1 jusqu'à m
     |   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)
     |   |
     |   |   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

transcrit comme suit:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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.

et son utilisation:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
 
procedure TdlgTestsUnitaires.Button6Click(Sender: TObject);
var
  MtxA, MtxB, MtxC: TDGCMatrix;
  m, n, i, j: Integer;
begin
  MtxA.Redim('A', 6);
  MtxA.FillRandom();
 
  MtxA.Inverse();
end;

**BeanzMaster** · 25/04/2022, 21h21

Salut JP ca ne serai pas l'inverse ?:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
for j := 1 to self.getNbRows() do
  begin
    AfficherMessageErreur('(0001');
    k := IdxMaxOfAij(j, r+1, self.getNbCols());

Je fais quelques recherche, mais je ne trouve rien de très concret pour le moment

Essayes de posté ta question ici : https://www.developpez.net/forums/f2...mathematiques/ dans un des sous-forums, Y a pas mal de matheux

A+
Jérôme

**JP CASSOU** · 25/04/2022, 21h35

Envoyé par BeanzMaster

Salut JP ca ne serai pas l'inverse ?:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
for j := 1 to self.getNbRows() do
  begin
    AfficherMessageErreur('(0001');
    k := IdxMaxOfAij(j, r+1, self.getNbCols());

Ben non. IdxMaxOfAij() est censé retourner l'indice de pivot max pour la passe j. D'ailleurs, même en 'déroulant' cette fonction, je n'arrive à rien.

**BeanzMaster** · 25/04/2022, 21h46

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
for j := 1 to self.getNbCols() do
  begin
    r := 0;

r:=0 ne devrait-il pas être initialiser avant la boucle et non dans la boucle ?

**JP CASSOU** · 25/04/2022, 21h53

Envoyé par BeanzMaster

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
for j := 1 to self.getNbCols() do
  begin
    r := 0;

r:=0 ne devrait-il pas être initialiser avant la boucle et non dans la boucle ?

Dans l'algo Wikipedia, r est initialisé avant la boucle mais j'obtiens un segfault si je ne remets pas à zéro dans la boucle. r va de 1 à n (n = nombre de lignes)

**BeanzMaster** · 25/04/2022, 22h28

Une erreur d'indice donc ca ne serait pas plutôt for j := 1 to self.getNbCols()-1 do dans toutes tes boucles ?

Les tableaux vont de 0 à n-1

**JP CASSOU** · 25/04/2022, 22h34

Envoyé par BeanzMaster

Une erreur d'indice donc ca ne serait pas plutôt for j := 1 to self.getNbCols()-1 do dans toutes tes boucles ?

Les tableaux vont de 0 à n-1

Les tableaux contenant les matrices ne sont accessibles que par getValue() et setValue(), qui prennent des indices entre 1 et n:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
 
 
procedure TDGCMatrixSquared.setValue(const ARow, ACol: integer; const V: double);
var
  R: TDGCMatrixRow; // TDGCMatrixRow est un 'array of double'
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;

**BeanzMaster** · 26/04/2022, 12h06

Salut JP

Envoyé par JP CASSOU

Les tableaux contenant les matrices ne sont accessibles que par getValue() et setValue(), qui prennent des indices entre 1 et n:

Autant pour moi

J'ai trouvé ça https://www.geeksforgeeks.org/findin...jordan-method/, cela pourra peut-être t'aidé

A+
Jérôme

**anapurna** · 26/04/2022, 18h40

Salut

Dans la méthode ProduitBy

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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 ri := 1 to self.NbRows do
    for rj := 1 to self.NbRows do
    begin
      s := 0.00;
      for ck := 1 to self.NbCols do 
	    s += self.getValue(ri, ck) * M.getValue(ck, rj);
      Result.setValue(i, j, s);
    end;
end;

il n'y aurait pas comme un petit souci ? Tu risques un debordement vu que tu n'utilises pas les élément de M.
De plus si M a la même taille que self... Tu as inversé les colonnes et les lignes.

**wheel** · 26/04/2022, 21h40

Il y a une implémentation dans la bibliothèque TPMath dans l'unité ugausjor.pas

http://www.unilim.fr/pages_perso/jea...ath/tpmath.htm

**Guesset** · 27/04/2022, 17h02

Bonjour,

Quelques interrogations :

Pourquoi avoir un nombre de colonnes et un nombre de lignes dans une matrice déclarée carrée ? Une seule, Dim par exemple, suffirait.
Dans redim la taille d'une ligne est SetLength(FRows[i], 2 * NbRows); soit deux fois plus large que haute. Pourquoi ? La matrice n'est plus carrée.
Pourquoi mettre un Extended (80 bits) dans un Double (64 bits) dans la fonction ProduitBy et la fonction Inverse ? Ca peut faire mal (ou être inutile comme dans Inverse).

Par ailleurs je me demande quel intérêt il y a à avoirs deux tableaux plutôt qu'un seul à deux dimensions.

L'utilisation d'une propriété Dim appelant éventuellement un GetDim simplifierait les selfs.getNbRows et autres selfs.getNbCols.

Salutations

**BeanzMaster** · 30/04/2022, 08h46

Hello

Envoyé par Guesset

Bonjour,

Quelques interrogations :

Pourquoi avoir un nombre de colonnes et un nombre de lignes dans une matrice déclarée carrée ? Une seule, Dim par exemple, suffirait.

C'est pour la résolution de l'algo. On à besoin d'une deuxième matrice. Mais il pourrait se passer de cette extension de taille en créant une variable contenant cette deuxième matrice.

Envoyé par Guesset

Dans redim la taille d'une ligne est SetLength(FRows[i], 2 * NbRows); soit deux fois plus large que haute. Pourquoi ? La matrice n'est plus carrée.
Pourquoi mettre un Extended (80 bits) dans un Double (64 bits) dans la fonction ProduitBy et la fonction Inverse ? Ca peut faire mal (ou être inutile comme dans Inverse).

L'utilisation d'une propriété Dim appelant éventuellement un GetDim simplifierait les selfs.getNbRows et autres selfs.getNbCols.

Je suis d'accord avec toi sur ces points.

**JP CASSOU** · 04/05/2022, 14h40

Après vérification complète du code, ce n'était pas celui-ci qui était en cause mais le jeu de données de test ... :anvil:

L'algo fonctionne.

Et maintenant, servez-vous

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
 
unit TDUnitMatrixAsVector;
 
{$mode delphi}
 
interface
 
uses
  //StructuresDonnees,
  math,
  Classes, SysUtils;
type
 
{ TTDIntVector }
 
 TTDIntVector = record
  private
    FListe: array of integer;
    function GetInternalIdx(const ARow: integer): integer; inline;
  public
    procedure Redim(const NbRows: integer; const DoFillWithZeros: boolean = false);
    function getNbElements(): integer;
    procedure setFrom(const ARow: integer; const Value: integer);
    function  GetFrom(const ARow: integer): integer;
 
end;
type
 
{ TTDIntVector }
 
 { TTDFloatVector }
 
 TTDFloatVector = record
  private
    FListe: array of double;
    function GetInternalIdx(const ARow: integer): integer; inline;
  public
    procedure Redim(const NbRows: integer; const DoFillWithZeros: boolean = false);
    function  getNbElements(): integer;
    procedure setFrom(const ARow: integer; const Value: double);
    function  GetFrom(const ARow: integer): double;
 
end;
type
 
{ TTDMatrixAsVector }
 
 TTDMatrixAsVector = record
  private
    FNbRows , FNbCols: integer;
    FListeElements: array of double;
    function  GetInternalIdx(const ARow, ACol: integer): integer; inline;
  public
    property  NbRows: integer read FNbRows;
    property  NbCols: integer read FNbCols;
 
    procedure Redim(const NbRows, NbCols: integer; const DoFillWithZeros: boolean = false);
    procedure FillRandom();
    function  getNbElements(): integer;
    procedure SetFrom(const ARow, ACol: integer; const V: double); inline;
    function  GetFrom(const ARow, ACol: integer): double; inline;
 
    function  GetFromVector(const idx: integer): double; inline;
    procedure SwapRows(const R1, R2: integer);
    procedure SwapCols(const C1, C2: integer);
 
    function  Transpose(): TTDMatrixAsVector;
    function  ProduitBy(const B: TTDMatrixAsVector): TTDMatrixAsVector;
    // Inversion de matrice (en place), retourne le déterminant
    function  Inverse(): double;
end;
 
 
 
 
implementation
uses
  Common, DGCDummyUnit;
 
{ TTDFloatVector }
 
function TTDFloatVector.GetInternalIdx(const ARow: integer): integer;
begin
  Result := ARow-1;
end;
 
 
 
procedure TTDFloatVector.Redim(const NbRows: integer; const DoFillWithZeros: boolean);
var
  i: Integer;
begin
  setLength(FListe, NbRows);
  if (DoFillWithZeros) then for i := 0 to High(FListe) do FListe[i] := 0.00;
 
end;
 
function TTDFloatVector.getNbElements(): integer;
begin
  result := length(FListe);
end;
 
procedure TTDFloatVector.setFrom(const ARow: integer; const Value: double);
begin
  FListe[GetInternalIdx(ARow)] := Value;
end;
 
function TTDFloatVector.GetFrom(const ARow: integer): double;
begin
  Result := FListe[GetInternalIdx(ARow)];
end;
 
{ TTDIntVector }
 
function TTDIntVector.GetInternalIdx(const ARow: integer): integer;
begin
  Result := ARow-1;
end;
 
procedure TTDIntVector.Redim(const NbRows: integer; const DoFillWithZeros: boolean);
var
  i: Integer;
begin
  setLength(FListe, NbRows);
  if (DoFillWithZeros) then for i := 0 to High(FListe) do FListe[i] := 0;
end;
 
function TTDIntVector.getNbElements(): integer;
begin
  result := length(FListe);
end;
 
procedure TTDIntVector.setFrom(const ARow: integer; const Value: integer);
begin
  FListe[GetInternalIdx(ARow)] := Value;
end;
 
function TTDIntVector.GetFrom(const ARow: integer): integer;
begin
  Result := FListe[GetInternalIdx(ARow)];
end;
 
{ TTDMatrixAsVector }
 
function TTDMatrixAsVector.GetInternalIdx(const ARow, ACol: integer): integer;
begin
  Result := FNbCols * (ARow-1) + ACol - 1;
end;
 
procedure TTDMatrixAsVector.SwapRows(const R1, R2: integer);
var
  V1, V2: Double;
  j: Integer;
begin
  for j := 1 to FNbCols do
  begin
    V1 := GetFrom(R1, j);
    V2 := GetFrom(R2, j);
    SetFrom(R1, j, V2);
    SetFrom(R2, j, V1);
  end;
end;
 
procedure TTDMatrixAsVector.SwapCols(const C1, C2: integer);
var
  i: Integer;
  V1, V2: double;
begin
  for i := 1 to FNbRows do
  begin
    V1 := GetFrom(i, C1);
    V2 := GetFrom(i, C2);
    SetFrom(i, C1, V2);
    SetFrom(i, C2, V1);
  end;
end;
 
function TTDMatrixAsVector.Transpose(): TTDMatrixAsVector;
var
  i, j: Integer;
begin
  Result.Redim(FNbCols, FNbRows);
  for i := 1 to FNbRows do
    for j := 1 to FNbCols do
      Result.setFrom(j, i, self.GetFrom(i, j));
end;
 
function TTDMatrixAsVector.ProduitBy(const B: TTDMatrixAsVector): TTDMatrixAsVector;
var
  i, j, k: Integer;
  S: double;
begin
  try
    if (B.NbRows <> FNbCols) then raise Exception.CreateFmt('Produit impossible Cols: %d != Rows: %d', [B.NbRows, FNbCols]);
    Result.Redim(FNbRows, B.NbCols);
    for i := 1 to Result.NbRows do
      for j := 1 to Result.NbCols do
      begin
        S := 0.00;
        for k := 1 to FNbCols do S += self.GetFrom(i, k) * B.GetFrom(k, j);
        Result.setFrom(i, j, S);
      end;
 
  except
    on E: Exception do AfficherMessageErreur(E.Message);
  end;
end;
 
procedure TTDMatrixAsVector.Redim(const NbRows, NbCols: integer; const DoFillWithZeros: boolean);
var
  i: Integer;
begin
  FNbRows := NbRows;
  FNbCols := NbCols;
  SetLength(FListeElements, FNbRows * FNbCols);
  if (DoFillWithZeros) then for i := 0 to High(FListeElements) do FListeElements[i] := 0.00;
end;
//******************************************************************************
function TTDMatrixAsVector.Inverse(): double;
var
  Pvt        : Float;       { Pivot }
  Ik, Jk     : Integer;     { Pivot's row and column }
  I, J, K    , Lb: Integer;     { Loop variables }
  T          : Float;       { Temporary variable }
  PRow, PCol : TTDIntVector;  { Stores pivot's row and column }
  MCol       : TTDFloatVector;  { Stores a column of matrix A }
  foo: Double;
begin
  result := 0.00;
  Lb := 1;
  if (FNbCols <> FNbRows) then exit;
  MCol.Redim(FNbRows, True);
  PRow.redim(FNbRows, True);
  PCol.Redim(FNbRows, true);
  Result := 1.0;
  K := Lb;
  while K <= FNbRows do                                                         //  while K <= Ub1 do
  begin
    { Search for largest pivot in submatrix A[K..NbRows, K..NbRows] }
    Pvt := self.GetFrom(K,K);                                                        //    Pvt := A[K,K];
    Ik := K;
    Jk := K;
    for I := K to FNbRows do begin
      for J := K to FNbRows do begin
        if (Abs(self.getFrom(I,J)) > Abs(Pvt)) then begin // if Abs(A[I,J]) > Abs(Pvt) then :
          Pvt := getFrom(I,J);                            //   Pvt := A[I,J];
          Ik := I;                                        //   Ik := I;
          Jk := J;                                        //   Jk := J;
        end;
      end;
    end;
    { Store pivot's position }
    PRow.setFrom(k, Ik);           //  PRow[K] := Ik;
    PCol.setFrom(K, Jk);           //  PCol[K] := Jk;
    { Update determinant }
    Result := Result * Pvt;
    if (Ik <> K) then Result := -Result;
    if (Jk <> K) then Result := -Result;
    { Too weak pivot ==> quasi-singular matrix }
    if Abs(Pvt) < 1E-15 then exit(0.000);
    { Exchange current row (K) with pivot row (Ik) }
    if (ik <> k) then self.SwapRows(Ik, k);                   // if Ik <> K then  for J := Lb to Ub2 do FSwap(A[Ik,J], A[K,J]);
    { Exchange current column (K) with pivot column (Jk) }
    if (Jk <> K) then self.SwapCols(Jk, K);                   //if Jk <> K then                  //  for I := Lb to FNbRows do FSwap(I,Jk,  I,K);
    { Store column K of matrix A into MCol and set this column to zero }
    for I := Lb to NbRows do begin                     // for I := Lb to Ub1 do
      if (I <> K) then  begin                             // if I <> K then:
        MCol.setFrom(i, self.GetFrom(i, k));              //   MCol[I] := A[I,K];
        self.SetFrom(I, K, 0.00);                         //   A[I,K] := 0.0;
      end else begin                                      // else:
        MCol.setFrom(I, 0.00);                            //   MCol[I] := 0.0;
        self.SetFrom(i, k, 1.00);                         //   A[I,K] := 1.0;
      end;
    end;
    { Transform pivot row }                               // { Transform pivot row }
    T := 1.0 / Pvt;                                       // T := 1.0 / Pvt;
    for J := Lb to FNbRows do                             // for J := Lb to Ub2 do  A[K,J] := T * A[K,J];
      self.setFrom(K, J,   T * getFrom(K, J));                   //A[K,J] := T * A[K,J];
 
    { Transform other rows }
    for I := Lb to FNbRows do begin                       // for I := Lb to Ub1 do:
      if (I <> K) then begin                              //   if I <> K then:
        T := MCol.GetFrom(I);                             //     T := MCol[I];
        for J := Lb to FNbRows do begin                   //     for J := Lb to Ub2 do  A[I,J] := A[I,J] - T * A[K,J];
          foo := self.getFrom(I,J);
          foo := foo -  T * self.GetFrom(K,J);
          self.setFrom(i, j,  foo); //A[I,J] := A[I,J] - T * A[K,J];
        end;
      end;
    end;
    K += 1;
  end; //  while K <= FNbRows do
  // Exchange lines of inverse matrix }
  for I := FNbRows downto Lb do begin             // for I := Ub1 downto Lb do:
    Ik := PCol.GetFrom(I);                        //   Ik := PCol[I];
    if (Ik <> I) then self.SwapRows(i, ik);       //   if Ik <> I then:
  end;                                            //     for J := Lb to Ub2 do FSwap(A[I,J], A[Ik,J]);
  { Exchange columns of inverse matrix }
  for J := FNbRows downto Lb do begin             // for J := Ub1 downto Lb do:
    Jk := PRow.GetFrom(J);                        //   Jk := PRow[J];
    if (Jk <> J) then self.SwapCols(J, Jk);       //   if Jk <> J then:
  end;                                            //     for I := Lb to Ub1 do FSwap(A[I,J], A[I,Jk]);
end;
//******************************************************************************
procedure TTDMatrixAsVector.FillRandom();
var
  i: Integer;
begin
  for i := 0 to High(FListeElements) do FListeElements[i] := Random;
end;
 
function TTDMatrixAsVector.getNbElements(): integer;
begin
  result := Length(FListeElements);
end;
 
procedure TTDMatrixAsVector.SetFrom(const ARow, ACol: integer; const V: double);
begin
  FListeElements[GetInternalIdx(ARow, ACol)] := V;
end;
 
function TTDMatrixAsVector.GetFrom(const ARow, ACol: integer): double;
begin
  Result := FListeElements[GetInternalIdx(ARow, ACol)];
end;
 
function TTDMatrixAsVector.GetFromVector(const idx: integer): double;
begin
  Result := FListeElements[Idx];
end;
 
end.

Algo Pivot de Gauss : résultats aberrants [Lazarus]

Lazarus Pascal

Discussions similaires

Partager

Partager