Gauss-Jordan matrix inversion
You are encouraged to solve this task according to the task description, using any language you may know.
- Task
Invert matrix A using Gauss-Jordan method.
A being an n × n matrix.
11l
V Eps = 1e-10
F transformToRref(&mat)
V lead = 0
L(r) 0 .< mat.len
I lead >= mat[0].len
R
V i = r
L mat[i][lead] == 0
i++
I i == mat.len
i = r
lead++
I lead == mat[0].len
R
swap(&mat[i], &mat[r])
V d = mat[r][lead]
I abs(d) > Eps
L(&item) mat[r]
item /= d
L(i) 0 .< mat.len
I i != r
V m = mat[i][lead]
L(c) 0 .< mat[0].len
mat[i][c] -= mat[r][c] * m
lead++
F inverse(mat)
V augmat = [[0.0] * (2 * mat.len)] * mat.len
L(i) 0 .< mat.len
L(j) 0 .< mat.len
augmat[i][j] = mat[i][j]
augmat[i][mat.len + i] = 1
transformToRref(&augmat)
V result = [[0.0] * mat.len] * mat.len
L(i) 0 .< mat.len
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
F print_mat(mat)
L(row) mat
V line = ‘’
L(val) row
I !line.empty
line ‘’= ‘ ’
line ‘’= ‘#3.5’.format(val)
print(line)
F runTest(mat)
print(‘Matrix:’)
print_mat(mat)
print()
print(‘Inverse:’)
print_mat(inverse(mat))
print()
print()
V m1 = [[Float(1), 2, 3],
[Float(4), 1, 6],
[Float(7), 8, 9]]
V m2 = [[Float( 2), -1, 0],
[Float(-1), 2, -1],
[Float( 0), -1, 2]]
V m3 = [[Float(-1), -2, 3, 2],
[Float(-4), -1, 6, 2],
[Float( 7), -8, 9, 1],
[Float( 1), -2, 1, 3]]
runTest(m1)
runTest(m2)
runTest(m3)- Output:
Matrix: 1.00000 2.00000 3.00000 4.00000 1.00000 6.00000 7.00000 8.00000 9.00000 Inverse: -0.81250 0.12500 0.18750 0.12500 -0.25000 0.12500 0.52083 0.12500 -0.14583 Matrix: 2.00000 -1.00000 0.00000 -1.00000 2.00000 -1.00000 0.00000 -1.00000 2.00000 Inverse: 0.75000 0.50000 0.25000 0.50000 1.00000 0.50000 0.25000 0.50000 0.75000 Matrix: -1.00000 -2.00000 3.00000 2.00000 -4.00000 -1.00000 6.00000 2.00000 7.00000 -8.00000 9.00000 1.00000 1.00000 -2.00000 1.00000 3.00000 Inverse: -0.91304 0.24638 0.09420 0.41304 -1.65217 0.65217 0.04348 0.65217 -0.69565 0.36232 0.07971 0.19565 -0.56522 0.23188 -0.02899 0.56522
360 Assembly
The COPY file of FORMAT, to format a floating point number, can be found at: Include files 360 Assembly.
* Gauss-Jordan matrix inversion 17/01/2021
GAUSSJOR CSECT
USING GAUSSJOR,R13 base register
B 72(R15) skip savearea
DC 17F'0' savearea
SAVE (14,12) save previous context
ST R13,4(R15) link backward
ST R15,8(R13) link forward
LR R13,R15 set addressability
LA R6,1 i=1
DO WHILE=(C,R6,LE,N) do i=1 to n
LA R7,1 j=1
DO WHILE=(C,R7,LE,N) do j=1 to n
LR R1,R6 i
BCTR R1,0 -1
MH R1,HN *n
LR R0,R7 j
BCTR R0,0 -1
AR R1,R0 (i-1)*n+j-1
SLA R1,2 *4
LE F0,AA(R1) aa(i,j)
LR R1,R6 i
BCTR R1,0 -1
MH R1,HN2 *n*2
LR R0,R7 j
BCTR R0,0 -1
AR R1,R0 (i-1)*n*2+j-1
SLA R1,2 *4
STE F0,TMP(R1) tmp(i,j)=aa(i,j)
LA R7,1(R7) j++
ENDDO , enddo j
LA R7,1 j=1
A R7,N j=n+1
DO WHILE=(C,R7,LE,N2) do j=n+1 to 2*n
LR R1,R6 i
BCTR R1,0 -1
MH R1,HN2 *n*2
LR R0,R7 j
BCTR R0,0 -1
AR R1,R0 (i-1)*n*2+j-1
SLA R1,2 *4
LE F0,=E'0' 0
STE F0,TMP(R1) tmp(i,j)=0
LA R7,1(R7) j++
ENDDO , enddo j
LR R2,R6 i
A R2,N i+n
LR R1,R6 i
BCTR R1,0 -1
MH R1,HN2 *n*2
BCTR R2,0 -1
AR R1,R2 (i+n-1)*n*2+i-1
SLA R1,2 *4
LE F0,=E'1' 1
STE F0,TMP(R1) tmp(i,i+n)=1
LA R6,1(R6) i++
ENDDO , enddo i
LA R6,1 i=1
DO WHILE=(C,R6,LE,N) do r=1 to n
LR R1,R6 r
BCTR R1,0 -1
MH R1,HN2 *n*2
LR R0,R6 r
BCTR R0,0 -1
AR R1,R0 (r-1)*n*2+r-1
SLA R1,2 *4
LE F0,TMP(R1) tmp(r,r)
STE F0,T t=tmp(r,r)
LA R7,1 c=1
DO WHILE=(C,R7,LE,N2) do c=1 to n*2
LR R1,R6 r
BCTR R1,0
MH R1,HN2
LR R0,R7 c
BCTR R0,0
AR R1,R0
SLA R1,2
LE F0,TMP(R1)
LTER F0,F0
BZ *+8
DE F0,T tmp(r,c)/t
LR R1,R6 r
BCTR R1,0
MH R1,HN2
LR R0,R7 c
BCTR R0,0
AR R1,R0
SLA R1,2
STE F0,TMP(R1) tmp(r,c)=tmp(r,c)/t
LA R7,1(R7) c++
ENDDO , enddo c
LA R8,1 i=1
DO WHILE=(C,R8,LE,N) do i=1 to n
IF CR,R8,NE,R6 THEN if i^=r then
LR R1,R8 i
BCTR R1,0
MH R1,HN2
LR R0,R6 r
BCTR R0,0
AR R1,R0
SLA R1,2
LE F0,TMP(R1) tmp(i,r)
STE F0,T t=tmp(i,r)
LA R7,1 c=1
DO WHILE=(C,R7,LE,N2) do c=1 to n*2
LR R2,R8 i
BCTR R2,0
MH R2,HN2
LR R0,R7 c
BCTR R0,0
AR R2,R0
SLA R2,2
LE F0,TMP(R2) f0=tmp(i,c)
LR R1,R6 r
BCTR R1,0
MH R1,HN2
LR R0,R7 c
BCTR R0,0
AR R1,R0
SLA R1,2
LE F2,TMP(R1) f2=tmp(r,c)
LE F4,T t
MER F4,F2 f4=t*tmp(r,c)
SER F0,F4 f0=tmp(i,c)-t*tmp(r,c)
STE F0,TMP(R2) tmp(i,c)=f0
LA R7,1(R7) c++
ENDDO , enddo c
ENDIF , endif
LA R8,1(R8) i++
ENDDO , enddo i
LA R6,1(R6) r++
ENDDO , enddo r
LA R6,1 i=1
DO WHILE=(C,R6,LE,N) do i=1 to n
L R7,N n
LA R7,1(R7) j=n+1
DO WHILE=(C,R7,LE,N2) do j=n+1 to 2*n
LR R2,R7 j
S R2,N -n
LR R1,R6 i
BCTR R1,0
MH R1,HN2 *2*n
LR R0,R7 j
BCTR R0,0
AR R1,R0
SLA R1,2
LE F0,TMP(R1) tmp(i,j)
LR R1,R6 i
BCTR R1,0
MH R1,HN *n
BCTR R2,0
AR R1,R2
SLA R1,2
STE F0,BB(R1) bb(i,j-n)=tmp(i,j)
LA R7,1(R7) j++
ENDDO , enddo j
LA R6,1(R6) i++
ENDDO , enddo i
LA R6,1 i=1
DO WHILE=(C,R6,LE,N) do i=1 to n
LA R3,PG @pg
LA R7,1 j=1
DO WHILE=(C,R7,LE,N) do j=1 to n
LR R1,R6 i
BCTR R1,0
MH R1,HN *n
LR R0,R7 j
BCTR R0,0
AR R1,R0
SLA R1,2
LE F0,BB(R1) bb(i,j)
LA R0,5 number of decimals
BAL R14,FORMATF FormatF(bb(i,j))
MVC 0(13,R3),0(R1) edit
LA R3,13(R3) @pg+=13
LA R7,1(R7) j++
ENDDO , enddo j
XPRNT PG,L'PG print buffer
LA R6,1(R6) i++
ENDDO , enddo i
L R13,4(0,R13) restore previous savearea pointer
RETURN (14,12),RC=0 restore registers from calling save
COPY plig\$_FORMATF.MLC format a float
NN EQU 5
N DC A(NN)
N2 DC A(NN*2)
AA DC E'3.0',E'1.0',E'5.0',E'9.0',E'7.0'
DC E'8.0',E'2.0',E'8.0',E'0.0',E'1.0'
DC E'1.0',E'7.0',E'2.0',E'0.0',E'3.0'
DC E'0.0',E'1.0',E'7.0',E'0.0',E'9.0'
DC E'3.0',E'5.0',E'6.0',E'1.0',E'1.0'
BB DS (NN*NN)E
TMP DS (NN*NN*2)E
T DS E
HN DC AL2(NN)
HN2 DC AL2(NN*2)
PG DC CL80' ' buffer
REGEQU
END GAUSSJOR- Output:
0.02863 0.19429 0.13303 -0.05956 -0.25778
-0.00580 -0.03255 0.12109 -0.03803 0.05226
-0.03020 -0.06824 -0.17903 0.06054 0.27187
0.10021 -0.06733 -0.05617 -0.06263 0.09806
0.02413 0.05669 0.12579 0.06824 -0.21726
ALGOL 60
begin
comment Gauss-Jordan matrix inversion - 22/01/2021;
integer n;
n:=4;
begin
procedure rref(m); real array m;
begin
integer r, c, i;
real d, w;
for r := 1 step 1 until n do begin
d := m[r,r];
if d notequal 0 then
for c := 1 step 1 until n*2 do
m[r,c] := m[r,c] / d
else
outstring(1,"inversion impossible!\n");
for i := 1 step 1 until n do
if i notequal r then begin
w := m[i,r];
for c := 1 step 1 until n*2 do
m[i,c] := m[i,c] - w * m[r,c]
end
end
end rref;
procedure inverse(mat,inv); real array mat, inv;
begin
integer i, j;
real array aug[1:n,1:n*2];
for i := 1 step 1 until n do begin
for j := 1 step 1 until n do
aug[i,j] := mat[i,j];
for j := 1+n step 1 until n*2 do
aug[i,j] := 0;
aug[i,i+n] := 1
end;
rref(aug);
for i := 1 step 1 until n do
for j := n+1 step 1 until 2*n do
inv[i,j-n] := aug[i,j]
end inverse;
procedure show(c, m); string c; real array m;
begin
integer i, j;
outstring(1,"matrix "); outstring(1,c); outstring(1,"\n");
for i := 1 step 1 until n do begin
for j := 1 step 1 until n do
outreal(1,m[i,j]);
outstring(1,"\n")
end
end show;
integer i;
real array a[1:n,1:n], b[1:n,1:n], c[1:n,1:n];
i:=1; a[i,1]:= 2; a[i,2]:= 1; a[i,3]:= 1; a[i,4]:= 4;
i:=2; a[i,1]:= 0; a[i,2]:=-1; a[i,3]:= 0; a[i,4]:=-1;
i:=3; a[i,1]:= 1; a[i,2]:= 0; a[i,3]:=-2; a[i,4]:= 4;
i:=4; a[i,1]:= 4; a[i,2]:= 1; a[i,3]:= 2; a[i,4]:= 2;
show("a",a);
inverse(a,b);
show("b",b);
inverse(b,c);
show("c",c)
end
end- Output:
matrix a 2 1 1 4 0 -1 0 -1 1 0 -2 4 4 1 2 2 matrix b -0.4 0 0.2 0.4 -0.4 -1.2 0 0.2 0.6 0.4 -0.4 -0.2 0.4 0.2 0 -0.2 matrix c 2 1 1 4 0 -1 0 -1 1 0 -2 4 4 1 2 2
ALGOL 68
BEGIN # Gauss-Jordan matrix inversion #
PROC rref = ( REF[,]REAL m )VOID:
BEGIN
INT n = 1 UPB m;
FOR r TO n DO
REAL d = m[ r, r ];
IF d /= 0 THEN
FOR c TO n * 2 DO m[ r, c ] := m[ r, c ] / d OD
ELSE
print( ( "inversion impossible!", newline ) )
FI;
FOR i TO n DO
IF i /= r THEN
REAL w = m[ i, r ];
FOR c TO n * 2 DO m[ i, c ] := m[ i, c ] - w * m[ r, c ] OD
FI
OD
OD
END # rref # ;
PROC inverse = ( REF[,]REAL mat )[,]REAL:
BEGIN
INT n = 1 UPB mat;
[ 1 : n, 1 : n ]REAL inv;
[ 1 : n, 1 : n * 2 ]REAL aug;
FOR i TO n DO
FOR j TO n DO aug[ i, j ] := mat[ i, j ] OD;
FOR j FROM 1 + n TO n * 2 DO aug[ i, j ] := 0 OD;
aug[ i, i + n ] := 1
OD;
rref( aug );
FOR i TO n DO
FOR j FROM n + 1 TO 2 * n DO inv[ i, j - n ] := aug[ i, j ] OD
OD;
inv
END # inverse # ;
PROC show = ( STRING c, [,]REAL m )VOID:
BEGIN
INT n = 1 UPB m;
print( ( "matrix ", c, newline ) );
FOR i TO n DO
FOR j TO n DO print( ( " ", fixed( m[ i, j ], -8, 3 ) ) ) OD;
print( ( newline ) )
OD
END # show #;
INT n = 4;
[ 1 : n, 1 : n ]REAL a, b, c;
a := [,]REAL( ( 2, 1, 1, 4 )
, ( 0, -1, 0, -1 )
, ( 1, 0, -2, 4 )
, ( 4, 1, 2, 2 )
);
show( "a", a );
b := inverse( a );
show( "b", b );
c := inverse( b );
show( "c", c )
END- Output:
matrix a
2.000 1.000 1.000 4.000
0.000 -1.000 0.000 -1.000
1.000 0.000 -2.000 4.000
4.000 1.000 2.000 2.000
matrix b
-0.400 0.000 0.200 0.400
-0.400 -1.200 0.000 0.200
0.600 0.400 -0.400 -0.200
0.400 0.200 0.000 -0.200
matrix c
2.000 1.000 1.000 4.000
0.000 -1.000 0.000 -1.000
1.000 0.000 -2.000 4.000
4.000 1.000 2.000 2.000
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.
(* 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
(*------------------------------------------------------------------*)- Output:
(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.)
$ 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 ]
BASIC
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!
#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
- Output:
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
VBA
Uses ToReducedRowEchelonForm() from Reduced_row_echelon_form#VBA
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
- Output:
0,75 0,5 0,25 0,5 1 0,5 0,25 0,5 0,75
VBScript
To run in console mode with cscript.
' 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"
- Output:
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
C
/*----------------------------------------------------------------------
gjinv - Invert a matrix, Gauss-Jordan algorithm
A is destroyed. Returns 1 for a singular matrix.
___Name_____Type______In/Out____Description_____________________________
a[n*n] double* In An N by N matrix
n int In Order of matrix
b[n*n] double* Out Inverse of A
----------------------------------------------------------------------*/
#include <math.h>
int gjinv (double *a, int n, double *b)
{
int i, j, k, p;
double f, g, tol;
if (n < 1) return -1; /* Function Body */
f = 0.; /* Frobenius norm of a */
for (i = 0; i < n; ++i) {
for (j = 0; j < n; ++j) {
g = a[j+i*n];
f += g * g;
}
}
f = sqrt(f);
tol = f * 2.2204460492503131e-016;
for (i = 0; i < n; ++i) { /* Set b to identity matrix. */
for (j = 0; j < n; ++j) {
b[j+i*n] = (i == j) ? 1. : 0.;
}
}
for (k = 0; k < n; ++k) { /* Main loop */
f = fabs(a[k+k*n]); /* Find pivot. */
p = k;
for (i = k+1; i < n; ++i) {
g = fabs(a[k+i*n]);
if (g > f) {
f = g;
p = i;
}
}
if (f < tol) return 1; /* Matrix is singular. */
if (p != k) { /* Swap rows. */
for (j = k; j < n; ++j) {
f = a[j+k*n];
a[j+k*n] = a[j+p*n];
a[j+p*n] = f;
}
for (j = 0; j < n; ++j) {
f = b[j+k*n];
b[j+k*n] = b[j+p*n];
b[j+p*n] = f;
}
}
f = 1. / a[k+k*n]; /* Scale row so pivot is 1. */
for (j = k; j < n; ++j) a[j+k*n] *= f;
for (j = 0; j < n; ++j) b[j+k*n] *= f;
for (i = 0; i < n; ++i) { /* Subtract to get zeros. */
if (i == k) continue;
f = a[k+i*n];
for (j = k; j < n; ++j) a[j+i*n] -= a[j+k*n] * f;
for (j = 0; j < n; ++j) b[j+i*n] -= b[j+k*n] * f;
}
}
return 0;
} /* end of gjinv */
Test program:
/* Test matrix inversion */
#include <stdio.h>
int main (void)
{
static double aorig[16] = { -1.,-2.,3.,2.,-4.,
-1.,6.,2.,7.,-8.,9.,1.,1.,-2.,1.,3. };
double a[16], b[16], c[16];
int n = 4;
int i, j, k, ierr;
for (i = 0; i < n; ++i) {
for (j = 0; j < n; ++j) {
a[j+i*n] = aorig[j+i*n];
}
}
ierr = gjinv (a, n, b);
printf("gjinv returns #%i\n\n", ierr);
printf("matrix:\n");
for (i = 0; i < n; ++i) {
for (j = 0; j < n; ++j) {
printf("%8.3f", aorig[j+i*n]);
}
printf("\n");
}
printf("\ninverse:\n");
for (i = 0; i < n; ++i) {
for (j = 0; j < n; ++j) {
printf("%8.3f", b[j+i*n]);
}
printf("\n");
}
for (j = 0; j < n; ++j) {
for (k = 0; k < n; ++k) {
c[k+j*n] = 0.;
for (i = 0; i < n; ++i) {
c[k+j*n] += aorig[i+j*n] * b[k+i*n];
}
}
}
printf("\nmatrix @ inverse:\n");
for (i = 0; i < n; ++i) {
for (j = 0; j < n; ++j) {
printf("%8.3f", c[j+i*n]);
}
printf("\n");
}
ierr = gjinv (b, n, a);
printf("\ngjinv returns #%i\n", ierr);
printf("\ninverse of inverse:\n");
for (i = 0; i < n; ++i) {
for (j = 0; j < n; ++j) {
printf("%8.3f", a[j+i*n]);
}
printf("\n");
}
return 0;
} /* end of test program */
Output is the same as the Fortran example.
C#
using System;
namespace Rosetta
{
internal class Vector
{
private double[] b;
internal readonly int rows;
internal Vector(int rows)
{
this.rows = rows;
b = new double[rows];
}
internal Vector(double[] initArray)
{
b = (double[])initArray.Clone();
rows = b.Length;
}
internal Vector Clone()
{
Vector v = new Vector(b);
return v;
}
internal double this[int row]
{
get { return b[row]; }
set { b[row] = value; }
}
internal void SwapRows(int r1, int r2)
{
if (r1 == r2) return;
double tmp = b[r1];
b[r1] = b[r2];
b[r2] = tmp;
}
internal double norm(double[] weights)
{
double sum = 0;
for (int i = 0; i < rows; i++)
{
double d = b[i] * weights[i];
sum += d*d;
}
return Math.Sqrt(sum);
}
internal void print()
{
for (int i = 0; i < rows; i++)
Console.WriteLine(b[i]);
Console.WriteLine();
}
public static Vector operator-(Vector lhs, Vector rhs)
{
Vector v = new Vector(lhs.rows);
for (int i = 0; i < lhs.rows; i++)
v[i] = lhs[i] - rhs[i];
return v;
}
}
class Matrix
{
private double[] b;
internal readonly int rows, cols;
internal Matrix(int rows, int cols)
{
this.rows = rows;
this.cols = cols;
b = new double[rows * cols];
}
internal Matrix(int size)
{
this.rows = size;
this.cols = size;
b = new double[rows * cols];
for (int i = 0; i < size; i++)
this[i, i] = 1;
}
internal Matrix(int rows, int cols, double[] initArray)
{
this.rows = rows;
this.cols = cols;
b = (double[])initArray.Clone();
if (b.Length != rows * cols) throw new Exception("bad init array");
}
internal double this[int row, int col]
{
get { return b[row * cols + col]; }
set { b[row * cols + col] = value; }
}
public static Vector operator*(Matrix lhs, Vector rhs)
{
if (lhs.cols != rhs.rows) throw new Exception("I can't multiply matrix by vector");
Vector v = new Vector(lhs.rows);
for (int i = 0; i < lhs.rows; i++)
{
double sum = 0;
for (int j = 0; j < rhs.rows; j++)
sum += lhs[i,j]*rhs[j];
v[i] = sum;
}
return v;
}
internal void SwapRows(int r1, int r2)
{
if (r1 == r2) return;
int firstR1 = r1 * cols;
int firstR2 = r2 * cols;
for (int i = 0; i < cols; i++)
{
double tmp = b[firstR1 + i];
b[firstR1 + i] = b[firstR2 + i];
b[firstR2 + i] = tmp;
}
}
//with partial pivot
internal bool InvPartial()
{
const double Eps = 1e-12;
if (rows != cols) throw new Exception("rows != cols for Inv");
Matrix M = new Matrix(rows); //unitary
for (int diag = 0; diag < rows; diag++)
{
int max_row = diag;
double max_val = Math.Abs(this[diag, diag]);
double d;
for (int row = diag + 1; row < rows; row++)
if ((d = Math.Abs(this[row, diag])) > max_val)
{
max_row = row;
max_val = d;
}
if (max_val <= Eps) return false;
SwapRows(diag, max_row);
M.SwapRows(diag, max_row);
double invd = 1 / this[diag, diag];
for (int col = diag; col < cols; col++)
{
this[diag, col] *= invd;
}
for (int col = 0; col < cols; col++)
{
M[diag, col] *= invd;
}
for (int row = 0; row < rows; row++)
{
d = this[row, diag];
if (row != diag)
{
for (int col = diag; col < this.cols; col++)
{
this[row, col] -= d * this[diag, col];
}
for (int col = 0; col < this.cols; col++)
{
M[row, col] -= d * M[diag, col];
}
}
}
}
b = M.b;
return true;
}
internal void print()
{
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < cols; j++)
Console.Write(this[i,j].ToString()+" ");
Console.WriteLine();
}
}
}
}
using System;
namespace Rosetta
{
class Program
{
static void Main(string[] args)
{
Matrix M = new Matrix(4, 4, new double[] { -1, -2, 3, 2, -4, -1, 6, 2, 7, -8, 9, 1, 1, -2, 1, 3 });
M.InvPartial();
M.print();
}
}
}
- Output:
-0.913043478260869 0.246376811594203 0.0942028985507246 0.413043478260869 -1.65217391304348 0.652173913043478 0.0434782608695652 0.652173913043478 -0.695652173913043 0.36231884057971 0.0797101449275362 0.195652173913043 -0.565217391304348 0.231884057971014 -0.0289855072463768 0.565217391304348
C++
#include <algorithm>
#include <cassert>
#include <iomanip>
#include <iostream>
#include <vector>
template <typename scalar_type> class matrix {
public:
matrix(size_t rows, size_t columns)
: rows_(rows), columns_(columns), elements_(rows * columns) {}
matrix(size_t rows, size_t columns, scalar_type value)
: rows_(rows), columns_(columns), elements_(rows * columns, value) {}
matrix(size_t rows, size_t columns, std::initializer_list<scalar_type> values)
: rows_(rows), columns_(columns), elements_(rows * columns) {
assert(values.size() <= rows_ * columns_);
std::copy(values.begin(), values.end(), elements_.begin());
}
size_t rows() const { return rows_; }
size_t columns() const { return columns_; }
scalar_type& operator()(size_t row, size_t column) {
assert(row < rows_);
assert(column < columns_);
return elements_[row * columns_ + column];
}
const scalar_type& operator()(size_t row, size_t column) const {
assert(row < rows_);
assert(column < columns_);
return elements_[row * columns_ + column];
}
private:
size_t rows_;
size_t columns_;
std::vector<scalar_type> elements_;
};
template <typename scalar_type>
matrix<scalar_type> product(const matrix<scalar_type>& a,
const matrix<scalar_type>& b) {
assert(a.columns() == b.rows());
size_t arows = a.rows();
size_t bcolumns = b.columns();
size_t n = a.columns();
matrix<scalar_type> c(arows, bcolumns);
for (size_t i = 0; i < arows; ++i) {
for (size_t j = 0; j < n; ++j) {
for (size_t k = 0; k < bcolumns; ++k)
c(i, k) += a(i, j) * b(j, k);
}
}
return c;
}
template <typename scalar_type>
void swap_rows(matrix<scalar_type>& m, size_t i, size_t j) {
size_t columns = m.columns();
for (size_t column = 0; column < columns; ++column)
std::swap(m(i, column), m(j, column));
}
// Convert matrix to reduced row echelon form
template <typename scalar_type>
void rref(matrix<scalar_type>& m) {
size_t rows = m.rows();
size_t columns = m.columns();
for (size_t row = 0, lead = 0; row < rows && lead < columns; ++row, ++lead) {
size_t i = row;
while (m(i, lead) == 0) {
if (++i == rows) {
i = row;
if (++lead == columns)
return;
}
}
swap_rows(m, i, row);
if (m(row, lead) != 0) {
scalar_type f = m(row, lead);
for (size_t column = 0; column < columns; ++column)
m(row, column) /= f;
}
for (size_t j = 0; j < rows; ++j) {
if (j == row)
continue;
scalar_type f = m(j, lead);
for (size_t column = 0; column < columns; ++column)
m(j, column) -= f * m(row, column);
}
}
}
template <typename scalar_type>
matrix<scalar_type> inverse(const matrix<scalar_type>& m) {
assert(m.rows() == m.columns());
size_t rows = m.rows();
matrix<scalar_type> tmp(rows, 2 * rows, 0);
for (size_t row = 0; row < rows; ++row) {
for (size_t column = 0; column < rows; ++column)
tmp(row, column) = m(row, column);
tmp(row, row + rows) = 1;
}
rref(tmp);
matrix<scalar_type> inv(rows, rows);
for (size_t row = 0; row < rows; ++row) {
for (size_t column = 0; column < rows; ++column)
inv(row, column) = tmp(row, column + rows);
}
return inv;
}
template <typename scalar_type>
void print(std::ostream& out, const matrix<scalar_type>& m) {
size_t rows = m.rows(), columns = m.columns();
out << std::fixed << std::setprecision(4);
for (size_t row = 0; row < rows; ++row) {
for (size_t column = 0; column < columns; ++column) {
if (column > 0)
out << ' ';
out << std::setw(7) << m(row, column);
}
out << '\n';
}
}
int main() {
matrix<double> m(3, 3, {2, -1, 0, -1, 2, -1, 0, -1, 2});
std::cout << "Matrix:\n";
print(std::cout, m);
auto inv(inverse(m));
std::cout << "Inverse:\n";
print(std::cout, inv);
std::cout << "Product of matrix and inverse:\n";
print(std::cout, product(m, inv));
std::cout << "Inverse of inverse:\n";
print(std::cout, inverse(inv));
return 0;
}
- Output:
Matrix: 2.0000 -1.0000 0.0000 -1.0000 2.0000 -1.0000 0.0000 -1.0000 2.0000 Inverse: 0.7500 0.5000 0.2500 0.5000 1.0000 0.5000 0.2500 0.5000 0.7500 Product of matrix and inverse: 1.0000 0.0000 0.0000 -0.0000 1.0000 -0.0000 0.0000 0.0000 1.0000 Inverse of inverse: 2.0000 -1.0000 0.0000 -1.0000 2.0000 -1.0000 0.0000 -1.0000 2.0000
EasyLang
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[][]
Factor
Normally you would call recip to calculate the inverse of a matrix, but it uses a different method than Gauss-Jordan, so here's Gauss-Jordan.
USING: kernel math.matrices math.matrices.elimination
prettyprint sequences ;
! Augment a matrix with its identity. E.g.
!
! 1 2 3 1 2 3 1 0 0
! 4 5 6 augment-identity -> 4 5 6 0 1 0
! 7 8 9 7 8 9 0 0 1
: augment-identity ( matrix -- new-matrix )
dup first length <identity-matrix>
[ flip ] bi@ append flip ;
! Note: the 'solution' word finds the reduced row echelon form
! of a matrix.
: gauss-jordan-invert ( matrix -- inverted )
dup square-matrix? [ "Matrix must be square." throw ] unless
augment-identity solution
! now remove the vestigial identity portion of the matrix
flip halves nip flip ;
{
{ -1 -2 3 2 }
{ -4 -1 6 2 }
{ 7 -8 9 1 }
{ 1 -2 1 3 }
} gauss-jordan-invert simple-table.
- Output:
-21/23 17/69 13/138 19/46 -1-15/23 15/23 1/23 15/23 -16/23 25/69 11/138 9/46 -13/23 16/69 -2/69 13/23
Fortran
Note that the Crout algorithm is a more efficient way to invert a matrix.
!-----------------------------------------------------------------------
! gjinv - Invert a matrix, Gauss-Jordan algorithm
! A is destroyed.
!
!___Name_______Type_______________In/Out____Description_________________
! A(LDA,N) Real In An N by N matrix
! LDA Integer In Row bound of A
! N Integer In Order of matrix
! B(LDB,N) Real Out Inverse of A
! LDB Integer In Row bound of B
! IERR Integer Out 0 = no errors
! 1 = singular matrix
!-----------------------------------------------------------------------
SUBROUTINE GJINV (A, LDA, N, B, LDB, IERR)
IMPLICIT NONE
INTEGER LDA, N, LDB, IERR
REAL A(LDA,N), B(LDB,N)
REAL EPS ! machine constant
PARAMETER (EPS = 1.1920929E-07)
INTEGER I, J, K, P ! local variables
REAL F, TOL
!-----------------------------------------------------------------------
! Begin.
!-----------------------------------------------------------------------
IF (N < 1) THEN ! Validate.
IERR = -1
RETURN
ELSE IF (N > LDA .OR. N > LDB) THEN
IERR = -2
RETURN
END IF
IERR = 0
F = 0. ! Frobenius norm of A
DO J = 1, N
DO I = 1, N
F = F + A(I,J)**2
END DO
END DO
F = SQRT(F)
TOL = F * EPS
DO J = 1, N ! Set B to identity matrix.
DO I = 1, N
IF (I .EQ. J) THEN
B(I,J) = 1.
ELSE
B(I,J) = 0.
END IF
END DO
END DO
! Main loop
DO K = 1, N
F = ABS(A(K,K)) ! Find pivot.
P = K
DO I = K+1, N
IF (ABS(A(I,K)) > F) THEN
F = ABS(A(I,K))
P = I
END IF
END DO
IF (F < TOL) THEN ! Matrix is singular.
IERR = 1
RETURN
END IF
IF (P .NE. K) THEN ! Swap rows.
DO J = K, N
F = A(K,J)
A(K,J) = A(P,J)
A(P,J) = F
END DO
DO J = 1, N
F = B(K,J)
B(K,J) = B(P,J)
B(P,J) = F
END DO
END IF
F = 1. / A(K,K) ! Scale row so pivot is 1.
DO J = K, N
A(K,J) = A(K,J) * F
END DO
DO J = 1, N
B(K,J) = B(K,J) * F
END DO
DO 10 I = 1, N ! Subtract to get zeros.
IF (I .EQ. K) GO TO 10
F = A(I,K)
DO J = K, N
A(I,J) = A(I,J) - A(K,J) * F
END DO
DO J = 1, N
B(I,J) = B(I,J) - B(K,J) * F
END DO
10 CONTINUE
END DO
RETURN
END ! of gjinv
Test program:
! Test matrix inversion
PROGRAM TINV
IMPLICIT NONE
INTEGER I, J, K, N, IERR
PARAMETER (N = 4)
REAL AORIG(N,N), A(N,N), B(N,N), C(N,N)
DATA AORIG / -1., -4., 7., 1.,
$ -2., -1., -8., -2.,
$ 3., 6., 9., 1.,
$ 2., 2., 1., 3. /
20 FORMAT (F8.3, $)
30 FORMAT ()
DO J = 1, N
DO I = 1, N
A(I,J) = AORIG(I,J)
END DO
END DO
CALL GJINV (A, N, N, B, N, IERR)
PRINT *, 'gjinv returns #', IERR
PRINT 30
PRINT *, 'matrix:'
DO I = 1, N
DO J = 1, N
PRINT 20, AORIG(I,J)
END DO
PRINT 30
END DO
PRINT 30
PRINT *, 'inverse:'
DO I = 1, N
DO J = 1, N
PRINT 20, B(I,J)
END DO
PRINT 30
END DO
PRINT 30
DO K = 1, N
DO J = 1, N
C(J,K) = 0.
DO I = 1, N
C(J,K) = C(J,K) + AORIG(J,I) * B(I,K)
END DO
END DO
END DO
PRINT *, 'matrix @ inverse:'
DO I = 1, N
DO J = 1, N
PRINT 20, C(I,J)
END DO
PRINT 30
END DO
PRINT 30
CALL GJINV (B, N, N, A, N, IERR)
PRINT *, 'gjinv returns #', IERR
PRINT 30
PRINT *, 'inverse of inverse:'
DO I = 1, N
DO J = 1, N
PRINT 20, A(I,J)
END DO
PRINT 30
END DO
END ! of test program
- Output:
gjinv returns # 0 matrix: -1.000 -2.000 3.000 2.000 -4.000 -1.000 6.000 2.000 7.000 -8.000 9.000 1.000 1.000 -2.000 1.000 3.000 inverse: -0.913 0.246 0.094 0.413 -1.652 0.652 0.043 0.652 -0.696 0.362 0.080 0.196 -0.565 0.232 -0.029 0.565 matrix @ inverse: 1.000 0.000 -0.000 0.000 0.000 1.000 -0.000 0.000 -0.000 -0.000 1.000 0.000 -0.000 -0.000 0.000 1.000 gjinv returns # 0 inverse of inverse: -1.000 -2.000 3.000 2.000 -4.000 -1.000 6.000 2.000 7.000 -8.000 9.000 1.000 1.000 -2.000 1.000 3.000
Generic
// The Generic Language is a database compiler. This code is compiled into database then executed out of database.
generic coordinaat
{
ecs;
uuii;
coordinaat() { ecs=+a;uuii=+a;}
coordinaat(ecs_set,uuii_set) { ecs = ecs_set; uuii=uuii_set;}
operaator<(c)
{
iph ecs < c.ecs return troo;
iph c.ecs < ecs return phals;
iph uuii < c.uuii return troo;
return phals;
}
operaator==(connpair) // eecuuols and not eecuuols deriiu phronn operaator<
{
iph this < connpair return phals;
iph connpair < this return phals;
return troo;
}
operaator!=(connpair)
{
iph this < connpair return troo;
iph connpair < this return troo;
return phals;
}
too_string()
{
return "(" + ecs.too_string() + "," + uuii.too_string() + ")";
}
print()
{
str = too_string();
str.print();
}
println()
{
str = too_string();
str.println();
}
}
generic nnaatrics
{
s; // this is a set of coordinaat/ualioo pairs.
iteraator; // this field holds an iteraator phor the nnaatrics.
nnaatrics() // no parameters required phor nnaatrics construction.
{
s = nioo set(); // create a nioo set of coordinaat/ualioo pairs.
iteraator = nul; // the iteraator is initially set to nul.
}
nnaatrics(copee) // copee the nnaatrics.
{
s = nioo set(); // create a nioo set of coordinaat/ualioo pairs.
iteraator = nul; // the iteraator is initially set to nul.
r = copee.rouus;
c = copee.cols;
i = 0;
uuiil i < r
{
j = 0;
uuiil j < c
{
this[i,j] = copee[i,j];
j++;
}
i++;
}
}
begin { get { return s.begin; } } // property: used to commence manual iteraashon.
end { get { return s.end; } } // property: used to dephiin the end itenn of iteraashon
operaator<(a) // les than operaator is corld bii the avl tree algorithnns
{ // this operaator innpliis phor instance that you could potenshalee hav sets ou nnaatricss.
iph cees < a.cees // connpair the cee sets phurst.
return troo;
els iph a.cees < cees
return phals;
els // the cee sets are eecuuol thairphor connpair nnaatrics elennents.
{
phurst1 = begin;
lahst1 = end;
phurst2 = a.begin;
lahst2 = a.end;
uuiil phurst1 != lahst1 && phurst2 != lahst2
{
iph phurst1.daata.ualioo < phurst2.daata.ualioo
return troo;
els
{
iph phurst2.daata.ualioo < phurst1.daata.ualioo
return phals;
els
{
phurst1 = phurst1.necst;
phurst2 = phurst2.necst;
}
}
}
return phals;
}
}
operaator==(connpair) // eecuuols and not eecuuols deriiu phronn operaator<
{
iph this < connpair return phals;
iph connpair < this return phals;
return troo;
}
operaator!=(connpair)
{
iph this < connpair return troo;
iph connpair < this return troo;
return phals;
}
operaator[](cee_a,cee_b) // this is the nnaatrics indexer.
{
set
{
trii { s >> nioo cee_ualioo(new coordinaat(cee_a,cee_b)); } catch {}
s << nioo cee_ualioo(new coordinaat(nioo integer(cee_a),nioo integer(cee_b)),ualioo);
}
get
{
d = s.get(nioo cee_ualioo(new coordinaat(cee_a,cee_b)));
return d.ualioo;
}
}
operaator>>(coordinaat) // this operaator reennoous an elennent phronn the nnaatrics.
{
s >> nioo cee_ualioo(coordinaat);
return this;
}
iteraat() // and this is how to iterate on the nnaatrics.
{
iph iteraator.nul()
{
iteraator = s.lepht_nnohst;
iph iteraator == s.heder
return nioo iteraator(phals,nioo nun());
els
return nioo iteraator(troo,iteraator.daata.ualioo);
}
els
{
iteraator = iteraator.necst;
iph iteraator == s.heder
{
iteraator = nul;
return nioo iteraator(phals,nioo nun());
}
els
return nioo iteraator(troo,iteraator.daata.ualioo);
}
}
couunt // this property returns a couunt ou elennents in the nnaatrics.
{
get
{
return s.couunt;
}
}
ennptee // is the nnaatrics ennptee?
{
get
{
return s.ennptee;
}
}
lahst // returns the ualioo of the lahst elennent in the nnaatrics.
{
get
{
iph ennptee
throuu "ennptee nnaatrics";
els
return s.lahst.ualioo;
}
}
too_string() // conuerts the nnaatrics too aa string
{
return s.too_string();
}
print() // prints the nnaatrics to the consohl.
{
out = too_string();
out.print();
}
println() // prints the nnaatrics as a liin too the consohl.
{
out = too_string();
out.println();
}
cees // return the set ou cees ou the nnaatrics (a set of coordinaats).
{
get
{
k = nioo set();
phor e : s k << e.cee;
return k;
}
}
operaator+(p)
{
ouut = nioo nnaatrics();
phurst1 = begin;
lahst1 = end;
phurst2 = p.begin;
lahst2 = p.end;
uuiil phurst1 != lahst1 && phurst2 != lahst2
{
ouut[phurst1.daata.cee.ecs,phurst1.daata.cee.uuii] = phurst1.daata.ualioo + phurst2.daata.ualioo;
phurst1 = phurst1.necst;
phurst2 = phurst2.necst;
}
return ouut;
}
operaator-(p)
{
ouut = nioo nnaatrics();
phurst1 = begin;
lahst1 = end;
phurst2 = p.begin;
lahst2 = p.end;
uuiil phurst1 != lahst1 && phurst2 != lahst2
{
ouut[phurst1.daata.cee.ecs,phurst1.daata.cee.uuii] = phurst1.daata.ualioo - phurst2.daata.ualioo;
phurst1 = phurst1.necst;
phurst2 = phurst2.necst;
}
return ouut;
}
rouus
{
get
{
r = +a;
phurst1 = begin;
lahst1 = end;
uuiil phurst1 != lahst1
{
iph r < phurst1.daata.cee.ecs r = phurst1.daata.cee.ecs;
phurst1 = phurst1.necst;
}
return r + +b;
}
}
cols
{
get
{
c = +a;
phurst1 = begin;
lahst1 = end;
uuiil phurst1 != lahst1
{
iph c < phurst1.daata.cee.uuii c = phurst1.daata.cee.uuii;
phurst1 = phurst1.necst;
}
return c + +b;
}
}
operaator*(o)
{
iph cols != o.rouus throw "rouus-cols nnisnnatch";
reesult = nioo nnaatrics();
rouu_couunt = rouus;
colunn_couunt = o.cols;
loop = cols;
i = +a;
uuiil i < rouu_couunt
{
g = +a;
uuiil g < colunn_couunt
{
sunn = +a.a;
h = +a;
uuiil h < loop
{
a = this[i, h];
b = o[h, g];
nn = a * b;
sunn = sunn + nn;
h++;
}
reesult[i, g] = sunn;
g++;
}
i++;
}
return reesult;
}
suuop_rouus(a, b)
{
c = cols;
i = 0;
uuiil u < cols
{
suuop = this[a, i];
this[a, i] = this[b, i];
this[b, i] = suuop;
i++;
}
}
suuop_colunns(a, b)
{
r = rouus;
i = 0;
uuiil i < rouus
{
suuopp = this[i, a];
this[i, a] = this[i, b];
this[i, b] = suuop;
i++;
}
}
transpohs
{
get
{
reesult = new nnaatrics();
r = rouus;
c = cols;
i=0;
uuiil i < r
{
g = 0;
uuiil g < c
{
reesult[g, i] = this[i, g];
g++;
}
i++;
}
return reesult;
}
}
deternninant
{
get
{
rouu_couunt = rouus;
colunn_count = cols;
if rouu_couunt != colunn_count
throw "not a scuuair nnaatrics";
if rouu_couunt == 0
throw "the nnaatrics is ennptee";
if rouu_couunt == 1
return this[0, 0];
if rouu_couunt == 2
return this[0, 0] * this[1, 1] -
this[0, 1] * this[1, 0];
temp = nioo nnaatrics();
det = 0.0;
parity = 1.0;
j = 0;
uuiil j < rouu_couunt
{
k = 0;
uuiil k < rouu_couunt-1
{
skip_col = phals;
l = 0;
uuiil l < rouu_couunt-1
{
if l == j skip_col = troo;
if skip_col
n = l + 1;
els
n = l;
temp[k, l] = this[k + 1, n];
l++;
}
k++;
}
det = det + parity * this[0, j] * temp.deeternninant;
parity = 0.0 - parity;
j++;
}
return det;
}
}
ad_rouu(a, b)
{
c = cols;
i = 0;
uuiil i < c
{
this[a, i] = this[a, i] + this[b, i];
i++;
}
}
ad_colunn(a, b)
{
c = rouus;
i = 0;
uuiil i < c
{
this[i, a] = this[i, a] + this[i, b];
i++;
}
}
subtract_rouu(a, b)
{
c = cols;
i = 0;
uuiil i < c
{
this[a, i] = this[a, i] - this[b, i];
i++;
}
}
subtract_colunn(a, b)
{
c = rouus;
i = 0;
uuiil i < c
{
this[i, a] = this[i, a] - this[i, b];
i++;
}
}
nnultiplii_rouu(rouu, scalar)
{
c = cols;
i = 0;
uuiil i < c
{
this[rouu, i] = this[rouu, i] * scalar;
i++;
}
}
nnultiplii_colunn(colunn, scalar)
{
r = rouus;
i = 0;
uuiil i < r
{
this[i, colunn] = this[i, colunn] * scalar;
i++;
}
}
diuiid_rouu(rouu, scalar)
{
c = cols;
i = 0;
uuiil i < c
{
this[rouu, i] = this[rouu, i] / scalar;
i++;
}
}
diuiid_colunn(colunn, scalar)
{
r = rouus;
i = 0;
uuiil i < r
{
this[i, colunn] = this[i, colunn] / scalar;
i++;
}
}
connbiin_rouus_ad(a,b,phactor)
{
c = cols;
i = 0;
uuiil i < c
{
this[a, i] = this[a, i] + phactor * this[b, i];
i++;
}
}
connbiin_rouus_subtract(a,b,phactor)
{
c = cols;
i = 0;
uuiil i < c
{
this[a, i] = this[a, i] - phactor * this[b, i];
i++;
}
}
connbiin_colunns_ad(a,b,phactor)
{
r = rouus;
i = 0;
uuiil i < r
{
this[i, a] = this[i, a] + phactor * this[i, b];
i++;
}
}
connbiin_colunns_subtract(a,b,phactor)
{
r = rouus;
i = 0;
uuiil i < r
{
this[i, a] = this[i, a] - phactor * this[i, b];
i++;
}
}
inuers
{
get
{
rouu_couunt = rouus;
colunn_couunt = cols;
iph rouu_couunt != colunn_couunt
throw "nnatrics not scuuair";
els iph rouu_couunt == 0
throw "ennptee nnatrics";
els iph rouu_couunt == 1
{
r = nioo nnaatrics();
r[0, 0] = 1.0 / this[0, 0];
return r;
}
gauss = nioo nnaatrics(this);
i = 0;
uuiil i < rouu_couunt
{
j = 0;
uuiil j < rouu_couunt
{
iph i == j
gauss[i, j + rouu_couunt] = 1.0;
els
gauss[i, j + rouu_couunt] = 0.0;
j++;
}
i++;
}
j = 0;
uuiil j < rouu_couunt
{
iph gauss[j, j] == 0.0
{
k = j + 1;
uuiil k < rouu_couunt
{
if gauss[k, j] != 0.0 {gauss.nnaat.suuop_rouus(j, k); break; }
k++;
}
if k == rouu_couunt throw "nnatrics is singioolar";
}
phactor = gauss[j, j];
iph phactor != 1.0 gauss.diuiid_rouu(j, phactor);
i = j+1;
uuiil i < rouu_couunt
{
gauss.connbiin_rouus_subtract(i, j, gauss[i, j]);
i++;
}
j++;
}
i = rouu_couunt - 1;
uuiil i > 0
{
k = i - 1;
uuiil k >= 0
{
gauss.connbiin_rouus_subtract(k, i, gauss[k, i]);
k--;
}
i--;
}
reesult = nioo nnaatrics();
i = 0;
uuiil i < rouu_couunt
{
j = 0;
uuiil j < rouu_couunt
{
reesult[i, j] = gauss[i, j + rouu_couunt];
j++;
}
i++;
}
return reesult;
}
}
}
Go
package main
import "fmt"
type vector = []float64
type matrix []vector
func (m matrix) inverse() matrix {
le := len(m)
for _, v := range m {
if len(v) != le {
panic("Not a square matrix")
}
}
aug := make(matrix, le)
for i := 0; i < le; i++ {
aug[i] = make(vector, 2*le)
copy(aug[i], m[i])
// augment by identity matrix to right
aug[i][i+le] = 1
}
aug.toReducedRowEchelonForm()
inv := make(matrix, le)
// remove identity matrix to left
for i := 0; i < le; i++ {
inv[i] = make(vector, le)
copy(inv[i], aug[i][le:])
}
return inv
}
// note: this mutates the matrix in place
func (m matrix) toReducedRowEchelonForm() {
lead := 0
rowCount, colCount := len(m), len(m[0])
for r := 0; r < rowCount; r++ {
if colCount <= lead {
return
}
i := r
for m[i][lead] == 0 {
i++
if rowCount == i {
i = r
lead++
if colCount == lead {
return
}
}
}
m[i], m[r] = m[r], m[i]
if div := m[r][lead]; div != 0 {
for j := 0; j < colCount; j++ {
m[r][j] /= div
}
}
for k := 0; k < rowCount; k++ {
if k != r {
mult := m[k][lead]
for j := 0; j < colCount; j++ {
m[k][j] -= m[r][j] * mult
}
}
}
lead++
}
}
func (m matrix) print(title string) {
fmt.Println(title)
for _, v := range m {
fmt.Printf("% f\n", v)
}
fmt.Println()
}
func main() {
a := matrix{{1, 2, 3}, {4, 1, 6}, {7, 8, 9}}
a.inverse().print("Inverse of A is:\n")
b := matrix{{2, -1, 0}, {-1, 2, -1}, {0, -1, 2}}
b.inverse().print("Inverse of B is:\n")
}
- Output:
Inverse of A is: [-0.812500 0.125000 0.187500] [ 0.125000 -0.250000 0.125000] [ 0.520833 0.125000 -0.145833] Inverse of B is: [ 0.750000 0.500000 0.250000] [ 0.500000 1.000000 0.500000] [ 0.250000 0.500000 0.750000]
Alternative
package main
import (
"fmt"
"math"
)
type vector = []float64
type matrix []vector
func (m matrix) print(title string) {
fmt.Println(title)
for _, v := range m {
fmt.Printf("% f\n", v)
}
}
func I(n int) matrix {
b := make(matrix, n)
for i := 0; i < n; i++ {
b[i] = make(vector, n)
b[i][i] = 1
}
return b
}
func (m matrix) inverse() matrix {
n := len(m)
for _, v := range m {
if n != len(v) {
panic("Not a square matrix")
}
}
b := I(n)
a := make(matrix, n)
for i, v := range m {
a[i] = make(vector, n)
copy(a[i], v)
}
for k := range a {
iMax := 0
max := -1.
for i := k; i < n; i++ {
row := a[i]
// compute scale factor s = max abs in row
s := -1.
for j := k; j < n; j++ {
x := math.Abs(row[j])
if x > s {
s = x
}
}
if s == 0 {
panic("Irregular matrix")
}
// scale the abs used to pick the pivot.
if abs := math.Abs(row[k]) / s; abs > max {
iMax = i
max = abs
}
}
if k != iMax {
a[k], a[iMax] = a[iMax], a[k]
b[k], b[iMax] = b[iMax], b[k]
}
akk := a[k][k]
for j := 0; j < n; j++ {
a[k][j] /= akk
b[k][j] /= akk
}
for i := 0; i < n; i++ {
if (i != k) {
aik := a[i][k]
for j := 0; j < n; j++ {
a[i][j] -= a[k][j] * aik
b[i][j] -= b[k][j] * aik
}
}
}
}
return b
}
func main() {
a := matrix{{1, 2, 3}, {4, 1, 6}, {7, 8, 9}}
a.inverse().print("Inverse of A is:\n")
b := matrix{{2, -1, 0}, {-1, 2, -1}, {0, -1, 2}}
b.inverse().print("Inverse of B is:\n")
}
- Output:
Same output than above
Haskell
isMatrix xs = null xs || all ((== (length.head $ xs)).length) xs
isSquareMatrix xs = null xs || all ((== (length xs)).length) xs
mult:: Num a => [[a]] -> [[a]] -> [[a]]
mult uss vss = map ((\xs -> if null xs then [] else foldl1 (zipWith (+)) xs). zipWith (\vs u -> map (u*) vs) vss) uss
matI::(Num a) => Int -> [[a]]
matI n = [ [fromIntegral.fromEnum $ i == j | j <- [1..n]] | i <- [1..n]]
inversion xs = gauss xs (matI $ length xs)
gauss::[[Double]] -> [[Double]] -> [[Double]]
gauss xs bs = map (map fromRational) $ solveGauss (toR xs) (toR bs)
where toR = map $ map toRational
solveGauss:: (Fractional a, Ord a) => [[a]] -> [[a]] -> [[a]]
solveGauss xs bs | null xs || null bs || length xs /= length bs || (not $ isSquareMatrix xs) || (not $ isMatrix bs) = []
| otherwise = uncurry solveTriangle $ triangle xs bs
solveTriangle::(Fractional a,Eq a) => [[a]] -> [[a]] -> [[a]]
solveTriangle us _ | not.null.dropWhile ((/= 0).head) $ us = []
solveTriangle ([c]:as) (b:bs) = go as bs [map (/c) b]
where
val us vs ws = let u = head us in map (/u) $ zipWith (-) vs (head $ mult [tail us] ws)
go [] _ zs = zs
go _ [] zs = zs
go (x:xs) (y:ys) zs = go xs ys $ (val x y zs):zs
triangle::(Num a, Ord a) => [[a]] -> [[a]] -> ([[a]],[[a]])
triangle xs bs = triang ([],[]) (xs,bs)
where
triang ts (_,[]) = ts
triang ts ([],_) = ts
triang (os,ps) zs = triang (us:os,cs:ps).unzip $ [(fun tus vs, fun cs es) | (v:vs,es) <- zip uss css,let fun = zipWith (\x y -> v*x - u*y)]
where ((us@(u:tus)):uss,cs:css) = bubble zs
bubble::(Num a, Ord a) => ([[a]],[[a]]) -> ([[a]],[[a]])
bubble (xs,bs) = (go xs, go bs)
where
idmax = snd.maximum.flip zip [0..].map (abs.head) $ xs
go ys = let (us,vs) = splitAt idmax ys in vs ++ us
main = do
let a = [[1, 2, 3], [4, 1, 6], [7, 8, 9]]
let b = [[2, -1, 0], [-1, 2, -1], [0, -1, 2]]
putStrLn "inversion a ="
mapM_ print $ inversion a
putStrLn "\ninversion b ="
mapM_ print $ inversion b
- Output:
inversion a = [-0.8125,0.125,0.1875] [0.125,-0.25,0.125] [0.5208333333333334,0.125,-0.14583333333333334] inversion b = [0.75,0.5,0.25] [0.5,1.0,0.5] [0.25,0.5,0.75]
J
Solution:
Uses Gauss-Jordan implementation (as described in Reduced_row_echelon_form#J) to find reduced row echelon form of the matrix after augmenting with an identity matrix.
require 'math/misc/linear'
augmentR_I1=: ,. e.@i.@# NB. augment matrix on the right with its Identity matrix
matrix_invGJ=: # }."1 [: gauss_jordan@augmentR_I1
Usage:
]A =: 1 2 3, 4 1 6,: 7 8 9
1 2 3
4 1 6
7 8 9
matrix_invGJ A
_0.8125 0.125 0.1875
0.125 _0.25 0.125
0.520833 0.125 _0.145833
Java
// GaussJordan.java
import java.util.Random;
public class GaussJordan {
public static void main(String[] args) {
int rows = 5;
Matrix m = new Matrix(rows, rows);
Random r = new Random();
for (int row = 0; row < rows; ++row) {
for (int column = 0; column < rows; ++column)
m.set(row, column, r.nextDouble());
}
System.out.println("Matrix:");
m.print();
System.out.println("Inverse:");
Matrix inv = m.inverse();
inv.print();
System.out.println("Product of matrix and inverse:");
Matrix.product(m, inv).print();
}
}
// Matrix.java
public class Matrix {
private int rows;
private int columns;
private double[][] elements;
public Matrix(int rows, int columns) {
this.rows = rows;
this.columns = columns;
elements = new double[rows][columns];
}
public void set(int row, int column, double value) {
elements[row][column] = value;
}
public double get(int row, int column) {
return elements[row][column];
}
public void print() {
for (int row = 0; row < rows; ++row) {
for (int column = 0; column < columns; ++column) {
if (column > 0)
System.out.print(' ');
System.out.printf("%7.3f", elements[row][column]);
}
System.out.println();
}
}
// Returns the inverse of this matrix
public Matrix inverse() {
assert(rows == columns);
// Augment by identity matrix
Matrix tmp = new Matrix(rows, columns * 2);
for (int row = 0; row < rows; ++row) {
System.arraycopy(elements[row], 0, tmp.elements[row], 0, columns);
tmp.elements[row][row + columns] = 1;
}
tmp.toReducedRowEchelonForm();
Matrix inv = new Matrix(rows, columns);
for (int row = 0; row < rows; ++row)
System.arraycopy(tmp.elements[row], columns, inv.elements[row], 0, columns);
return inv;
}
// Converts this matrix into reduced row echelon form
public void toReducedRowEchelonForm() {
for (int row = 0, lead = 0; row < rows && lead < columns; ++row, ++lead) {
int i = row;
while (elements[i][lead] == 0) {
if (++i == rows) {
i = row;
if (++lead == columns)
return;
}
}
swapRows(i, row);
if (elements[row][lead] != 0) {
double f = elements[row][lead];
for (int column = 0; column < columns; ++column)
elements[row][column] /= f;
}
for (int j = 0; j < rows; ++j) {
if (j == row)
continue;
double f = elements[j][lead];
for (int column = 0; column < columns; ++column)
elements[j][column] -= f * elements[row][column];
}
}
}
// Returns the matrix product of a and b
public static Matrix product(Matrix a, Matrix b) {
assert(a.columns == b.rows);
Matrix result = new Matrix(a.rows, b.columns);
for (int i = 0; i < a.rows; ++i) {
double[] resultRow = result.elements[i];
double[] aRow = a.elements[i];
for (int j = 0; j < a.columns; ++j) {
double[] bRow = b.elements[j];
for (int k = 0; k < b.columns; ++k)
resultRow[k] += aRow[j] * bRow[k];
}
}
return result;
}
private void swapRows(int i, int j) {
double[] tmp = elements[i];
elements[i] = elements[j];
elements[j] = tmp;
}
}
- Output:
Matrix: 0.913 0.438 0.002 0.053 0.588 0.739 0.043 0.593 0.599 0.115 0.542 0.734 0.273 0.889 0.921 0.503 0.960 0.707 0.088 0.921 0.777 0.829 0.190 0.673 0.728 Inverse: 0.873 0.502 -0.705 -0.272 0.452 -1.361 -0.578 -2.118 0.400 3.368 -0.539 0.883 -0.106 0.910 -0.722 -0.720 0.291 0.625 -0.628 0.539 1.425 -0.377 2.616 0.178 -3.257 Product of matrix and inverse: 1.000 0.000 0.000 -0.000 0.000 -0.000 1.000 -0.000 -0.000 0.000 0.000 0.000 1.000 -0.000 -0.000 -0.000 0.000 0.000 1.000 -0.000 -0.000 0.000 0.000 -0.000 1.000
jq
Adapted from Wren-matrix
Works with gojq, the Go implementation of jq
# Create an m x n matrix,
# it being understood that:
# matrix(0; _; _) evaluates to []
def matrix(m; n; init):
if m == 0 then []
elif m == 1 then [[range(0;n) | init]]
elif m > 0 then
[range(0;n) | init] as $row
| [range(0;m) | $row ]
else error("matrix\(m);_;_) invalid")
end;
def mprint($dec):
def max(s): reduce s as $x (null; if . == null or $x > . then $x else . end);
def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;
pow(10; $dec) as $power
| def p: (. * $power | round) / $power;
(max(.[][] | p | tostring | length) + 1) as $w
| . as $in
| range(0; length) as $i
| reduce range(0; .[$i]|length) as $j ("|"; . + ($in[$i][$j]|p|lpad($w)))
| . + " |" ;
# Swaps two rows of the input matrix using IO==0
def swapRows(rowNum1; rowNum2):
if (rowNum1 == rowNum2)
then .
else .[rowNum1] as $t
| .[rowNum1] = .[rowNum2]
| .[rowNum2] = $t
end;
def toReducedRowEchelonForm:
. as $in
| length as $nr
| (.[0]|length) as $nc
| {a: $in, lead: 0 }
| label $out
| last(foreach range(0; $nr) as $r (.;
if $nc <= .lead then ., break $out else . end
| .i = $r
| until( .a[.i][.lead] != 0 or .lead == $nc - 1;
.i += 1
| if $nr == .i
then .i = $r
| .lead += 1
else .
end
)
| swapRows(.i; $r)
| if .a[$r][.lead] != 0
then .a[$r][.lead] as $div
| reduce range(0; $nc) as $j (.; .a[$r][$j] /= $div)
else .
end
| reduce range(0; $nr) as $k (.;
if $k != $r
then .a[$k][.lead] as $mult
| reduce range(0; $nc) as $j (.; .a[$k][$j] -= .a[$r][$j] * $mult)
else .
end)
| .lead += 1
))
| .a ;
# Assumption: the input is a square matrix with an inverse
# Uses the Gauss-Jordan method.
def inverse:
. as $a
| length as $nr
| reduce range(0; $nr) as $i (
matrix($nr; 2 * $nr; 0);
reduce range( 0; $nr) as $j (.;
.[$i][$j] = $a[$i][$j]
| .[$i][$i + $nr] = 1 ))
| last(toReducedRowEchelonForm)
| . as $ary
| reduce range(0; $nr) as $i ( matrix($nr; $nr; 0);
reduce range($nr; 2 *$nr) as $j (.;
.[$i][$j - $nr] = $ary[$i][$j] )) ;The Task
def arrays: [
[ [ 1, 2, 3],
[ 4, 1, 6],
[ 7, 8, 9] ],
[ [ 2, -1, 0 ],
[-1, 2, -1 ],
[ 0, -1, 2 ] ],
[ [ -1, -2, 3, 2 ],
[ -4, -1, 6, 2 ],
[ 7, -8, 9, 1 ],
[ 1, -2, 1, 3 ] ]
];
def task:
arrays[]
| "Original:",
mprint(1),
"\nInverse:",
(inverse|mprint(6)),
""
;
task- Output:
Original: | 1 2 3 | | 4 1 6 | | 7 8 9 | Inverse: | -0.8125 0.125 0.1875 | | 0.125 -0.25 0.125 | | 0.520833 0.125 -0.145833 | Original: | 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 | Original: | -1 -2 3 2 | | -4 -1 6 2 | | 7 -8 9 1 | | 1 -2 1 3 | Inverse: | -0.913043 0.246377 0.094203 0.413043 | | -1.652174 0.652174 0.043478 0.652174 | | -0.695652 0.362319 0.07971 0.195652 | | -0.565217 0.231884 -0.028986 0.565217 |
Julia
Built-in LAPACK-based linear solver uses partial-pivoted Gauss elimination):
A = [1 2 3; 4 1 6; 7 8 9]
@show I / A
@show inv(A)
Native implementation:
function gaussjordan(A::Matrix)
size(A, 1) == size(A, 2) || throw(ArgumentError("A must be squared"))
n = size(A, 1)
M = [convert(Matrix{float(eltype(A))}, A) I]
i = 1
local tmp = Vector{eltype(M)}(2n)
# forward
while i ≤ n
if M[i, i] ≈ 0.0
local j = i + 1
while j ≤ n && M[j, i] ≈ 0.0
j += 1
end
if j ≤ n
tmp .= M[i, :]
M[i, :] .= M[j, :]
M[j, :] .= tmp
else
throw(ArgumentError("matrix is singular, cannot compute the inverse"))
end
end
for j in (i + 1):n
M[j, :] .-= M[j, i] / M[i, i] .* M[i, :]
end
i += 1
end
i = n
# backward
while i ≥ 1
if M[i, i] ≈ 0.0
local j = i - 1
while j ≥ 1 && M[j, i] ≈ 0.0
j -= 1
end
if j ≥ 1
tmp .= M[i, :]
M[i, :] .= M[j, :]
M[j, :] .= tmp
else
throw(ArgumentError("matrix is singular, cannot compute the inverse"))
end
end
for j in (i - 1):-1:1
M[j, :] .-= M[j, i] / M[i, i] .* M[i, :]
end
i -= 1
end
M ./= diag(M) # normalize
return M[:, n+1:2n]
end
@show gaussjordan(A)
@assert gaussjordan(A) ≈ inv(A)
A = rand(10, 10)
@assert gaussjordan(A) ≈ inv(A)
- Output:
I / A = [-0.8125 0.125 0.1875; 0.125 -0.25 0.125; 0.520833 0.125 -0.145833] inv(A) = [-0.8125 0.125 0.1875; 0.125 -0.25 0.125; 0.520833 0.125 -0.145833] gaussjordan(A) = [-0.8125 0.125 0.1875; 0.125 -0.25 0.125; 0.520833 0.125 -0.145833]
Kotlin
This follows the description of Gauss-Jordan elimination in Wikipedia whereby the original square matrix is first augmented to the right by its identity matrix, its reduced row echelon form is then found and finally the identity matrix to the left is removed to leave the inverse of the original square matrix.
// version 1.2.21
typealias Matrix = Array<DoubleArray>
fun Matrix.inverse(): Matrix {
val len = this.size
require(this.all { it.size == len }) { "Not a square matrix" }
val aug = Array(len) { DoubleArray(2 * len) }
for (i in 0 until len) {
for (j in 0 until len) aug[i][j] = this[i][j]
// augment by identity matrix to right
aug[i][i + len] = 1.0
}
aug.toReducedRowEchelonForm()
val inv = Array(len) { DoubleArray(len) }
// remove identity matrix to left
for (i in 0 until len) {
for (j in len until 2 * len) inv[i][j - len] = aug[i][j]
}
return inv
}
fun Matrix.toReducedRowEchelonForm() {
var lead = 0
val rowCount = this.size
val colCount = this[0].size
for (r in 0 until rowCount) {
if (colCount <= lead) return
var i = r
while (this[i][lead] == 0.0) {
i++
if (rowCount == i) {
i = r
lead++
if (colCount == lead) return
}
}
val temp = this[i]
this[i] = this[r]
this[r] = temp
if (this[r][lead] != 0.0) {
val div = this[r][lead]
for (j in 0 until colCount) this[r][j] /= div
}
for (k in 0 until rowCount) {
if (k != r) {
val mult = this[k][lead]
for (j in 0 until colCount) this[k][j] -= this[r][j] * mult
}
}
lead++
}
}
fun Matrix.printf(title: String) {
println(title)
val rowCount = this.size
val colCount = this[0].size
for (r in 0 until rowCount) {
for (c in 0 until colCount) {
if (this[r][c] == -0.0) this[r][c] = 0.0 // get rid of negative zeros
print("${"% 10.6f".format(this[r][c])} ")
}
println()
}
println()
}
fun main(args: Array<String>) {
val a = arrayOf(
doubleArrayOf(1.0, 2.0, 3.0),
doubleArrayOf(4.0, 1.0, 6.0),
doubleArrayOf(7.0, 8.0, 9.0)
)
a.inverse().printf("Inverse of A is :\n")
val b = arrayOf(
doubleArrayOf( 2.0, -1.0, 0.0),
doubleArrayOf(-1.0, 2.0, -1.0),
doubleArrayOf( 0.0, -1.0, 2.0)
)
b.inverse().printf("Inverse of B is :\n")
}
- Output:
Inverse of A is : -0.812500 0.125000 0.187500 0.125000 -0.250000 0.125000 0.520833 0.125000 -0.145833 Inverse of B is : 0.750000 0.500000 0.250000 0.500000 1.000000 0.500000 0.250000 0.500000 0.750000
Lambdatalk
{require lib_matrix}
{M.gaussJordan
{M.new [[1,2,3],
[4,5,6],
[7,8,-9]]}}
->
[[-1.722222222222222,0.7777777777777777,-0.05555555555555555],
[1.4444444444444444,-0.5555555555555556,0.1111111111111111],
[-0.05555555555555555,0.1111111111111111,-0.05555555555555555]]
M2000 Interpreter
We can use supported types like Decimal, Currency, Double, Single. Also matrix may have different base on each dimension.
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"- Output:
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
Mathematica / Wolfram Language
a = N@{{1, 7, 5, 0}, {5, 8, 6, 9}, {2, 1, 6, 4}, {8, 1, 2, 4}};
elimmat = RowReduce[Join[a, IdentityMatrix[Length[a]], 2]];
inv = elimmat[[All, -Length[a] ;;]]
- Output:
{{0.0463243, -0.0563948, -0.0367573, 0.163646}, {0.0866062, 0.0684794, -0.133938, -0.020141}, {0.0694864, -0.0845921, 0.194864, -0.00453172}, {-0.149043, 0.137966, 0.00956697, -0.0699899}}
MATLAB
function GaussInverse(M)
original = M;
[n,~] = size(M);
I = eye(n);
for j = 1:n
for i = j:n
if ~(M(i,j) == 0)
for k = 1:n
q = M(j,k); M(j,k) = M(i,k); M(i,k) = q;
q = I(j,k); I(j,k) = I(i,k); I(i,k) = q;
end
p = 1/M(j,j);
for k = 1:n
M(j,k) = p*M(j,k);
I(j,k) = p*I(j,k);
end
for L = 1:n
if ~(L == j)
p = -M(L,j);
for k = 1:n
M(L,k) = M(L,k) + p*M(j,k);
I(L,k) = I(L,k) + p*I(j,k);
end
end
end
end
end
inverted = I;
end
fprintf("Matrix:\n")
disp(original)
fprintf("Inverted matrix:\n")
disp(inverted)
fprintf("Inverted matrix calculated by built-in function:\n")
disp(inv(original))
fprintf("Product of matrix and inverse:\n")
disp(original*inverted)
end
A = [ 2, -1, 0;
-1, 2, -1;
0, -1, 2 ];
GaussInverse(A)
- Output:
Matrix:
2 -1 0
-1 2 -1
0 -1 2
Inverted matrix:
0.7500 0.5000 0.2500
0.5000 1.0000 0.5000
0.2500 0.5000 0.7500
Inverted matrix calculated by built-in function:
0.7500 0.5000 0.2500
0.5000 1.0000 0.5000
0.2500 0.5000 0.7500
Product of matrix and inverse:
1.0000 0 0
-0.0000 1.0000 -0.0000
-0.0000 -0.0000 1.0000
Nim
We reuse the algorithm of task https://rosettacode.org/wiki/Reduced_row_echelon_form (version using floats).
import strformat, strutils
const Eps = 1e-10
type
Matrix[M, N: static Positive] = array[M, array[N, float]]
SquareMatrix[N: static Positive] = Matrix[N, N]
func toSquareMatrix[N: static Positive](a: array[N, array[N, int]]): SquareMatrix[N] =
## Convert a square matrix of integers to a square matrix of floats.
for i in 0..<N:
for j in 0..<N:
result[i][j] = a[i][j].toFloat
func transformToRref(mat: var Matrix) =
## Transform a matrix to reduced row echelon form.
var lead = 0
for r in 0..<mat.M:
if lead >= mat.N: return
var i = r
while mat[i][lead] == 0:
inc i
if i == mat.M:
i = r
inc lead
if lead == mat.N: return
swap mat[i], mat[r]
let d = mat[r][lead]
if abs(d) > Eps: # Checking "d != 0" will give wrong results in some cases.
for item in mat[r].mitems:
item /= d
for i in 0..<mat.M:
if i != r:
let m = mat[i][lead]
for c in 0..<mat.N:
mat[i][c] -= mat[r][c] * m
inc lead
func inverse(mat: SquareMatrix): SquareMatrix[mat.N] =
## Return the inverse of a matrix.
# Build augmented matrix.
var augmat: Matrix[mat.N, 2 * mat.N]
for i in 0..<mat.N:
augmat[i][0..<mat.N] = mat[i]
augmat[i][mat.N + i] = 1
# Transform it to reduced row echelon form.
augmat.transformToRref()
# Check if the first half is the identity matrix and extract second half.
for i in 0..<mat.N:
for j in 0..<mat.N:
if augmat[i][j] != float(i == j):
raise newException(ValueError, "matrix is singular")
result[i][j] = augmat[i][mat.N + j]
proc `$`(mat: Matrix): string =
## Display a matrix (which may be a square matrix).
for row in mat:
var line = ""
for val in row:
line.addSep(" ", 0)
line.add &"{val:9.5f}"
echo line
#———————————————————————————————————————————————————————————————————————————————————————————————————
template runTest(mat: SquareMatrix) =
## Run a test using square matrix "mat".
echo "Matrix:"
echo $mat
echo "Inverse:"
echo mat.inverse
echo ""
let m1 = [[1, 2, 3],
[4, 1, 6],
[7, 8, 9]].toSquareMatrix()
let m2 = [[ 2, -1, 0],
[-1, 2, -1],
[ 0, -1, 2]].toSquareMatrix()
let m3 = [[ -1, -2, 3, 2],
[ -4, -1, 6, 2],
[ 7, -8, 9, 1],
[ 1, -2, 1, 3]].toSquareMatrix()
runTest(m1)
runTest(m2)
runTest(m3)
- Output:
Matrix: 1.00000 2.00000 3.00000 4.00000 1.00000 6.00000 7.00000 8.00000 9.00000 Inverse: -0.81250 0.12500 0.18750 0.12500 -0.25000 0.12500 0.52083 0.12500 -0.14583 Matrix: 2.00000 -1.00000 0.00000 -1.00000 2.00000 -1.00000 0.00000 -1.00000 2.00000 Inverse: 0.75000 0.50000 0.25000 0.50000 1.00000 0.50000 0.25000 0.50000 0.75000 Matrix: -1.00000 -2.00000 3.00000 2.00000 -4.00000 -1.00000 6.00000 2.00000 7.00000 -8.00000 9.00000 1.00000 1.00000 -2.00000 1.00000 3.00000 Inverse: -0.91304 0.24638 0.09420 0.41304 -1.65217 0.65217 0.04348 0.65217 -0.69565 0.36232 0.07971 0.19565 -0.56522 0.23188 -0.02899 0.56522
Pascal
Free 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.
- Output:
== 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
PascalABC.NET
uses NumLibABC;
begin
var A := new Matrix(3,3,1,2,3,4,1,6,7,8,9);
'Original:'.Println;
A.Println(10,5);
'Inverse:'.Println;
A.Inv.Println(10,5);
end.
- Output:
Original: 1.00000 2.00000 3.00000 4.00000 1.00000 6.00000 7.00000 8.00000 9.00000 Inverse: -0.81250 0.12500 0.18750 0.12500 -0.25000 0.12500 0.52083 0.12500 -0.14583
Perl
Included code from Reduced row echelon form task.
sub rref {
our @m; local *m = shift;
@m or return;
my ($lead, $rows, $cols) = (0, scalar(@m), scalar(@{$m[0]}));
foreach my $r (0 .. $rows - 1) {
$lead < $cols or return;
my $i = $r;
until ($m[$i][$lead])
{++$i == $rows or next;
$i = $r;
++$lead == $cols and return;}
@m[$i, $r] = @m[$r, $i];
my $lv = $m[$r][$lead];
$_ /= $lv foreach @{ $m[$r] };
my @mr = @{ $m[$r] };
foreach my $i (0 .. $rows - 1)
{$i == $r and next;
($lv, my $n) = ($m[$i][$lead], -1);
$_ -= $lv * $mr[++$n] foreach @{ $m[$i] };}
++$lead;}
}
sub display { join("\n" => map join(" " => map(sprintf("%6.2f", $_), @$_)), @{+shift})."\n" }
sub gauss_jordan_invert {
my(@m) = @_;
my $rows = @m;
my @i = identity(scalar @m);
push @{$m[$_]}, @{$i[$_]} for 0..$rows-1;
rref(\@m);
map { splice @$_, 0, $rows } @m;
@m;
}
sub identity {
my($n) = @_;
map { [ (0) x $_, 1, (0) x ($n-1 - $_) ] } 0..$n-1
}
my @tests = (
[
[ 2, -1, 0 ],
[-1, 2, -1 ],
[ 0, -1, 2 ]
],
[
[ -1, -2, 3, 2 ],
[ -4, -1, 6, 2 ],
[ 7, -8, 9, 1 ],
[ 1, -2, 1, 3 ]
],
);
for my $matrix (@tests) {
print "Original Matrix:\n" . display(\@$matrix) . "\n";
my @gj = gauss_jordan_invert( @$matrix );
print "Gauss-Jordan Inverted Matrix:\n" . display(\@gj) . "\n";
my @rt = gauss_jordan_invert( @gj );
print "After round-trip:\n" . display(\@rt) . "\n";} . "\n"
}
- Output:
Original Matrix: 2.00 -1.00 0.00 -1.00 2.00 -1.00 0.00 -1.00 2.00 Gauss-Jordan Inverted Matrix: 0.75 0.50 0.25 0.50 1.00 0.50 0.25 0.50 0.75 After round-trip: 2.00 -1.00 0.00 -1.00 2.00 -1.00 0.00 -1.00 2.00 Original Matrix: -1.00 -2.00 3.00 2.00 -4.00 -1.00 6.00 2.00 7.00 -8.00 9.00 1.00 1.00 -2.00 1.00 3.00 Gauss-Jordan Inverted Matrix: -0.91 0.25 0.09 0.41 -1.65 0.65 0.04 0.65 -0.70 0.36 0.08 0.20 -0.57 0.23 -0.03 0.57 After round-trip: -1.00 -2.00 3.00 2.00 -4.00 -1.00 6.00 2.00 7.00 -8.00 9.00 1.00 1.00 -2.00 1.00 3.00
Phix
uses ToReducedRowEchelonForm() from Reduced_row_echelon_form#Phix
with javascript_semantics function ToReducedRowEchelonForm(sequence M) integer lead = 1, rowCount = length(M), columnCount = length(M[1]), i for r=1 to rowCount do if lead>=columnCount then exit end if i = r while M[i][lead]=0 do i += 1 if i=rowCount then i = r lead += 1 if lead=columnCount then exit end if end if end while object mr = sq_div(M[i],M[i][lead]) M[i] = M[r] M[r] = mr for j=1 to rowCount do if j!=r then M[j] = sq_sub(M[j],sq_mul(M[j][lead],M[r])) end if end for lead += 1 end for return M end function --</end of copy> function inverse(sequence mat) integer len = length(mat) sequence aug = repeat(repeat(0,2*len),len) for i=1 to len do if length(mat[i])!=len then ?9/0 end if -- "Not a square matrix" for j=1 to len do aug[i][j] = mat[i][j] end for -- augment by identity matrix to right aug[i][i + len] = 1 end for aug = ToReducedRowEchelonForm(aug) sequence inv = repeat(repeat(0,len),len) -- remove identity matrix to left for i=1 to len do for j=len+1 to 2*len do inv[i][j-len] = aug[i][j] end for end for return inv end function constant test = {{ 2, -1, 0}, {-1, 2, -1}, { 0, -1, 2}} pp(inverse(test),{pp_Nest,1})
- Output:
{{0.75,0.5,0.25},
{0.5,1,0.5},
{0.25,0.5,0.75}}
alternate
with javascript_semantics function gauss_jordan_inv(sequence a) integer n = length(a) sequence b = repeat(repeat(0,n),n) for i=1 to n do b[i,i] = 1 end for for k=1 to n do integer lmax = k atom amax = abs(a[k][k]) for l=k+1 to n do atom tmp = abs(a[l][k]) if amax<tmp then {amax,lmax} = {tmp,l} end if end for if k!=lmax then {a[k], a[lmax]} = {a[lmax], a[k]} {b[k], b[lmax]} = {b[lmax], b[k]} end if atom akk = a[k][k] if akk=0 then throw("Irregular matrix") end if for j=1 to n do a[k][j] /= akk b[k][j] /= akk end for for i=1 to n do if i!=k then atom aik = a[i][k] for j=1 to n do a[i][j] -= a[k][j]*aik b[i][j] -= b[k][j]*aik end for end if end for end for return b end function sequence a = {{ 2, -1, 0}, {-1, 2, -1}, { 0, -1, 2}}, inva = gauss_jordan_inv(deep_copy(a)) puts(1,"a = ") pp(a,{pp_Nest,1,pp_Indent,4,pp_IntFmt,"%2d"}) puts(1,"inv(a) = ") pp(inva,{pp_Nest,1,pp_Indent,9,pp_FltFmt,"%4.2f"})
- Output:
a = {{ 2,-1, 0},
{-1, 2,-1},
{ 0,-1, 2}}
inv(a) = {{0.75,0.50,0.25},
{0.50,1.00,0.50},
{0.25,0.50,0.75}}
PL/I
/* Gauss-Jordan matrix inversion */
G_J: procedure options (main); /* 4 November 2020 */
declare t float;
declare (i, j, k, n) fixed binary;
open file (sysin) title ('/GAUSSJOR.DAT');
get (n); /* Read in the order of the matrix. */
put skip data (n);
begin;
declare a(n,n)float, aux(n,n) float;
aux = 0;
do i = 1 to n; aux(i,i) = 1; end;
get (a); /* Read in the matrix. */
put skip list ('The matrix to be inverted is:');
put edit (a) ( skip, (n) F(10,4)); /* Print the matrix. */
do k = 1 to n;
/* Divide row k by a(k,k) */
t = a(k,k); a(k,*) = a(k,*) / t; aux(k,*) = aux(k,*) / t;
do i = 1 to k-1, k+1 to n; /* Work down the rows. */
t = a(i,k);
a(i,*) = a(i,*) - t*a(k,*); aux(i,*) = aux(i,*) - t*aux(k,*);
end;
end;
put skip (2) list ('The inverse is:');
put edit (aux) ( skip, (n) F(10,4));
end; /* of BEGIN block */
end G_J;
Output:
N= 4;
The matrix to be inverted is:
8.0000 3.0000 7.0000 5.0000
9.0000 12.0000 10.0000 11.0000
6.0000 2.0000 4.0000 13.0000
14.0000 15.0000 16.0000 17.0000
The inverse is:
0.3590 0.5385 0.0769 -0.5128
-0.0190 0.3419 -0.0199 -0.2004
-0.1833 -0.7009 -0.1425 0.6163
-0.1064 -0.0855 0.0883 0.0779
PowerShell
function gauss-jordan-inv([double[][]]$a) {
$n = $a.count
[double[][]]$b = 0..($n-1) | foreach{[double[]]$row = @(0) * $n; $row[$_] = 1; ,$row}
for ($k = 0; $k -lt $n; $k++) {
$lmax, $max = $k, [Math]::Abs($a[$k][$k])
for ($l = $k+1; $l -lt $n; $l++) {
$tmp = [Math]::Abs($a[$l][$k])
if($max -lt $tmp) {
$max, $lmax = $tmp, $l
}
}
if ($k -ne $lmax) {
$a[$k], $a[$lmax] = $a[$lmax], $a[$k]
$b[$k], $b[$lmax] = $b[$lmax], $b[$k]
}
$akk = $a[$k][$k]
if (0 -eq $akk) {throw "Irregular matrix"}
for ($j = 0; $j -lt $n; $j++) {
$a[$k][$j] /= $akk
$b[$k][$j] /= $akk
}
for ($i = 0; $i -lt $n; $i++){
if ($i -ne $k) {
$aik = $a[$i][$k]
for ($j = 0; $j -lt $n; $j++) {
$a[$i][$j] -= $a[$k][$j]*$aik
$b[$i][$j] -= $b[$k][$j]*$aik
}
}
}
}
$b
}
function show($a) { $a | foreach{ "$_"} }
$a = @(@(@(1, 2, 3), @(4, 1, 6), @(7, 8, 9)))
$inva = gauss-jordan-inv $a
"a ="
show $a
""
"inv(a) ="
show $inva
""
$b = @(@(2, -1, 0), @(-1, 2, -1), @(0, -1, 2))
"b ="
show $b
""
$invb = gauss-jordan-inv $b
"inv(b) ="
show $invb
Output:
a = 1 2 3 4 1 6 7 8 9 inv(a) = -0.8125 0.125 0.1875 0.125 -0.25 0.125 0.520833333333333 0.125 -0.145833333333333 b = 2 -1 0 -1 2 -1 0 -1 2 inv(b) = 0.75 0.5 0.25 0.5 1 0.5 0.25 0.5 0.75
Python
Use numpy matrix inversion
import numpy as np
from numpy.linalg import inv
a = np.array([[1., 2., 3.], [4., 1., 6.],[ 7., 8., 9.]])
ainv = inv(a)
print(a)
print(ainv)
- Output:
[[1. 2. 3.] [4. 1. 6.] [7. 8. 9.]] [[-0.8125 0.125 0.1875 ] [ 0.125 -0.25 0.125 ] [ 0.52083333 0.125 -0.14583333]]
Racket
Using the matrix library that comes with default Racket installation
#lang racket
(require math/matrix
math/array)
(define (inverse M)
(define dim (square-matrix-size M))
(define MI (matrix-augment (list M (identity-matrix dim))))
(submatrix (matrix-row-echelon MI #t #t) (::) (:: dim #f)))
(define A (matrix [[1 2 3] [4 1 6] [7 8 9]]))
(inverse A)
(matrix-inverse A)
- Output:
(array #[#[-13/16 1/8 3/16] #[1/8 -1/4 1/8] #[25/48 1/8 -7/48]]) (array #[#[-13/16 1/8 3/16] #[1/8 -1/4 1/8] #[25/48 1/8 -7/48]])
Raku
(formerly Perl 6)
Uses bits and pieces from other tasks, Reduced row echelon form primarily.
sub gauss-jordan-invert (@m where &is-square) {
^@m .map: { @m[$_].append: identity(+@m)[$_] };
@m.&rref[*]»[+@m .. *];
}
sub is-square (@m) { so @m == all @m[*] }
sub identity ($n) { [ 1, |(0 xx $n-1) ], *.rotate(-1).Array ... *.tail }
# reduced row echelon form (from 'Gauss-Jordan elimination' task)
sub rref (@m) {
my ($lead, $rows, $cols) = 0, @m, @m[0];
for ^$rows -> $r {
$lead < $cols or return @m;
my $i = $r;
until @m[$i;$lead] {
++$i == $rows or next;
$i = $r;
++$lead == $cols and return @m;
}
@m[$i, $r] = @m[$r, $i] if $r != $i;
@m[$r] »/=» $ = @m[$r;$lead];
for ^$rows -> $n {
next if $n == $r;
@m[$n] »-=» @m[$r] »×» (@m[$n;$lead] // 0);
}
++$lead;
}
@m
}
sub rat-or-int ($num) {
return $num unless $num ~~ Rat|FatRat;
return $num.narrow if $num.narrow ~~ Int;
$num.nude.join: '/';
}
sub say_it ($message, @array) {
my $max;
@array.map: {$max max= max $_».&rat-or-int.comb(/\S+/)».chars};
say "\n$message";
$_».&rat-or-int.fmt(" %{$max}s").put for @array;
}
multi to-matrix ($str) { [$str.split(';').map(*.words.Array)] }
multi to-matrix (@array) { @array }
sub hilbert-matrix ($h) {
[ (1..$h).map( -> $n { [ ($n ..^ $n + $h).map: { (1/$_).FatRat } ] } ) ]
}
my @tests =
'1 2 3; 4 1 6; 7 8 9',
'2 -1 0; -1 2 -1; 0 -1 2',
'-1 -2 3 2; -4 -1 6 2; 7 -8 9 1; 1 -2 1 3',
'1 2 3 4; 5 6 7 8; 9 33 11 12; 13 14 15 17',
'3 1 8 9 6; 6 2 8 10 1; 5 7 2 10 3; 3 2 7 7 9; 3 5 6 1 1',
# Test with a Hilbert matrix
hilbert-matrix 10;
@tests.map: {
my @matrix = .&to-matrix;
say_it( " {'=' x 20} Original Matrix: {'=' x 20}", @matrix );
say_it( ' Gauss-Jordan Inverted:', my @invert = gauss-jordan-invert @matrix );
say_it( ' Re-inverted:', gauss-jordan-invert @invert».Array );
}
- Output:
==================== Original Matrix: ====================
1 2 3
4 1 6
7 8 9
Gauss-Jordan Inverted:
-13/16 1/8 3/16
1/8 -1/4 1/8
25/48 1/8 -7/48
Re-inverted:
1 2 3
4 1 6
7 8 9
==================== Original Matrix: ====================
2 -1 0
-1 2 -1
0 -1 2
Gauss-Jordan Inverted:
3/4 1/2 1/4
1/2 1 1/2
1/4 1/2 3/4
Re-inverted:
2 -1 0
-1 2 -1
0 -1 2
==================== Original Matrix: ====================
-1 -2 3 2
-4 -1 6 2
7 -8 9 1
1 -2 1 3
Gauss-Jordan Inverted:
-21/23 17/69 13/138 19/46
-38/23 15/23 1/23 15/23
-16/23 25/69 11/138 9/46
-13/23 16/69 -2/69 13/23
Re-inverted:
-1 -2 3 2
-4 -1 6 2
7 -8 9 1
1 -2 1 3
==================== Original Matrix: ====================
1 2 3 4
5 6 7 8
9 33 11 12
13 14 15 17
Gauss-Jordan Inverted:
19/184 -199/184 -1/46 1/2
1/23 -2/23 1/23 0
-441/184 813/184 -1/46 -3/2
2 -3 0 1
Re-inverted:
1 2 3 4
5 6 7 8
9 33 11 12
13 14 15 17
==================== Original Matrix: ====================
3 1 8 9 6
6 2 8 10 1
5 7 2 10 3
3 2 7 7 9
3 5 6 1 1
Gauss-Jordan Inverted:
-4525/6238 2529/6238 -233/3119 1481/3119 -639/6238
1033/6238 -1075/6238 342/3119 -447/3119 871/6238
1299/6238 -289/6238 -204/3119 -390/3119 739/6238
782/3119 -222/3119 237/3119 -556/3119 -177/3119
-474/3119 -17/3119 -24/3119 688/3119 -140/3119
Re-inverted:
3 1 8 9 6
6 2 8 10 1
5 7 2 10 3
3 2 7 7 9
3 5 6 1 1
==================== Original Matrix: ====================
1 1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10
1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11
1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12
1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13
1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14
1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15
1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16
1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17
1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18
1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19
Gauss-Jordan Inverted:
100 -4950 79200 -600600 2522520 -6306300 9609600 -8751600 4375800 -923780
-4950 326700 -5880600 47567520 -208107900 535134600 -832431600 770140800 -389883780 83140200
79200 -5880600 112907520 -951350400 4281076800 -11237826600 17758540800 -16635041280 8506555200 -1829084400
-600600 47567520 -951350400 8245036800 -37875637800 101001700800 -161602721280 152907955200 -78843164400 17071454400
2522520 -208107900 4281076800 -37875637800 176752976400 -477233036280 771285715200 -735869534400 382086104400 -83223340200
-6306300 535134600 -11237826600 101001700800 -477233036280 1301544644400 -2121035716800 2037792556800 -1064382719400 233025352560
9609600 -832431600 17758540800 -161602721280 771285715200 -2121035716800 3480673996800 -3363975014400 1766086882560 -388375587600
-8751600 770140800 -16635041280 152907955200 -735869534400 2037792556800 -3363975014400 3267861442560 -1723286307600 380449555200
4375800 -389883780 8506555200 -78843164400 382086104400 -1064382719400 1766086882560 -1723286307600 912328045200 -202113826200
-923780 83140200 -1829084400 17071454400 -83223340200 233025352560 -388375587600 380449555200 -202113826200 44914183600
Re-inverted:
1 1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10
1/2 1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11
1/3 1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12
1/4 1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13
1/5 1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14
1/6 1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15
1/7 1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16
1/8 1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17
1/9 1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18
1/10 1/11 1/12 1/13 1/14 1/15 1/16 1/17 1/18 1/19REXX
version 1
/* REXX */
Parse Arg seed nn
If seed='' Then
seed=23345
If nn='' Then nn=5
If seed='?' Then Do
Say 'rexx gjmi seed n computes a random matrix with n rows and columns'
Say 'Default is 23345 5'
Exit
End
Numeric Digits 50
Call random 1,2,seed
a=''
Do i=1 To nn**2
a=a random(9)+1
End
n2=words(a)
Do n=2 To n2/2
If n**2=n2 Then
Leave
End
If n>n2/2 Then
Call exit 'Not a square matrix:' a '('n2 'elements).'
det=determinante(a,n)
If det=0 Then
Call exit 'Determinant is 0'
Do j=1 To n
Do i=1 To n
Parse Var A a.i.j a
aa.i.j=a.i.j
End
Do ii=1 To n
z=(ii=j)
iii=ii+n
a.iii.j=z
End
End
Call show 1,'The given matrix'
Do m=1 To n-1
If a.m.m=0 Then Do
Do j=m+1 To n
If a.m.j<>0 Then Leave
End
If j>n Then Do
Say 'No pivot>0 found in column' m
Exit
End
Do i=1 To n*2
temp=a.i.m
a.i.m=a.i.j
a.i.j=temp
End
End
Do j=m+1 To n
If a.m.j<>0 Then Do
jj=m
fact=divide(a.m.m,a.m.j)
Do i=1 To n*2
a.i.j=subtract(multiply(a.i.j,fact),a.i.jj)
End
End
End
Call show 2 m
End
Say 'Lower part has all zeros'
Say ''
Do j=1 To n
If denom(a.j.j)<0 Then Do
Do i=1 To 2*n
a.i.j=subtract(0,a.i.j)
End
End
End
Call show 3
Do m=n To 2 By -1
Do j=1 To m-1
jj=m
fact=divide(a.m.j,a.m.jj)
Do i=1 To n*2
a.i.j=subtract(a.i.j,multiply(a.i.jj,fact))
End
End
Call show 4 m
End
Say 'Upper half has all zeros'
Say ''
Do j=1 To n
If decimal(a.j.j)<>1 Then Do
z=a.j.j
Do i=1 To 2*n
a.i.j=divide(a.i.j,z)
End
End
End
Call show 5
Say 'Main diagonal has all ones'
Say ''
Do j=1 To n
Do i=1 To n
z=i+n
a.i.j=a.z.j
End
End
Call show 6,'The inverse matrix'
do i = 1 to n
do j = 1 to n
sum=0
Do k=1 To n
sum=add(sum,multiply(aa.i.k,a.k.j))
End
c.i.j = sum
end
End
Call showc 7,'The product of input and inverse matrix'
Exit
show:
Parse Arg num,text
Say 'show' arg(1) text
If arg(1)<>6 Then rows=n*2
Else rows=n
len=0
Do j=1 To n
Do i=1 To rows
len=max(len,length(a.i.j))
End
End
Do j=1 To n
ol=''
Do i=1 To rows
ol=ol||right(a.i.j,len+1)
End
Say ol
End
Say ''
Return
showc:
Parse Arg num,text
Say text
clen=0
Do j=1 To n
Do i=1 To n
clen=max(clen,length(c.i.j))
End
End
Do j=1 To n
ol=''
Do i=1 To n
ol=ol||right(c.i.j,clen+1)
End
Say ol
End
Say ''
Return
denom: Procedure
/* Return the denominator */
Parse Arg d '/' n
Return d
decimal: Procedure
/* compute the fraction's value */
Parse Arg a
If pos('/',a)=0 Then a=a'/1'; Parse Var a ad '/' an
Return ad/an
gcd: procedure
/**********************************************************************
* Greatest commn divisor
**********************************************************************/
Parse Arg a,b
If b = 0 Then Return abs(a)
Return gcd(b,a//b)
add: Procedure
Parse Arg a,b
If pos('/',a)=0 Then a=a'/1'; Parse Var a ad '/' an
If pos('/',b)=0 Then b=b'/1'; Parse Var b bd '/' bn
sum=divide(ad*bn+bd*an,an*bn)
Return sum
multiply: Procedure
Parse Arg a,b
If pos('/',a)=0 Then a=a'/1'; Parse Var a ad '/' an
If pos('/',b)=0 Then b=b'/1'; Parse Var b bd '/' bn
prd=divide(ad*bd,an*bn)
Return prd
subtract: Procedure
Parse Arg a,b
If pos('/',a)=0 Then a=a'/1'; Parse Var a ad '/' an
If pos('/',b)=0 Then b=b'/1'; Parse Var b bd '/' bn
div=divide(ad*bn-bd*an,an*bn)
Return div
divide: Procedure
Parse Arg a,b
If pos('/',a)=0 Then a=a'/1'; Parse Var a ad '/' an
If pos('/',b)=0 Then b=b'/1'; Parse Var b bd '/' bn
sd=ad*bn
sn=an*bd
g=gcd(sd,sn)
Select
When sd=0 Then res='0'
When abs(sn/g)=1 Then res=(sd/g)*sign(sn/g)
Otherwise Do
den=sd/g
nom=sn/g
If nom<0 Then Do
If den<0 Then den=abs(den)
Else den=-den
nom=abs(nom)
End
res=den'/'nom
End
End
Return res
determinante: Procedure
/* REXX ***************************************************************
* determinant.rex
* compute the determinant of the given square matrix
* Input: as: the representation of the matrix as vector (n**2 elements)
* 21.05.2013 Walter Pachl
**********************************************************************/
Parse Arg as,n
Do i=1 To n
Do j=1 To n
Parse Var as a.i.j as
End
End
Select
When n=2 Then det=subtract(multiply(a.1.1,a.2.2),multiply(a.1.2,a.2.1))
When n=3 Then Do
det=multiply(multiply(a.1.1,a.2.2),a.3.3)
det=add(det,multiply(multiply(a.1.2,a.2.3),a.3.1))
det=add(det,multiply(multiply(a.1.3,a.2.1),a.3.2))
det=subtract(det,multiply(multiply(a.1.3,a.2.2),a.3.1))
det=subtract(det,multiply(multiply(a.1.2,a.2.1),a.3.3))
det=subtract(det,multiply(multiply(a.1.1,a.2.3),a.3.2))
End
Otherwise Do
det=0
Do k=1 To n
sign=((-1)**(k+1))
If sign=1 Then
det=add(det,multiply(a.1.k,determinante(subm(k),n-1)))
Else
det=subtract(det,multiply(a.1.k,determinante(subm(k),n-1)))
End
End
End
Return det
subm: Procedure Expose a. n
/**********************************************************************
* compute the submatrix resulting when row 1 and column k are removed
* Input: a.*.*, k
* Output: bs the representation of the submatrix as vector
**********************************************************************/
Parse Arg k
bs=''
do i=2 To n
Do j=1 To n
If j=k Then Iterate
bs=bs a.i.j
End
End
Return bs
Exit: Say arg(1)
- Output:
Using the defaults for seed and n
show 1 The given matrix
10 3 8 6 3 1 0 0 0 0
5 7 8 8 2 0 1 0 0 0
4 10 5 4 7 0 0 1 0 0
9 4 5 3 3 0 0 0 1 0
6 3 3 3 7 0 0 0 0 1
show 2 1
10 3 8 6 3 1 0 0 0 0
0 11 8 10 1 -1 2 0 0 0
0 22 9/2 4 29/2 -1 0 5/2 0 0
0 13/9 -22/9 -8/3 1/3 -1 0 0 10/9 0
0 2 -3 -1 26/3 -1 0 0 0 5/3
show 2 2
10 3 8 6 3 1 0 0 0 0
0 11 8 10 1 -1 2 0 0 0
0 0 -23/4 -8 25/4 1/2 -2 5/4 0 0
0 0 -346/13 -394/13 20/13 -86/13 -2 0 110/13 0
0 0 -49/2 -31/2 140/3 -9/2 -2 0 0 55/6
show 2 3
10 3 8 6 3 1 0 0 0 0
0 11 8 10 1 -1 2 0 0 0
0 0 -23/4 -8 25/4 1/2 -2 5/4 0 0
0 0 0 1005/692 -4095/692 -1335/692 1085/692 -5/4 1265/692 0
0 0 0 855/196 395/84 -305/196 75/49 -5/4 0 1265/588
show 2 4
10 3 8 6 3 1 0 0 0 0
0 11 8 10 1 -1 2 0 0 0
0 0 -23/4 -8 25/4 1/2 -2 5/4 0 0
0 0 0 1005/692 -4095/692 -1335/692 1085/692 -5/4 1265/692 0
0 0 0 0 221375/29583 13915/9861 -13915/13148 16445/19722 -1265/692 84755/118332
Lower part has all zeros
show 3
10 3 8 6 3 1 0 0 0 0
0 11 8 10 1 -1 2 0 0 0
0 0 23/4 8 -25/4 -1/2 2 -5/4 0 0
0 0 0 1005/692 -4095/692 -1335/692 1085/692 -5/4 1265/692 0
0 0 0 0 221375/29583 13915/9861 -13915/13148 16445/19722 -1265/692 84755/118332
show 4 5
10 3 8 6 0 76/175 297/700 -117/350 513/700 -201/700
0 11 8 10 0 -208/175 1499/700 -39/350 171/700 -67/700
0 0 23/4 8 0 19/28 125/112 -31/56 -171/112 67/112
0 0 0 1005/692 0 -1407/1730 10117/13840 -4087/6920 5293/13840 7839/13840
0 0 0 0 221375/29583 13915/9861 -13915/13148 16445/19722 -1265/692 84755/118332
show 4 4
10 3 8 0 0 664/175 -1817/700 737/350 -593/700 -1839/700
0 11 8 0 0 772/175 -6073/2100 4153/1050 -5017/2100 -2797/700
0 0 23/4 0 0 3611/700 -24449/8400 11339/4200 -30521/8400 -7061/2800
0 0 0 1005/692 0 -1407/1730 10117/13840 -4087/6920 5293/13840 7839/13840
0 0 0 0 221375/29583 13915/9861 -13915/13148 16445/19722 -1265/692 84755/118332
show 4 3
10 3 0 0 0 -592/175 3053/2100 -1733/1050 8837/2100 617/700
0 11 0 0 0 -484/175 2431/2100 209/1050 5599/2100 -341/700
0 0 23/4 0 0 3611/700 -24449/8400 11339/4200 -30521/8400 -7061/2800
0 0 0 1005/692 0 -1407/1730 10117/13840 -4087/6920 5293/13840 7839/13840
0 0 0 0 221375/29583 13915/9861 -13915/13148 16445/19722 -1265/692 84755/118332
show 4 2
10 0 0 0 0 -92/35 239/210 -179/105 731/210 71/70
0 11 0 0 0 -484/175 2431/2100 209/1050 5599/2100 -341/700
0 0 23/4 0 0 3611/700 -24449/8400 11339/4200 -30521/8400 -7061/2800
0 0 0 1005/692 0 -1407/1730 10117/13840 -4087/6920 5293/13840 7839/13840
0 0 0 0 221375/29583 13915/9861 -13915/13148 16445/19722 -1265/692 84755/118332
Upper half has all zeros
show 5
1 0 0 0 0 -46/175 239/2100 -179/1050 731/2100 71/700
0 1 0 0 0 -44/175 221/2100 19/1050 509/2100 -31/700
0 0 1 0 0 157/175 -1063/2100 493/1050 -1327/2100 -307/700
0 0 0 1 0 -14/25 151/300 -61/150 79/300 39/100
0 0 0 0 1 33/175 -99/700 39/350 -171/700 67/700
Main diagonal has all ones
show 6 The inverse matrix
-46/175 239/2100 -179/1050 731/2100 71/700
-44/175 221/2100 19/1050 509/2100 -31/700
157/175 -1063/2100 493/1050 -1327/2100 -307/700
-14/25 151/300 -61/150 79/300 39/100
33/175 -99/700 39/350 -171/700 67/700
The product of input and inverse matrix
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
version 2
/*REXX program performs a (square) matrix inversion using the Gauss─Jordan method. */
data= 8 3 7 5 9 12 10 11 6 2 4 13 14 15 16 17 /*the matrix element values. */
call build 4 /*assign data elements to the matrix. */
call show '@', 'The matrix of order ' n " is:" /*display the (square) matrix. */
call aux /*define the auxiliary (identity) array*/
call invert /*invert the matrix, store result in X.*/
call show 'X', "The inverted matrix is:" /*display (inverted) auxiliary matrix. */
exit 0 /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
aux: x.= 0; do i=1 for n; x.i.i= 1; end; return
/*──────────────────────────────────────────────────────────────────────────────────────*/
build: arg n; #=0; w=0; do r=1 for n /*read a row of the matrix.*/
do c=1 for n; #= # + 1 /* " " col " " " */
@.r.c= word(data, #); w= max(w, length(@.r.c) )
end /*c*/ /*W: max width of a number*/
end /*r*/; return
/*──────────────────────────────────────────────────────────────────────────────────────*/
invert: do k=1 for n; t= @.k.k /*divide each matrix row by T. */
do c=1 for n; @.k.c= @.k.c / t /*process row of original matrix.*/
x.k.c= x.k.c / t /* " " " auxiliary " */
end /*c*/
do r=1 for n; if r==k then iterate /*skip if R is the same row as K.*/
t= @.r.k
do c=1 for n; @.r.c= @.r.c - t*@.k.c /*process row of original matrix.*/
x.r.c= x.r.c - t*x.k.c /* " " " auxiliary " */
end /*c*/
end /*r*/
end /*k*/; return
/*──────────────────────────────────────────────────────────────────────────────────────*/
show: parse arg ?, title; say; say title; f= 4 /*F: fractional digits precision.*/
do r=1 for n; _=
do c=1 for n; if ?=='@' then _= _ right( @.r.c, w)
else _= _ right(format(x.r.c, w, f), w + f + length(.))
end /*c*/; say _
end /*r*/; return
- output when using the internal default input:
The matrix of order 4 is: 8 3 7 5 9 12 10 11 6 2 4 13 14 15 16 17 The inverted matrix is: 0.3590 0.5385 0.0769 -0.5128 -0.0190 0.3419 -0.0199 -0.2004 -0.1833 -0.7009 -0.1425 0.6163 -0.1064 -0.0855 0.0883 0.0779
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
- Output:
[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)]]
Rust
fn main() {
let mut a: Vec<Vec<f64>> = vec![vec![1.0, 2.0, 3.0],
vec![4.0, 1.0, 6.0],
vec![7.0, 8.0, 9.0]
];
let mut b: Vec<Vec<f64>> = vec![vec![2.0, -1.0, 0.0],
vec![-1.0, 2.0, -1.0],
vec![0.0, -1.0, 2.0]
];
let mut ref_a = &mut a;
let rref_a = &mut ref_a;
let mut ref_b = &mut b;
let rref_b = &mut ref_b;
println!("Matrix A:\n");
print_matrix(rref_a);
println!("\nInverse of Matrix A:\n");
print_matrix(&mut matrix_inverse(rref_a));
println!("\n\nMatrix B:\n");
print_matrix(rref_b);
println!("\nInverse of Matrix B:\n");
print_matrix(&mut matrix_inverse(rref_b));
}
//Begin Matrix Inversion
fn matrix_inverse(matrix: &mut Vec<Vec<f64>>) -> Vec<Vec<f64>>{
let len = matrix.len();
let mut aug = zero_matrix(len, len * 2);
for i in 0..len {
for j in 0.. len {
aug[i][j] = matrix[i][j];
}
aug[i][i + len] = 1.0;
}
gauss_jordan_general(&mut aug);
let mut unaug = zero_matrix(len, len);
for i in 0..len {
for j in 0..len {
unaug[i][j] = aug[i][j+len];
}
}
unaug
}
//End Matrix Inversion
//Begin Generalised Reduced Row Echelon Form
fn gauss_jordan_general(matrix: &mut Vec<Vec<f64>>) {
let mut lead = 0;
let row_count = matrix.len();
let col_count = matrix[0].len();
for r in 0..row_count {
if col_count <= lead {
break;
}
let mut i = r;
while matrix[i][lead] == 0.0 {
i = i + 1;
if row_count == i {
i = r;
lead = lead + 1;
if col_count == lead {
break;
}
}
}
let temp = matrix[i].to_owned();
matrix[i] = matrix[r].to_owned();
matrix[r] = temp.to_owned();
if matrix[r][lead] != 0.0 {
let div = matrix[r][lead];
for j in 0..col_count {
matrix[r][j] = matrix[r][j] / div;
}
}
for k in 0..row_count {
if k != r {
let mult = matrix[k][lead];
for j in 0..col_count {
matrix[k][j] = matrix[k][j] - matrix[r][j] * mult;
}
}
}
lead = lead + 1;
}
//matrix.to_owned()
}
fn zero_matrix(rows: usize, cols: usize) -> Vec<Vec<f64>> {
let mut matrix = Vec::with_capacity(cols);
for _ in 0..rows {
let mut col: Vec<f64> = Vec::with_capacity(rows);
for _ in 0..cols {
col.push(0.0);
}
matrix.push(col);
}
matrix
}
fn print_matrix(mat: &mut Vec<Vec<f64>>) {
for row in 0..mat.len(){
println!("{:?}", mat[row]);
}
}
- Output:
Matrix A: [1.0, 2.0, 3.0] [4.0, 1.0, 6.0] [7.0, 8.0, 9.0] Inverse of Matrix A: [-0.8125, 0.125, 0.1875] [0.12499999999999994, -0.24999999999999997, 0.12499999999999997] [0.5208333333333334, 0.12499999999999999, -0.14583333333333331] Matrix B: [2.0, -1.0, 0.0] [-1.0, 2.0, -1.0] [0.0, -1.0, 2.0] Inverse of Matrix B: [0.75, 0.49999999999999994, 0.24999999999999994] [0.49999999999999994, 0.9999999999999999, 0.4999999999999999] [0.24999999999999997, 0.49999999999999994, 0.7499999999999999]
Sidef
Uses the rref(M) function from Reduced row echelon form.
func gauss_jordan_invert (M) {
var I = M.len.of {|i|
M.len.of {|j|
i == j ? 1 : 0
}
}
var A = gather {
^M -> each {|i| take(M[i] + I[i]) }
}
rref(A).map { .last(M.len) }
}
var A = [
[-1, -2, 3, 2],
[-4, -1, 6, 2],
[ 7, -8, 9, 1],
[ 1, -2, 1, 3],
]
say gauss_jordan_invert(A).map {
.map { "%6s" % .as_rat }.join(" ")
}.join("\n")
- Output:
-21/23 17/69 13/138 19/46 -38/23 15/23 1/23 15/23 -16/23 25/69 11/138 9/46 -13/23 16/69 -2/69 13/23
Wren
The above module does in fact use the Gauss-Jordan method for matrix inversion.
import "./matrix" for Matrix
import "./fmt" for Fmt
var arrays = [
[ [ 1, 2, 3],
[ 4, 1, 6],
[ 7, 8, 9] ],
[ [ 2, -1, 0 ],
[-1, 2, -1 ],
[ 0, -1, 2 ] ],
[ [ -1, -2, 3, 2 ],
[ -4, -1, 6, 2 ],
[ 7, -8, 9, 1 ],
[ 1, -2, 1, 3 ] ]
]
for (array in arrays) {
System.print("Original:\n")
var m = Matrix.new(array)
Fmt.mprint(m, 2, 0)
System.print("\nInverse:\n")
Fmt.mprint(m.inverse, 9, 6)
System.print()
}
- Output:
Original: | 1 2 3| | 4 1 6| | 7 8 9| Inverse: |-0.812500 0.125000 0.187500| | 0.125000 -0.250000 0.125000| | 0.520833 0.125000 -0.145833| Original: | 2 -1 0| |-1 2 -1| | 0 -1 2| Inverse: | 0.750000 0.500000 0.250000| | 0.500000 1.000000 0.500000| | 0.250000 0.500000 0.750000| Original: |-1 -2 3 2| |-4 -1 6 2| | 7 -8 9 1| | 1 -2 1 3| Inverse: |-0.913043 0.246377 0.094203 0.413043| |-1.652174 0.652174 0.043478 0.652174| |-0.695652 0.362319 0.079710 0.195652| |-0.565217 0.231884 -0.028986 0.565217|
zkl
This uses GSL to invert a matrix via LU decomposition, not Gauss-Jordan.
var [const] GSL=Import.lib("zklGSL"); // libGSL (GNU Scientific Library)
m:=GSL.Matrix(3,3).set(1,2,3, 4,1,6, 7,8,9);
i:=m.invert();
i.format(10,4).println("\n");
(m*i).format(10,4).println();- Output:
-0.8125, 0.1250, 0.1875
0.1250, -0.2500, 0.1250
0.5208, 0.1250, -0.1458
1.0000, 0.0000, 0.0000
-0.0000, 1.0000, 0.0000
-0.0000, 0.0000, 1.0000
m:=GSL.Matrix(3,3).set(2,-1,0, -1,2,-1, 0,-1,2);
m.invert().format(10,4).println("\n");- Output:
0.7500, 0.5000, 0.2500
0.5000, 1.0000, 0.5000
0.2500, 0.5000, 0.7500
Categories:
- Programming Tasks
- Solutions by Programming Task
- Matrices
- 11l
- 360 Assembly
- ALGOL 60
- ALGOL 68
- ATS
- BASIC
- FreeBASIC
- VBA
- VBScript
- C
- C sharp
- C++
- EasyLang
- Factor
- Fortran
- Generic
- Go
- Haskell
- J
- Java
- Jq
- Julia
- Kotlin
- Lambdatalk
- M2000 Interpreter
- Mathematica
- Wolfram Language
- MATLAB
- Nim
- Pascal
- Free Pascal
- PascalABC.NET
- Perl
- Phix
- PL/I
- PowerShell
- Python
- Racket
- Raku
- REXX
- Ruby
- Rust
- Sidef
- Wren
- Wren-fmt
- Wren-matrix
- Zkl
- PL/0/Omit