Gauss-Jordan matrix inversion: Difference between revisions

From Rosetta Code
Content added Content deleted
(→‎{{header|Go}}: Trans PowerShell)
(Put Generic at the right place)
Line 480: Line 480:
</pre>
</pre>


=={{header|Go}}==
{{trans|Kotlin}}
<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")
}</lang>
{{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|PoweShell}}
<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++ {
a := make(vector, n)
a[i] = 1
b[i] = a
}
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 {
row := make([]float64, n)
copy(row, v)
a[i] = row
}
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
}
}
// scale the abs used to pick the pivot.
if abs := math.Abs(row[k]) / s; abs > max {
iMax = i
max = abs
}
}
if a[iMax][k] == 0 {
panic("Irregular matrix")
}
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")
}</lang>
{{out}}
<pre>Same output than above</pre>

=={{header|Haskell}}==
<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</lang>
{{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|Generic}}==
=={{header|Generic}}==
<lang cpp>
<lang cpp>
Line 1,416: Line 1,154:
}
}
</lang>
</lang>
=={{header|Go}}==
{{trans|Kotlin}}
<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")
}</lang>
{{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|PoweShell}}
<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++ {
a := make(vector, n)
a[i] = 1
b[i] = a
}
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 {
row := make([]float64, n)
copy(row, v)
a[i] = row
}
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
}
}
// scale the abs used to pick the pivot.
if abs := math.Abs(row[k]) / s; abs > max {
iMax = i
max = abs
}
}
if a[iMax][k] == 0 {
panic("Irregular matrix")
}
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")
}</lang>
{{out}}
<pre>Same output than above</pre>

=={{header|Haskell}}==
<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</lang>
{{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}}==
=={{header|J}}==

Revision as of 13:30, 15 January 2021

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

Invert matrix   A   using Gauss-Jordan method.

A   being an   n × n   matrix.

C++

<lang cpp>#include <algorithm>

  1. include <cassert>
  2. include <iomanip>
  3. include <iostream>
  4. 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;

}</lang>

Output:
Matrix:
 2.0000 -1.0000  0.0000
-1.0000  2.0000 -1.0000
 0.0000 -1.0000  2.0000
Inverse:
 0.7500  0.5000  0.2500
 0.5000  1.0000  0.5000
 0.2500  0.5000  0.7500
Product of matrix and inverse:
 1.0000  0.0000  0.0000
-0.0000  1.0000 -0.0000
 0.0000  0.0000  1.0000
Inverse of inverse:
 2.0000 -1.0000  0.0000
-1.0000  2.0000 -1.0000
 0.0000 -1.0000  2.0000

C#

<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();
           }
       }
   }

} </lang> <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();
       }
   }

} </lang>

Output:


-0.913043478260869  0.246376811594203  0.0942028985507246  0.413043478260869
-1.65217391304348  0.652173913043478  0.0434782608695652  0.652173913043478
-0.695652173913043  0.36231884057971  0.0797101449275362  0.195652173913043
-0.565217391304348  0.231884057971014  -0.0289855072463768  0.565217391304348

Factor

Normally you would call recip to calculate the inverse of a matrix, but it uses a different method than Gauss-Jordan, so here's Gauss-Jordan.

Works with: Factor version 0.99 2020-01-23

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

Output:
-21/23   17/69 13/138 19/46
-1-15/23 15/23 1/23   15/23
-16/23   25/69 11/138 9/46
-13/23   16/69 -2/69  13/23

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!

<lang freebasic>#include once "matmult.bas"

  1. 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</lang>

Output:
 1             5             7            
 3             5             9            
 5             9             11           

-0.5416666666666667          0.1666666666666667          0.2083333333333333         
 0.2499999999999999         -0.5           0.25         
 0.0416666666666667          0.3333333333333333         -0.2083333333333333

Generic

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

} </lang>

Go

Translation of: Kotlin

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

}</lang>

Output:
Inverse of A is:

[-0.812500  0.125000  0.187500]
[ 0.125000 -0.250000  0.125000]
[ 0.520833  0.125000 -0.145833]

Inverse of B is:

[ 0.750000  0.500000  0.250000]
[ 0.500000  1.000000  0.500000]
[ 0.250000  0.500000  0.750000]

Alternative

Translation of: PoweShell

