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
|
{ ******************************************************************
Solution of a system of linear equations by Gauss-Jordan method
****************************************************************** }
unit ugausjor;
interface
uses
utypes, uminmax;
procedure GaussJordan(A : PMatrix;
Lb, Ub1, Ub2 : Integer;
var Det : Float);
{ ------------------------------------------------------------------
Transforms a matrix according to the Gauss-Jordan method
------------------------------------------------------------------
Input parameters : A = system matrix
Lb = lower matrix bound in both dim.
Ub1, Ub2 = upper matrix bounds
------------------------------------------------------------------
Output parameters: A = transformed matrix
Det = determinant of A
------------------------------------------------------------------
Possible results : MatOk : No error
MatErrDim : Non-compatible dimensions
MatSing : Singular matrix
------------------------------------------------------------------ }
implementation
procedure GaussJordan(A : PMatrix;
Lb, Ub1, Ub2 : Integer;
var Det : Float);
var
Pvt : Float; { Pivot }
Ik, Jk : Integer; { Pivot's row and column }
I, J, K : Integer; { Loop variables }
T : Float; { Temporary variable }
PRow, PCol : PIntVector; { Stores pivot's row and column }
MCol : PVector; { Stores a column of matrix A }
procedure Terminate(ErrCode : Integer);
{ Set error code and deallocate arrays }
begin
DelIntVector(PRow, Ub1);
DelIntVector(PCol, Ub1);
DelVector(MCol, Ub1);
SetErrCode(ErrCode);
end;
begin
if Ub1 > Ub2 then
begin
SetErrCode(MatErrDim);
Exit
end;
DimIntVector(PRow, Ub1);
DimIntVector(PCol, Ub1);
DimVector(MCol, Ub1);
Det := 1.0;
K := Lb;
while K <= Ub1 do
begin
{ Search for largest pivot in submatrix A[K..Ub1, K..Ub1] }
Pvt := A^[K]^[K];
Ik := K;
Jk := K;
for I := K to Ub1 do
for J := K to Ub1 do
if Abs(A^[I]^[J]) > Abs(Pvt) then
begin
Pvt := A^[I]^[J];
Ik := I;
Jk := J;
end;
{ Store pivot's position }
PRow^[K] := Ik;
PCol^[K] := Jk;
{ Update determinant }
Det := Det * Pvt;
if Ik <> K then Det := - Det;
if Jk <> K then Det := - Det;
{ Too weak pivot ==> quasi-singular matrix }
if Abs(Pvt) < MachEp then
begin
Terminate(MatSing);
Exit
end;
{ Exchange current row (K) with pivot row (Ik) }
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
for I := Lb to Ub1 do
FSwap(A^[I]^[Jk], A^[I]^[K]);
{ Store column K of matrix A into MCol
and set this column to zero }
for I := Lb to Ub1 do
if I <> K then
begin
MCol^[I] := A^[I]^[K];
A^[I]^[K] := 0.0;
end
else
begin
MCol^[I] := 0.0;
A^[I]^[K] := 1.0;
end;
{ Transform pivot row }
T := 1.0 / Pvt;
for J := Lb to Ub2 do
A^[K]^[J] := T * A^[K]^[J];
{ Transform other rows }
for I := Lb to Ub1 do
if I <> K then
begin
T := MCol^[I];
for J := Lb to Ub2 do
A^[I]^[J] := A^[I]^[J] - T * A^[K]^[J];
end;
Inc(K);
end;
{ Exchange lines of inverse matrix }
for I := Ub1 downto Lb do
begin
Ik := PCol^[I];
if Ik <> I then
for J := Lb to Ub2 do
FSwap(A^[I]^[J], A^[Ik]^[J]);
end;
{ Exchange columns of inverse matrix }
for J := Ub1 downto Lb do
begin
Jk := PRow^[J];
if Jk <> J then
for I := Lb to Ub1 do
FSwap(A^[I]^[J], A^[I]^[Jk]);
end;
Terminate(MatOk);
end;
end.
|