View source for Gauss-Jordan matrix inversion

{{task}}[[Category:Matrices]]

;Task:
Invert matrix &nbsp; '''A''' &nbsp; using Gauss-Jordan method.

'''A''' &nbsp; being an &nbsp; '''n''' &times; '''n''' &nbsp; matrix.
<br><br>

=={{header|11l}}==
{{trans|Nim}}

<syntaxhighlight lang="11l">V Eps = 1e-10

F transformToRref(&mat)
   V lead = 0

   L(r) 0 .< mat.len
      I lead >= mat[0].len
         R

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

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

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

      lead++

F inverse(mat)
   V augmat = [[0.0] * (2 * mat.len)] * mat.len
   L(i) 0 .< mat.len
      L(j) 0 .< mat.len
         augmat[i][j] = mat[i][j]
      augmat[i][mat.len + i] = 1

   transformToRref(&augmat)

   V result = [[0.0] * mat.len] * mat.len
   L(i) 0 .< mat.len
      L(j) 0 .< mat.len
         I augmat[i][j] != Float(i == j)
            X.throw ValueError(‘matrix is singular’)
         result[i][j] = augmat[i][mat.len + j]
   R result

F print_mat(mat)
   L(row) mat
      V line = ‘’
      L(val) row
         I !line.empty
            line ‘’= ‘ ’
         line ‘’= ‘#3.5’.format(val)
      print(line)

F runTest(mat)
   print(‘Matrix:’)
   print_mat(mat)
   print()
   print(‘Inverse:’)
   print_mat(inverse(mat))
   print()
   print()

V m1 = [[Float(1), 2, 3],
        [Float(4), 1, 6],
        [Float(7), 8, 9]]

V m2 = [[Float( 2), -1,  0],
        [Float(-1),  2, -1],
        [Float( 0), -1,  2]]

V m3 = [[Float(-1), -2,  3,  2],
        [Float(-4), -1,  6,  2],
        [Float( 7), -8,  9,  1],
        [Float( 1), -2,  1,  3]]

runTest(m1)
runTest(m2)
runTest(m3)</syntaxhighlight>

{{out}}
<pre>
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


</pre>

=={{header|360 Assembly}}==
The COPY file of FORMAT, to format a floating point number, can be found at:
[[360_Assembly_include|Include files 360 Assembly]].
<syntaxhighlight lang="360asm">*        Gauss-Jordan matrix inversion 17/01/2021
GAUSSJOR CSECT
         USING  GAUSSJOR,R13       base register
         B      72(R15)            skip savearea
         DC     17F'0'             savearea
         SAVE   (14,12)            save previous context
         ST     R13,4(R15)         link backward
         ST     R15,8(R13)         link forward
         LR     R13,R15            set addressability
         LA     R6,1               i=1
       DO  WHILE=(C,R6,LE,N)       do i=1 to n
         LA     R7,1                 j=1
       DO  WHILE=(C,R7,LE,N)         do j=1 to n
         LR     R1,R6                  i
         BCTR   R1,0                   -1
         MH     R1,HN                  *n
         LR     R0,R7                   j
         BCTR   R0,0                   -1
         AR     R1,R0                  (i-1)*n+j-1
         SLA    R1,2                   *4
         LE     F0,AA(R1)              aa(i,j)
         LR     R1,R6                  i
         BCTR   R1,0                   -1
         MH     R1,HN2                 *n*2
         LR     R0,R7                  j
         BCTR   R0,0                   -1
         AR     R1,R0                  (i-1)*n*2+j-1
         SLA    R1,2                   *4
         STE    F0,TMP(R1)             tmp(i,j)=aa(i,j)
         LA     R7,1(R7)               j++
       ENDDO    ,                    enddo j
         LA     R7,1                 j=1
         A      R7,N                 j=n+1
       DO  WHILE=(C,R7,LE,N2)        do j=n+1 to 2*n
         LR     R1,R6                  i
         BCTR   R1,0                   -1
         MH     R1,HN2                 *n*2
         LR     R0,R7                  j
         BCTR   R0,0                   -1
         AR     R1,R0                  (i-1)*n*2+j-1
         SLA    R1,2                   *4
         LE     F0,=E'0'               0
         STE    F0,TMP(R1)             tmp(i,j)=0
         LA     R7,1(R7)               j++
       ENDDO    ,                    enddo j
         LR     R2,R6                i
         A      R2,N                 i+n
         LR     R1,R6                i
         BCTR   R1,0                 -1
         MH     R1,HN2               *n*2
         BCTR   R2,0                 -1
         AR     R1,R2                (i+n-1)*n*2+i-1
         SLA    R1,2                 *4
         LE     F0,=E'1'             1
         STE    F0,TMP(R1)           tmp(i,i+n)=1
         LA     R6,1(R6)             i++
       ENDDO    ,                  enddo i
         LA     R6,1               i=1
       DO  WHILE=(C,R6,LE,N)       do r=1 to n
         LR     R1,R6                r
         BCTR   R1,0                 -1
         MH     R1,HN2               *n*2
         LR     R0,R6                r
         BCTR   R0,0                 -1
         AR     R1,R0                (r-1)*n*2+r-1
         SLA    R1,2                 *4
         LE     F0,TMP(R1)           tmp(r,r)
         STE    F0,T                 t=tmp(r,r)
         LA     R7,1                 c=1 
       DO WHILE=(C,R7,LE,N2)         do c=1 to n*2
         LR     R1,R6                  r
         BCTR   R1,0
         MH     R1,HN2
         LR     R0,R7                  c
         BCTR   R0,0
         AR     R1,R0
         SLA    R1,2
         LE     F0,TMP(R1)
         LTER   F0,F0
         BZ     *+8
         DE     F0,T                   tmp(r,c)/t
         LR     R1,R6                  r
         BCTR   R1,0
         MH     R1,HN2
         LR     R0,R7                  c
         BCTR   R0,0
         AR     R1,R0
         SLA    R1,2
         STE    F0,TMP(R1)             tmp(r,c)=tmp(r,c)/t
         LA     R7,1(R7)               c++
       ENDDO    ,                    enddo c
         LA     R8,1                   i=1 
       DO WHILE=(C,R8,LE,N)            do i=1 to n
       IF    CR,R8,NE,R6 THEN            if i^=r then
         LR     R1,R8                      i                    
         BCTR   R1,0
         MH     R1,HN2
         LR     R0,R6                      r
         BCTR   R0,0
         AR     R1,R0
         SLA    R1,2
         LE     F0,TMP(R1)                 tmp(i,r)
         STE    F0,T                       t=tmp(i,r)
         LA     R7,1                       c=1 
       DO WHILE=(C,R7,LE,N2)               do c=1 to n*2
         LR     R2,R8                        i
         BCTR   R2,0
         MH     R2,HN2
         LR     R0,R7                        c
         BCTR   R0,0
         AR     R2,R0
         SLA    R2,2
         LE     F0,TMP(R2)                   f0=tmp(i,c)
         LR     R1,R6                        r
         BCTR   R1,0
         MH     R1,HN2
         LR     R0,R7                        c
         BCTR   R0,0
         AR     R1,R0
         SLA    R1,2
         LE     F2,TMP(R1)                   f2=tmp(r,c)
         LE     F4,T                         t
         MER    F4,F2                        f4=t*tmp(r,c)
         SER    F0,F4                        f0=tmp(i,c)-t*tmp(r,c)
         STE    F0,TMP(R2)                   tmp(i,c)=f0
         LA     R7,1(R7)                     c++
       ENDDO    ,                          enddo c
       ENDIF    ,                        endif
         LA     R8,1(R8)                 i++
       ENDDO    ,                      enddo i
         LA     R6,1(R6)             r++
       ENDDO    ,                  enddo r
         LA     R6,1               i=1 
       DO WHILE=(C,R6,LE,N)        do i=1 to n
         L      R7,N                 n
         LA     R7,1(R7)             j=n+1 
       DO WHILE=(C,R7,LE,N2)       do j=n+1 to 2*n
         LR     R2,R7                  j
         S      R2,N                   -n
         LR     R1,R6                  i
         BCTR   R1,0
         MH     R1,HN2                 *2*n
         LR     R0,R7                  j
         BCTR   R0,0
         AR     R1,R0
         SLA    R1,2
         LE     F0,TMP(R1)             tmp(i,j)
         LR     R1,R6                  i
         BCTR   R1,0
         MH     R1,HN                  *n
         BCTR   R2,0
         AR     R1,R2
         SLA    R1,2
         STE    F0,BB(R1)              bb(i,j-n)=tmp(i,j)
         LA     R7,1(R7)               j++
       ENDDO    ,                    enddo j
         LA     R6,1(R6)             i++ 
       ENDDO    ,                  enddo i
         LA     R6,1               i=1 
       DO WHILE=(C,R6,LE,N)        do i=1 to n
         LA     R3,PG                @pg
         LA     R7,1                 j=1 
       DO WHILE=(C,R7,LE,N)          do j=1 to n
         LR     R1,R6                  i
         BCTR   R1,0
         MH     R1,HN                  *n
         LR     R0,R7                  j
         BCTR   R0,0
         AR     R1,R0
         SLA    R1,2
         LE     F0,BB(R1)              bb(i,j) 
         LA     R0,5                   number of decimals
         BAL    R14,FORMATF            FormatF(bb(i,j))
         MVC    0(13,R3),0(R1)         edit
         LA     R3,13(R3)              @pg+=13
         LA     R7,1(R7)               j++
       ENDDO    ,                    enddo j
         XPRNT  PG,L'PG              print buffer
         LA     R6,1(R6)             i++ 
       ENDDO    ,                  enddo i
         L      R13,4(0,R13)       restore previous savearea pointer
         RETURN (14,12),RC=0       restore registers from calling save
         COPY   plig\$_FORMATF.MLC format a float
NN       EQU    5
N        DC     A(NN)
N2       DC     A(NN*2)
AA       DC     E'3.0',E'1.0',E'5.0',E'9.0',E'7.0'
         DC     E'8.0',E'2.0',E'8.0',E'0.0',E'1.0'
         DC     E'1.0',E'7.0',E'2.0',E'0.0',E'3.0'
         DC     E'0.0',E'1.0',E'7.0',E'0.0',E'9.0'
         DC     E'3.0',E'5.0',E'6.0',E'1.0',E'1.0'
BB       DS     (NN*NN)E
TMP      DS     (NN*NN*2)E
T        DS     E
HN       DC     AL2(NN)
HN2      DC     AL2(NN*2)
PG       DC     CL80' '            buffer
         REGEQU
         END    GAUSSJOR </syntaxhighlight>
{{out}}
<pre>
      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
</pre>

=={{header|ALGOL 60}}==
{{works with|A60}}
<syntaxhighlight lang="algol60">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 </syntaxhighlight>
{{out}}
<pre>
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
</pre>

=={{header|ALGOL 68}}==
{{Trans|ALGOL 60}}
<syntaxhighlight lang="algol68">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</syntaxhighlight>
{{out}}
<pre>
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
</pre>

=={{header|ATS}}==

Not all the parts of the "little matrix library" included in the code are actually used, although some of those not used were (in fact) used during debugging.

<syntaxhighlight lang="ats">
(* Set to 1 for debugging: to fill in ones and zeros and other values
   that are not actually used but are part of the theory of the
   Gauss-Jordan algorithm. *)
#define DO_THINGS_THAT_DO_NOT_NEED_TO_BE_DONE 0

(* Setting this to 1 may cause rounding to change, and the change in
   rounding is not unlikely to cause detection of singularity of a
   matrix to change. (To invert a matrix that might be nearly
   singular, the SVD seems a popular method.)  The
   -fexpensive-optimizations option to GCC also may cause the same
   rounding changes (due to fused-multiply-and-add instructions being
   generated). *)
#define USE_MULTIPLY_AND_ADD 0

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

#include "share/atspre_staload.hats"

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

#if USE_MULTIPLY_AND_ADD #then
(* "fma" from the C math library, although your system may have it as
   a built-in. *)
extern fn {tk : tkind} g0float_fma : (g0float tk, g0float tk, g0float tk) -<> g0float tk
implement g0float_fma<fltknd> (x, y, z) = $extfcall (float, "fmaf", x, y, z)
implement g0float_fma<dblknd> (x, y, z) = $extfcall (double, "fma", x, y, z)
implement g0float_fma<ldblknd> (x, y, z) = $extfcall (ldouble, "fmal", x, y, z)
overload fma with g0float_fma
macdef multiply_and_add = fma
#else
macdef multiply_and_add (x, y, z) = (,(x) * ,(y)) + ,(z)
#endif

(*------------------------------------------------------------------*)
(* A "little matrix library"                                        *)

typedef Matrix_Index_Map (m1 : int, n1 : int, m0 : int, n0 : int) =
  {i1, j1 : pos | i1 <= m1; j1 <= n1}
  (int i1, int j1) -<cloref0>
  [i0, j0 : pos | i0 <= m0; j0 <= n0]
    @(int i0, int j0)

datatype Real_Matrix (tk : tkind,
                      m1 : int, n1 : int,
                      m0 : int, n0 : int) =
| Real_Matrix of (matrixref (g0float tk, m0, n0),
                  int m1, int n1, int m0, int n0,
                  Matrix_Index_Map (m1, n1, m0, n0))
typedef Real_Matrix (tk : tkind, m1 : int, n1 : int) =
  [m0, n0 : pos] Real_Matrix (tk, m1, n1, m0, n0)
typedef Real_Vector (tk : tkind, m1 : int, n1 : int) =
  [m1 == 1 || n1 == 1] Real_Matrix (tk, m1, n1)
typedef Real_Row (tk : tkind, n1 : int) = Real_Vector (tk, 1, n1)
typedef Real_Column (tk : tkind, m1 : int) = Real_Vector (tk, m1, 1)

extern fn {tk : tkind}
Real_Matrix_make_elt :
  {m0, n0 : pos}
  (int m0, int n0, g0float tk) -< !wrt >
    Real_Matrix (tk, m0, n0, m0, n0)

extern fn {tk : tkind}
Real_Matrix_copy :
  {m1, n1 : pos}
  Real_Matrix (tk, m1, n1) -< !refwrt > Real_Matrix (tk, m1, n1)

extern fn {tk : tkind}
Real_Matrix_copy_to :
  {m1, n1 : pos}
  (Real_Matrix (tk, m1, n1),    (* destination *)
   Real_Matrix (tk, m1, n1)) -< !refwrt >
    void

extern fn {tk : tkind}
Real_Matrix_fill_with_elt :
  {m1, n1 : pos}
  (Real_Matrix (tk, m1, n1), g0float tk) -< !refwrt > void

extern fn {}
Real_Matrix_dimension :
  {tk : tkind}
  {m1, n1 : pos}
  Real_Matrix (tk, m1, n1) -<> @(int m1, int n1)

extern fn {tk : tkind}
Real_Matrix_get_at :
  {m1, n1 : pos}
  {i1, j1 : pos | i1 <= m1; j1 <= n1}
  (Real_Matrix (tk, m1, n1), int i1, int j1) -< !ref > g0float tk

extern fn {tk : tkind}
Real_Matrix_set_at :
  {m1, n1 : pos}
  {i1, j1 : pos | i1 <= m1; j1 <= n1}
  (Real_Matrix (tk, m1, n1), int i1, int j1, g0float tk) -< !refwrt >
    void

extern fn {}
Real_Matrix_apply_index_map :
  {tk : tkind}
  {m1, n1 : pos}
  {m0, n0 : pos}
  (Real_Matrix (tk, m0, n0), int m1, int n1,
   Matrix_Index_Map (m1, n1, m0, n0)) -<>
    Real_Matrix (tk, m1, n1)

extern fn {}
Real_Matrix_transpose :
  (* This is transposed INDEXING. It does NOT copy the data. *)
  {tk : tkind}
  {m1, n1 : pos}
  {m0, n0 : pos}
  Real_Matrix (tk, m1, n1, m0, n0) -<>
    Real_Matrix (tk, n1, m1, m0, n0)

extern fn {}
Real_Matrix_block :
  (* This is block (submatrix) INDEXING. It does NOT copy the data. *)
  {tk : tkind}
  {p0, p1 : pos | p0 <= p1}
  {q0, q1 : pos | q0 <= q1}
  {m1, n1 : pos | p1 <= m1; q1 <= n1}
  {m0, n0 : pos}
  (Real_Matrix (tk, m1, n1, m0, n0),
   int p0, int p1, int q0, int q1) -<>
    Real_Matrix (tk, p1 - p0 + 1, q1 - q0 + 1, m0, n0)