<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++ {
       a := make(vector, n)
       a[i] = 1
       b[i] = a
   }
   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 {
       row := make([]float64, n)
       copy(row, v)
       a[i] = row
  }
  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
              }
          }
          // scale the abs used to pick the pivot.
          if abs := math.Abs(row[k]) / s; abs > max {
              iMax = i
              max = abs
          }
      }
      if a[iMax][k] == 0 {
          panic("Irregular matrix")
      }
      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")

}</lang>

Output:
Same output than above

Haskell

<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</lang>
Output:
inversion a =
[-0.8125,0.125,0.1875]
[0.125,-0.25,0.125]
[0.5208333333333334,0.125,-0.14583333333333334]

inversion b =
[0.75,0.5,0.25]
[0.5,1.0,0.5]
[0.25,0.5,0.75]


J

Solution:

Uses Gauss-Jordan implementation (as described in Reduced_row_echelon_form#J) to find reduced row echelon form of the matrix after augmenting with an identity matrix. <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</lang>

Usage: <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</lang>

Java

<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();
   }

}</lang>

<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.print(String.format("%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;
   }

}</lang>

Output:
Matrix:
  0.913   0.438   0.002   0.053   0.588
  0.739   0.043   0.593   0.599   0.115
  0.542   0.734   0.273   0.889   0.921
  0.503   0.960   0.707   0.088   0.921
  0.777   0.829   0.190   0.673   0.728
Inverse:
  0.873   0.502  -0.705  -0.272   0.452
 -1.361  -0.578  -2.118   0.400   3.368
 -0.539   0.883  -0.106   0.910  -0.722
 -0.720   0.291   0.625  -0.628   0.539
  1.425  -0.377   2.616   0.178  -3.257
Product of matrix and inverse:
  1.000   0.000   0.000  -0.000   0.000
 -0.000   1.000  -0.000  -0.000   0.000
  0.000   0.000   1.000  -0.000  -0.000
 -0.000   0.000   0.000   1.000  -0.000
 -0.000   0.000   0.000  -0.000   1.000

Julia

Works with: Julia version 0.6

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

Native implementation: <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)</lang>

Output:
I / A = [-0.8125 0.125 0.1875; 0.125 -0.25 0.125; 0.520833 0.125 -0.145833]
inv(A) = [-0.8125 0.125 0.1875; 0.125 -0.25 0.125; 0.520833 0.125 -0.145833]
gaussjordan(A) = [-0.8125 0.125 0.1875; 0.125 -0.25 0.125; 0.520833 0.125 -0.145833]

Kotlin

This follows the description of Gauss-Jordan elimination in Wikipedia whereby the original square matrix is first augmented to the right by its identity matrix, its reduced row echelon form is then found and finally the identity matrix to the left is removed to leave the inverse of the original square matrix. <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")

}</lang>

Output:
Inverse of A is :

 -0.812500    0.125000    0.187500  
  0.125000   -0.250000    0.125000  
  0.520833    0.125000   -0.145833  

Inverse of B is :

  0.750000    0.500000    0.250000  
  0.500000    1.000000    0.500000  
  0.250000    0.500000    0.750000

Lambdatalk

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

</lang>

MATLAB

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

Output:
Matrix:
     2    -1     0
    -1     2    -1
     0    -1     2

Inverted matrix:
    0.7500    0.5000    0.2500
    0.5000    1.0000    0.5000
    0.2500    0.5000    0.7500

Inverted matrix calculated by built-in function:
    0.7500    0.5000    0.2500
    0.5000    1.0000    0.5000
    0.2500    0.5000    0.7500

Product of matrix and inverse:
    1.0000         0         0
   -0.0000    1.0000   -0.0000
   -0.0000   -0.0000    1.0000

Nim

We reuse the algorithm of task https://rosettacode.org/wiki/Reduced_row_echelon_form (version using floats). <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


  1. ———————————————————————————————————————————————————————————————————————————————————————————————————

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

Output:
Matrix:
  1.00000   2.00000   3.00000
  4.00000   1.00000   6.00000
  7.00000   8.00000   9.00000

Inverse:
 -0.81250   0.12500   0.18750
  0.12500  -0.25000   0.12500
  0.52083   0.12500  -0.14583


Matrix:
  2.00000  -1.00000   0.00000
 -1.00000   2.00000  -1.00000
  0.00000  -1.00000   2.00000

Inverse:
  0.75000   0.50000   0.25000
  0.50000   1.00000   0.50000
  0.25000   0.50000   0.75000


