Gauss-Jordan matrix inversion: Difference between revisions

=={{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; Order: Integer);

var

i, j: Integer;

begin

for i := 1 to Order do

begin

for j := 1 to Order do

Write(A[i, j]:8:4, ' ');

WriteLn;

end;

end;

procedure PrintInverse(const A: Matrix; Order: Integer);

var

i, j: Integer;

begin

for i := 1 to Order do

begin

for j := Order + 1 to 2 * Order do

Write(A[i, j]:8:4, ' ');

WriteLn;

end;

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;

procedure GaussJordan(var A: Matrix; Order: Integer);

var

i, j: Integer;

Pivot, Multiplier: Real;

begin

// 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;

end;

procedure DoGaussJordan(B: Matrix);

var

A: Matrix;

Order: Integer;

begin

Order := Length(B);

A := PopulateMatrix(B, Order);

WriteLn('=== Original Matrix ===');

PrintMatrix(A, Order);

GaussJordan(A, Order);

WriteLn('=== Inverse Matrix ===');

PrintInverse(A, Order);

end;

begin

DoGaussJordan(MATRIX_1);

WriteLn;

DoGaussJordan(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>