extern fn {tk : tkind}
Real_Matrix_unit_matrix :
  {m : pos}
  int m -< !refwrt > Real_Matrix (tk, m, m)

extern fn {tk : tkind}
Real_Matrix_unit_matrix_to :
  {m : pos}
  Real_Matrix (tk, m, m) -< !refwrt > void

extern fn {tk : tkind}
Real_Matrix_matrix_sum :
  {m, n : pos}
  (Real_Matrix (tk, m, n), Real_Matrix (tk, m, n)) -< !refwrt >
    Real_Matrix (tk, m, n)

extern fn {tk : tkind}
Real_Matrix_matrix_sum_to :
  {m, n : pos}
  (Real_Matrix (tk, m, n),      (* destination*)
   Real_Matrix (tk, m, n),
   Real_Matrix (tk, m, n)) -< !refwrt >
    void

extern fn {tk : tkind}
Real_Matrix_matrix_difference :
  {m, n : pos}
  (Real_Matrix (tk, m, n), Real_Matrix (tk, m, n)) -< !refwrt >
    Real_Matrix (tk, m, n)

extern fn {tk : tkind}
Real_Matrix_matrix_difference_to :
  {m, n : pos}
  (Real_Matrix (tk, m, n),      (* destination*)
   Real_Matrix (tk, m, n),
   Real_Matrix (tk, m, n)) -< !refwrt >
    void

extern fn {tk : tkind}
Real_Matrix_matrix_product :
  {m, n, p : pos}
  (Real_Matrix (tk, m, n), Real_Matrix (tk, n, p)) -< !refwrt >
    Real_Matrix (tk, m, p)

extern fn {tk : tkind}
Real_Matrix_matrix_product_to :
  {m, n, p : pos}
  (Real_Matrix (tk, m, p),      (* destination*)
   Real_Matrix (tk, m, n),
   Real_Matrix (tk, n, p)) -< !refwrt >
    void

extern fn {tk : tkind}
Real_Matrix_scalar_product :
  {m, n : pos}
  (Real_Matrix (tk, m, n), g0float tk) -< !refwrt >
    Real_Matrix (tk, m, n)

extern fn {tk : tkind}
Real_Matrix_scalar_product_2 :
  {m, n : pos}
  (g0float tk, Real_Matrix (tk, m, n)) -< !refwrt >
    Real_Matrix (tk, m, n)

extern fn {tk : tkind}
Real_Matrix_scalar_product_to :
  {m, n : pos}
  (Real_Matrix (tk, m, n),      (* destination*)
   Real_Matrix (tk, m, n), g0float tk) -< !refwrt > void

extern fn {tk : tkind}          (* Useful for debugging. *)
Real_Matrix_fprint :
  {m, n : pos}
  (FILEref, Real_Matrix (tk, m, n)) -<1> void

overload copy with Real_Matrix_copy
overload copy_to with Real_Matrix_copy_to
overload fill_with_elt with Real_Matrix_fill_with_elt
overload dimension with Real_Matrix_dimension
overload [] with Real_Matrix_get_at
overload [] with Real_Matrix_set_at
overload apply_index_map with Real_Matrix_apply_index_map
overload transpose with Real_Matrix_transpose
overload block with Real_Matrix_block
overload unit_matrix with Real_Matrix_unit_matrix
overload unit_matrix_to with Real_Matrix_unit_matrix_to
overload matrix_sum with Real_Matrix_matrix_sum
overload matrix_sum_to with Real_Matrix_matrix_sum_to
overload matrix_difference with Real_Matrix_matrix_difference
overload matrix_difference_to with Real_Matrix_matrix_difference_to
overload matrix_product with Real_Matrix_matrix_product
overload matrix_product_to with Real_Matrix_matrix_product_to
overload scalar_product with Real_Matrix_scalar_product
overload scalar_product with Real_Matrix_scalar_product_2
overload scalar_product_to with Real_Matrix_scalar_product_to
overload + with matrix_sum
overload - with matrix_difference
overload * with matrix_product
overload * with scalar_product

(*------------------------------------------------------------------*)
(* Implementation of the "little matrix library"                    *)

implement {tk}
Real_Matrix_make_elt (m0, n0, elt) =
  Real_Matrix (matrixref_make_elt<g0float tk> (i2sz m0, i2sz n0, elt),
               m0, n0, m0, n0, lam (i1, j1) => @(i1, j1))

implement {}
Real_Matrix_dimension A =
  case+ A of Real_Matrix (_, m1, n1, _, _, _) => @(m1, n1)

implement {tk}
Real_Matrix_get_at (A, i1, j1) =
  let
    val+ Real_Matrix (storage, _, _, _, n0, index_map) = A
    val @(i0, j0) = index_map (i1, j1)
  in
    matrixref_get_at<g0float tk> (storage, pred i0, n0, pred j0)
  end

implement {tk}
Real_Matrix_set_at (A, i1, j1, x) =
  let
    val+ Real_Matrix (storage, _, _, _, n0, index_map) = A
    val @(i0, j0) = index_map (i1, j1)
  in
    matrixref_set_at<g0float tk> (storage, pred i0, n0, pred j0, x)
  end

implement {}
Real_Matrix_apply_index_map (A, m1, n1, index_map) =
  (* This is not the most efficient way to acquire new indexing, but
     it will work. It requires three closures, instead of the two
     needed by our implementations of "transpose" and "block". *)
  let
    val+ Real_Matrix (storage, m1a, n1a, m0, n0, index_map_1a) = A
  in
    Real_Matrix (storage, m1, n1, m0, n0,
                 lam (i1, j1) =>
                   index_map_1a (i1a, j1a) where
                     { val @(i1a, j1a) = index_map (i1, j1) })
  end

implement {}
Real_Matrix_transpose A =
  let
    val+ Real_Matrix (storage, m1, n1, m0, n0, index_map) = A
  in
    Real_Matrix (storage, n1, m1, m0, n0,
                 lam (i1, j1) => index_map (j1, i1))
  end

implement {}
Real_Matrix_block (A, p0, p1, q0, q1) =
  let
    val+ Real_Matrix (storage, m1, n1, m0, n0, index_map) = A
  in
    Real_Matrix (storage, succ (p1 - p0), succ (q1 - q0), m0, n0,
                 lam (i1, j1) =>
                  index_map (p0 + pred i1, q0 + pred j1))
  end

implement {tk}
Real_Matrix_copy A =
  let
    val @(m1, n1) = dimension A
    val C = Real_Matrix_make_elt<tk> (m1, n1, A[1, 1])
    val () = copy_to<tk> (C, A)
  in
    C
  end

implement {tk}
Real_Matrix_copy_to (Dst, Src) =
  let
    val @(m1, n1) = dimension Src
    prval [m1 : int] EQINT () = eqint_make_gint m1
    prval [n1 : int] EQINT () = eqint_make_gint n1

    var i : intGte 1
  in
    for* {i : pos | i <= m1 + 1} .<(m1 + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m1; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= n1 + 1} .<(n1 + 1) - j>.
               (j : int j) =>
            (j := 1; j <> succ n1; j := succ j)
              Dst[i, j] := Src[i, j]
        end
  end

implement {tk}
Real_Matrix_fill_with_elt (A, elt) =
  let
    val @(m1, n1) = dimension A
    prval [m1 : int] EQINT () = eqint_make_gint m1
    prval [n1 : int] EQINT () = eqint_make_gint n1

    var i : intGte 1
  in
    for* {i : pos | i <= m1 + 1} .<(m1 + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m1; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= n1 + 1} .<(n1 + 1) - j>.
               (j : int j) =>
            (j := 1; j <> succ n1; j := succ j)
              A[i, j] := elt
        end
  end

implement {tk}
Real_Matrix_unit_matrix {m} m =
  let
    val A = Real_Matrix_make_elt<tk> (m, m, Zero)
    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        A[i, i] := One;
    A
  end

implement {tk}
Real_Matrix_unit_matrix_to A =
  let
    val @(m, _) = dimension A
    prval [m : int] EQINT () = eqint_make_gint m

    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= m + 1} .<(m + 1) - j>.
               (j : int j) =>
               (j := 1; j <> succ m; j := succ j)
            A[i, j] := (if i = j then One else Zero)
        end
  end

implement {tk}
Real_Matrix_matrix_sum (A, B) =
  let
    val @(m, n) = dimension A
    val C = Real_Matrix_make_elt<tk> (m, n, NAN)
    val () = matrix_sum_to<tk> (C, A, B)
  in
    C
  end

implement {tk}
Real_Matrix_matrix_sum_to (C, A, B) =
  let
    val @(m, n) = dimension A
    prval [m : int] EQINT () = eqint_make_gint m
    prval [n : int] EQINT () = eqint_make_gint n

    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= n + 1} .<(n + 1) - j>.
               (j : int j) =>
            (j := 1; j <> succ n; j := succ j)
              C[i, j] := A[i, j] + B[i, j]
        end
  end

implement {tk}
Real_Matrix_matrix_difference (A, B) =
  let
    val @(m, n) = dimension A
    val C = Real_Matrix_make_elt<tk> (m, n, NAN)
    val () = matrix_difference_to<tk> (C, A, B)
  in
    C
  end

implement {tk}
Real_Matrix_matrix_difference_to (C, A, B) =
  let
    val @(m, n) = dimension A
    prval [m : int] EQINT () = eqint_make_gint m
    prval [n : int] EQINT () = eqint_make_gint n

    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= n + 1} .<(n + 1) - j>.
               (j : int j) =>
            (j := 1; j <> succ n; j := succ j)
              C[i, j] := A[i, j] - B[i, j]
        end
  end

implement {tk}
Real_Matrix_matrix_product (A, B) =
  let
    val @(m, n) = dimension A and @(_, p) = dimension B
    val C = Real_Matrix_make_elt<tk> (m, p, NAN)
    val () = matrix_product_to<tk> (C, A, B)
  in
    C
  end

implement {tk}
Real_Matrix_matrix_product_to (C, A, B) =
  let
    val @(m, n) = dimension A and @(_, p) = dimension B
    prval [m : int] EQINT () = eqint_make_gint m
    prval [n : int] EQINT () = eqint_make_gint n
    prval [p : int] EQINT () = eqint_make_gint p

    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        let
          var k : intGte 1
        in
          for* {k : pos | k <= p + 1} .<(p + 1) - k>.
               (k : int k) =>
            (k := 1; k <> succ p; k := succ k)
              let
                var j : intGte 1
              in
                C[i, k] := A[i, 1] * B[1, k];
                for* {j : pos | j <= n + 1} .<(n + 1) - j>.
                     (j : int j) =>
                  (j := 2; j <> succ n; j := succ j)
                    C[i, k] :=
                      multiply_and_add (A[i, j], B[j, k], C[i, k])
              end
        end
  end

implement {tk}
Real_Matrix_scalar_product (A, r) =
  let
    val @(m, n) = dimension A
    val C = Real_Matrix_make_elt<tk> (m, n, NAN)
    val () = scalar_product_to<tk> (C, A, r)
  in
    C
  end

implement {tk}
Real_Matrix_scalar_product_2 (r, A) =
  Real_Matrix_scalar_product<tk> (A, r)

implement {tk}
Real_Matrix_scalar_product_to (C, A, r) =
  let
    val @(m, n) = dimension A
    prval [m : int] EQINT () = eqint_make_gint m
    prval [n : int] EQINT () = eqint_make_gint n

    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= n + 1} .<(n + 1) - j>.
               (j : int j) =>
            (j := 1; j <> succ n; j := succ j)
              C[i, j] := A[i, j] * r
        end
  end

implement {tk}
Real_Matrix_fprint {m, n} (outf, A) =
  let
    val @(m, n) = dimension A
    var i : intGte 1
  in
    for* {i : pos | i <= m + 1} .<(m + 1) - i>.
         (i : int i) =>
      (i := 1; i <> succ m; i := succ i)
        let
          var j : intGte 1
        in
          for* {j : pos | j <= n + 1} .<(n + 1) - j>.
               (j : int j) =>
            (j := 1; j <> succ n; j := succ j)
              let
                typedef FILEstar = $extype"FILE *"
                extern castfn FILEref2star : FILEref -<> FILEstar
                val _ = $extfcall (int, "fprintf", FILEref2star outf,
                                   "%16.6g", A[i, j])
              in
              end;
          fprintln! (outf)
        end
  end

(*------------------------------------------------------------------*)
(* Gauss-Jordan inversion                                           *)

extern fn {tk : tkind}
Real_Matrix_inverse_by_gauss_jordan :
  {n : pos}
  Real_Matrix (tk, n, n) -< !refwrt > Option (Real_Matrix (tk, n, n))

#if DO_THINGS_THAT_DO_NOT_NEED_TO_BE_DONE #then
  macdef do_needless_things = true
#else
  macdef do_needless_things = false
#endif

fn {tk : tkind}
needlessly_set_to_value
          {n    : pos}
          {i, j : pos | i <= n; j <= n}
          (A : Real_Matrix (tk, n, n),
           i : int i,
           j : int j,
           x : g0float tk) :<!refwrt> void =
  if do_needless_things then A[i, j] := x