Matrix:
 -1.00000  -2.00000   3.00000   2.00000
 -4.00000  -1.00000   6.00000   2.00000
  7.00000  -8.00000   9.00000   1.00000
  1.00000  -2.00000   1.00000   3.00000

Inverse:
 -0.91304   0.24638   0.09420   0.41304
 -1.65217   0.65217   0.04348   0.65217
 -0.69565   0.36232   0.07971   0.19565
 -0.56522   0.23188  -0.02899   0.56522

Perl

Included code from Reduced row echelon form task. <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"

}</lang>

Output:
Original Matrix:
  2.00  -1.00   0.00
 -1.00   2.00  -1.00
  0.00  -1.00   2.00

Gauss-Jordan Inverted Matrix:
  0.75   0.50   0.25
  0.50   1.00   0.50
  0.25   0.50   0.75

After round-trip:
  2.00  -1.00   0.00
 -1.00   2.00  -1.00
  0.00  -1.00   2.00

Original Matrix:
 -1.00  -2.00   3.00   2.00
 -4.00  -1.00   6.00   2.00
  7.00  -8.00   9.00   1.00
  1.00  -2.00   1.00   3.00

Gauss-Jordan Inverted Matrix:
 -0.91   0.25   0.09   0.41
 -1.65   0.65   0.04   0.65
 -0.70   0.36   0.08   0.20
 -0.57   0.23  -0.03   0.57

After round-trip:
 -1.00  -2.00   3.00   2.00
 -4.00  -1.00   6.00   2.00
  7.00  -8.00   9.00   1.00
  1.00  -2.00   1.00   3.00

Phix

Translation of: Kotlin

uses ToReducedRowEchelonForm() from Reduced_row_echelon_form#Phix <lang Phix>function inverse(sequence mat) 

   integer len = length(mat)
   sequence aug = repeat(repeat(0,2*len),len)
   for i=1 to len do
       if length(mat[i])!=len then ?9/0 end if -- "Not a square matrix"
       for j=1 to len do
           aug[i][j] = mat[i][j]
       end for
       -- augment by identity matrix to right
       aug[i][i + len] = 1
   end for
   aug = ToReducedRowEchelonForm(aug)
   sequence inv = repeat(repeat(0,len),len)
   -- remove identity matrix to left
   for i=1 to len do
       for j=len+1 to 2*len do
           inv[i][j-len] = aug[i][j]
       end for
   end for
   return inv

end function

constant test = {{ 2, -1, 0}, 

                {-1,  2, -1},
                { 0, -1,  2}}

pp(inverse(test),{pp_Nest,1})</lang>

Output:
{{0.75,0.5,0.25},
 {0.5,1,0.5},
 {0.25,0.5,0.75}}

alternate

Translation of: PowerShell

<lang Phix>function gauss_jordan_inv(sequence a) 

   integer n = length(a)
   sequence b = repeat(repeat(0,n),n)
   for i=1 to n do b[i,i] = 1 end for
   for k=1 to n do
       integer lmax = k
       atom amax = abs(a[k][k])
       for l=k+1 to n do
           atom tmp = abs(a[l][k])
           if amax<tmp then
               {amax,lmax} = {tmp,l}
           end if
       end for
       if k!=lmax then
           {a[k], a[lmax]} = {a[lmax], a[k]}
           {b[k], b[lmax]} = {b[lmax], b[k]}
       end if
       atom akk = a[k][k]
       if akk=0 then throw("Irregular matrix") end if
       for j=1 to n do
           a[k][j] /= akk
           b[k][j] /= akk
       end for
       for i=1 to n do
           if i!=k then
               atom aik = a[i][k]
               for j=1 to n do
                   a[i][j] -= a[k][j]*aik
                   b[i][j] -= b[k][j]*aik
               end for
           end if
       end for
   end for
   return b

end function

sequence a = {{ 2, -1, 0}, 

             {-1,  2, -1},
             { 0, -1,  2}},
        inva = gauss_jordan_inv(a)

puts(1,"a = ") pp(a,{pp_Nest,1,pp_Indent,4,pp_IntFmt,"%2d"}) puts(1,"inv(a) = ") pp(inva,{pp_Nest,1,pp_Indent,9,pp_FltFmt,"%4.2f"})</lang>

Output:
a = {{ 2,-1, 0},
     {-1, 2,-1},
     { 0,-1, 2}}
