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

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

Lazarus Pascal

Vue hybride

Discussions similaires

Partager

Partager