implement {tk}
Real_Matrix_inverse_by_gauss_jordan {n} A =
  let
    val @(n, _) = dimension A
    typedef one_to_n = intBtwe (1, n)

    (* Partial pivoting. *)
    implement
    array_tabulate$fopr<one_to_n> i =
      let
        val i = g1ofg0 (sz2i (succ i))
        val () = assertloc ((1 <= i) * (i <= n))
      in
        i
      end
    val rows_permutation =
      $effmask_all arrayref_tabulate<one_to_n> (i2sz n)
    fn
    index_map : Matrix_Index_Map (n, n, n, n) =
      lam (i1, j1) => $effmask_ref
        (@(i0, j1) where { val i0 = rows_permutation[i1 - 1] })

    val A = apply_index_map (copy<tk> A, n, n, index_map)
    and B = apply_index_map (unit_matrix<tk> n, n, n, index_map)

    fn {}
    exchange_rows (i1 : one_to_n,
                   i2 : one_to_n) :<!refwrt> void =
      if i1 <> i2 then
        let
          val k1 = rows_permutation[pred i1]
          and k2 = rows_permutation[pred i2]
        in
          rows_permutation[pred i1] := k2;
          rows_permutation[pred i2] := k1
        end

    fn {}
    normalize_pivot_row (j : one_to_n) :<!refwrt> void =
      let
        prval [j : int] EQINT () = eqint_make_gint j
        val pivot_val = A[j, j]
        var k : intGte 1
      in
        needlessly_set_to_value (A, j, j, One);
        for* {k : int | j + 1 <= k; k <= n + 1} .<(n + 1) - k>.
             (k : int k) =>
          (k := succ j; k <> succ n; k := succ k)
            A[j, k] := A[j, k] / pivot_val;
        for* {k : int | 1 <= k; k <= n + 1} .<(n + 1) - k>.
             (k : int k) =>
          (k := 1; k <> succ n; k := succ k)
            B[j, k] := B[j, k] / pivot_val;
      end

    fn
    subtract_normalized_pivot_row
              (i : one_to_n, j : one_to_n) :<!refwrt> void =
      let
        prval [j : int] EQINT () = eqint_make_gint j
        val lead_val = A[i, j]
      in
        if lead_val <> Zero then
          let
            val factor = ~lead_val
            var k : intGte 1
          in
            needlessly_set_to_value (A, i, j, Zero);
            for* {k : int | j + 1 <= k; k <= n + 1} .<(n + 1) - k>.
                 (k : int k) =>
              (k := succ j; k <> succ n; k := succ k)
                A[i, k] :=
                  multiply_and_add (A[j, k], factor, A[i, k]);
            for* {k : int | 1 <= k; k <= n + 1} .<(n + 1) - k>.
                 (k : int k) =>
              (k := 1; k <> succ n; k := succ k)
                B[i, k] :=
                  multiply_and_add (B[j, k], factor, B[i, k])
          end
      end

    fun
    main_loop {j       : pos | j <= n + 1} .<(n + 1) - j>.
              (j       : int j,
               success : &bool? >> bool) :<!refwrt> void =
      if j = succ n then
        success := true
      else
        let
          fun
          select_pivot {i : int | j <= i; i <= n + 1}
                       .<(n + 1) - i>.
                       (i         : int i,
                        max_abs   : g0float tk,
                        i_max_abs : intBtwe (j - 1, n))
              :<!ref> intBtwe (j - 1, n) =
            if i = succ n then
              i_max_abs
            else
              let
                val abs_aij = abs A[i, j]
              in
                if abs_aij > max_abs then
                  select_pivot (succ i, abs_aij, i)
                else
                  select_pivot (succ i, max_abs, i_max_abs)
              end

          val i_pivot = select_pivot (j, Zero, pred j)
          prval [i_pivot : int] EQINT () = eqint_make_gint i_pivot
        in
          if i_pivot = pred j then
            success := false
          else
            let
              var i : intGte 1
            in
              exchange_rows (i_pivot, j);
              normalize_pivot_row (j);
              for* {i : int | 1 <= i; i <= j}
                   .<j - i>.
                   (i : int i) =>
                (i := 1; i <> j; i := succ i)
                  subtract_normalized_pivot_row (i, j);
              for* {i : int | j + 1 <= i; i <= n + 1}
                   .<(n + 1) - i>.
                   (i : int i) =>
                (i := succ j; i <> succ n; i := succ i)
                  subtract_normalized_pivot_row (i, j);
              main_loop (succ j, success)
            end
        end

    var success : bool
    val () = main_loop (1, success)
  in
    (* The returned B will "contain" the rows_permutation array and
       some extra closures. If you want a "clean" B, use "val B = copy
       B" or some such. *)
    if success then Some B else None ()
  end

overload inverse_by_gauss_jordan with
  Real_Matrix_inverse_by_gauss_jordan

(*------------------------------------------------------------------*)

fn {tk : tkind}
fprint_matrix_and_its_inverse
          {n    : pos}
          (outf : FILEref,
           A    : Real_Matrix (tk, n, n)) : void =
  let
    typedef FILEstar = $extype"FILE *"
    extern castfn FILEref2star : FILEref -<> FILEstar

    macdef fmt = "%8.4f"
    macdef left_bracket = "["
    macdef right_bracket = "  ]"
    macdef product = "  ✕  "
    macdef equals = "  =  "
    macdef spacing = "     "
    macdef msg_for_singular = "  appears to be singular"

    macdef print_num (x) =
      ignoret ($extfcall (int, "fprintf", FILEref2star outf,
                          fmt, ,(x)))

    val @(n, _) = dimension A
  in
    case+ inverse_by_gauss_jordan<tk> (A) of
    | None () =>
      let
        var i : intGte 1
      in
        for* {i : pos | i <= n + 1} .<(n + 1) - i>.
             (i : int i) =>
          (i := 1; i <> succ n; i := succ i)
            let
              var j : intGte 1
            in
              fprint! (outf, left_bracket);
              for* {j : pos | j <= n + 1} .<(n + 1) - j>.
                   (j : int j) =>
                (j := 1; j <> succ n; j := succ j)
                  print_num A[i, j];
              fprint! (outf, right_bracket);
              if pred i = half n then
                fprint! (outf, msg_for_singular);
              fprintln! (outf)
            end
      end
    | Some B =>
      let
        val AB = A * B
        var i : intGte 1
      in
        for* {i : pos | i <= n + 1} .<(n + 1) - i>.
             (i : int i) =>
          (i := 1; i <> succ n; i := succ i)
            let
              var j : intGte 1
            in
              fprint! (outf, left_bracket);
              for* {j : pos | j <= n + 1} .<(n + 1) - j>.
                   (j : int j) =>
                (j := 1; j <> succ n; j := succ j)
                  print_num A[i, j];
              fprint! (outf, right_bracket);
              if pred i = half n then
                fprint! (outf, product)
              else
                fprint! (outf, spacing);
              fprint! (outf, left_bracket);
              for* {j : pos | j <= n + 1} .<(n + 1) - j>.
                   (j : int j) =>
                (j := 1; j <> succ n; j := succ j)
                  print_num B[i, j];
              fprint! (outf, right_bracket);
              if pred i = half n then
                fprint! (outf, equals)
              else
                fprint! (outf, spacing);
              fprint! (outf, left_bracket);
              for* {j : pos | j <= n + 1} .<(n + 1) - j>.
                   (j : int j) =>
                (j := 1; j <> succ n; j := succ j)
                  print_num AB[i, j];
              fprint! (outf, right_bracket);
              fprintln! (outf)
            end
      end
  end

fn {tk : tkind}           (* We like Fortran, so COLUMN major here. *)
column_major_list_to_square_matrix
          {n   : pos}
          (n   : int n,
           lst : list (g0float tk, n * n))
    : Real_Matrix (tk, n, n) =
  let
    #define :: list_cons
    prval () = mul_gte_gte_gte {n, n} ()
    val A = Real_Matrix_make_elt (n, n, NAN)
    val lstref : ref (List0 (g0float tk)) = ref lst
    var j : intGte 1
  in
    for* {j : pos | j <= n + 1} .<(n + 1) - j>.
         (j : int j) =>
      (j := 1; j <> succ n; j := succ j)
        let
          var i : intGte 1
        in
          for* {i : pos | i <= n + 1} .<(n + 1) - i>.
               (i : int i) =>
            (i := 1; i <> succ n; i := succ i)
              case- !lstref of
              | hd :: tl =>
                begin
                  A[i, j] := hd;
                  !lstref := tl
                end
        end;
    A
  end

macdef separator = "\n"

fn
print_example
          {n   : pos}
          (n   : int n,
           lst : list (double, n * n)) : void =
  let
    val A = column_major_list_to_square_matrix (n, lst)
  in
    print! separator;
    fprint_matrix_and_its_inverse<dblknd> (stdout_ref, A)
  end

implement
main0 () =
  begin
    println! ("\n(The examples are printed here after ",
              "rounding to 4 decimal places.)");

    print_example (1, $list (0.0));
    print_example (1, $list (4.0));
    print_example (2, $list (1.0, 1.0, 1.0, 1.0));
    print_example (2, $list (2.0, 0.0, 0.0, 2.0));
    print_example (2, $list (1.0, 2.0, 3.0, 4.0));

    println! ("\nAn orthonormal matrix of rank 5 ",
              "(its inverse will equal its transpose):");
    print_example
      (5, $list (~0.08006407690254358, ~0.1884825001438613, 
                 ~0.8954592676839166, 0.2715238629598956,
                 ~0.2872134789517761, ~0.5604485383178048,
                 ~0.4587939881550582, ~0.01946650581921516, 
                 0.1609525690123523, 0.6701647842208114, 
                 ~0.2401922307076307, 0.1517054269450593,
                 ~0.3178281430868086, ~0.8979459113320678,
                 0.1094146586482966, ~0.3202563076101743,
                 ~0.6104994151001173, 0.2983616372675922, 
                 ~0.1997416922923375, ~0.6291342872277008, 
                 ~0.7205766921228921, 0.5985468663105067,
                 0.08797363206760908, 0.2327347396799102,
                 ~0.2461829819586655));

    (* The following matrix may or may not be detected as singular,
       depending on how rounding is done. *)
    (*
    print_example
      (5, $list (~0.08006407690254358, ~0.1884825001438613, 
                 ~0.8954592676839166, 0.2715238629598956,
                 ~0.2872134789517761, ~0.5604485383178048,
                 ~0.4587939881550582, ~0.01946650581921516, 
                 0.1609525690123523, 0.6701647842208114, 
                 ~0.2401922307076307, 0.1517054269450593,
                 ~0.3178281430868086, ~0.8979459113320678,
                 0.1094146586482966, ~0.3202563076101743,
                 ~0.6104994151001173, 0.2983616372675922, 
                 ~0.1997416922923375, ~0.6291342872277008, 
                 5.0 * ~0.08006407690254358,
                 5.0 * ~0.1884825001438613, 
                 5.0 * ~0.8954592676839166,
                 5.0 * 0.2715238629598956,
                 5.0 * ~0.2872134789517761));
    *)

    println! ("\nA 5✕5 singular matrix:");
    print_example
      (5, $list (1.0, 1.0, 2.0, 2.0, 2.0,
                 5.0, 1.0, 2.0, 2.0, 2.0,
                 3.0, 2.0, 3.0, 5.0, 7.0,
                 1.0, 1.0, 2.0, 2.0, 2.0,
                 4.0, 1.0, 2.0, 2.0, 2.0));

    print! separator
  end

(*------------------------------------------------------------------*)
</syntaxhighlight>

{{out}}

(The following settings on my computer ''will'' result in GCC generating fused-multiply-and-add instructions [and thus the ''-lm'' flag actually is not needed]. See comments in the program for a discussion of the matter.)

<pre style="font-size:90%">$ patscc -std=gnu2x -g -O3 -march=native -DATS_MEMALLOC_GCBDW gauss_jordan_task.dats -lgc -lm && ./a.out

(The examples are printed here after rounding to 4 decimal places.)

[  0.0000  ]  appears to be singular

[  4.0000  ]  ✕  [  0.2500  ]  =  [  1.0000  ]

[  1.0000  1.0000  ]
[  1.0000  1.0000  ]  appears to be singular

[  2.0000  0.0000  ]     [  0.5000  0.0000  ]     [  1.0000  0.0000  ]
[  0.0000  2.0000  ]  ✕  [  0.0000  0.5000  ]  =  [  0.0000  1.0000  ]

[  1.0000  3.0000  ]     [ -2.0000  1.5000  ]     [  1.0000  0.0000  ]
[  2.0000  4.0000  ]  ✕  [  1.0000 -0.5000  ]  =  [  0.0000  1.0000  ]

An orthonormal matrix of rank 5 (its inverse will equal its transpose):

[ -0.0801 -0.5604 -0.2402 -0.3203 -0.7206  ]     [ -0.0801 -0.1885 -0.8955  0.2715 -0.2872  ]     [  1.0000 -0.0000 -0.0000  0.0000  0.0000  ]
[ -0.1885 -0.4588  0.1517 -0.6105  0.5985  ]     [ -0.5604 -0.4588 -0.0195  0.1610  0.6702  ]     [  0.0000  1.0000 -0.0000  0.0000  0.0000  ]
[ -0.8955 -0.0195 -0.3178  0.2984  0.0880  ]  ✕  [ -0.2402  0.1517 -0.3178 -0.8979  0.1094  ]  =  [  0.0000 -0.0000  1.0000  0.0000  0.0000  ]
[  0.2715  0.1610 -0.8979 -0.1997  0.2327  ]     [ -0.3203 -0.6105  0.2984 -0.1997 -0.6291  ]     [ -0.0000  0.0000 -0.0000  1.0000 -0.0000  ]
[ -0.2872  0.6702  0.1094 -0.6291 -0.2462  ]     [ -0.7206  0.5985  0.0880  0.2327 -0.2462  ]     [ -0.0000  0.0000 -0.0000  0.0000  1.0000  ]

A 5✕5 singular matrix:

[  1.0000  5.0000  3.0000  1.0000  4.0000  ]
[  1.0000  1.0000  2.0000  1.0000  1.0000  ]
[  2.0000  2.0000  3.0000  2.0000  2.0000  ]  appears to be singular
[  2.0000  2.0000  5.0000  2.0000  2.0000  ]
[  2.0000  2.0000  7.0000  2.0000  2.0000  ]

</pre>

=={{header|BASIC}}==
==={{header|FreeBASIC}}===
The include statements use code from [[Matrix multiplication#FreeBASIC]], which contains the Matrix type used here, and [[Reduced row echelon form#FreeBASIC]] which contains the function for reducing a matrix to row-echelon form. Make sure to take out all the print statements first!

<syntaxhighlight lang="freebasic">#include once "matmult.bas"
#include once "rowech.bas"

function matinv( A as Matrix ) as Matrix
    dim ret as Matrix, working as Matrix
    dim as uinteger i, j
    if not ubound( A.m, 1 ) = ubound( A.m, 2 ) then return ret
    dim as integer n = ubound(A.m, 1)
    redim ret.m( n, n )
    working = Matrix( n+1, 2*(n+1) )
    for i = 0 to n
        for j = 0 to n
            working.m(i, j) = A.m(i, j)
        next j
        working.m(i, i+n+1) = 1
    next i
    working = rowech(working)
    for i = 0 to n
        for j = 0 to n
            ret.m(i, j) = working.m(i, j+n+1)
        next j
    next i
    return ret
end function

dim as integer i, j
dim as Matrix M = Matrix(3,3)
for i = 0 to 2
    for j = 0 to 2
        M.m(i, j) = 1 + i*i + 3*j + (i+j) mod 2   'just some arbitrary values
        print M.m(i, j),
    next j
    print
next i
print
M = matinv(M)
for i = 0 to 2
    for j = 0 to 2
        print M.m(i, j),
    next j
    print
next i</syntaxhighlight>
{{out}}
<pre>
 1             5             7            
 3             5             9            
 5             9             11           

-0.5416666666666667          0.1666666666666667          0.2083333333333333         
 0.2499999999999999         -0.5           0.25         
 0.0416666666666667          0.3333333333333333         -0.2083333333333333
</pre>

==={{header|VBA}}===
{{trans|Phix}}
Uses ToReducedRowEchelonForm() from [[Reduced_row_echelon_form#VBA]]<syntaxhighlight lang="vb">Private Function inverse(mat As Variant) As Variant
    Dim len_ As Integer: len_ = UBound(mat)
    Dim tmp() As Variant
    ReDim tmp(2 * len_ + 1)
    Dim aug As Variant
    ReDim aug(len_)
    For i = 0 To len_
        If UBound(mat(i)) <> len_ Then Debug.Print 9 / 0 '-- "Not a square matrix"
        aug(i) = tmp
        For j = 0 To len_
            aug(i)(j) = mat(i)(j)
        Next j
        '-- augment by identity matrix to right
        aug(i)(i + len_ + 1) = 1
    Next i
    aug = ToReducedRowEchelonForm(aug)
    Dim inv As Variant
    inv = mat
    '-- remove identity matrix to left
    For i = 0 To len_
        For j = len_ + 1 To 2 * len_ + 1
            inv(i)(j - len_ - 1) = aug(i)(j)
        Next j
    Next i
    inverse = inv
End Function
 
Public Sub main()
    Dim test As Variant
    test = inverse(Array( _
        Array(2, -1, 0), _
        Array(-1, 2, -1), _
        Array(0, -1, 2)))
    For i = LBound(test) To UBound(test)
        For j = LBound(test(0)) To UBound(test(0))
            Debug.Print test(i)(j),
        Next j
    Debug.Print
    Next i
End Sub</syntaxhighlight>{{out}}
<pre> 0,75          0,5           0,25         
 0,5           1             0,5          
 0,25          0,5           0,75 </pre>

==={{header|VBScript}}===
{{trans|Rexx}}
To run in console mode with cscript.
<syntaxhighlight lang="vb">' Gauss-Jordan matrix inversion - VBScript - 22/01/2021
Option Explicit

Function rref(m)
    Dim r, c, i, n, div, wrk
    n=UBound(m)
    For r = 1 To n  'row
        div = m(r,r)
        If div <> 0 Then
            For c = 1 To n*2  'col
                m(r,c) = m(r,c) / div
            Next 'c
        Else
            WScript.Echo "inversion impossible!"
            WScript.Quit
        End If
        For i = 1 To n  'row
            If i <> r Then
                wrk = m(i,r)
                For c = 1 To n*2
                    m(i,c) = m(i,c) - wrk * m(r,c)
                Next' c
            End If
        Next 'i
    Next 'r
    rref = m
End Function 'rref

Function inverse(mat)
    Dim i, j, aug, inv, n
    n = UBound(mat)
    ReDim inv(n,n), aug(n,2*n)
    For i = 1 To n
        For j = 1 To n
            aug(i,j) = mat(i,j)
        Next 'j
        aug(i,i+n) = 1
    Next 'i
    aug = rref(aug)
    For i = 1 To n
        For j = n+1 To 2*n
            inv(i,j-n) = aug(i,j)
        Next 'j
    Next 'i
    inverse = inv
End Function 'inverse

Sub wload(m)
    Dim i, j, k
    k = -1
    For i = 1 To n
        For j = 1 To n
            k = k + 1
            m(i,j) = w(k)
        Next 'j
    Next 'i
End Sub 'wload 

Sub show(c, m, t)
    Dim i, j, buf
    buf = "Matrix " & c &"=" & vbCrlf & vbCrlf
    For i = 1 To n
        For j = 1 To n
            If t="fixed" Then
                buf = buf & FormatNumber(m(i,j),6,,,0) & "  "
            Else
                buf = buf & m(i,j) & "  "
            End If
        Next 'j
        buf=buf & vbCrlf
    Next 'i
    WScript.Echo buf
End Sub 'show 

Dim n, a, b, c, w
w = Array( _
	2,  1,  1,  4, _
	0, -1,  0, -1, _
	1,  0, -2,  4, _
	4,  1,  2,  2)
n=Sqr(UBound(w)+1)
ReDim a(n,n), b(n,n), c(n,n)
wload a
show "a", a, "simple"
b = inverse(a)
show "b", b, "fixed"
c = inverse(b)
show "c", c, "fixed"
</syntaxhighlight>
{{out}}
<pre>
Matrix a=

2  1  1  4
0  -1  0  -1
1  0  -2  4
4  1  2  2

Matrix b=

-0,400000  0,000000  0,200000  0,400000
-0,400000  -1,200000  0,000000  0,200000
0,600000  0,400000  -0,400000  -0,200000
0,400000  0,200000  -0,000000  -0,200000

Matrix c=

2,000000  1,000000  1,000000  4,000000
0,000000  -1,000000  0,000000  -1,000000
1,000000  0,000000  -2,000000  4,000000
4,000000  1,000000  2,000000  2,000000

</pre>

=={{header|C}}==
{{trans|Fortran}}
<syntaxhighlight lang="c">/*----------------------------------------------------------------------
gjinv - Invert a matrix, Gauss-Jordan algorithm
A is destroyed.  Returns 1 for a singular matrix.

___Name_____Type______In/Out____Description_____________________________
   a[n*n]   double*   In        An N by N matrix
   n        int       In        Order of matrix
   b[n*n]   double*   Out       Inverse of A
----------------------------------------------------------------------*/
#include <math.h>
int gjinv (double *a, int n, double *b)
{
	int i, j, k, p;
	double f, g, tol;
	if (n < 1) return -1;  /* Function Body */
	f = 0.;  /* Frobenius norm of a */
	for (i = 0; i < n; ++i) {
		for (j = 0; j < n; ++j) {
			g = a[j+i*n];
			f += g * g;
		}
	}
	f = sqrt(f);
	tol = f * 2.2204460492503131e-016;
	for (i = 0; i < n; ++i) {  /* Set b to identity matrix. */
		for (j = 0; j < n; ++j) {
			b[j+i*n] = (i == j) ? 1. : 0.;
		}
	}
	for (k = 0; k < n; ++k) {  /* Main loop */
		f = fabs(a[k+k*n]);  /* Find pivot. */
		p = k;
		for (i = k+1; i < n; ++i) {
			g = fabs(a[k+i*n]);
			if (g > f) {
				f = g;
				p = i;
			}
		}
		if (f < tol) return 1;  /* Matrix is singular. */
		if (p != k) {  /* Swap rows. */
			for (j = k; j < n; ++j) {
				f = a[j+k*n];
				a[j+k*n] = a[j+p*n];
				a[j+p*n] = f;
			}
			for (j = 0; j < n; ++j) {
				f = b[j+k*n];
				b[j+k*n] = b[j+p*n];
				b[j+p*n] = f;
			}
		}
		f = 1. / a[k+k*n];  /* Scale row so pivot is 1. */
		for (j = k; j < n; ++j) a[j+k*n] *= f;
		for (j = 0; j < n; ++j) b[j+k*n] *= f;
		for (i = 0; i < n; ++i) {  /* Subtract to get zeros. */
			if (i == k) continue;
			f = a[k+i*n];
			for (j = k; j < n; ++j) a[j+i*n] -= a[j+k*n] * f;
			for (j = 0; j < n; ++j) b[j+i*n] -= b[j+k*n] * f;
		}
	}
	return 0;
} /* end of gjinv */</syntaxhighlight>
Test program:
<syntaxhighlight lang="c">/* 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 */</syntaxhighlight>
Output is the same as the Fortran example.

=={{header|C sharp|C#}}==
<syntaxhighlight lang="csharp">
using System;

namespace Rosetta
{
    internal class Vector
    {
        private double[] b;
        internal readonly int rows;

        internal Vector(int rows)
        {
            this.rows = rows;
            b = new double[rows];
        }

        internal Vector(double[] initArray)
        {
            b = (double[])initArray.Clone();
            rows = b.Length;
        }

        internal Vector Clone()
        {
            Vector v = new Vector(b);
            return v;
        }

        internal double this[int row]
        {
            get { return b[row]; }
            set { b[row] = value; }
        }

        internal void SwapRows(int r1, int r2)
        {
            if (r1 == r2) return;
            double tmp = b[r1];
            b[r1] = b[r2];
            b[r2] = tmp;
        }

        internal double norm(double[] weights)
        {
            double sum = 0;
            for (int i = 0; i < rows; i++)
            {
                double d = b[i] * weights[i];
                sum +=  d*d;
            }
            return Math.Sqrt(sum);
        }

        internal void print()
        {
            for (int i = 0; i < rows; i++)
                Console.WriteLine(b[i]);
            Console.WriteLine();
        }

        public static Vector operator-(Vector lhs, Vector rhs)
        {
            Vector v = new Vector(lhs.rows);
            for (int i = 0; i < lhs.rows; i++)
                v[i] = lhs[i] - rhs[i];
            return v;
        }
    }

    class Matrix
    {
        private double[] b;
        internal readonly int rows, cols;

        internal Matrix(int rows, int cols)
        {
            this.rows = rows;
            this.cols = cols;
            b = new double[rows * cols];            
        }

        internal Matrix(int size)
        {
            this.rows = size;
            this.cols = size;
            b = new double[rows * cols];
            for (int i = 0; i < size; i++)
                this[i, i] = 1;
        }

        internal Matrix(int rows, int cols, double[] initArray)
        {
            this.rows = rows;
            this.cols = cols;
            b = (double[])initArray.Clone();
            if (b.Length != rows * cols) throw new Exception("bad init array");
        }

        internal double this[int row, int col]
        {
            get { return b[row * cols + col]; }
            set { b[row * cols + col] = value; }
        }        
        
        public static Vector operator*(Matrix lhs, Vector rhs)
        {
            if (lhs.cols != rhs.rows) throw new Exception("I can't multiply matrix by vector");
            Vector v = new Vector(lhs.rows);
            for (int i = 0; i < lhs.rows; i++)
            {
                double sum = 0;
                for (int j = 0; j < rhs.rows; j++)
                    sum += lhs[i,j]*rhs[j];
                v[i] = sum;
            }
            return v;
        }

        internal void SwapRows(int r1, int r2)
        {
            if (r1 == r2) return;
            int firstR1 = r1 * cols;
            int firstR2 = r2 * cols;
            for (int i = 0; i < cols; i++)
            {
                double tmp = b[firstR1 + i];
                b[firstR1 + i] = b[firstR2 + i];
                b[firstR2 + i] = tmp;
            }
        }

        //with partial pivot
        internal bool InvPartial()
        {
            const double Eps = 1e-12;
            if (rows != cols) throw new Exception("rows != cols for Inv");
            Matrix M = new Matrix(rows); //unitary
            for (int diag = 0; diag < rows; diag++)
            {
                int max_row = diag;
                double max_val = Math.Abs(this[diag, diag]);
                double d;
                for (int row = diag + 1; row < rows; row++)
                    if ((d = Math.Abs(this[row, diag])) > max_val)
                    {
                        max_row = row;
                        max_val = d;
                    }
                if (max_val <= Eps) return false;
                SwapRows(diag, max_row);
                M.SwapRows(diag, max_row);
                double invd = 1 / this[diag, diag];
                for (int col = diag; col < cols; col++)
                {
                    this[diag, col] *= invd;
                }
                for (int col = 0; col < cols; col++)
                {
                    M[diag, col] *= invd;
                }
                for (int row = 0; row < rows; row++)
                {
                    d = this[row, diag];
                    if (row != diag)
                    {
                        for (int col = diag; col < this.cols; col++)
                        {
                            this[row, col] -= d * this[diag, col];
                        }
                        for (int col = 0; col < this.cols; col++)
                        {
                            M[row, col] -= d * M[diag, col];
                        }
                    }
                }
            }
            b = M.b;
            return true;
        }

        internal void print()
        {
            for (int i = 0; i < rows; i++)
            {
                for (int j = 0; j < cols; j++)
                    Console.Write(this[i,j].ToString()+"  ");
                Console.WriteLine();
            }
        }
    }
}
</syntaxhighlight>
<syntaxhighlight lang="csharp">
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();
        }
    }
}
</syntaxhighlight>
{{out}}<pre>
-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
</pre>

=={{header|C++}}==
<syntaxhighlight lang="cpp">#include <algorithm>
#include <cassert>
#include <iomanip>
#include <iostream>
#include <vector>

template <typename scalar_type> class matrix {
public:
    matrix(size_t rows, size_t columns)
        : rows_(rows), columns_(columns), elements_(rows * columns) {}

    matrix(size_t rows, size_t columns, scalar_type value)
        : rows_(rows), columns_(columns), elements_(rows * columns, value) {}

    matrix(size_t rows, size_t columns, std::initializer_list<scalar_type> values)
        : rows_(rows), columns_(columns), elements_(rows * columns) {
        assert(values.size() <= rows_ * columns_);
        std::copy(values.begin(), values.end(), elements_.begin());
    }

    size_t rows() const { return rows_; }
    size_t columns() const { return columns_; }

    scalar_type& operator()(size_t row, size_t column) {
        assert(row < rows_);
        assert(column < columns_);
        return elements_[row * columns_ + column];
    }
    const scalar_type& operator()(size_t row, size_t column) const {
        assert(row < rows_);
        assert(column < columns_);
        return elements_[row * columns_ + column];
    }
private:
    size_t rows_;
    size_t columns_;
    std::vector<scalar_type> elements_;
};

template <typename scalar_type>
matrix<scalar_type> product(const matrix<scalar_type>& a,
                            const matrix<scalar_type>& b) {
    assert(a.columns() == b.rows());
    size_t arows = a.rows();
    size_t bcolumns = b.columns();
    size_t n = a.columns();
    matrix<scalar_type> c(arows, bcolumns);
    for (size_t i = 0; i < arows; ++i) {
        for (size_t j = 0; j < n; ++j) {
            for (size_t k = 0; k < bcolumns; ++k)
                c(i, k) += a(i, j) * b(j, k);
        }
    }
    return c;
}

template <typename scalar_type>
void swap_rows(matrix<scalar_type>& m, size_t i, size_t j) {
    size_t columns = m.columns();
    for (size_t column = 0; column < columns; ++column)
        std::swap(m(i, column), m(j, column));
}

// Convert matrix to reduced row echelon form
template <typename scalar_type>
void rref(matrix<scalar_type>& m) {
    size_t rows = m.rows();
    size_t columns = m.columns();
    for (size_t row = 0, lead = 0; row < rows && lead < columns; ++row, ++lead) {
        size_t i = row;
        while (m(i, lead) == 0) {
            if (++i == rows) {
                i = row;
                if (++lead == columns)
                    return;
            }
        }
        swap_rows(m, i, row);
        if (m(row, lead) != 0) {
            scalar_type f = m(row, lead);
            for (size_t column = 0; column < columns; ++column)
                m(row, column) /= f;
        }
        for (size_t j = 0; j < rows; ++j) {
            if (j == row)
                continue;
            scalar_type f = m(j, lead);
            for (size_t column = 0; column < columns; ++column)
                m(j, column) -= f * m(row, column);
        }
    }
}

template <typename scalar_type>
matrix<scalar_type> inverse(const matrix<scalar_type>& m) {
    assert(m.rows() == m.columns());
    size_t rows = m.rows();
    matrix<scalar_type> tmp(rows, 2 * rows, 0);
    for (size_t row = 0; row < rows; ++row) {
        for (size_t column = 0; column < rows; ++column)
            tmp(row, column) = m(row, column);
        tmp(row, row + rows) = 1;
    }
    rref(tmp);
    matrix<scalar_type> inv(rows, rows);
    for (size_t row = 0; row < rows; ++row) {
        for (size_t column = 0; column < rows; ++column)
            inv(row, column) = tmp(row, column + rows);
    }
    return inv;
}

template <typename scalar_type>
void print(std::ostream& out, const matrix<scalar_type>& m) {
    size_t rows = m.rows(), columns = m.columns();
    out << std::fixed << std::setprecision(4);
    for (size_t row = 0; row < rows; ++row) {
        for (size_t column = 0; column < columns; ++column) {
            if (column > 0)
                out << ' ';
            out << std::setw(7) << m(row, column);
        }
        out << '\n';
    }
}

int main() {
    matrix<double> m(3, 3, {2, -1, 0, -1, 2, -1, 0, -1, 2});
    std::cout << "Matrix:\n";
    print(std::cout, m);
    auto inv(inverse(m));
    std::cout << "Inverse:\n";
    print(std::cout, inv);
    std::cout << "Product of matrix and inverse:\n";
    print(std::cout, product(m, inv));
    std::cout << "Inverse of inverse:\n";
    print(std::cout, inverse(inv));
    return 0;
}</syntaxhighlight>

{{out}}
<pre>
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
</pre>

=={{header|EasyLang}}==
{{trans|Go}}
<syntaxhighlight>
proc rref . m[][] .
   nrow = len m[][]
   ncol = len m[1][]
   lead = 1
   for r to nrow
      if lead > ncol
         return
      .
      i = r
      while m[i][lead] = 0
         i += 1
         if i > nrow
            i = r
            lead += 1
            if lead > ncol
               return
            .
         .
      .
      swap m[i][] m[r][]
      m = m[r][lead]
      for k to ncol
         m[r][k] /= m
      .
      for i to nrow
         if i <> r
            m = m[i][lead]
            for k to ncol
               m[i][k] -= m * m[r][k]
            .
         .
      .
      lead += 1
   .
.
proc inverse . mat[][] inv[][] .
   inv[][] = [ ]
   ln = len mat[][]
   for i to ln
      if len mat[i][] <> ln
         # not a square matrix
         return
      .
      aug[][] &= [ ]
      len aug[i][] 2 * ln
      for j to ln
         aug[i][j] = mat[i][j]
      .
      aug[i][ln + i] = 1
   .
   rref aug[][]
   for i to ln
      inv[][] &= [ ]
      for j = ln + 1 to 2 * ln
         inv[i][] &= aug[i][j]
      .
   .
.
test[][] = [ [ 1 2 3 ] [ 4 1 6 ] [ 7 8 9 ] ]
inverse test[][] inv[][]
print inv[][]
</syntaxhighlight>

=={{header|Factor}}==
Normally you would call <code>recip</code> to calculate the inverse of a matrix, but it uses a different method than Gauss-Jordan, so here's Gauss-Jordan.
{{works with|Factor|0.99 2020-01-23}}
<syntaxhighlight lang="factor">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.</syntaxhighlight>
{{out}}
<pre>
-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
</pre>

=={{header|Fortran}}==
Note that the Crout algorithm is a more efficient way to invert a matrix.
<syntaxhighlight lang="fortran">!-----------------------------------------------------------------------
! 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</syntaxhighlight>
Test program:
{{works with|gfortran|8.3.0}}
<syntaxhighlight lang="fortran">!             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</syntaxhighlight>
{{out}}
<pre>
 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
</pre>

=={{header|Generic}}==
<syntaxhighlight lang="cpp">
// The Generic Language is a database compiler. This code is compiled into database then executed out of database.

generic coordinaat
{
 ecs;
 uuii;

 coordinaat() { ecs=+a;uuii=+a;}

 coordinaat(ecs_set,uuii_set) { ecs = ecs_set; uuii=uuii_set;}

 operaator<(c)
 {
   iph ecs < c.ecs return troo;
   iph c.ecs < ecs return phals;
   iph uuii < c.uuii return troo;
   return phals;
 }

 operaator==(connpair) // eecuuols and not eecuuols deriiu phronn operaator<
 {
   iph this < connpair return phals;
   iph connpair < this return phals;
   return troo;
 }

 operaator!=(connpair)
 {
   iph this < connpair return troo;
   iph connpair < this return troo;
   return phals;
 }

 too_string()
 {
     return "(" + ecs.too_string() + "," + uuii.too_string() + ")";
 }

 print()
 {
   str = too_string();
   str.print();
 }

 println()
 {
   str = too_string();
   str.println();
 }
}

generic nnaatrics
{
  s;         // this is a set of coordinaat/ualioo pairs.
  iteraator; // this field holds an iteraator phor the nnaatrics.

  nnaatrics()   // no parameters required phor nnaatrics construction.
  {
   s = nioo set();   // create a nioo set of coordinaat/ualioo pairs.
   iteraator = nul; // the iteraator is initially set to nul.
  }

  nnaatrics(copee)   // copee the nnaatrics.
  {
   s = nioo set();   // create a nioo set of coordinaat/ualioo pairs.
   iteraator = nul; // the iteraator is initially set to nul.

   r = copee.rouus;
   c = copee.cols;
   i = 0;
   uuiil i < r
   {
     j = 0;
     uuiil j < c
     {
        this[i,j] = copee[i,j];
       j++;
     }
     i++;
   }
  }

  begin { get { return s.begin; } } // property: used to commence manual iteraashon.

  end { get { return s.end; } } // property: used to dephiin the end itenn of iteraashon

  operaator<(a) // les than operaator is corld bii the avl tree algorithnns
  {             // this operaator innpliis phor instance that you could potenshalee hav sets ou nnaatricss.
    iph cees < a.cees  // connpair the cee sets phurst.
       return troo;
    els iph a.cees < cees
       return phals;
    els                // the cee sets are eecuuol thairphor connpair nnaatrics elennents.
    {
     phurst1 = begin;
     lahst1 = end;
     phurst2 = a.begin;
     lahst2 = a.end;

     uuiil phurst1 != lahst1 && phurst2 != lahst2
     {
       iph phurst1.daata.ualioo < phurst2.daata.ualioo
        return troo;
       els
       {
         iph phurst2.daata.ualioo < phurst1.daata.ualioo
          return phals;
         els
         {
            phurst1 = phurst1.necst;
            phurst2 = phurst2.necst;
         }
       }
     }

     return phals;
    }
  }

    operaator==(connpair) // eecuuols and not eecuuols deriiu phronn operaator<
    {
     iph this < connpair return phals;
     iph connpair < this return phals;
     return troo;
    }

    operaator!=(connpair)
    {
     iph this < connpair return troo;
     iph connpair < this return troo;
     return phals;
    }


   operaator[](cee_a,cee_b) // this is the nnaatrics indexer.
   {
       set
       {
           trii { s >> nioo cee_ualioo(new coordinaat(cee_a,cee_b)); } catch {}
           s << nioo cee_ualioo(new coordinaat(nioo integer(cee_a),nioo integer(cee_b)),ualioo);
       }
       get
       {
          d = s.get(nioo cee_ualioo(new coordinaat(cee_a,cee_b)));
          return d.ualioo;
       }
   }

  operaator>>(coordinaat) // this operaator reennoous an elennent phronn the nnaatrics.
  {
   s >> nioo cee_ualioo(coordinaat);
   return this;
  }

  iteraat() // and this is how to iterate on the nnaatrics.
  {
      iph iteraator.nul()
      {
         iteraator = s.lepht_nnohst;
         iph iteraator == s.heder
             return nioo iteraator(phals,nioo nun());
         els
             return nioo iteraator(troo,iteraator.daata.ualioo);
      }
      els
      {
         iteraator = iteraator.necst;
  
         iph iteraator == s.heder
         {
             iteraator = nul;
             return nioo iteraator(phals,nioo nun());
         }
         els
             return nioo iteraator(troo,iteraator.daata.ualioo);
      }
   }

   couunt // this property returns a couunt ou elennents in the nnaatrics.
   {
      get
      {
         return s.couunt;
      }
   }

   ennptee // is the nnaatrics ennptee?
   {
       get
       {
           return s.ennptee;
       }
   }


   lahst // returns the ualioo of the lahst elennent in the nnaatrics.
   {
       get
       {
           iph ennptee
                 throuu "ennptee nnaatrics";
           els
               return s.lahst.ualioo;
       }
   }

    too_string() // conuerts the nnaatrics too aa string
    {
       return s.too_string();
    }

    print() // prints the nnaatrics to the consohl.
    {
        out = too_string();
        out.print();
    }

   println() // prints the nnaatrics as a liin too the consohl.
    {
        out = too_string();
        out.println();
    } 

   cees // return the set ou cees ou the nnaatrics (a set of coordinaats).
   {
      get
      {
          k = nioo set();
          phor e : s k << e.cee;     
          return k;
      }
   }

  operaator+(p)
  {
     ouut = nioo nnaatrics();
     phurst1 = begin;
     lahst1 = end;
     phurst2 = p.begin;
     lahst2 = p.end;
     uuiil phurst1 != lahst1 && phurst2 != lahst2
     {
        ouut[phurst1.daata.cee.ecs,phurst1.daata.cee.uuii] = phurst1.daata.ualioo + phurst2.daata.ualioo;
        phurst1 = phurst1.necst;
        phurst2 = phurst2.necst;
     }
     return ouut;
  }
  
  operaator-(p)
  {
     ouut = nioo nnaatrics();
     phurst1 = begin;
     lahst1 = end;
     phurst2 = p.begin;
     lahst2 = p.end;
     uuiil phurst1 != lahst1 && phurst2 != lahst2
     {
        ouut[phurst1.daata.cee.ecs,phurst1.daata.cee.uuii] = phurst1.daata.ualioo - phurst2.daata.ualioo;
        phurst1 = phurst1.necst;
        phurst2 = phurst2.necst;
     }
     return ouut;
  }

  rouus
  {
     get
     {
        r = +a;
        phurst1 = begin;
        lahst1 = end;
        uuiil phurst1 != lahst1
        {
            iph r < phurst1.daata.cee.ecs r = phurst1.daata.cee.ecs;
            phurst1 = phurst1.necst;
        }
        return r + +b;
     }
  }

  cols
  {
     get
     {
        c = +a;
        phurst1 = begin;
        lahst1 = end;
        uuiil phurst1 != lahst1
        {
            iph c < phurst1.daata.cee.uuii c = phurst1.daata.cee.uuii;
            phurst1 = phurst1.necst;
        }
        return c + +b;
     }
  }

  operaator*(o)
  {
      iph cols != o.rouus throw "rouus-cols nnisnnatch";
      reesult = nioo nnaatrics();
      rouu_couunt = rouus;
      colunn_couunt = o.cols;
      loop = cols;
      i = +a;
      uuiil i < rouu_couunt
      {
        g = +a;
        uuiil g < colunn_couunt
        {
           sunn = +a.a;
           h = +a;
           uuiil h < loop
           {
             a = this[i, h];

             b = o[h, g];
             nn = a * b;
             sunn = sunn +  nn;
             h++;
           }

            reesult[i, g] = sunn;

            g++;
         }
        i++;
     }
     return reesult;
 }

  suuop_rouus(a, b)
  {
    c = cols;
    i = 0;
    uuiil u < cols
    {
      suuop = this[a, i];
      this[a, i] = this[b, i];
      this[b, i] = suuop;
      i++;
    }
  }

  suuop_colunns(a, b)
  {
    r = rouus;
    i = 0;
    uuiil i < rouus
    {
     suuopp = this[i, a];
     this[i, a] = this[i, b];
     this[i, b] = suuop;
     i++;
    }
  }

  transpohs
  {
     get
     {
         reesult = new nnaatrics();

         r = rouus;
         c = cols;
         i=0;
         uuiil i < r
         {
              g = 0;
              uuiil g < c
              {
                 reesult[g, i] = this[i, g];
                 g++;
              }
             i++;
         }

         return reesult;
      }
   }

   deternninant
   {
        get
        {
            rouu_couunt = rouus;
            colunn_count = cols;

            if rouu_couunt != colunn_count
                throw "not a scuuair nnaatrics";

            if rouu_couunt == 0
               throw "the nnaatrics is ennptee";

            if rouu_couunt == 1
               return this[0, 0];

            if rouu_couunt == 2
                return this[0, 0] * this[1, 1] -
                       this[0, 1] * this[1, 0];

            temp = nioo nnaatrics();

            det = 0.0;
            parity = 1.0;

            j = 0;
            uuiil j < rouu_couunt
            {
                 k = 0;
                 uuiil k < rouu_couunt-1
                 {
                      skip_col = phals;

                      l = 0;
                      uuiil l < rouu_couunt-1
                      {
                           if l == j skip_col = troo;

                           if skip_col
                               n = l + 1;
                           els
                               n = l;

                            temp[k, l] = this[k + 1, n];
                            l++;
                       }
                       k++;
                  }

                  det = det + parity * this[0, j] * temp.deeternninant;

                  parity = 0.0 - parity;
                  j++;
            }

        return det;
      }
  }

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

   ad_colunn(a, b)
   {
       c = rouus;
       i = 0;
       uuiil i < c
       {
         this[i, a] = this[i, a] + this[i, b];
         i++;
        }
   }

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

   subtract_colunn(a, b)
   {
       c = rouus;
       i = 0;
       uuiil i < c
       {
         this[i, a] = this[i, a] - this[i, b];
         i++;
        }
   }

   nnultiplii_rouu(rouu, scalar)
   {
       c = cols;
       i = 0;
       uuiil i < c
       {  
          this[rouu, i] = this[rouu, i] * scalar;
           i++;
       }
   }

   nnultiplii_colunn(colunn, scalar)
   {
       r = rouus;
       i = 0;
       uuiil i < r
       {  
          this[i, colunn] = this[i, colunn] * scalar;
           i++;
       }
   }

   diuiid_rouu(rouu, scalar)
   {
       c = cols;
       i = 0;
       uuiil i < c
       {  
          this[rouu, i] = this[rouu, i] / scalar;
           i++;
       }
   }

   diuiid_colunn(colunn, scalar)
   {
       r = rouus;
       i = 0;
       uuiil i < r
       {  
          this[i, colunn] = this[i, colunn] / scalar;
           i++;
       }
   }

   connbiin_rouus_ad(a,b,phactor)
   {
       c = cols;
       i = 0;
       uuiil i < c
       {
          this[a, i] = this[a, i] + phactor * this[b, i];
          i++;
       }
    }

   connbiin_rouus_subtract(a,b,phactor)
   {
       c = cols;
       i = 0;
       uuiil i < c
       {
          this[a, i] = this[a, i] - phactor * this[b, i];
          i++;
       }
    }

   connbiin_colunns_ad(a,b,phactor)
   {
       r = rouus;
       i = 0;
       uuiil i < r
       {
          this[i, a] = this[i, a] + phactor * this[i, b];
          i++;
       }
    }

   connbiin_colunns_subtract(a,b,phactor)
   {
       r = rouus;
       i = 0;
       uuiil i < r
       {
          this[i, a] = this[i, a] - phactor * this[i, b];
          i++;
       }
    }

    inuers
    {
        get
        {
            rouu_couunt = rouus;
            colunn_couunt = cols;

            iph rouu_couunt != colunn_couunt
                throw "nnatrics not scuuair";

            els iph rouu_couunt == 0
                throw "ennptee nnatrics";

            els iph rouu_couunt == 1
            {
                r = nioo nnaatrics();
                r[0, 0] = 1.0 / this[0, 0];
                return r;
            }

            gauss = nioo nnaatrics(this);

            i = 0;
            uuiil i < rouu_couunt
            {
                 j = 0;
                 uuiil j < rouu_couunt
                 {
                       iph i == j 

                            gauss[i, j + rouu_couunt] = 1.0;
                        els
                            gauss[i, j + rouu_couunt] = 0.0;
                      j++;
                 }

                 i++;
             }

              j = 0;
              uuiil j < rouu_couunt
              {
                 iph gauss[j, j] == 0.0
                 {
                      k = j + 1;

                      uuiil k < rouu_couunt
                      {
                          if gauss[k, j] != 0.0 {gauss.nnaat.suuop_rouus(j, k); break; }
                           k++;
                       }

                       if k == rouu_couunt throw "nnatrics is singioolar";
                 }

                 phactor = gauss[j, j];
                 iph phactor != 1.0 gauss.diuiid_rouu(j, phactor);

                 i = j+1;
                 uuiil i < rouu_couunt
                 {
                     gauss.connbiin_rouus_subtract(i, j, gauss[i, j]);
                     i++;
                  }

                 j++;
              }

              i = rouu_couunt - 1;
              uuiil i > 0
              {
                  k = i - 1;
                  uuiil k >= 0
                  {
                      gauss.connbiin_rouus_subtract(k, i, gauss[k, i]);
                      k--;
                  }
                  i--;
              }

               reesult = nioo nnaatrics();

                i = 0;
                uuiil i < rouu_couunt
                {
                     j = 0;
                     uuiil  j < rouu_couunt
                     {
                         reesult[i, j] = gauss[i, j + rouu_couunt];
                         j++;
                     }
                    i++;
                }

              return reesult;
            }
        }

}
</syntaxhighlight>
=={{header|Go}}==
{{trans|Kotlin}}
<syntaxhighlight lang="go">package main

import "fmt"

type vector = []float64
type matrix []vector

func (m matrix) inverse() matrix {
    le := len(m)
    for _, v := range m {
        if len(v) != le {
            panic("Not a square matrix")
        }
    }
    aug := make(matrix, le)
    for i := 0; i < le; i++ {
        aug[i] = make(vector, 2*le)
        copy(aug[i], m[i])
        // augment by identity matrix to right
        aug[i][i+le] = 1
    }
    aug.toReducedRowEchelonForm()
    inv := make(matrix, le)
    // remove identity matrix to left
    for i := 0; i < le; i++ {
        inv[i] = make(vector, le)
        copy(inv[i], aug[i][le:])
    }
    return inv
}

// note: this mutates the matrix in place
func (m matrix) toReducedRowEchelonForm() {
    lead := 0
    rowCount, colCount := len(m), len(m[0])
    for r := 0; r < rowCount; r++ {
        if colCount <= lead {
            return
        }
        i := r

        for m[i][lead] == 0 {
            i++
            if rowCount == i {
                i = r
                lead++
                if colCount == lead {
                    return
                }
            }
        }

        m[i], m[r] = m[r], m[i]
        if div := m[r][lead]; div != 0 {
            for j := 0; j < colCount; j++ {
                m[r][j] /= div
            }
        }

        for k := 0; k < rowCount; k++ {
            if k != r {
                mult := m[k][lead]
                for j := 0; j < colCount; j++ {
                    m[k][j] -= m[r][j] * mult
                }
            }
        }
        lead++
    }
}

func (m matrix) print(title string) {
    fmt.Println(title)
    for _, v := range m {
        fmt.Printf("% f\n", v)
    }
    fmt.Println()
}

func main() {
    a := matrix{{1, 2, 3}, {4, 1, 6}, {7, 8, 9}}
    a.inverse().print("Inverse of A is:\n")

    b := matrix{{2, -1, 0}, {-1, 2, -1}, {0, -1, 2}}
    b.inverse().print("Inverse of B is:\n")
}</syntaxhighlight>
{{out}}
<pre>
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]
</pre>
===Alternative===
{{trans|PowerShell}}
<syntaxhighlight lang="go">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")
}</syntaxhighlight>
{{out}}
<pre>Same output than above</pre>

=={{header|Haskell}}==
<syntaxhighlight lang="haskell">isMatrix xs = null xs || all ((== (length.head $ xs)).length) xs

isSquareMatrix xs = null xs || all ((== (length xs)).length) xs

mult:: Num a => [[a]] -> [[a]] -> [[a]]
mult uss vss = map ((\xs -> if null xs then [] else foldl1 (zipWith (+)) xs). zipWith (\vs u -> map (u*) vs) vss) uss

matI::(Num a) => Int -> [[a]]
matI n = [ [fromIntegral.fromEnum $ i == j | j <- [1..n]] | i <- [1..n]]

inversion xs = gauss xs (matI $ length xs)

gauss::[[Double]] -> [[Double]] -> [[Double]]
gauss xs bs = map (map fromRational) $ solveGauss (toR xs) (toR bs)
    where toR = map $ map toRational

solveGauss:: (Fractional a, Ord a) => [[a]] -> [[a]] -> [[a]]
solveGauss xs bs | null xs || null bs || length xs /= length bs || (not $ isSquareMatrix xs) || (not $ isMatrix bs) = []
                 | otherwise = uncurry solveTriangle $ triangle xs bs

solveTriangle::(Fractional a,Eq a) => [[a]] -> [[a]] -> [[a]]
solveTriangle us _ | not.null.dropWhile ((/= 0).head) $ us = []
solveTriangle ([c]:as) (b:bs) = go as bs [map (/c) b]
  where
  val us vs ws = let u = head us in map (/u) $ zipWith (-) vs (head $ mult [tail us] ws)
  go [] _ zs          = zs
  go _ [] zs          = zs
  go (x:xs) (y:ys) zs = go xs ys $ (val x y zs):zs

triangle::(Num a, Ord a) => [[a]] -> [[a]] -> ([[a]],[[a]])
triangle xs bs = triang ([],[]) (xs,bs)
    where
    triang ts (_,[]) = ts
    triang ts ([],_) = ts
    triang (os,ps) zs = triang (us:os,cs:ps).unzip $ [(fun tus vs, fun cs es) | (v:vs,es) <- zip uss css,let fun = zipWith (\x y -> v*x - u*y)]
        where ((us@(u:tus)):uss,cs:css) = bubble zs

bubble::(Num a, Ord a) => ([[a]],[[a]]) -> ([[a]],[[a]])
bubble (xs,bs) = (go xs, go bs)
    where
    idmax = snd.maximum.flip zip [0..].map (abs.head) $ xs
    go ys = let (us,vs) = splitAt idmax ys in vs ++ us
 
main = do
  let a = [[1, 2, 3], [4, 1, 6], [7, 8, 9]]
  let b = [[2, -1, 0], [-1, 2, -1], [0, -1, 2]]
  putStrLn "inversion a ="
  mapM_ print $ inversion a
  putStrLn "\ninversion b ="
  mapM_ print $ inversion b</syntaxhighlight>
{{out}}
<pre>
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]
</pre>


=={{header|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.
<syntaxhighlight lang="j">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</syntaxhighlight>

'''Usage:'''
<syntaxhighlight lang="j">   ]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</syntaxhighlight>

=={{header|Java}}==
<syntaxhighlight lang="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();
    }
}</syntaxhighlight>