inv(a) = {{0.75,0.50,0.25},
          {0.50,1.00,0.50},
          {0.25,0.50,0.75}}

PL/I

/* Gauss-Jordan matrix inversion */

G_J: procedure options (main);       /* 4 November 2020 */

   declare t float;
   declare (i, j, k, n) fixed binary;

   open file (sysin) title ('/GAUSSJOR.DAT');

   get (n);                     /* Read in the order of the matrix. */
   put skip data (n);
   begin;
      declare a(n,n)float, aux(n,n) float;

      aux = 0;
      do i = 1 to n; aux(i,i) = 1; end;
      get (a);                  /* Read in the matrix. */
      put skip list ('The matrix to be inverted is:');
      put edit (a) ( skip, (n) F(10,4)); /* Print  the matrix. */

      do k = 1 to n;
         /* Divide row k by a(k,k) */
         t = a(k,k); a(k,*) = a(k,*) / t; aux(k,*) = aux(k,*) / t;

         do i = 1 to k-1, k+1 to n; /* Work down the rows. */
            t = a(i,k);
            a(i,*) = a(i,*) - t*a(k,*); aux(i,*) = aux(i,*)  - t*aux(k,*);
         end;
      end;
      put skip (2) list ('The inverse is:');
      put edit (aux) ( skip, (n) F(10,4));
   end; /* of BEGIN block */
end G_J;

Output: 

N=        4;
The matrix to be inverted is: 
    8.0000    3.0000    7.0000    5.0000
    9.0000   12.0000   10.0000   11.0000
    6.0000    2.0000    4.0000   13.0000
   14.0000   15.0000   16.0000   17.0000

The inverse is: 
    0.3590    0.5385    0.0769   -0.5128
   -0.0190    0.3419   -0.0199   -0.2004
   -0.1833   -0.7009   -0.1425    0.6163
   -0.1064   -0.0855    0.0883    0.0779

PowerShell

<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

</lang> Output: 

a =
1 2 3
4 1 6
7 8 9

inv(a) =
-0.8125 0.125 0.1875
0.125 -0.25 0.125
0.520833333333333 0.125 -0.145833333333333

b =
2 -1 0
-1 2 -1
0 -1 2

inv(b) =
0.75 0.5 0.25
0.5 1 0.5
0.25 0.5 0.75

Python

Use numpy matrix inversion

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

Output:
[[1. 2. 3.]
 [4. 1. 6.]
 [7. 8. 9.]]
[[-0.8125      0.125       0.1875    ]
 [ 0.125      -0.25        0.125     ]
 [ 0.52083333  0.125      -0.14583333]]

Racket

Using the matrix library that comes with default Racket installation

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

Output:
(array #[#[-13/16 1/8 3/16] #[1/8 -1/4 1/8] #[25/48 1/8 -7/48]])
(array #[#[-13/16 1/8 3/16] #[1/8 -1/4 1/8] #[25/48 1/8 -7/48]])

Raku

(formerly Perl 6)

Works with: Rakudo version 2019.03.1

Uses bits and pieces from other tasks, Reduced row echelon form primarily.

<lang perl6>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 }

  1. reduced row echelon form (Gauss-Jordan elimination)

sub rref (@m) { 

   return unless @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;
       my $lv = @m[$r;$lead];
       @m[$r] »/=» $lv;
       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.WHAT ~~ 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 );

}</lang>

