Gauss-Jordan matrix inversion: Difference between revisions
m
→{{header|11l}}: X.throw
Alextretyak (talk | contribs) m (→{{header|11l}}: X.throw) |
|||
| (6 intermediate revisions by 3 users not shown) | |||
Line 55:
L(j) 0 .< mat.len
I augmat[i][j] != Float(i == j)
X.throw ValueError(‘matrix is singular’)
result[i][j] = augmat[i][mat.len + j]
R result
Line 2,195:
0.0000 -1.0000 2.0000
</pre>
=={{header|EasyLang}}==
{{trans|Go}}
<syntaxhighlight>
proc rref . m[][] .
nrow = len m[][]
ncol = len m[1][]
lead = 1
for r to nrow
if lead > ncol
return
.
i = r
while m[i][lead] = 0
i += 1
if i > nrow
i = r
lead += 1
if lead > ncol
return
.
.
.
swap m[i][] m[r][]
m = m[r][lead]
for k to ncol
m[r][k] /= m
.
for i to nrow
if i <> r
m = m[i][lead]
for k to ncol
m[i][k] -= m * m[r][k]
.
.
.
lead += 1
.
.
proc inverse . mat[][] inv[][] .
inv[][] = [ ]
ln = len mat[][]
for i to ln
if len mat[i][] <> ln
# not a square matrix
return
.
aug[][] &= [ ]
len aug[i][] 2 * ln
for j to ln
aug[i][j] = mat[i][j]
.
aug[i][ln + i] = 1
.
rref aug[][]
for i to ln
inv[][] &= [ ]
for j = ln + 1 to 2 * ln
inv[i][] &= aug[i][j]
.
.
.
test[][] = [ [ 1 2 3 ] [ 4 1 6 ] [ 7 8 9 ] ]
inverse test[][] inv[][]
print inv[][]
</syntaxhighlight>
=={{header|Factor}}==
Line 2,235 ⟶ 2,301:
-13/23 16/69 -2/69 13/23
</pre>
=={{header|Fortran}}==
Note that the Crout algorithm is a more efficient way to invert a matrix.
Line 4,240 ⟶ 4,307:
-0.69565 0.36232 0.07971 0.19565
-0.56522 0.23188 -0.02899 0.56522</pre>
=={{header|Pascal}}==
==={{header|Free Pascal}}===
<syntaxhighlight lang="pascal">
program MatrixInverse;
{$mode ObjFPC}{$H+}
const
MATRIX_1: array of array of Real = ((1, 2, 3),
(4, 1, 6),
(7, 8, 9));
MATRIX_2: array of array of Real = ((3.0, 1.0, 5.0, 9.0, 7.0),
(8.0, 2.0, 8.0, 0.0, 1.0),
(1.0, 7.0, 2.0, 0.0, 3.0),
(0.0, 1.0, 7.0, 0.0, 9.0),
(3.0, 5.0, 6.0, 1.0, 1.0));
type
Matrix = array of array of Real;
function PopulateMatrix(A: Matrix; Order: Integer): Matrix;
var
i, j: Integer;
begin
SetLength(Result, Succ(Order), Succ(Order * 2));
for i := 0 to Pred(Order) do
for j := 0 to Pred(Order) do
Result[i + 1, j + 1] := A[i, j];
end;
procedure PrintMatrix(const A: Matrix);
var
i, j, Order: Integer;
begin
Order := Length(A);
for i := 0 to Pred(Order) do
begin
for j := 0 to Pred(Order) do
Write(A[i, j]:8:4);
WriteLn;
end;
WriteLn;
end;
procedure InterchangeRows(var A: Matrix; Order: Integer; Row1, Row2: Integer);
var
j: Integer;
temp: Real;
begin
for j := 1 to 2 * Order do
begin
temp := A[Row1, j];
A[Row1, j] := A[Row2, j];
A[Row2, j] := temp;
end;
end;
procedure DivideRow(var A: Matrix; Order: Integer; Row: Integer; Divisor: Real);
var
j: Integer;
begin
for j := 1 to 2 * Order do
A[Row, j] := A[Row, j] / Divisor;
end;
procedure SubtractRows(var A: Matrix; Order: Integer; Row1, Row2: Integer; Multiplier: Real);
var
j: Integer;
begin
for j := 1 to 2 * Order do
A[Row1, j] := A[Row1, j] - Multiplier * A[Row2, j];
end;
function GaussJordan(B: Matrix): Matrix;
var
i, j, Order: Integer;
Pivot, Multiplier: Real;
A: Matrix;
begin
Order := Length(B);
SetLength(Result, Order, Order);
A := PopulateMatrix(B, Order);
// Create the augmented matrix
for i := 1 to Order do
for j := Order + 1 to 2 * Order do
if j = (i + Order) then
A[i, j] := 1;
// Interchange rows if needed
for i := Order downto 2 do
if A[i - 1, 1] < A[i, 1] then
InterchangeRows(A, Order, i - 1, i);
// Perform Gauss-Jordan elimination
for i := 1 to Order do
begin
Pivot := A[i, i];
DivideRow(A, Order, i, Pivot);
for j := 1 to Order do
if j <> i then
begin
Multiplier := A[j, i];
SubtractRows(A, Order, j, i, Multiplier);
end;
end;
for i := 0 to Pred(Order) do
for j := 0 to Pred(Order) do
Result[i, j] := A[Succ(i), Succ(j) + Order];
end;
begin
writeln('== Original Matrix ==');
PrintMatrix(MATRIX_1);
writeln('== Inverse Matrix ==');
PrintMatrix(GaussJordan(MATRIX_1));
writeln('== Original Matrix ==');
PrintMatrix(MATRIX_2);
writeln('== Inverse Matrix ==');
PrintMatrix(GaussJordan(MATRIX_2));
end.
</syntaxhighlight>
{{out}}
<pre>
== Original Matrix ==
1.0000 2.0000 3.0000
4.0000 1.0000 6.0000
7.0000 8.0000 9.0000
== Inverse Matrix ==
-0.8125 0.1250 0.1875
0.1250 -0.2500 0.1250
0.5208 0.1250 -0.1458
== Original Matrix ==
3.0000 1.0000 5.0000 9.0000 7.0000
8.0000 2.0000 8.0000 0.0000 1.0000
1.0000 7.0000 2.0000 0.0000 3.0000
0.0000 1.0000 7.0000 0.0000 9.0000
3.0000 5.0000 6.0000 1.0000 1.0000
== Inverse Matrix ==
0.0286 0.1943 0.1330 -0.0596 -0.2578
-0.0058 -0.0326 0.1211 -0.0380 0.0523
-0.0302 -0.0682 -0.1790 0.0605 0.2719
0.1002 -0.0673 -0.0562 -0.0626 0.0981
0.0241 0.0567 0.1258 0.0682 -0.2173
</pre>
=={{header|Perl}}==
Line 5,476 ⟶ 5,693:
{{libheader|Wren-matrix}}
The above module does in fact use the Gauss-Jordan method for matrix inversion.
<syntaxhighlight lang="
import "./fmt" for Fmt
var arrays = [
| |||