<syntaxhighlight lang="java">// Matrix.java

public class Matrix {
    private int rows;
    private int columns;
    private double[][] elements;

    public Matrix(int rows, int columns) {
        this.rows = rows;
        this.columns = columns;
        elements = new double[rows][columns];
    }

    public void set(int row, int column, double value) {
        elements[row][column] = value;
    }

    public double get(int row, int column) {
        return elements[row][column];
    }

    public void print() {
        for (int row = 0; row < rows; ++row) {
            for (int column = 0; column < columns; ++column) {
                if (column > 0)
                    System.out.print(' ');
                System.out.printf("%7.3f", elements[row][column]);
            }
            System.out.println();
        }
    }

    // Returns the inverse of this matrix
    public Matrix inverse() {
        assert(rows == columns);
        // Augment by identity matrix
        Matrix tmp = new Matrix(rows, columns * 2);
        for (int row = 0; row < rows; ++row) {
            System.arraycopy(elements[row], 0, tmp.elements[row], 0, columns);
            tmp.elements[row][row + columns] = 1;
        }
        tmp.toReducedRowEchelonForm();
        Matrix inv = new Matrix(rows, columns);
        for (int row = 0; row < rows; ++row)
            System.arraycopy(tmp.elements[row], columns, inv.elements[row], 0, columns);
        return inv;
    }