Output:
 ==================== Original Matrix: ====================
 1  2  3
 4  1  6
 7  8  9

 Gauss-Jordan Inverted:
 -13/16     1/8    3/16
    1/8    -1/4     1/8
  25/48     1/8   -7/48

 Re-inverted:
 1  2  3
 4  1  6
 7  8  9

 ==================== Original Matrix: ====================
  2  -1   0
 -1   2  -1
  0  -1   2

 Gauss-Jordan Inverted:
 3/4  1/2  1/4
 1/2    1  1/2
 1/4  1/2  3/4

 Re-inverted:
  2  -1   0
 -1   2  -1
  0  -1   2

 ==================== Original Matrix: ====================
 -1  -2   3   2
 -4  -1   6   2
  7  -8   9   1
  1  -2   1   3

 Gauss-Jordan Inverted:
 -21/23   17/69  13/138   19/46
 -38/23   15/23    1/23   15/23
 -16/23   25/69  11/138    9/46
 -13/23   16/69   -2/69   13/23

 Re-inverted:
 -1  -2   3   2
 -4  -1   6   2
  7  -8   9   1
  1  -2   1   3

 ==================== Original Matrix: ====================
  1   2   3   4
  5   6   7   8
  9  33  11  12
 13  14  15  17

 Gauss-Jordan Inverted:
   19/184  -199/184     -1/46       1/2
     1/23     -2/23      1/23         0
 -441/184   813/184     -1/46      -3/2
        2        -3         0         1

 Re-inverted:
  1   2   3   4
  5   6   7   8
  9  33  11  12
 13  14  15  17

 ==================== Original Matrix: ====================
  3   1   8   9   6
  6   2   8  10   1
  5   7   2  10   3
  3   2   7   7   9
  3   5   6   1   1

 Gauss-Jordan Inverted:
 -4525/6238   2529/6238   -233/3119   1481/3119   -639/6238
  1033/6238  -1075/6238    342/3119   -447/3119    871/6238
  1299/6238   -289/6238   -204/3119   -390/3119    739/6238
   782/3119   -222/3119    237/3119   -556/3119   -177/3119
  -474/3119    -17/3119    -24/3119    688/3119   -140/3119

 Re-inverted:
  3   1   8   9   6
  6   2   8  10   1
  5   7   2  10   3
  3   2   7   7   9
  3   5   6   1   1

 ==================== Original Matrix: ====================
    1   1/2   1/3   1/4   1/5   1/6   1/7   1/8   1/9  1/10
  1/2   1/3   1/4   1/5   1/6   1/7   1/8   1/9  1/10  1/11
  1/3   1/4   1/5   1/6   1/7   1/8   1/9  1/10  1/11  1/12
  1/4   1/5   1/6   1/7   1/8   1/9  1/10  1/11  1/12  1/13
  1/5   1/6   1/7   1/8   1/9  1/10  1/11  1/12  1/13  1/14
  1/6   1/7   1/8   1/9  1/10  1/11  1/12  1/13  1/14  1/15
  1/7   1/8   1/9  1/10  1/11  1/12  1/13  1/14  1/15  1/16
  1/8   1/9  1/10  1/11  1/12  1/13  1/14  1/15  1/16  1/17
  1/9  1/10  1/11  1/12  1/13  1/14  1/15  1/16  1/17  1/18
 1/10  1/11  1/12  1/13  1/14  1/15  1/16  1/17  1/18  1/19

 Gauss-Jordan Inverted:
            100           -4950           79200         -600600         2522520        -6306300         9609600        -8751600         4375800         -923780
          -4950          326700        -5880600        47567520      -208107900       535134600      -832431600       770140800      -389883780        83140200
          79200        -5880600       112907520      -951350400      4281076800    -11237826600     17758540800    -16635041280      8506555200     -1829084400
        -600600        47567520      -951350400      8245036800    -37875637800    101001700800   -161602721280    152907955200    -78843164400     17071454400
        2522520      -208107900      4281076800    -37875637800    176752976400   -477233036280    771285715200   -735869534400    382086104400    -83223340200
       -6306300       535134600    -11237826600    101001700800   -477233036280   1301544644400  -2121035716800   2037792556800  -1064382719400    233025352560
        9609600      -832431600     17758540800   -161602721280    771285715200  -2121035716800   3480673996800  -3363975014400   1766086882560   -388375587600
       -8751600       770140800    -16635041280    152907955200   -735869534400   2037792556800  -3363975014400   3267861442560  -1723286307600    380449555200
        4375800      -389883780      8506555200    -78843164400    382086104400  -1064382719400   1766086882560  -1723286307600    912328045200   -202113826200
        -923780        83140200     -1829084400     17071454400    -83223340200    233025352560   -388375587600    380449555200   -202113826200     44914183600

 Re-inverted:
    1   1/2   1/3   1/4   1/5   1/6   1/7   1/8   1/9  1/10
  1/2   1/3   1/4   1/5   1/6   1/7   1/8   1/9  1/10  1/11
  1/3   1/4   1/5   1/6   1/7   1/8   1/9  1/10  1/11  1/12
  1/4   1/5   1/6   1/7   1/8   1/9  1/10  1/11  1/12  1/13
  1/5   1/6   1/7   1/8   1/9  1/10  1/11  1/12  1/13  1/14
  1/6   1/7   1/8   1/9  1/10  1/11  1/12  1/13  1/14  1/15
  1/7   1/8   1/9  1/10  1/11  1/12  1/13  1/14  1/15  1/16
  1/8   1/9  1/10  1/11  1/12  1/13  1/14  1/15  1/16  1/17
  1/9  1/10  1/11  1/12  1/13  1/14  1/15  1/16  1/17  1/18
 1/10  1/11  1/12  1/13  1/14  1/15  1/16  1/17  1/18  1/19

