# Gauss-Jordan matrix inversion

Gauss-Jordan matrix inversion
You are encouraged to solve this task according to the task description, using any language you may know.

Invert matrix   A   using Gauss-Jordan method.

A   being an   n × n   matrix.

## 11l

Translation of: Nim
```V Eps = 1e-10

F transformToRref(&mat)

L(r) 0 .< mat.len
R

V i = r
i++
I i == mat.len
i = r
R
swap(&mat[i], &mat[r])

I abs(d) > Eps
L(&item) mat[r]
item /= d

L(i) 0 .< mat.len
I i != r
L(c) 0 .< mat[0].len
mat[i][c] -= mat[r][c] * m

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

Works with: A60
```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

Translation of: ALGOL 60
```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). *)

%{^
#include <math.h>
#include <float.h>
%}

macdef NAN = g0f2f (\$extval (float, "NAN"))
macdef Zero = g0i2f 0
macdef One = g0i2f 1
macdef Two = g0i2f 2

(* "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)
#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

(*------------------------------------------------------------------*)
(* 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 =
fn
index_map : Matrix_Index_Map (n, n, n, n) =
(@(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
in
let
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

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

Translation of: Phix

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

Translation of: Rexx

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

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

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)
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

Translation of: Fortran
```/*----------------------------------------------------------------------
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 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 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;
return;
}
}
swap_rows(m, i, row);
if (m(row, lead) != 0) {
for (size_t column = 0; column < columns; ++column)
m(row, column) /= f;
}
for (size_t j = 0; j < rows; ++j) {
if (j == row)
continue;
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

Translation of: Go
```proc rref . m[][] .
nrow = len m[][]
ncol = len m[1][]
for r to nrow
return
.
i = r
i += 1
if i > nrow
i = r
return
.
.
.
swap m[i][] m[r][]
for k to ncol
m[r][k] /= m
.
for i to nrow
if i <> r
for k to ncol
m[i][k] -= m * m[r][k]
.
.
.
.
.
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.

Works with: Factor version 0.99 2020-01-23
```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:

Works with: gfortran version 8.3.0
```!             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;
}
}

{
c = cols;
i = 0;
uuiil i < c
{
this[a, i] = this[a, i] + this[b, i];
i++;
}
}

{
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++;
}
}

{
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++;
}
}

{
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

Translation of: Kotlin
```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() {
rowCount, colCount := len(m), len(m[0])
for r := 0; r < rowCount; r++ {
return
}
i := r

i++
if rowCount == i {
i = r
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 {
for j := 0; j < colCount; j++ {
m[k][j] -= m[r][j] * mult
}
}
}
}
}

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

Translation of: PowerShell
```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`

```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;
if (++i == rows) {
i = row;
return;
}
}
swapRows(i, row);
for (int column = 0; column < columns; ++column)
elements[row][column] /= f;
}
for (int j = 0; j < rows; ++j) {
if (j == row)
continue;
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

Works with: jq

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
.i += 1
| if \$nr == .i
then .i = \$r
else .
end
)
| swapRows(.i; \$r)
| reduce range(0; \$nc) as \$j (.; .a[\$r][\$j] /= \$div)
else .
end
| reduce range(0; \$nr) as \$k (.;
if \$k != \$r
| reduce range(0; \$nc) as \$j (.; .a[\$k][\$j] -= .a[\$r][\$j] * \$mult)
else .
end)
))
| .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] )) ;```

```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 ] ]
];

arrays[]
| "Original:",
mprint(1),
"\nInverse:",
(inverse|mprint(6)),
""
;

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

Works with: Julia version 0.6

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() {
val rowCount = this.size
val colCount = this[0].size
for (r in 0 until rowCount) {
var i = r

i++
if (rowCount == i) {
i = r
}
}

val temp = this[i]
this[i] = this[r]
this[r] = temp

for (j in 0 until colCount) this[r][j] /= div
}

for (k in 0 until rowCount) {
if (k != r) {
for (j in 0 until colCount) this[k][j] -= this[r][j] * mult
}
}

}
}

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 {
i=r
i++
}

for c =bj to columncount-1+bj {swap A(i, c), A(r, c)}

For c =bj to columncount-1+bj {A(r, c)/=div1}
end if
for i=bi to rowcount-1+bi {
if i<>r then {
for j=bj to columncount-1+bj {A(i,j)-=A(r,j)*mult}
}
}
}
=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"
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.

for r in 0..<mat.M:

var i = r
inc i
if i == mat.M:
i = r
swap mat[i], mat[r]

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:
for c in 0..<mat.N:
mat[i][c] -= mat[r][c] * m

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:
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) {
my \$i = \$r;

{++\$i == \$rows or next;
\$i = \$r;

@m[\$i, \$r] = @m[\$r, \$i];
\$_ /= \$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] };}

}

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

Translation of: Kotlin

uses ToReducedRowEchelonForm() from Reduced_row_echelon_form#Phix

```with javascript_semantics
function ToReducedRowEchelonForm(sequence M)
rowCount = length(M),
columnCount = length(M[1]),
i
for r=1 to rowCount do
if lead>=columnCount then exit end if
i = r
i += 1
if i=rowCount then
i = r
if lead=columnCount then exit end if
end if
end while
M[i] = M[r]
M[r] = mr
for j=1 to rowCount do
if j!=r then
end if
end for
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

Translation of: PowerShell
```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;
++\$i == \$rows or next;
\$i = \$r;
++\$lead == \$cols and return @m;
}
@m[\$i, \$r] = @m[\$r, \$i] if \$r != \$i;
for ^\$rows -> \$n {
next if \$n == \$r;
@m[\$n] »-=» @m[\$r] »×» (@m[\$n;\$lead] // 0);
}
}
@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/19```

## REXX

### 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
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

gcd: procedure
/**********************************************************************
* Greatest commn divisor
**********************************************************************/
Parse Arg a,b
If b = 0 Then Return abs(a)
Return gcd(b,a//b)

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
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
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
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
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=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
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

Translation of: PL/I
```/*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 row_count = matrix.len();
let col_count = matrix[0].len();

for r in 0..row_count {
break;
}
let mut i = r;
i = i + 1;
if row_count == i {
i = r;
break;
}
}
}

let temp = matrix[i].to_owned();
matrix[i] = matrix[r].to_owned();
matrix[r] = temp.to_owned();

for j in 0..col_count {
matrix[r][j] = matrix[r][j] / div;
}
}

for k in 0..row_count {
if k != r {
for j in 0..col_count {
matrix[k][j] = matrix[k][j] - matrix[r][j] * mult;
}
}
}

}
//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.

Translation of: Raku
```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

Library: Wren-fmt
Library: Wren-matrix

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
```