    // Converts this matrix into reduced row echelon form
    public void toReducedRowEchelonForm() {
        for (int row = 0, lead = 0; row < rows && lead < columns; ++row, ++lead) {
            int i = row;
            while (elements[i][lead] == 0) {
                if (++i == rows) {
                    i = row;
                    if (++lead == columns)
                        return;
                }
            }
            swapRows(i, row);
            if (elements[row][lead] != 0) {
                double f = elements[row][lead];
                for (int column = 0; column < columns; ++column)
                    elements[row][column] /= f;
            }
            for (int j = 0; j < rows; ++j) {
                if (j == row)
                    continue;
                double f = elements[j][lead];
                for (int column = 0; column < columns; ++column)
                    elements[j][column] -= f * elements[row][column];
            }
        }
    }

    // Returns the matrix product of a and b
    public static Matrix product(Matrix a, Matrix b) {
        assert(a.columns == b.rows);
        Matrix result = new Matrix(a.rows, b.columns);
        for (int i = 0; i < a.rows; ++i) {
            double[] resultRow = result.elements[i];
            double[] aRow = a.elements[i];
            for (int j = 0; j < a.columns; ++j) {
                double[] bRow = b.elements[j];
                for (int k = 0; k < b.columns; ++k)
                    resultRow[k] += aRow[j] * bRow[k];
            }
        }
        return result;
    }