REXX

version 1

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

Output:

Using the defaults for seed and n 

show 1 The given matrix
 10  3  8  6  3  1  0  0  0  0
  5  7  8  8  2  0  1  0  0  0
  4 10  5  4  7  0  0  1  0  0
  9  4  5  3  3  0  0  0  1  0
  6  3  3  3  7  0  0  0  0  1

show 2 1
    10     3     8     6     3     1     0     0     0     0
     0    11     8    10     1    -1     2     0     0     0
     0    22   9/2     4  29/2    -1     0   5/2     0     0
     0  13/9 -22/9  -8/3   1/3    -1     0     0  10/9     0
     0     2    -3    -1  26/3    -1     0     0     0   5/3

show 2 2
      10       3       8       6       3       1       0       0       0       0
       0      11       8      10       1      -1       2       0       0       0
       0       0   -23/4      -8    25/4     1/2      -2     5/4       0       0
       0       0 -346/13 -394/13   20/13  -86/13      -2       0  110/13       0
       0       0   -49/2   -31/2   140/3    -9/2      -2       0       0    55/6

show 2 3
        10         3         8         6         3         1         0         0         0         0
         0        11         8        10         1        -1         2         0         0         0
         0         0     -23/4        -8      25/4       1/2        -2       5/4         0         0
         0         0         0  1005/692 -4095/692 -1335/692  1085/692      -5/4  1265/692         0
         0         0         0   855/196    395/84  -305/196     75/49      -5/4         0  1265/588

show 2 4
           10            3            8            6            3            1            0            0            0            0
            0           11            8           10            1           -1            2            0            0            0
            0            0        -23/4           -8         25/4          1/2           -2          5/4            0            0
            0            0            0     1005/692    -4095/692    -1335/692     1085/692         -5/4     1265/692            0
            0            0            0            0 221375/29583   13915/9861 -13915/13148  16445/19722    -1265/692 84755/118332

Lower part has all zeros

show 3
           10            3            8            6            3            1            0            0            0            0
            0           11            8           10            1           -1            2            0            0            0
            0            0         23/4            8        -25/4         -1/2            2         -5/4            0            0
            0            0            0     1005/692    -4095/692    -1335/692     1085/692         -5/4     1265/692            0
            0            0            0            0 221375/29583   13915/9861 -13915/13148  16445/19722    -1265/692 84755/118332

show 4 5
           10            3            8            6            0       76/175      297/700     -117/350      513/700     -201/700
            0           11            8           10            0     -208/175     1499/700      -39/350      171/700      -67/700
            0            0         23/4            8            0        19/28      125/112       -31/56     -171/112       67/112
            0            0            0     1005/692            0   -1407/1730  10117/13840   -4087/6920   5293/13840   7839/13840
            0            0            0            0 221375/29583   13915/9861 -13915/13148  16445/19722    -1265/692 84755/118332

show 4 4
           10            3            8            0            0      664/175    -1817/700      737/350     -593/700    -1839/700
            0           11            8            0            0      772/175   -6073/2100    4153/1050   -5017/2100    -2797/700
            0            0         23/4            0            0     3611/700  -24449/8400   11339/4200  -30521/8400   -7061/2800
            0            0            0     1005/692            0   -1407/1730  10117/13840   -4087/6920   5293/13840   7839/13840
            0            0            0            0 221375/29583   13915/9861 -13915/13148  16445/19722    -1265/692 84755/118332

show 4 3
           10            3            0            0            0     -592/175    3053/2100   -1733/1050    8837/2100      617/700
            0           11            0            0            0     -484/175    2431/2100     209/1050    5599/2100     -341/700
            0            0         23/4            0            0     3611/700  -24449/8400   11339/4200  -30521/8400   -7061/2800
            0            0            0     1005/692            0   -1407/1730  10117/13840   -4087/6920   5293/13840   7839/13840
            0            0            0            0 221375/29583   13915/9861 -13915/13148  16445/19722    -1265/692 84755/118332

