Gauss-Jordan matrix inversion: Difference between revisions
m
→{{header|11l}}: X.throw
(→{{header|jq}}: fix bug that did not affect output (add: .lead == $nc - 1)) |
Alextretyak (talk | contribs) m (→{{header|11l}}: X.throw) |
||
| (17 intermediate revisions by 8 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 508:
1.000 0.000 -2.000 4.000
4.000 1.000 2.000 2.000
</pre>
=={{header|ATS}}==
Not all the parts of the "little matrix library" included in the code are actually used, although some of those not used were (in fact) used during debugging.
<syntaxhighlight lang="ats">
(* Set to 1 for debugging: to fill in ones and zeros and other values
that are not actually used but are part of the theory of the
Gauss-Jordan algorithm. *)
#define DO_THINGS_THAT_DO_NOT_NEED_TO_BE_DONE 0
(* Setting this to 1 may cause rounding to change, and the change in
rounding is not unlikely to cause detection of singularity of a
matrix to change. (To invert a matrix that might be nearly
singular, the SVD seems a popular method.) The
-fexpensive-optimizations option to GCC also may cause the same
rounding changes (due to fused-multiply-and-add instructions being
generated). *)
#define USE_MULTIPLY_AND_ADD 0
%{^
#include <math.h>
#include <float.h>
%}
#include "share/atspre_staload.hats"
macdef NAN = g0f2f ($extval (float, "NAN"))
macdef Zero = g0i2f 0
macdef One = g0i2f 1
macdef Two = g0i2f 2
#if USE_MULTIPLY_AND_ADD #then
(* "fma" from the C math library, although your system may have it as
a built-in. *)
extern fn {tk : tkind} g0float_fma : (g0float tk, g0float tk, g0float tk) -<> g0float tk
implement g0float_fma<fltknd> (x, y, z) = $extfcall (float, "fmaf", x, y, z)
implement g0float_fma<dblknd> (x, y, z) = $extfcall (double, "fma", x, y, z)
implement g0float_fma<ldblknd> (x, y, z) = $extfcall (ldouble, "fmal", x, y, z)
overload fma with g0float_fma
macdef multiply_and_add = fma
#else
macdef multiply_and_add (x, y, z) = (,(x) * ,(y)) + ,(z)
#endif
(*------------------------------------------------------------------*)
(* A "little matrix library" *)
typedef Matrix_Index_Map (m1 : int, n1 : int, m0 : int, n0 : int) =
{i1, j1 : pos | i1 <= m1; j1 <= n1}
(int i1, int j1) -<cloref0>
[i0, j0 : pos | i0 <= m0; j0 <= n0]
@(int i0, int j0)
datatype Real_Matrix (tk : tkind,
m1 : int, n1 : int,
m0 : int, n0 : int) =
| Real_Matrix of (matrixref (g0float tk, m0, n0),
int m1, int n1, int m0, int n0,
Matrix_Index_Map (m1, n1, m0, n0))
typedef Real_Matrix (tk : tkind, m1 : int, n1 : int) =
[m0, n0 : pos] Real_Matrix (tk, m1, n1, m0, n0)
typedef Real_Vector (tk : tkind, m1 : int, n1 : int) =
[m1 == 1 || n1 == 1] Real_Matrix (tk, m1, n1)
typedef Real_Row (tk : tkind, n1 : int) = Real_Vector (tk, 1, n1)
typedef Real_Column (tk : tkind, m1 : int) = Real_Vector (tk, m1, 1)
extern fn {tk : tkind}
Real_Matrix_make_elt :
{m0, n0 : pos}
(int m0, int n0, g0float tk) -< !wrt >
Real_Matrix (tk, m0, n0, m0, n0)
extern fn {tk : tkind}
Real_Matrix_copy :
{m1, n1 : pos}
Real_Matrix (tk, m1, n1) -< !refwrt > Real_Matrix (tk, m1, n1)
extern fn {tk : tkind}
Real_Matrix_copy_to :
{m1, n1 : pos}
(Real_Matrix (tk, m1, n1), (* destination *)
Real_Matrix (tk, m1, n1)) -< !refwrt >
void
extern fn {tk : tkind}
Real_Matrix_fill_with_elt :
{m1, n1 : pos}
(Real_Matrix (tk, m1, n1), g0float tk) -< !refwrt > void
extern fn {}
Real_Matrix_dimension :
{tk : tkind}
{m1, n1 : pos}
Real_Matrix (tk, m1, n1) -<> @(int m1, int n1)
extern fn {tk : tkind}
Real_Matrix_get_at :
{m1, n1 : pos}
{i1, j1 : pos | i1 <= m1; j1 <= n1}
(Real_Matrix (tk, m1, n1), int i1, int j1) -< !ref > g0float tk
extern fn {tk : tkind}
Real_Matrix_set_at :
{m1, n1 : pos}
{i1, j1 : pos | i1 <= m1; j1 <= n1}
(Real_Matrix (tk, m1, n1), int i1, int j1, g0float tk) -< !refwrt >
void
extern fn {}
Real_Matrix_apply_index_map :
{tk : tkind}
{m1, n1 : pos}
{m0, n0 : pos}
(Real_Matrix (tk, m0, n0), int m1, int n1,
Matrix_Index_Map (m1, n1, m0, n0)) -<>
Real_Matrix (tk, m1, n1)
extern fn {}
Real_Matrix_transpose :
(* This is transposed INDEXING. It does NOT copy the data. *)
{tk : tkind}
{m1, n1 : pos}
{m0, n0 : pos}
Real_Matrix (tk, m1, n1, m0, n0) -<>
Real_Matrix (tk, n1, m1, m0, n0)
extern fn {}
Real_Matrix_block :
(* This is block (submatrix) INDEXING. It does NOT copy the data. *)
{tk : tkind}
{p0, p1 : pos | p0 <= p1}
{q0, q1 : pos | q0 <= q1}
{m1, n1 : pos | p1 <= m1; q1 <= n1}
{m0, n0 : pos}
(Real_Matrix (tk, m1, n1, m0, n0),
int p0, int p1, int q0, int q1) -<>
Real_Matrix (tk, p1 - p0 + 1, q1 - q0 + 1, m0, n0)
extern fn {tk : tkind}
Real_Matrix_unit_matrix :
{m : pos}
int m -< !refwrt > Real_Matrix (tk, m, m)
extern fn {tk : tkind}
Real_Matrix_unit_matrix_to :
{m : pos}
Real_Matrix (tk, m, m) -< !refwrt > void
extern fn {tk : tkind}
Real_Matrix_matrix_sum :
{m, n : pos}
(Real_Matrix (tk, m, n), Real_Matrix (tk, m, n)) -< !refwrt >
Real_Matrix (tk, m, n)
extern fn {tk : tkind}
Real_Matrix_matrix_sum_to :
{m, n : pos}
(Real_Matrix (tk, m, n), (* destination*)
Real_Matrix (tk, m, n),
Real_Matrix (tk, m, n)) -< !refwrt >
void
extern fn {tk : tkind}
Real_Matrix_matrix_difference :
{m, n : pos}
(Real_Matrix (tk, m, n), Real_Matrix (tk, m, n)) -< !refwrt >
Real_Matrix (tk, m, n)
extern fn {tk : tkind}
Real_Matrix_matrix_difference_to :
{m, n : pos}
(Real_Matrix (tk, m, n), (* destination*)
Real_Matrix (tk, m, n),
Real_Matrix (tk, m, n)) -< !refwrt >
void
extern fn {tk : tkind}
Real_Matrix_matrix_product :
{m, n, p : pos}
(Real_Matrix (tk, m, n), Real_Matrix (tk, n, p)) -< !refwrt >
Real_Matrix (tk, m, p)
extern fn {tk : tkind}
Real_Matrix_matrix_product_to :
{m, n, p : pos}
(Real_Matrix (tk, m, p), (* destination*)
Real_Matrix (tk, m, n),
Real_Matrix (tk, n, p)) -< !refwrt >
void
extern fn {tk : tkind}
Real_Matrix_scalar_product :
{m, n : pos}
(Real_Matrix (tk, m, n), g0float tk) -< !refwrt >
Real_Matrix (tk, m, n)
extern fn {tk : tkind}
Real_Matrix_scalar_product_2 :
{m, n : pos}
(g0float tk, Real_Matrix (tk, m, n)) -< !refwrt >
Real_Matrix (tk, m, n)
extern fn {tk : tkind}
Real_Matrix_scalar_product_to :
{m, n : pos}
(Real_Matrix (tk, m, n), (* destination*)
Real_Matrix (tk, m, n), g0float tk) -< !refwrt > void
extern fn {tk : tkind} (* Useful for debugging. *)
Real_Matrix_fprint :
{m, n : pos}
(FILEref, Real_Matrix (tk, m, n)) -<1> void
overload copy with Real_Matrix_copy
overload copy_to with Real_Matrix_copy_to
overload fill_with_elt with Real_Matrix_fill_with_elt
overload dimension with Real_Matrix_dimension
overload [] with Real_Matrix_get_at
overload [] with Real_Matrix_set_at
overload apply_index_map with Real_Matrix_apply_index_map
overload transpose with Real_Matrix_transpose
overload block with Real_Matrix_block
overload unit_matrix with Real_Matrix_unit_matrix
overload unit_matrix_to with Real_Matrix_unit_matrix_to
overload matrix_sum with Real_Matrix_matrix_sum
overload matrix_sum_to with Real_Matrix_matrix_sum_to
overload matrix_difference with Real_Matrix_matrix_difference
overload matrix_difference_to with Real_Matrix_matrix_difference_to
overload matrix_product with Real_Matrix_matrix_product
overload matrix_product_to with Real_Matrix_matrix_product_to
overload scalar_product with Real_Matrix_scalar_product
overload scalar_product with Real_Matrix_scalar_product_2
overload scalar_product_to with Real_Matrix_scalar_product_to
overload + with matrix_sum
overload - with matrix_difference
overload * with matrix_product
overload * with scalar_product
(*------------------------------------------------------------------*)
(* Implementation of the "little matrix library" *)
implement {tk}
Real_Matrix_make_elt (m0, n0, elt) =
Real_Matrix (matrixref_make_elt<g0float tk> (i2sz m0, i2sz n0, elt),
m0, n0, m0, n0, lam (i1, j1) => @(i1, j1))
implement {}
Real_Matrix_dimension A =
case+ A of Real_Matrix (_, m1, n1, _, _, _) => @(m1, n1)
implement {tk}
Real_Matrix_get_at (A, i1, j1) =
let
val+ Real_Matrix (storage, _, _, _, n0, index_map) = A
val @(i0, j0) = index_map (i1, j1)
in
matrixref_get_at<g0float tk> (storage, pred i0, n0, pred j0)
end
implement {tk}
Real_Matrix_set_at (A, i1, j1, x) =
let
val+ Real_Matrix (storage, _, _, _, n0, index_map) = A
val @(i0, j0) = index_map (i1, j1)
in
matrixref_set_at<g0float tk> (storage, pred i0, n0, pred j0, x)
end
implement {}
Real_Matrix_apply_index_map (A, m1, n1, index_map) =
(* This is not the most efficient way to acquire new indexing, but
it will work. It requires three closures, instead of the two
needed by our implementations of "transpose" and "block". *)
let
val+ Real_Matrix (storage, m1a, n1a, m0, n0, index_map_1a) = A
in
Real_Matrix (storage, m1, n1, m0, n0,
lam (i1, j1) =>
index_map_1a (i1a, j1a) where
{ val @(i1a, j1a) = index_map (i1, j1) })
end
implement {}
Real_Matrix_transpose A =
let
val+ Real_Matrix (storage, m1, n1, m0, n0, index_map) = A
in
Real_Matrix (storage, n1, m1, m0, n0,
lam (i1, j1) => index_map (j1, i1))
end
implement {}
Real_Matrix_block (A, p0, p1, q0, q1) =
let
val+ Real_Matrix (storage, m1, n1, m0, n0, index_map) = A
in
Real_Matrix (storage, succ (p1 - p0), succ (q1 - q0), m0, n0,
lam (i1, j1) =>
index_map (p0 + pred i1, q0 + pred j1))
end
implement {tk}
Real_Matrix_copy A =
let
val @(m1, n1) = dimension A
val C = Real_Matrix_make_elt<tk> (m1, n1, A[1, 1])
val () = copy_to<tk> (C, A)
in
C
end
implement {tk}
Real_Matrix_copy_to (Dst, Src) =
let
val @(m1, n1) = dimension Src
prval [m1 : int] EQINT () = eqint_make_gint m1
prval [n1 : int] EQINT () = eqint_make_gint n1
var i : intGte 1
in
for* {i : pos | i <= m1 + 1} .<(m1 + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m1; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= n1 + 1} .<(n1 + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n1; j := succ j)
Dst[i, j] := Src[i, j]
end
end
implement {tk}
Real_Matrix_fill_with_elt (A, elt) =
let
val @(m1, n1) = dimension A
prval [m1 : int] EQINT () = eqint_make_gint m1
prval [n1 : int] EQINT () = eqint_make_gint n1
var i : intGte 1
in
for* {i : pos | i <= m1 + 1} .<(m1 + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m1; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= n1 + 1} .<(n1 + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n1; j := succ j)
A[i, j] := elt
end
end
implement {tk}
Real_Matrix_unit_matrix {m} m =
let
val A = Real_Matrix_make_elt<tk> (m, m, Zero)
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
A[i, i] := One;
A
end
implement {tk}
Real_Matrix_unit_matrix_to A =
let
val @(m, _) = dimension A
prval [m : int] EQINT () = eqint_make_gint m
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= m + 1} .<(m + 1) - j>.
(j : int j) =>
(j := 1; j <> succ m; j := succ j)
A[i, j] := (if i = j then One else Zero)
end
end
implement {tk}
Real_Matrix_matrix_sum (A, B) =
let
val @(m, n) = dimension A
val C = Real_Matrix_make_elt<tk> (m, n, NAN)
val () = matrix_sum_to<tk> (C, A, B)
in
C
end
implement {tk}
Real_Matrix_matrix_sum_to (C, A, B) =
let
val @(m, n) = dimension A
prval [m : int] EQINT () = eqint_make_gint m
prval [n : int] EQINT () = eqint_make_gint n
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
C[i, j] := A[i, j] + B[i, j]
end
end
implement {tk}
Real_Matrix_matrix_difference (A, B) =
let
val @(m, n) = dimension A
val C = Real_Matrix_make_elt<tk> (m, n, NAN)
val () = matrix_difference_to<tk> (C, A, B)
in
C
end
implement {tk}
Real_Matrix_matrix_difference_to (C, A, B) =
let
val @(m, n) = dimension A
prval [m : int] EQINT () = eqint_make_gint m
prval [n : int] EQINT () = eqint_make_gint n
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
C[i, j] := A[i, j] - B[i, j]
end
end
implement {tk}
Real_Matrix_matrix_product (A, B) =
let
val @(m, n) = dimension A and @(_, p) = dimension B
val C = Real_Matrix_make_elt<tk> (m, p, NAN)
val () = matrix_product_to<tk> (C, A, B)
in
C
end
implement {tk}
Real_Matrix_matrix_product_to (C, A, B) =
let
val @(m, n) = dimension A and @(_, p) = dimension B
prval [m : int] EQINT () = eqint_make_gint m
prval [n : int] EQINT () = eqint_make_gint n
prval [p : int] EQINT () = eqint_make_gint p
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
let
var k : intGte 1
in
for* {k : pos | k <= p + 1} .<(p + 1) - k>.
(k : int k) =>
(k := 1; k <> succ p; k := succ k)
let
var j : intGte 1
in
C[i, k] := A[i, 1] * B[1, k];
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 2; j <> succ n; j := succ j)
C[i, k] :=
multiply_and_add (A[i, j], B[j, k], C[i, k])
end
end
end
implement {tk}
Real_Matrix_scalar_product (A, r) =
let
val @(m, n) = dimension A
val C = Real_Matrix_make_elt<tk> (m, n, NAN)
val () = scalar_product_to<tk> (C, A, r)
in
C
end
implement {tk}
Real_Matrix_scalar_product_2 (r, A) =
Real_Matrix_scalar_product<tk> (A, r)
implement {tk}
Real_Matrix_scalar_product_to (C, A, r) =
let
val @(m, n) = dimension A
prval [m : int] EQINT () = eqint_make_gint m
prval [n : int] EQINT () = eqint_make_gint n
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
C[i, j] := A[i, j] * r
end
end
implement {tk}
Real_Matrix_fprint {m, n} (outf, A) =
let
val @(m, n) = dimension A
var i : intGte 1
in
for* {i : pos | i <= m + 1} .<(m + 1) - i>.
(i : int i) =>
(i := 1; i <> succ m; i := succ i)
let
var j : intGte 1
in
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
let
typedef FILEstar = $extype"FILE *"
extern castfn FILEref2star : FILEref -<> FILEstar
val _ = $extfcall (int, "fprintf", FILEref2star outf,
"%16.6g", A[i, j])
in
end;
fprintln! (outf)
end
end
(*------------------------------------------------------------------*)
(* Gauss-Jordan inversion *)
extern fn {tk : tkind}
Real_Matrix_inverse_by_gauss_jordan :
{n : pos}
Real_Matrix (tk, n, n) -< !refwrt > Option (Real_Matrix (tk, n, n))
#if DO_THINGS_THAT_DO_NOT_NEED_TO_BE_DONE #then
macdef do_needless_things = true
#else
macdef do_needless_things = false
#endif
fn {tk : tkind}
needlessly_set_to_value
{n : pos}
{i, j : pos | i <= n; j <= n}
(A : Real_Matrix (tk, n, n),
i : int i,
j : int j,
x : g0float tk) :<!refwrt> void =
if do_needless_things then A[i, j] := x
implement {tk}
Real_Matrix_inverse_by_gauss_jordan {n} A =
let
val @(n, _) = dimension A
typedef one_to_n = intBtwe (1, n)
(* Partial pivoting. *)
implement
array_tabulate$fopr<one_to_n> i =
let
val i = g1ofg0 (sz2i (succ i))
val () = assertloc ((1 <= i) * (i <= n))
in
i
end
val rows_permutation =
$effmask_all arrayref_tabulate<one_to_n> (i2sz n)
fn
index_map : Matrix_Index_Map (n, n, n, n) =
lam (i1, j1) => $effmask_ref
(@(i0, j1) where { val i0 = rows_permutation[i1 - 1] })
val A = apply_index_map (copy<tk> A, n, n, index_map)
and B = apply_index_map (unit_matrix<tk> n, n, n, index_map)
fn {}
exchange_rows (i1 : one_to_n,
i2 : one_to_n) :<!refwrt> void =
if i1 <> i2 then
let
val k1 = rows_permutation[pred i1]
and k2 = rows_permutation[pred i2]
in
rows_permutation[pred i1] := k2;
rows_permutation[pred i2] := k1
end
fn {}
normalize_pivot_row (j : one_to_n) :<!refwrt> void =
let
prval [j : int] EQINT () = eqint_make_gint j
val pivot_val = A[j, j]
var k : intGte 1
in
needlessly_set_to_value (A, j, j, One);
for* {k : int | j + 1 <= k; k <= n + 1} .<(n + 1) - k>.
(k : int k) =>
(k := succ j; k <> succ n; k := succ k)
A[j, k] := A[j, k] / pivot_val;
for* {k : int | 1 <= k; k <= n + 1} .<(n + 1) - k>.
(k : int k) =>
(k := 1; k <> succ n; k := succ k)
B[j, k] := B[j, k] / pivot_val;
end
fn
subtract_normalized_pivot_row
(i : one_to_n, j : one_to_n) :<!refwrt> void =
let
prval [j : int] EQINT () = eqint_make_gint j
val lead_val = A[i, j]
in
if lead_val <> Zero then
let
val factor = ~lead_val
var k : intGte 1
in
needlessly_set_to_value (A, i, j, Zero);
for* {k : int | j + 1 <= k; k <= n + 1} .<(n + 1) - k>.
(k : int k) =>
(k := succ j; k <> succ n; k := succ k)
A[i, k] :=
multiply_and_add (A[j, k], factor, A[i, k]);
for* {k : int | 1 <= k; k <= n + 1} .<(n + 1) - k>.
(k : int k) =>
(k := 1; k <> succ n; k := succ k)
B[i, k] :=
multiply_and_add (B[j, k], factor, B[i, k])
end
end
fun
main_loop {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j,
success : &bool? >> bool) :<!refwrt> void =
if j = succ n then
success := true
else
let
fun
select_pivot {i : int | j <= i; i <= n + 1}
.<(n + 1) - i>.
(i : int i,
max_abs : g0float tk,
i_max_abs : intBtwe (j - 1, n))
:<!ref> intBtwe (j - 1, n) =
if i = succ n then
i_max_abs
else
let
val abs_aij = abs A[i, j]
in
if abs_aij > max_abs then
select_pivot (succ i, abs_aij, i)
else
select_pivot (succ i, max_abs, i_max_abs)
end
val i_pivot = select_pivot (j, Zero, pred j)
prval [i_pivot : int] EQINT () = eqint_make_gint i_pivot
in
if i_pivot = pred j then
success := false
else
let
var i : intGte 1
in
exchange_rows (i_pivot, j);
normalize_pivot_row (j);
for* {i : int | 1 <= i; i <= j}
.<j - i>.
(i : int i) =>
(i := 1; i <> j; i := succ i)
subtract_normalized_pivot_row (i, j);
for* {i : int | j + 1 <= i; i <= n + 1}
.<(n + 1) - i>.
(i : int i) =>
(i := succ j; i <> succ n; i := succ i)
subtract_normalized_pivot_row (i, j);
main_loop (succ j, success)
end
end
var success : bool
val () = main_loop (1, success)
in
(* The returned B will "contain" the rows_permutation array and
some extra closures. If you want a "clean" B, use "val B = copy
B" or some such. *)
if success then Some B else None ()
end
overload inverse_by_gauss_jordan with
Real_Matrix_inverse_by_gauss_jordan
(*------------------------------------------------------------------*)
fn {tk : tkind}
fprint_matrix_and_its_inverse
{n : pos}
(outf : FILEref,
A : Real_Matrix (tk, n, n)) : void =
let
typedef FILEstar = $extype"FILE *"
extern castfn FILEref2star : FILEref -<> FILEstar
macdef fmt = "%8.4f"
macdef left_bracket = "["
macdef right_bracket = " ]"
macdef product = " ✕ "
macdef equals = " = "
macdef spacing = " "
macdef msg_for_singular = " appears to be singular"
macdef print_num (x) =
ignoret ($extfcall (int, "fprintf", FILEref2star outf,
fmt, ,(x)))
val @(n, _) = dimension A
in
case+ inverse_by_gauss_jordan<tk> (A) of
| None () =>
let
var i : intGte 1
in
for* {i : pos | i <= n + 1} .<(n + 1) - i>.
(i : int i) =>
(i := 1; i <> succ n; i := succ i)
let
var j : intGte 1
in
fprint! (outf, left_bracket);
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
print_num A[i, j];
fprint! (outf, right_bracket);
if pred i = half n then
fprint! (outf, msg_for_singular);
fprintln! (outf)
end
end
| Some B =>
let
val AB = A * B
var i : intGte 1
in
for* {i : pos | i <= n + 1} .<(n + 1) - i>.
(i : int i) =>
(i := 1; i <> succ n; i := succ i)
let
var j : intGte 1
in
fprint! (outf, left_bracket);
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
print_num A[i, j];
fprint! (outf, right_bracket);
if pred i = half n then
fprint! (outf, product)
else
fprint! (outf, spacing);
fprint! (outf, left_bracket);
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
print_num B[i, j];
fprint! (outf, right_bracket);
if pred i = half n then
fprint! (outf, equals)
else
fprint! (outf, spacing);
fprint! (outf, left_bracket);
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
print_num AB[i, j];
fprint! (outf, right_bracket);
fprintln! (outf)
end
end
end
fn {tk : tkind} (* We like Fortran, so COLUMN major here. *)
column_major_list_to_square_matrix
{n : pos}
(n : int n,
lst : list (g0float tk, n * n))
: Real_Matrix (tk, n, n) =
let
#define :: list_cons
prval () = mul_gte_gte_gte {n, n} ()
val A = Real_Matrix_make_elt (n, n, NAN)
val lstref : ref (List0 (g0float tk)) = ref lst
var j : intGte 1
in
for* {j : pos | j <= n + 1} .<(n + 1) - j>.
(j : int j) =>
(j := 1; j <> succ n; j := succ j)
let
var i : intGte 1
in
for* {i : pos | i <= n + 1} .<(n + 1) - i>.
(i : int i) =>
(i := 1; i <> succ n; i := succ i)
case- !lstref of
| hd :: tl =>
begin
A[i, j] := hd;
!lstref := tl
end
end;
A
end
macdef separator = "\n"
fn
print_example
{n : pos}
(n : int n,
lst : list (double, n * n)) : void =
let
val A = column_major_list_to_square_matrix (n, lst)
in
print! separator;
fprint_matrix_and_its_inverse<dblknd> (stdout_ref, A)
end
implement
main0 () =
begin
println! ("\n(The examples are printed here after ",
"rounding to 4 decimal places.)");
print_example (1, $list (0.0));
print_example (1, $list (4.0));
print_example (2, $list (1.0, 1.0, 1.0, 1.0));
print_example (2, $list (2.0, 0.0, 0.0, 2.0));
print_example (2, $list (1.0, 2.0, 3.0, 4.0));
println! ("\nAn orthonormal matrix of rank 5 ",
"(its inverse will equal its transpose):");
print_example
(5, $list (~0.08006407690254358, ~0.1884825001438613,
~0.8954592676839166, 0.2715238629598956,
~0.2872134789517761, ~0.5604485383178048,
~0.4587939881550582, ~0.01946650581921516,
0.1609525690123523, 0.6701647842208114,
~0.2401922307076307, 0.1517054269450593,
~0.3178281430868086, ~0.8979459113320678,
0.1094146586482966, ~0.3202563076101743,
~0.6104994151001173, 0.2983616372675922,
~0.1997416922923375, ~0.6291342872277008,
~0.7205766921228921, 0.5985468663105067,
0.08797363206760908, 0.2327347396799102,
~0.2461829819586655));
(* The following matrix may or may not be detected as singular,
depending on how rounding is done. *)
(*
print_example
(5, $list (~0.08006407690254358, ~0.1884825001438613,
~0.8954592676839166, 0.2715238629598956,
~0.2872134789517761, ~0.5604485383178048,
~0.4587939881550582, ~0.01946650581921516,
0.1609525690123523, 0.6701647842208114,
~0.2401922307076307, 0.1517054269450593,
~0.3178281430868086, ~0.8979459113320678,
0.1094146586482966, ~0.3202563076101743,
~0.6104994151001173, 0.2983616372675922,
~0.1997416922923375, ~0.6291342872277008,
5.0 * ~0.08006407690254358,
5.0 * ~0.1884825001438613,
5.0 * ~0.8954592676839166,
5.0 * 0.2715238629598956,
5.0 * ~0.2872134789517761));
*)
println! ("\nA 5✕5 singular matrix:");
print_example
(5, $list (1.0, 1.0, 2.0, 2.0, 2.0,
5.0, 1.0, 2.0, 2.0, 2.0,
3.0, 2.0, 3.0, 5.0, 7.0,
1.0, 1.0, 2.0, 2.0, 2.0,
4.0, 1.0, 2.0, 2.0, 2.0));
print! separator
end
(*------------------------------------------------------------------*)
</syntaxhighlight>
{{out}}
(The following settings on my computer ''will'' result in GCC generating fused-multiply-and-add instructions [and thus the ''-lm'' flag actually is not needed]. See comments in the program for a discussion of the matter.)
<pre style="font-size:90%">$ patscc -std=gnu2x -g -O3 -march=native -DATS_MEMALLOC_GCBDW gauss_jordan_task.dats -lgc -lm && ./a.out
(The examples are printed here after rounding to 4 decimal places.)
[ 0.0000 ] appears to be singular
[ 4.0000 ] ✕ [ 0.2500 ] = [ 1.0000 ]
[ 1.0000 1.0000 ]
[ 1.0000 1.0000 ] appears to be singular
[ 2.0000 0.0000 ] [ 0.5000 0.0000 ] [ 1.0000 0.0000 ]
[ 0.0000 2.0000 ] ✕ [ 0.0000 0.5000 ] = [ 0.0000 1.0000 ]
[ 1.0000 3.0000 ] [ -2.0000 1.5000 ] [ 1.0000 0.0000 ]
[ 2.0000 4.0000 ] ✕ [ 1.0000 -0.5000 ] = [ 0.0000 1.0000 ]
An orthonormal matrix of rank 5 (its inverse will equal its transpose):
[ -0.0801 -0.5604 -0.2402 -0.3203 -0.7206 ] [ -0.0801 -0.1885 -0.8955 0.2715 -0.2872 ] [ 1.0000 -0.0000 -0.0000 0.0000 0.0000 ]
[ -0.1885 -0.4588 0.1517 -0.6105 0.5985 ] [ -0.5604 -0.4588 -0.0195 0.1610 0.6702 ] [ 0.0000 1.0000 -0.0000 0.0000 0.0000 ]
[ -0.8955 -0.0195 -0.3178 0.2984 0.0880 ] ✕ [ -0.2402 0.1517 -0.3178 -0.8979 0.1094 ] = [ 0.0000 -0.0000 1.0000 0.0000 0.0000 ]
[ 0.2715 0.1610 -0.8979 -0.1997 0.2327 ] [ -0.3203 -0.6105 0.2984 -0.1997 -0.6291 ] [ -0.0000 0.0000 -0.0000 1.0000 -0.0000 ]
[ -0.2872 0.6702 0.1094 -0.6291 -0.2462 ] [ -0.7206 0.5985 0.0880 0.2327 -0.2462 ] [ -0.0000 0.0000 -0.0000 0.0000 1.0000 ]
A 5✕5 singular matrix:
[ 1.0000 5.0000 3.0000 1.0000 4.0000 ]
[ 1.0000 1.0000 2.0000 1.0000 1.0000 ]
[ 2.0000 2.0000 3.0000 2.0000 2.0000 ] appears to be singular
[ 2.0000 2.0000 5.0000 2.0000 2.0000 ]
[ 2.0000 2.0000 7.0000 2.0000 2.0000 ]
</pre>
=={{header|BASIC}}==
==={{header|FreeBASIC}}===
The include statements use code from [[Matrix multiplication#FreeBASIC]], which contains the Matrix type used here, and [[Reduced row echelon form#FreeBASIC]] which contains the function for reducing a matrix to row-echelon form. Make sure to take out all the print statements first!
<syntaxhighlight lang="freebasic">#include once "matmult.bas"
#include once "rowech.bas"
function matinv( A as Matrix ) as Matrix
dim ret as Matrix, working as Matrix
dim as uinteger i, j
if not ubound( A.m, 1 ) = ubound( A.m, 2 ) then return ret
dim as integer n = ubound(A.m, 1)
redim ret.m( n, n )
working = Matrix( n+1, 2*(n+1) )
for i = 0 to n
for j = 0 to n
working.m(i, j) = A.m(i, j)
next j
working.m(i, i+n+1) = 1
next i
working = rowech(working)
for i = 0 to n
for j = 0 to n
ret.m(i, j) = working.m(i, j+n+1)
next j
next i
return ret
end function
dim as integer i, j
dim as Matrix M = Matrix(3,3)
for i = 0 to 2
for j = 0 to 2
M.m(i, j) = 1 + i*i + 3*j + (i+j) mod 2 'just some arbitrary values
print M.m(i, j),
next j
print
next i
print
M = matinv(M)
for i = 0 to 2
for j = 0 to 2
print M.m(i, j),
next j
print
next i</syntaxhighlight>
{{out}}
<pre>
1 5 7
3 5 9
5 9 11
-0.5416666666666667 0.1666666666666667 0.2083333333333333
0.2499999999999999 -0.5 0.25
0.0416666666666667 0.3333333333333333 -0.2083333333333333
</pre>
==={{header|VBA}}===
{{trans|Phix}}
Uses ToReducedRowEchelonForm() from [[Reduced_row_echelon_form#VBA]]<syntaxhighlight lang="vb">Private Function inverse(mat As Variant) As Variant
Dim len_ As Integer: len_ = UBound(mat)
Dim tmp() As Variant
ReDim tmp(2 * len_ + 1)
Dim aug As Variant
ReDim aug(len_)
For i = 0 To len_
If UBound(mat(i)) <> len_ Then Debug.Print 9 / 0 '-- "Not a square matrix"
aug(i) = tmp
For j = 0 To len_
aug(i)(j) = mat(i)(j)
Next j
'-- augment by identity matrix to right
aug(i)(i + len_ + 1) = 1
Next i
aug = ToReducedRowEchelonForm(aug)
Dim inv As Variant
inv = mat
'-- remove identity matrix to left
For i = 0 To len_
For j = len_ + 1 To 2 * len_ + 1
inv(i)(j - len_ - 1) = aug(i)(j)
Next j
Next i
inverse = inv
End Function
Public Sub main()
Dim test As Variant
test = inverse(Array( _
Array(2, -1, 0), _
Array(-1, 2, -1), _
Array(0, -1, 2)))
For i = LBound(test) To UBound(test)
For j = LBound(test(0)) To UBound(test(0))
Debug.Print test(i)(j),
Next j
Debug.Print
Next i
End Sub</syntaxhighlight>{{out}}
<pre> 0,75 0,5 0,25
0,5 1 0,5
0,25 0,5 0,75 </pre>
==={{header|VBScript}}===
{{trans|Rexx}}
To run in console mode with cscript.
<syntaxhighlight lang="vb">' Gauss-Jordan matrix inversion - VBScript - 22/01/2021
Option Explicit
Function rref(m)
Dim r, c, i, n, div, wrk
n=UBound(m)
For r = 1 To n 'row
div = m(r,r)
If div <> 0 Then
For c = 1 To n*2 'col
m(r,c) = m(r,c) / div
Next 'c
Else
WScript.Echo "inversion impossible!"
WScript.Quit
End If
For i = 1 To n 'row
If i <> r Then
wrk = m(i,r)
For c = 1 To n*2
m(i,c) = m(i,c) - wrk * m(r,c)
Next' c
End If
Next 'i
Next 'r
rref = m
End Function 'rref
Function inverse(mat)
Dim i, j, aug, inv, n
n = UBound(mat)
ReDim inv(n,n), aug(n,2*n)
For i = 1 To n
For j = 1 To n
aug(i,j) = mat(i,j)
Next 'j
aug(i,i+n) = 1
Next 'i
aug = rref(aug)
For i = 1 To n
For j = n+1 To 2*n
inv(i,j-n) = aug(i,j)
Next 'j
Next 'i
inverse = inv
End Function 'inverse
Sub wload(m)
Dim i, j, k
k = -1
For i = 1 To n
For j = 1 To n
k = k + 1
m(i,j) = w(k)
Next 'j
Next 'i
End Sub 'wload
Sub show(c, m, t)
Dim i, j, buf
buf = "Matrix " & c &"=" & vbCrlf & vbCrlf
For i = 1 To n
For j = 1 To n
If t="fixed" Then
buf = buf & FormatNumber(m(i,j),6,,,0) & " "
Else
buf = buf & m(i,j) & " "
End If
Next 'j
buf=buf & vbCrlf
Next 'i
WScript.Echo buf
End Sub 'show
Dim n, a, b, c, w
w = Array( _
2, 1, 1, 4, _
0, -1, 0, -1, _
1, 0, -2, 4, _
4, 1, 2, 2)
n=Sqr(UBound(w)+1)
ReDim a(n,n), b(n,n), c(n,n)
wload a
show "a", a, "simple"
b = inverse(a)
show "b", b, "fixed"
c = inverse(b)
show "c", c, "fixed"
</syntaxhighlight>
{{out}}
<pre>
Matrix a=
2 1 1 4
0 -1 0 -1
1 0 -2 4
4 1 2 2
Matrix b=
-0,400000 0,000000 0,200000 0,400000
-0,400000 -1,200000 0,000000 0,200000
0,600000 0,400000 -0,400000 -0,200000
0,400000 0,200000 -0,000000 -0,200000
Matrix c=
2,000000 1,000000 1,000000 4,000000
0,000000 -1,000000 0,000000 -1,000000
1,000000 0,000000 -2,000000 4,000000
4,000000 1,000000 2,000000 2,000000
</pre>
Line 1,011 ⟶ 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 1,051 ⟶ 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 1,262 ⟶ 2,513:
7.000 -8.000 9.000 1.000
1.000 -2.000 1.000 3.000
</pre>
Line 2,471 ⟶ 3,666:
then .i = $r
| .lead += 1
else .
end
Line 2,760 ⟶ 3,954:
[-0.05555555555555555,0.1111111111111111,-0.05555555555555555]]
</syntaxhighlight>
=={{header|M2000 Interpreter}}==
We can use supported types like Decimal, Currency, Double, Single. Also matrix may have different base on each dimension.
<syntaxhighlight lang="m2000 interpreter">
Module Gauss_Jordan_matrix_inversion (there$="") {
open there$ for output as #f
dim matrix(1 to 3, 2 to 4) as Double
matrix(1,2)=2, -1, 0, -1, 2, -1, 0, -1, 2
print #f, "Matrix 3x3:"
disp(matrix())
print #f, "Inverse:"
disp(@inverse(matrix()))
dim matrix2(2 to 4, 3 to 5) as single
matrix2(2, 3)=1, 2, 3, 4, 1, 6, 7, 8, 9
print #f,"Matrix 3x3:"
disp(matrix2())
print #f,"Inverse:"
disp(@inverse(matrix2()))
dim matrix3(1 to 4,1 to 4) as Currency, matrix4()
matrix3(1,1)=-1,-2,3,2,-4,-1,6,2,7,-8,9,1,1,-2,1,3
print #f,"Matrix 4x4:"
disp(matrix3())
matrix4()=@inverse(matrix3())
Rem Print Type$(matrix4(1,1))="Currency"
print #f,"Inverse:"
disp(matrix4())
close #f
function inverse(a())
if dimension(a())<>2 then Error "Not 2 Dimensions Matrix"
if not dimension(a(),1)=dimension(a(),2) then Error "Not Square Matrix"
Local aug() // this is an empty array (type of variant)
aug()=a() // so we get the same type from a()
dim aug(dimension(a(),1), dimension(a(),1)*2)=0
Rem gosub disp(aug())
Local long bi=dimension(a(), 1, 0), bj=dimension(a(), 2, 0)
Local long bim=dimension(a(), 1)-1, bjm=dimension(a(), 2)-1
Local long i, j
for i=0 to bim
for j=0 to bjm
aug(i,j)=a(i+bi,j+bj)
next
aug(i, i+bjm+1)=1
next
aug()=@ToReducedRowEchelonForm(aug())
for i=0 to bim
for j=0 to bjm
a(i+bi,j+bj)=aug(i,j+bjm+1)
next
next
Rem : Print type$(aug(0,0))
=a()
end function
function ToReducedRowEchelonForm(a())
Local long bi=dimension(a(), 1, 0), bj=dimension(a(), 2, 0)
Local long rowcount=dimension(a(), 1), columncount=dimension(a(), 2)
Local long lead=bj, r, i, j
for r=bi to rowcount-1+bi {
if columncount<=lead then exit
i=r
while A(i,lead)=0 {
i++
if rowcount=i then i=r : lead++ : if columncount<lead then exit
}
for c =bj to columncount-1+bj {swap A(i, c), A(r, c)}
if A(r, lead)<>0 then
div1=A(r,lead)
For c =bj to columncount-1+bj {A(r, c)/=div1}
end if
for i=bi to rowcount-1+bi {
if i<>r then {
mult=A(i,lead)
for j=bj to columncount-1+bj {A(i,j)-=A(r,j)*mult}
}
}
lead++
}
=a()
End function
sub disp(a())
local long i, j, c=8 // column width plus 1 space
for i=dimension(a(), 1, 0) to dimension(a(), 1, 1)
for j=dimension(a(), 2, 0) to dimension(a(), 2, 1)
Print #f, format$("{0:-"+c+"} ",str$(round(A(i, j), 5),1033));
Next
// if pos>0 then print // not used here
Print #f, ""
Next
Print #f, ""
End sub
}
Gauss_Jordan_matrix_inversion // to screen
Gauss_Jordan_matrix_inversion "out.txt"
win "notepad", dir$+"out.txt"
</syntaxhighlight>
{{out}}
<pre>
Matrix 3x3:
2 -1 0
-1 2 -1
0 -1 2
Inverse:
0.75 0.5 0.25
0.5 1 0.5
0.25 0.5 0.75
Matrix 3x3:
1 2 3
4 1 6
7 8 9
Inverse:
-0.8125 0.125 0.1875
0.125 -0.25 0.125
0.52083 0.125 -0.14583
Matrix 4x4:
-1 -2 3 2
-4 -1 6 2
7 -8 9 1
1 -2 1 3
Inverse:
-0.9129 0.2463 0.0941 0.413
-1.6519 0.6522 0.0434 0.6521
-0.6954 0.3623 0.0797 0.1956
-0.5651 0.2319 -0.029 0.5652
</pre>
=={{header|Mathematica}} / {{header|Wolfram Language}}==
Line 2,981 ⟶ 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 4,014 ⟶ 5,490:
-0.1833 -0.7009 -0.1425 0.6163
-0.1064 -0.0855 0.0883 0.0779
</pre>
=={{header|Ruby}}==
<syntaxhighlight lang="ruby">require 'matrix'
m = Matrix[[-1, -2, 3, 2],
[-4, -1, 6, 2],
[ 7, -8, 9, 1],
[ 1, -2, 1, 3]]
pp m.inv.row_vectors
</syntaxhighlight>
{{out}}
<pre>[Vector[(-21/23), (17/69), (13/138), (19/46)],
Vector[(-38/23), (15/23), (1/23), (15/23)],
Vector[(-16/23), (25/69), (11/138), (9/46)],
Vector[(-13/23), (16/69), (-2/69), (13/23)]]
</pre>
Line 4,159 ⟶ 5,652:
[0.24999999999999997, 0.49999999999999994, 0.7499999999999999]
</pre>
=={{header|Sidef}}==
Uses the '''rref(M)''' function from [https://rosettacode.org/wiki/Reduced_row_echelon_form#Sidef Reduced row echelon form].
Line 4,194 ⟶ 5,688:
-13/23 16/69 -2/69 13/23
</pre>
=={{header|Wren}}==
Line 4,363 ⟶ 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 = [
Line 4,456 ⟶ 5,786:
0.2500, 0.5000, 0.7500
</pre>
{{omit from|PL/0}}
| |||