    private void swapRows(int i, int j) {
        double[] tmp = elements[i];
        elements[i] = elements[j];
        elements[j] = tmp;
    }
}</syntaxhighlight>

{{out}}
<pre>
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
</pre>

=={{header|jq}}==
'''Adapted from [https://rosettacode.org/wiki/Category:Wren-matrix Wren-matrix]'''
{{works with|jq}}
'''Works with gojq, the Go implementation of jq'''

<syntaxhighlight lang="jq"># Create an m x n matrix,
# it being understood that:
# matrix(0; _; _) evaluates to []
def matrix(m; n; init):
  if m == 0 then []
  elif m == 1 then [[range(0;n) | init]]
  elif m > 0 then
    [range(0;n) | init] as $row
    | [range(0;m) | $row ]
  else error("matrix\(m);_;_) invalid")
  end;

def mprint($dec):
   def max(s): reduce s as $x (null; if . == null or $x > . then $x else . end);
   def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;
    
   pow(10; $dec) as $power
   | def p: (. * $power | round) / $power;
   (max(.[][] | p | tostring | length) + 1) as $w
   | . as $in
   | range(0; length) as $i
   | reduce range(0; .[$i]|length) as $j ("|"; . + ($in[$i][$j]|p|lpad($w)))
   | . + " |" ;

# Swaps two rows of the input matrix using IO==0
def swapRows(rowNum1; rowNum2):
  if (rowNum1 == rowNum2)
  then .
  else .[rowNum1] as $t
  | .[rowNum1] = .[rowNum2]
  | .[rowNum2] = $t
  end;

def toReducedRowEchelonForm:
   . as $in
   | length as $nr
   | (.[0]|length) as $nc
   | {a: $in, lead: 0 }
   | label $out
   | last(foreach range(0; $nr) as $r (.;
        if $nc <= .lead then ., break $out else . end
        | .i = $r
        | until( .a[.i][.lead] != 0 or .lead == $nc - 1;
            .i += 1
            | if $nr == .i
              then .i = $r
              | .lead += 1
	      else .
	      end
           )
        | swapRows(.i; $r)
        | if .a[$r][.lead] != 0
          then .a[$r][.lead] as $div
          | reduce range(0; $nc) as $j (.; .a[$r][$j] /= $div)
	  else .
	  end
        | reduce range(0; $nr) as $k (.;
            if $k != $r
            then .a[$k][.lead] as $mult
            | reduce range(0; $nc) as $j (.; .a[$k][$j] -= .a[$r][$j] * $mult)
            else .
	    end)
        | .lead += 1
	))
  | .a ;


# Assumption: the input is a square matrix with an inverse
# Uses the Gauss-Jordan method.
def inverse:
  . as $a
  | length as $nr
  | reduce range(0; $nr) as $i (
      matrix($nr; 2 * $nr; 0);
      reduce range( 0; $nr) as $j (.;
         .[$i][$j] = $a[$i][$j]
         | .[$i][$i + $nr] = 1 ))
  | last(toReducedRowEchelonForm)
  | . as $ary
  | reduce range(0; $nr) as $i ( matrix($nr; $nr; 0);
      reduce range($nr; 2 *$nr) as $j (.;
        .[$i][$j - $nr] = $ary[$i][$j] )) ;</syntaxhighlight>
'''The Task'''
<syntaxhighlight lang="jq">def arrays: [
    [ [ 1,  2,  3],
      [ 4,  1,  6],
      [ 7,  8,  9] ],
 
    [ [ 2, -1,  0 ],
      [-1,  2, -1 ],
      [ 0, -1,  2 ] ],
 
    [ [ -1, -2, 3, 2 ],
      [ -4, -1, 6, 2 ],
      [  7, -8, 9, 1 ],
      [  1, -2, 1, 3 ] ]
];

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

task</syntaxhighlight>
{{out}}
<pre>
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 |
</pre>

=={{header|Julia}}==
{{works with|Julia|0.6}}

'''Built-in''' LAPACK-based linear solver uses partial-pivoted Gauss elimination):
<syntaxhighlight lang="julia">A = [1 2 3; 4 1 6; 7 8 9]
@show I / A
@show inv(A)</syntaxhighlight>

'''Native implementation''':
<syntaxhighlight lang="julia">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)</syntaxhighlight>

{{out}}
<pre>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]</pre>