show 4 2
           10            0            0            0            0       -92/35      239/210     -179/105      731/210        71/70
            0           11            0            0            0     -484/175    2431/2100     209/1050    5599/2100     -341/700
            0            0         23/4            0            0     3611/700  -24449/8400   11339/4200  -30521/8400   -7061/2800
            0            0            0     1005/692            0   -1407/1730  10117/13840   -4087/6920   5293/13840   7839/13840
            0            0            0            0 221375/29583   13915/9861 -13915/13148  16445/19722    -1265/692 84755/118332

Upper half has all zeros

show 5
          1          0          0          0          0    -46/175   239/2100  -179/1050   731/2100     71/700
          0          1          0          0          0    -44/175   221/2100    19/1050   509/2100    -31/700
          0          0          1          0          0    157/175 -1063/2100   493/1050 -1327/2100   -307/700
          0          0          0          1          0     -14/25    151/300    -61/150     79/300     39/100
          0          0          0          0          1     33/175    -99/700     39/350   -171/700     67/700

Main diagonal has all ones

show 6 The inverse matrix
    -46/175   239/2100  -179/1050   731/2100     71/700
    -44/175   221/2100    19/1050   509/2100    -31/700
    157/175 -1063/2100   493/1050 -1327/2100   -307/700
     -14/25    151/300    -61/150     79/300     39/100
     33/175    -99/700     39/350   -171/700     67/700

The product of input and inverse matrix
 1 0 0 0 0
 0 1 0 0 0
 0 0 1 0 0
 0 0 0 1 0
 0 0 0 0 1

version 2

Translation of: PL/I

<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 /*──────────────────────────────────────────────────────────────────────────────────────*/ 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</lang>
output   when using the internal default input:
The matrix of order  4  is:
  8  3  7  5
  9 12 10 11
  6  2  4 13
 14 15 16 17

The inverted matrix is:
   0.3590   0.5385   0.0769  -0.5128
  -0.0190   0.3419  -0.0199  -0.2004
  -0.1833  -0.7009  -0.1425   0.6163
  -0.1064  -0.0855   0.0883   0.0779

Sidef

Uses the rref(M) function from Reduced row echelon form.

Translation of: Raku

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

Output:
-21/23   17/69  13/138   19/46
-38/23   15/23    1/23   15/23
-16/23   25/69  11/138    9/46
-13/23   16/69   -2/69   13/23

VBA

Translation of: Phix

Uses ToReducedRowEchelonForm() from Reduced_row_echelon_form#VBA<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</lang>

Output:
 0,75          0,5           0,25         
 0,5           1             0,5          
 0,25          0,5           0,75

Wren

Library: Wren-fmt
Library: Wren-matrix

The above module does in fact use the Gauss-Jordan method for matrix inversion. <lang ecmascript>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()

}</lang>

Output:
Original:

| 1  2  3|
| 4  1  6|
| 7  8  9|

Inverse:

|-0.812500  0.125000  0.187500|
| 0.125000 -0.250000  0.125000|
| 0.520833  0.125000 -0.145833|

Original:

| 2 -1  0|
|-1  2 -1|
| 0 -1  2|

Inverse:

| 0.750000  0.500000  0.250000|
| 0.500000  1.000000  0.500000|
| 0.250000  0.500000  0.750000|

Original:

|-1 -2  3  2|
|-4 -1  6  2|
| 7 -8  9  1|
| 1 -2  1  3|

Inverse:

|-0.913043  0.246377  0.094203  0.413043|
|-1.652174  0.652174  0.043478  0.652174|
|-0.695652  0.362319  0.079710  0.195652|
|-0.565217  0.231884 -0.028986  0.565217|

zkl

This uses GSL to invert a matrix via LU decomposition, not Gauss-Jordan. <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();</lang>

Output:
   -0.8125,    0.1250,    0.1875
    0.1250,   -0.2500,    0.1250
    0.5208,    0.1250,   -0.1458

    1.0000,    0.0000,    0.0000
   -0.0000,    1.0000,    0.0000
   -0.0000,    0.0000,    1.0000

<lang zkl>m:=GSL.Matrix(3,3).set(2,-1,0, -1,2,-1, 0,-1,2); m.invert().format(10,4).println("\n");</lang>

Output:
    0.7500,    0.5000,    0.2500
    0.5000,    1.0000,    0.5000
    0.2500,    0.5000,    0.7500
Retrieved from "https://rosettacode.org/wiki/Gauss-Jordan_matrix_inversion?oldid=279825"