=={{header|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.
<syntaxhighlight lang="scala">// version 1.2.21

typealias Matrix = Array<DoubleArray>

fun Matrix.inverse(): Matrix {
    val len = this.size
    require(this.all { it.size == len }) { "Not a square matrix" }
    val aug = Array(len) { DoubleArray(2 * len) }
    for (i in 0 until len) {
        for (j in 0 until len) aug[i][j] = this[i][j]
        // augment by identity matrix to right
        aug[i][i + len] = 1.0
    }
    aug.toReducedRowEchelonForm()
    val inv = Array(len) { DoubleArray(len) }
    // remove identity matrix to left
    for (i in 0 until len) {
        for (j in len until 2 * len) inv[i][j - len] = aug[i][j]
    }
    return inv
}

fun Matrix.toReducedRowEchelonForm() {
    var lead = 0
    val rowCount = this.size
    val colCount = this[0].size
    for (r in 0 until rowCount) {
        if (colCount <= lead) return
        var i = r

        while (this[i][lead] == 0.0) {
            i++
            if (rowCount == i) {
                i = r
                lead++
                if (colCount == lead) return
            }
        }

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

        if (this[r][lead] != 0.0) {
           val div = this[r][lead]
           for (j in 0 until colCount) this[r][j] /= div
        }

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

        lead++
    }
}

fun Matrix.printf(title: String) {
    println(title)
    val rowCount = this.size
    val colCount = this[0].size

    for (r in 0 until rowCount) {
        for (c in 0 until colCount) {
            if (this[r][c] == -0.0) this[r][c] = 0.0  // get rid of negative zeros
            print("${"% 10.6f".format(this[r][c])}  ")
        }
        println()
    }

    println()
}

fun main(args: Array<String>) {
    val a = arrayOf(
        doubleArrayOf(1.0, 2.0, 3.0),
        doubleArrayOf(4.0, 1.0, 6.0),
        doubleArrayOf(7.0, 8.0, 9.0)
    )
    a.inverse().printf("Inverse of A is :\n")

    val b = arrayOf(
        doubleArrayOf( 2.0, -1.0,  0.0),
        doubleArrayOf(-1.0,  2.0, -1.0),
        doubleArrayOf( 0.0, -1.0,  2.0)
    )
    b.inverse().printf("Inverse of B is :\n")    
}</syntaxhighlight>

{{out}}
<pre>
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  
</pre>

=={{header|Lambdatalk}}==
<syntaxhighlight lang="scheme">
{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]]
</syntaxhighlight>
=={{header|M2000 Interpreter}}==
We can use supported types like Decimal, Currency, Double, Single. Also matrix may have different base on each dimension.

<syntaxhighlight lang="m2000 interpreter">
Module Gauss_Jordan_matrix_inversion (there$="") {
	open there$ for output as #f
	dim matrix(1 to 3, 2 to 4) as Double
	matrix(1,2)=2, -1, 0, -1, 2, -1, 0, -1, 2
	print #f, "Matrix 3x3:"
	disp(matrix())
	print #f, "Inverse:"
	disp(@inverse(matrix()))
	dim matrix2(2 to 4, 3 to 5) as single
	matrix2(2, 3)=1, 2, 3, 4, 1, 6, 7, 8, 9
	print #f,"Matrix 3x3:"
	disp(matrix2())
	print #f,"Inverse:"	
	disp(@inverse(matrix2()))
	dim matrix3(1 to 4,1 to 4) as Currency, matrix4()
	matrix3(1,1)=-1,-2,3,2,-4,-1,6,2,7,-8,9,1,1,-2,1,3
	print #f,"Matrix 4x4:"
	disp(matrix3())
	matrix4()=@inverse(matrix3())
	Rem Print Type$(matrix4(1,1))="Currency"
	print #f,"Inverse:"	
	disp(matrix4())
	close #f
	function inverse(a())
		if dimension(a())<>2 then Error "Not 2 Dimensions Matrix"
		if not dimension(a(),1)=dimension(a(),2) then Error "Not Square Matrix"
		Local aug()  // this is an empty array (type of variant)
		aug()=a()  // so we get the same type from a()
		dim aug(dimension(a(),1), dimension(a(),1)*2)=0
		Rem gosub disp(aug())
		Local long bi=dimension(a(), 1, 0), bj=dimension(a(), 2, 0)
		Local long bim=dimension(a(), 1)-1, bjm=dimension(a(), 2)-1
		Local long i, j
		for i=0 to bim
			for j=0 to bjm
				aug(i,j)=a(i+bi,j+bj)
			next
			aug(i, i+bjm+1)=1
		next
		aug()=@ToReducedRowEchelonForm(aug())
		for i=0 to bim
			for j=0 to bjm
				a(i+bi,j+bj)=aug(i,j+bjm+1)
			next
		next
		Rem : Print type$(aug(0,0))
		=a()
	end function
	function ToReducedRowEchelonForm(a())
		Local long bi=dimension(a(), 1, 0), bj=dimension(a(), 2, 0)
		Local long rowcount=dimension(a(), 1), columncount=dimension(a(), 2)
		Local long lead=bj, r, i, j
	
		for r=bi to rowcount-1+bi {
			if columncount<=lead then exit
			i=r
			while A(i,lead)=0 {
				i++
				if rowcount=i then i=r : lead++ : if columncount<lead then exit
			}
		
			for c =bj to columncount-1+bj {swap A(i, c), A(r, c)}
		
			if A(r, lead)<>0 then
				div1=A(r,lead)
				For c =bj to columncount-1+bj {A(r, c)/=div1} 
			end if
			for i=bi to rowcount-1+bi {
				if i<>r then {
					mult=A(i,lead)
					for j=bj to columncount-1+bj {A(i,j)-=A(r,j)*mult}
				}
			} 
			lead++
		}
		=a()
	End function
	sub disp(a())
		local long i, j, c=8 // column width plus 1 space
		for i=dimension(a(), 1, 0) to dimension(a(), 1, 1)
			for j=dimension(a(), 2, 0) to dimension(a(), 2, 1)
				Print #f, format$("{0:-"+c+"} ",str$(round(A(i, j), 5),1033));
			Next
			// if pos>0 then print  // not used here
			Print #f, ""
		Next
		Print #f, ""
	End sub
}
Gauss_Jordan_matrix_inversion  // to screen
Gauss_Jordan_matrix_inversion "out.txt"
win "notepad", dir$+"out.txt"
</syntaxhighlight>

{{out}}
<pre>
Matrix 3x3:
       2       -1        0 
      -1        2       -1 
       0       -1        2 

Inverse:
    0.75      0.5     0.25 
     0.5        1      0.5 
    0.25      0.5     0.75 

Matrix 3x3:
       1        2        3 
       4        1        6 
       7        8        9 

Inverse:
 -0.8125    0.125   0.1875 
   0.125    -0.25    0.125 
 0.52083    0.125 -0.14583 

Matrix 4x4:
      -1       -2        3        2 
      -4       -1        6        2 
       7       -8        9        1 
       1       -2        1        3 

Inverse:
 -0.9129   0.2463   0.0941    0.413 
 -1.6519   0.6522   0.0434   0.6521 
 -0.6954   0.3623   0.0797   0.1956 
 -0.5651   0.2319   -0.029   0.5652 
</pre>

=={{header|Mathematica}} / {{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">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] ;;]]</syntaxhighlight>
{{out}}
<pre>{{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}}</pre>

=={{header|MATLAB}}==
<syntaxhighlight lang="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)</syntaxhighlight>

{{out}}
<pre>
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
</pre>

=={{header|Nim}}==
We reuse the algorithm of task https://rosettacode.org/wiki/Reduced_row_echelon_form (version using floats).
<syntaxhighlight lang="nim">import strformat, strutils

const Eps = 1e-10

type
  Matrix[M, N: static Positive] = array[M, array[N, float]]
  SquareMatrix[N: static Positive] = Matrix[N, N]


func toSquareMatrix[N: static Positive](a: array[N, array[N, int]]): SquareMatrix[N] =
  ## Convert a square matrix of integers to a square matrix of floats.

  for i in 0..<N:
    for j in 0..<N:
      result[i][j] = a[i][j].toFloat


func transformToRref(mat: var Matrix) =
  ## Transform a matrix to reduced row echelon form.

  var lead = 0

  for r in 0..<mat.M:

    if lead >= mat.N: return

    var i = r
    while mat[i][lead] == 0:
      inc i
      if i == mat.M:
        i = r
        inc lead
        if lead == mat.N: return
    swap mat[i], mat[r]

    let d = mat[r][lead]
    if abs(d) > Eps:    # Checking "d != 0" will give wrong results in some cases.
      for item in mat[r].mitems:
        item /= d

    for i in 0..<mat.M:
      if i != r:
        let m = mat[i][lead]
        for c in 0..<mat.N:
          mat[i][c] -= mat[r][c] * m

    inc lead


func inverse(mat: SquareMatrix): SquareMatrix[mat.N] =
  ## Return the inverse of a matrix.

  # Build augmented matrix.
  var augmat: Matrix[mat.N, 2 * mat.N]
  for i in 0..<mat.N:
    augmat[i][0..<mat.N] = mat[i]
    augmat[i][mat.N + i] = 1

  # Transform it to reduced row echelon form.
  augmat.transformToRref()

  # Check if the first half is the identity matrix and extract second half.
  for i in 0..<mat.N:
    for j in 0..<mat.N:
      if augmat[i][j] != float(i == j):
        raise newException(ValueError, "matrix is singular")
      result[i][j] = augmat[i][mat.N + j]


proc `$`(mat: Matrix): string =
  ## Display a matrix (which may be a square matrix).

  for row in mat:
    var line = ""
    for val in row:
      line.addSep(" ", 0)
      line.add &"{val:9.5f}"
    echo line


#———————————————————————————————————————————————————————————————————————————————————————————————————

template runTest(mat: SquareMatrix) =
  ## Run a test using square matrix "mat".

  echo "Matrix:"
  echo $mat
  echo "Inverse:"
  echo mat.inverse
  echo ""

let m1 = [[1, 2, 3],
          [4, 1, 6],
          [7, 8, 9]].toSquareMatrix()

let m2 = [[ 2, -1,  0],
          [-1,  2, -1],
          [ 0, -1,  2]].toSquareMatrix()

let m3 = [[ -1, -2,  3,  2],
          [ -4, -1,  6,  2],
          [  7, -8,  9,  1],
          [  1, -2,  1,  3]].toSquareMatrix()

runTest(m1)
runTest(m2)
runTest(m3)</syntaxhighlight>

{{out}}
<pre>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</pre>
=={{header|Pascal}}==
==={{header|Free Pascal}}===
<syntaxhighlight lang="pascal">
program MatrixInverse;
{$mode ObjFPC}{$H+}

const
  MATRIX_1: array of array of Real = ((1, 2, 3), 
                                      (4, 1, 6), 
                                      (7, 8, 9));
                                      
  MATRIX_2: array of array of Real = ((3.0, 1.0, 5.0, 9.0, 7.0),
                                      (8.0, 2.0, 8.0, 0.0, 1.0),
                                      (1.0, 7.0, 2.0, 0.0, 3.0),
                                      (0.0, 1.0, 7.0, 0.0, 9.0),
                                      (3.0, 5.0, 6.0, 1.0, 1.0));

type
  Matrix = array of array of Real;

function PopulateMatrix(A: Matrix; Order: Integer): Matrix;
var
  i, j: Integer;
begin
  SetLength(Result, Succ(Order), Succ(Order * 2));
  for i := 0 to Pred(Order) do
    for j := 0 to Pred(Order) do
      Result[i + 1, j + 1] := A[i, j];
end;

procedure PrintMatrix(const A: Matrix);
var
  i, j, Order: Integer;
begin
  Order := Length(A);
  for i := 0 to Pred(Order) do
  begin
    for j := 0 to Pred(Order) do
      Write(A[i, j]:8:4);
    WriteLn;
  end;
  WriteLn;
end;

procedure InterchangeRows(var A: Matrix; Order: Integer; Row1, Row2: Integer);
var
  j: Integer;
  temp: Real;
begin
  for j := 1 to 2 * Order do
  begin
    temp := A[Row1, j];
    A[Row1, j] := A[Row2, j];
    A[Row2, j] := temp;
  end;
end;

procedure DivideRow(var A: Matrix; Order: Integer; Row: Integer; Divisor: Real);
var
  j: Integer;
begin
  for j := 1 to 2 * Order do
    A[Row, j] := A[Row, j] / Divisor;
end;

procedure SubtractRows(var A: Matrix; Order: Integer; Row1, Row2: Integer; Multiplier: Real);
var
  j: Integer;
begin
  for j := 1 to 2 * Order do
    A[Row1, j] := A[Row1, j] - Multiplier * A[Row2, j];
end;

function GaussJordan(B: Matrix): Matrix;
var
  i, j, Order: Integer;
  Pivot, Multiplier: Real;
  A: Matrix;
begin
  Order := Length(B);
  SetLength(Result, Order, Order);
  A := PopulateMatrix(B, Order);

  // Create the augmented matrix
  for i := 1 to Order do
    for j := Order + 1 to 2 * Order do
      if j = (i + Order) then
        A[i, j] := 1;

  // Interchange rows if needed
  for i := Order downto 2 do
    if A[i - 1, 1] < A[i, 1] then
      InterchangeRows(A, Order, i - 1, i);

  // Perform Gauss-Jordan elimination
  for i := 1 to Order do
  begin
    Pivot := A[i, i];
    DivideRow(A, Order, i, Pivot);
    for j := 1 to Order do
      if j <> i then
      begin
        Multiplier := A[j, i];
        SubtractRows(A, Order, j, i, Multiplier);
      end;
  end;

  for i := 0 to Pred(Order) do
    for j := 0 to Pred(Order) do
      Result[i, j] := A[Succ(i), Succ(j) + Order];
end;

begin
  writeln('== Original Matrix ==');
  PrintMatrix(MATRIX_1);
  writeln('== Inverse Matrix ==');
  PrintMatrix(GaussJordan(MATRIX_1));
  
  writeln('== Original Matrix ==');
  PrintMatrix(MATRIX_2);
  writeln('== Inverse Matrix ==');
  PrintMatrix(GaussJordan(MATRIX_2));
end.
</syntaxhighlight>
{{out}}
<pre>
== Original Matrix ==
  1.0000  2.0000  3.0000
  4.0000  1.0000  6.0000
  7.0000  8.0000  9.0000

== Inverse Matrix ==
 -0.8125  0.1250  0.1875
  0.1250 -0.2500  0.1250
  0.5208  0.1250 -0.1458

== Original Matrix ==
  3.0000  1.0000  5.0000  9.0000  7.0000
  8.0000  2.0000  8.0000  0.0000  1.0000
  1.0000  7.0000  2.0000  0.0000  3.0000
  0.0000  1.0000  7.0000  0.0000  9.0000
  3.0000  5.0000  6.0000  1.0000  1.0000

== Inverse Matrix ==
  0.0286  0.1943  0.1330 -0.0596 -0.2578
 -0.0058 -0.0326  0.1211 -0.0380  0.0523
 -0.0302 -0.0682 -0.1790  0.0605  0.2719
  0.1002 -0.0673 -0.0562 -0.0626  0.0981
  0.0241  0.0567  0.1258  0.0682 -0.2173
</pre>

=={{header|Perl}}==
Included code from [[Reduced_row_echelon_form#Perl|Reduced row echelon form]] task.
<syntaxhighlight lang="perl">sub rref {
  our @m; local *m = shift;
  @m or return;
  my ($lead, $rows, $cols) = (0, scalar(@m), scalar(@{$m[0]}));

  foreach my $r (0 .. $rows - 1) {
     $lead < $cols or return;
      my $i = $r;

      until ($m[$i][$lead])
         {++$i == $rows or next;
          $i = $r;
          ++$lead == $cols and return;}

      @m[$i, $r] = @m[$r, $i];
      my $lv = $m[$r][$lead];
      $_ /= $lv foreach @{ $m[$r] };

      my @mr = @{ $m[$r] };
      foreach my $i (0 .. $rows - 1)
         {$i == $r and next;
          ($lv, my $n) = ($m[$i][$lead], -1);
          $_ -= $lv * $mr[++$n] foreach @{ $m[$i] };}

      ++$lead;}
}

sub display { join("\n" => map join(" " => map(sprintf("%6.2f", $_), @$_)), @{+shift})."\n" }

sub gauss_jordan_invert {
    my(@m) = @_;
    my $rows = @m;
    my @i = identity(scalar @m);
    push @{$m[$_]}, @{$i[$_]} for 0..$rows-1;
    rref(\@m);
    map { splice @$_, 0, $rows } @m;
    @m;
}

sub identity {
    my($n) = @_;
    map { [ (0) x $_, 1, (0) x ($n-1 - $_) ] } 0..$n-1
}

my @tests = (
    [
      [ 2, -1,  0 ],
      [-1,  2, -1 ],
      [ 0, -1,  2 ]
    ],
    [
      [ -1, -2, 3, 2 ],
      [ -4, -1, 6, 2 ],
      [  7, -8, 9, 1 ],
      [  1, -2, 1, 3 ]
    ],
);

for my $matrix (@tests) {
    print "Original Matrix:\n" . display(\@$matrix) . "\n";
    my @gj = gauss_jordan_invert( @$matrix );
    print "Gauss-Jordan Inverted Matrix:\n" . display(\@gj) . "\n";
    my @rt = gauss_jordan_invert( @gj );
    print "After round-trip:\n" . display(\@rt) . "\n";} . "\n"
}</syntaxhighlight>
{{out}}
<pre>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</pre>

=={{header|Phix}}==
{{trans|Kotlin}}
uses ToReducedRowEchelonForm() from [[Reduced_row_echelon_form#Phix]]
<!--<syntaxhighlight lang="phix">(phixonline)-->
 <span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">ToReducedRowEchelonForm</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">M</span><span style="color: #0000FF;">)</span>
 <span style="color: #004080;">integer</span> <span style="color: #000000;">lead</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">,</span>
         <span style="color: #000000;">rowCount</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">M</span><span style="color: #0000FF;">),</span>
         <span style="color: #000000;">columnCount</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">M</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]),</span>
         <span style="color: #000000;">i</span>
     <span style="color: #008080;">for</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">rowCount</span> <span style="color: #008080;">do</span>
         <span style="color: #008080;">if</span> <span style="color: #000000;">lead</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">columnCount</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
         <span style="color: #000000;">i</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">r</span>
         <span style="color: #008080;">while</span> <span style="color: #000000;">M</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">lead</span><span style="color: #0000FF;">]=</span><span style="color: #000000;">0</span> <span style="color: #008080;">do</span>
             <span style="color: #000000;">i</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
             <span style="color: #008080;">if</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">rowCount</span> <span style="color: #008080;">then</span>
                 <span style="color: #000000;">i</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">r</span>
                 <span style="color: #000000;">lead</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
                 <span style="color: #008080;">if</span> <span style="color: #000000;">lead</span><span style="color: #0000FF;">=</span><span style="color: #000000;">columnCount</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
             <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
         <span style="color: #004080;">object</span> <span style="color: #000000;">mr</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_div</span><span style="color: #0000FF;">(</span><span style="color: #000000;">M</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #000000;">M</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">lead</span><span style="color: #0000FF;">])</span>
         <span style="color: #000000;">M</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">M</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">]</span>
         <span style="color: #000000;">M</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mr</span>
         <span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">rowCount</span> <span style="color: #008080;">do</span>
             <span style="color: #008080;">if</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">r</span> <span style="color: #008080;">then</span>
                 <span style="color: #000000;">M</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sq_sub</span><span style="color: #0000FF;">(</span><span style="color: #000000;">M</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">],</span><span style="color: #7060A8;">sq_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">M</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">][</span><span style="color: #000000;">lead</span><span style="color: #0000FF;">],</span><span style="color: #000000;">M</span><span style="color: #0000FF;">[</span><span style="color: #000000;">r</span><span style="color: #0000FF;">]))</span>
             <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
         <span style="color: #000000;">lead</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
     <span style="color: #008080;">return</span> <span style="color: #000000;">M</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 <span style="color: #000080;font-style:italic;">--&lt;/end of copy&gt;</span>
  
 <span style="color: #008080;">function</span> <span style="color: #000000;">inverse</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">mat</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">integer</span> <span style="color: #000000;">len</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mat</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">sequence</span> <span style="color: #000000;">aug</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">len</span><span style="color: #0000FF;">),</span><span style="color: #000000;">len</span><span style="color: #0000FF;">)</span>
     <span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">len</span> <span style="color: #008080;">do</span>
         <span style="color: #008080;">if</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mat</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])!=</span><span style="color: #000000;">len</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span> <span style="color: #000080;font-style:italic;">-- "Not a square matrix"</span>
         <span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">len</span> <span style="color: #008080;">do</span>
             <span style="color: #000000;">aug</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">mat</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
         <span style="color: #000080;font-style:italic;">-- augment by identity matrix to right</span>
         <span style="color: #000000;">aug</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">i</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">len</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
     <span style="color: #000000;">aug</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ToReducedRowEchelonForm</span><span style="color: #0000FF;">(</span><span style="color: #000000;">aug</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">sequence</span> <span style="color: #000000;">inv</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">len</span><span style="color: #0000FF;">),</span><span style="color: #000000;">len</span><span style="color: #0000FF;">)</span>
     <span style="color: #000080;font-style:italic;">-- remove identity matrix to left</span>
     <span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">len</span> <span style="color: #008080;">do</span>
         <span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">len</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #000000;">len</span> <span style="color: #008080;">do</span>
             <span style="color: #000000;">inv</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">-</span><span style="color: #000000;">len</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">aug</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
     <span style="color: #008080;">return</span> <span style="color: #000000;">inv</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
  
 <span style="color: #008080;">constant</span> <span style="color: #000000;">test</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span>  <span style="color: #000000;">0</span><span style="color: #0000FF;">},</span>
                  <span style="color: #0000FF;">{-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span>  <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
                  <span style="color: #0000FF;">{</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span>  <span style="color: #000000;">2</span><span style="color: #0000FF;">}}</span>
 <span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inverse</span><span style="color: #0000FF;">(</span><span style="color: #000000;">test</span><span style="color: #0000FF;">),{</span><span style="color: #004600;">pp_Nest</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">})</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
{{0.75,0.5,0.25},
 {0.5,1,0.5},
 {0.25,0.5,0.75}}
</pre>

===alternate===
{{trans|PowerShell}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
 <span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">gauss_jordan_inv</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">integer</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
     <span style="color: #004080;">sequence</span> <span style="color: #000000;">b</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">),</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
     <span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
     <span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
         <span style="color: #004080;">integer</span> <span style="color: #000000;">lmax</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">k</span>
         <span style="color: #004080;">atom</span> <span style="color: #000000;">amax</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">][</span><span style="color: #000000;">k</span><span style="color: #0000FF;">])</span>
         <span style="color: #008080;">for</span> <span style="color: #000000;">l</span><span style="color: #0000FF;">=</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
             <span style="color: #004080;">atom</span> <span style="color: #000000;">tmp</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">l</span><span style="color: #0000FF;">][</span><span style="color: #000000;">k</span><span style="color: #0000FF;">])</span>
             <span style="color: #008080;">if</span> <span style="color: #000000;">amax</span><span style="color: #0000FF;"><</span><span style="color: #000000;">tmp</span> <span style="color: #008080;">then</span>
                 <span style="color: #0000FF;">{</span><span style="color: #000000;">amax</span><span style="color: #0000FF;">,</span><span style="color: #000000;">lmax</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">tmp</span><span style="color: #0000FF;">,</span><span style="color: #000000;">l</span><span style="color: #0000FF;">}</span>
             <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
         <span style="color: #008080;">if</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">lmax</span> <span style="color: #008080;">then</span>
             <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">lmax</span><span style="color: #0000FF;">]}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">lmax</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]}</span>
             <span style="color: #0000FF;">{</span><span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">lmax</span><span style="color: #0000FF;">]}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">lmax</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]}</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
         <span style="color: #004080;">atom</span> <span style="color: #000000;">akk</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">][</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span>
         <span style="color: #008080;">if</span> <span style="color: #000000;">akk</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">throw</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"Irregular matrix"</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
         <span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
             <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">/=</span> <span style="color: #000000;">akk</span>
             <span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">/=</span> <span style="color: #000000;">akk</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
         <span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
             <span style="color: #008080;">if</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">k</span> <span style="color: #008080;">then</span>
                 <span style="color: #004080;">atom</span> <span style="color: #000000;">aik</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span>
                 <span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
                     <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">aik</span>
                     <span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">-=</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">aik</span>
                 <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
             <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
         <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
     <span style="color: #008080;">return</span> <span style="color: #000000;">b</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
  
 <span style="color: #004080;">sequence</span> <span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span>  <span style="color: #000000;">0</span><span style="color: #0000FF;">},</span>
               <span style="color: #0000FF;">{-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span>  <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">},</span>
               <span style="color: #0000FF;">{</span> <span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span>  <span style="color: #000000;">2</span><span style="color: #0000FF;">}},</span>
          <span style="color: #000000;">inva</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">gauss_jordan_inv</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">))</span>
 <span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"a = "</span><span style="color: #0000FF;">)</span>
 <span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,{</span><span style="color: #004600;">pp_Nest</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_Indent</span><span style="color: #0000FF;">,</span><span style="color: #000000;">4</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_IntFmt</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%2d"</span><span style="color: #0000FF;">})</span>
 <span style="color: #7060A8;">puts</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"inv(a) = "</span><span style="color: #0000FF;">)</span>
 <span style="color: #7060A8;">pp</span><span style="color: #0000FF;">(</span><span style="color: #000000;">inva</span><span style="color: #0000FF;">,{</span><span style="color: #004600;">pp_Nest</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_Indent</span><span style="color: #0000FF;">,</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,</span><span style="color: #004600;">pp_FltFmt</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%4.2f"</span><span style="color: #0000FF;">})</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
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}}
</pre>

=={{header|PL/I}}==
<pre>
/* 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;
</pre>
Output:
<pre>
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
</pre>

=={{header|PowerShell}}==
<syntaxhighlight lang="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

</syntaxhighlight>
<b>Output:</b>
<pre>
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
</pre>

=={{header|Python}}==

Use numpy matrix inversion

<syntaxhighlight lang="python">
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)
</syntaxhighlight>

{{out}}
<pre>
[[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]]
</pre>

=={{header|Racket}}==

Using the matrix library that comes with default Racket installation

<syntaxhighlight lang="racket">#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)</syntaxhighlight>

{{out}}
<pre>
(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]])
</pre>

=={{header|Raku}}==
(formerly Perl 6)

Uses bits and pieces from other tasks, [[Reduced_row_echelon_form#Raku|Reduced row echelon form]] primarily.

<syntaxhighlight lang="raku" line>sub gauss-jordan-invert (@m where &is-square) {
    ^@m .map: { @m[$_].append: identity(+@m)[$_] };
    @m.&rref[*]»[+@m .. *];
}

sub is-square (@m) { so @m == all @m[*] }

sub identity ($n) { [ 1, |(0 xx $n-1) ], *.rotate(-1).Array ... *.tail }

# reduced row echelon form (from 'Gauss-Jordan elimination' task)
sub rref (@m) {
    my ($lead, $rows, $cols) = 0, @m, @m[0];

    for ^$rows -> $r {
        $lead < $cols or return @m;
        my $i = $r;
        until @m[$i;$lead] {
            ++$i == $rows or next;
            $i = $r;
            ++$lead == $cols and return @m;
        }
        @m[$i, $r] = @m[$r, $i] if $r != $i;
        @m[$r] »/=» $ = @m[$r;$lead];
        for ^$rows -> $n {
            next if $n == $r;
            @m[$n] »-=» @m[$r] »×» (@m[$n;$lead] // 0);
        }
        ++$lead;
    }
    @m
}
sub rat-or-int ($num) {
    return $num unless $num ~~ Rat|FatRat;
    return $num.narrow if $num.narrow ~~ Int;
    $num.nude.join: '/';
}

sub say_it ($message, @array) {
    my $max;
    @array.map: {$max max= max $_».&rat-or-int.comb(/\S+/)».chars};
    say "\n$message";
    $_».&rat-or-int.fmt(" %{$max}s").put for @array;
}

multi to-matrix ($str) { [$str.split(';').map(*.words.Array)] }
multi to-matrix (@array) { @array }

sub hilbert-matrix ($h) {
    [ (1..$h).map( -> $n { [ ($n ..^ $n + $h).map: { (1/$_).FatRat } ] } ) ]
}

my @tests =
    '1 2 3; 4 1 6; 7 8 9',
    '2 -1 0; -1 2 -1; 0 -1 2',
    '-1 -2 3 2; -4 -1 6 2; 7 -8 9 1; 1 -2 1 3',
    '1 2 3 4; 5 6 7 8; 9 33 11 12; 13 14 15 17',
    '3 1 8 9 6; 6 2 8 10 1; 5 7 2 10 3; 3 2 7 7 9; 3 5 6 1 1',

    # Test with a Hilbert matrix
    hilbert-matrix 10;

@tests.map: {
    my @matrix = .&to-matrix;
    say_it( " {'=' x 20} Original Matrix: {'=' x 20}", @matrix );
    say_it( ' Gauss-Jordan Inverted:', my @invert = gauss-jordan-invert @matrix );
    say_it( ' Re-inverted:', gauss-jordan-invert @invert».Array );
}</syntaxhighlight>

{{out}}
<div style="font-size:75%">
<pre> ==================== 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</pre></div>

=={{header|REXX}}==
===version 1===
<syntaxhighlight lang="rexx">/* REXX */
Parse Arg seed nn
If seed='' Then
  seed=23345
If nn='' Then nn=5
If seed='?' Then Do
  Say 'rexx gjmi seed n computes a random matrix with n rows and columns'
  Say 'Default is 23345 5'
  Exit
  End
Numeric Digits 50
Call random 1,2,seed
a=''
Do i=1 To nn**2
  a=a random(9)+1
  End
n2=words(a)
Do n=2 To n2/2
  If n**2=n2 Then
    Leave
  End
If n>n2/2 Then
  Call exit 'Not a square matrix:' a '('n2 'elements).'
det=determinante(a,n)
If det=0 Then
  Call exit 'Determinant is 0'
Do j=1 To n
  Do i=1 To n
    Parse Var A a.i.j a
    aa.i.j=a.i.j
    End
  Do ii=1 To n
    z=(ii=j)
    iii=ii+n
    a.iii.j=z
    End
  End
Call show 1,'The given matrix'
Do m=1 To n-1
  If a.m.m=0 Then Do
    Do j=m+1 To n
      If a.m.j<>0 Then Leave
      End
    If j>n Then Do
      Say 'No pivot>0 found in column' m
      Exit
      End
    Do i=1 To n*2
      temp=a.i.m
      a.i.m=a.i.j
      a.i.j=temp
      End
    End
  Do j=m+1 To n
    If a.m.j<>0 Then Do
      jj=m
      fact=divide(a.m.m,a.m.j)
      Do i=1 To n*2
        a.i.j=subtract(multiply(a.i.j,fact),a.i.jj)
        End
      End
    End
  Call show 2 m
  End
Say 'Lower part has all zeros'
Say ''

Do j=1 To n
  If denom(a.j.j)<0 Then Do
    Do i=1 To 2*n
      a.i.j=subtract(0,a.i.j)
      End
    End
  End
Call show 3

Do m=n To 2 By -1
  Do j=1 To m-1
    jj=m
    fact=divide(a.m.j,a.m.jj)
    Do i=1 To n*2
      a.i.j=subtract(a.i.j,multiply(a.i.jj,fact))
      End
    End
  Call show 4 m
  End
Say 'Upper half has all zeros'
Say ''
Do j=1 To n
  If decimal(a.j.j)<>1 Then Do
    z=a.j.j
    Do i=1 To 2*n
      a.i.j=divide(a.i.j,z)
      End
    End
  End
Call show 5
Say 'Main diagonal has all ones'
Say ''

Do j=1 To n
  Do i=1 To n
    z=i+n
    a.i.j=a.z.j
    End
  End
Call show 6,'The inverse matrix'

do i = 1 to n
  do j = 1 to n
    sum=0
    Do k=1 To n
      sum=add(sum,multiply(aa.i.k,a.k.j))
      End
    c.i.j = sum
    end
  End
Call showc 7,'The product of input and inverse matrix'
Exit

show:
  Parse Arg num,text
  Say 'show' arg(1) text
  If arg(1)<>6 Then rows=n*2
               Else rows=n
  len=0
  Do j=1 To n
    Do i=1 To rows
      len=max(len,length(a.i.j))
      End
    End
  Do j=1 To n
    ol=''
    Do i=1 To rows
      ol=ol||right(a.i.j,len+1)
      End
    Say ol
    End
  Say ''
  Return

showc:
  Parse Arg num,text
  Say text
  clen=0
  Do j=1 To n
    Do i=1 To n
      clen=max(clen,length(c.i.j))
      End
    End
  Do j=1 To n
    ol=''
    Do i=1 To n
      ol=ol||right(c.i.j,clen+1)
      End
    Say ol
    End
  Say ''
  Return

denom: Procedure
  /* Return the denominator */
  Parse Arg d '/' n
  Return d

decimal: Procedure
  /* compute the fraction's value */
  Parse Arg a
  If pos('/',a)=0 Then a=a'/1'; Parse Var a ad '/' an
  Return ad/an

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

add: Procedure
  Parse Arg a,b
  If pos('/',a)=0 Then a=a'/1'; Parse Var a ad '/' an
  If pos('/',b)=0 Then b=b'/1'; Parse Var b bd '/' bn
  sum=divide(ad*bn+bd*an,an*bn)
  Return sum

multiply: Procedure
  Parse Arg a,b
  If pos('/',a)=0 Then a=a'/1'; Parse Var a ad '/' an
  If pos('/',b)=0 Then b=b'/1'; Parse Var b bd '/' bn
  prd=divide(ad*bd,an*bn)
  Return prd

subtract: Procedure
  Parse Arg a,b
  If pos('/',a)=0 Then a=a'/1'; Parse Var a ad '/' an
  If pos('/',b)=0 Then b=b'/1'; Parse Var b bd '/' bn
  div=divide(ad*bn-bd*an,an*bn)
  Return div

divide: Procedure
  Parse Arg a,b
  If pos('/',a)=0 Then a=a'/1'; Parse Var a ad '/' an
  If pos('/',b)=0 Then b=b'/1'; Parse Var b bd '/' bn
  sd=ad*bn
  sn=an*bd
  g=gcd(sd,sn)
  Select
    When sd=0 Then res='0'
    When abs(sn/g)=1 Then res=(sd/g)*sign(sn/g)
    Otherwise Do
      den=sd/g
      nom=sn/g
      If nom<0 Then Do
        If den<0 Then den=abs(den)
        Else den=-den
        nom=abs(nom)
        End
      res=den'/'nom
      End
    End
  Return res

determinante: Procedure
/* REXX ***************************************************************
* determinant.rex
* compute the determinant of the given square matrix
* Input: as: the representation of the matrix as vector (n**2 elements)
* 21.05.2013 Walter Pachl
**********************************************************************/
  Parse Arg as,n
  Do i=1 To n
    Do j=1 To n
      Parse Var as a.i.j as
      End
    End
  Select
    When n=2 Then det=subtract(multiply(a.1.1,a.2.2),multiply(a.1.2,a.2.1))
    When n=3 Then Do
      det=multiply(multiply(a.1.1,a.2.2),a.3.3)
      det=add(det,multiply(multiply(a.1.2,a.2.3),a.3.1))
      det=add(det,multiply(multiply(a.1.3,a.2.1),a.3.2))
      det=subtract(det,multiply(multiply(a.1.3,a.2.2),a.3.1))
      det=subtract(det,multiply(multiply(a.1.2,a.2.1),a.3.3))
      det=subtract(det,multiply(multiply(a.1.1,a.2.3),a.3.2))
      End
    Otherwise Do
      det=0
      Do k=1 To n
        sign=((-1)**(k+1))
        If sign=1 Then
          det=add(det,multiply(a.1.k,determinante(subm(k),n-1)))
        Else
          det=subtract(det,multiply(a.1.k,determinante(subm(k),n-1)))
        End
      End
    End
  Return det

subm: Procedure Expose a. n
/**********************************************************************
* compute the submatrix resulting when row 1 and column k are removed
* Input: a.*.*, k
* Output: bs the representation of the submatrix as vector
**********************************************************************/
  Parse Arg k
  bs=''
  do i=2 To n
    Do j=1 To n
      If j=k Then Iterate
      bs=bs a.i.j
      End
    End
  Return bs

Exit: Say arg(1)</syntaxhighlight>
{{out}} Using the defaults for seed and n
<pre>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</pre>

===version 2===
{{trans|PL/I}}
<syntaxhighlight lang="rexx">/*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</syntaxhighlight>
{{out|output|text=&nbsp; when using the internal default input:}}
<pre>
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
</pre>

=={{header|Ruby}}==
<syntaxhighlight lang="ruby">require 'matrix'

m = Matrix[[-1, -2, 3, 2],
           [-4, -1, 6, 2],
           [ 7, -8, 9, 1],
           [ 1, -2, 1, 3]]

pp m.inv.row_vectors
</syntaxhighlight>
{{out}}
<pre>[Vector[(-21/23), (17/69), (13/138), (19/46)],
 Vector[(-38/23), (15/23), (1/23), (15/23)],
 Vector[(-16/23), (25/69), (11/138), (9/46)],
 Vector[(-13/23), (16/69), (-2/69), (13/23)]]
</pre>

=={{header|Rust}}==
<syntaxhighlight lang="rust">
fn main() {
    let mut a: Vec<Vec<f64>> = vec![vec![1.0, 2.0, 3.0],
                                vec![4.0, 1.0, 6.0],
                                vec![7.0, 8.0, 9.0]
                                ];
    let mut b: Vec<Vec<f64>> = vec![vec![2.0, -1.0, 0.0],
                                vec![-1.0, 2.0, -1.0],
                                vec![0.0, -1.0, 2.0]
                                ];

    let mut ref_a = &mut a;
    let rref_a = &mut ref_a;
    let mut ref_b = &mut b;
    let rref_b = &mut ref_b;

    println!("Matrix A:\n");
    print_matrix(rref_a);
    println!("\nInverse of Matrix A:\n");
    print_matrix(&mut matrix_inverse(rref_a));
    println!("\n\nMatrix B:\n");
    print_matrix(rref_b);
    println!("\nInverse of Matrix B:\n");
    print_matrix(&mut matrix_inverse(rref_b));
}

//Begin Matrix Inversion
fn matrix_inverse(matrix: &mut Vec<Vec<f64>>) -> Vec<Vec<f64>>{
    let len = matrix.len();
    let mut aug = zero_matrix(len, len * 2);
    for i in 0..len {
        for j in 0.. len {
            aug[i][j] = matrix[i][j];
        }
        aug[i][i + len] = 1.0;
    }

    gauss_jordan_general(&mut aug);
    
    
    let mut unaug = zero_matrix(len, len);
    for i in 0..len {
        for j in 0..len {
            unaug[i][j] = aug[i][j+len];
        }
    }
    unaug
}
//End Matrix Inversion

//Begin Generalised Reduced Row Echelon Form
fn gauss_jordan_general(matrix: &mut Vec<Vec<f64>>) {
    let mut lead = 0;
    let row_count = matrix.len();
    let col_count = matrix[0].len();

    for r in 0..row_count {
        if col_count <= lead {
            break;
        }
        let mut i = r;
        while matrix[i][lead] == 0.0 {
            i = i + 1;
            if row_count == i {
                i = r;
                lead = lead + 1;
                if col_count == lead {
                    break;
                }
            }
        }

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

        if matrix[r][lead] != 0.0 {
            let div = matrix[r][lead];
            for j in 0..col_count {
                matrix[r][j] = matrix[r][j] / div;
            }
        }

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

    }
    //matrix.to_owned()
}

fn zero_matrix(rows: usize, cols: usize) -> Vec<Vec<f64>> {
    let mut matrix = Vec::with_capacity(cols);
    for _ in 0..rows {
        let mut col: Vec<f64> = Vec::with_capacity(rows);
        for _ in 0..cols {
            col.push(0.0);
        }
        matrix.push(col);
    }
    matrix
}

fn print_matrix(mat: &mut Vec<Vec<f64>>) {
    for row in 0..mat.len(){
        println!("{:?}", mat[row]);
    }
}
</syntaxhighlight>
{{out}}
<pre>
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]
</pre>

=={{header|Sidef}}==
Uses the '''rref(M)''' function from [https://rosettacode.org/wiki/Reduced_row_echelon_form#Sidef Reduced row echelon form].
{{trans|Raku}}
<syntaxhighlight lang="ruby">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")</syntaxhighlight>
{{out}}
<pre>
-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
</pre>

=={{header|Wren}}==
{{libheader|Wren-fmt}}
{{libheader|Wren-matrix}}
The above module does in fact use the Gauss-Jordan method for matrix inversion.
<syntaxhighlight lang="wren">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()
}</syntaxhighlight>

{{out}}
<pre>
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|
</pre>

=={{header|zkl}}==
This uses GSL to invert a matrix via LU decomposition, not Gauss-Jordan.
<syntaxhighlight lang="zkl">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();</syntaxhighlight>
{{out}}
<pre>
   -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
</pre>
<syntaxhighlight lang="zkl">m:=GSL.Matrix(3,3).set(2,-1,0, -1,2,-1, 0,-1,2);
m.invert().format(10,4).println("\n");</syntaxhighlight>
{{out}}
<pre>
    0.7500,    0.5000,    0.2500
    0.5000,    1.0000,    0.5000
    0.2500,    0.5000,    0.7500
</pre>

{{omit from|PL/0}}