Gauss-Jordan matrix inversion
Invert matrix A using Gauss-Jordan method.
You are encouraged to solve this task according to the task description, using any language you may know.
- Task
A being an n × n matrix.
360 Assembly
The COPY file of FORMAT, to format a floating point number, can be found at: Include files 360 Assembly. <lang 360asm>* Gauss-Jordan matrix inversion 17/01/2021 GAUSSJOR CSECT
USING GAUSSJOR,R13 base register B 72(R15) skip savearea DC 17F'0' savearea SAVE (14,12) save previous context ST R13,4(R15) link backward ST R15,8(R13) link forward LR R13,R15 set addressability LA R6,1 i=1 DO WHILE=(C,R6,LE,N) do i=1 to n LA R7,1 j=1 DO WHILE=(C,R7,LE,N) do j=1 to n LR R1,R6 i BCTR R1,0 -1 MH R1,HN *n LR R0,R7 j BCTR R0,0 -1 AR R1,R0 (i-1)*n+j-1 SLA R1,2 *4 LE F0,AA(R1) aa(i,j) LR R1,R6 i BCTR R1,0 -1 MH R1,HN2 *n*2 LR R0,R7 j BCTR R0,0 -1 AR R1,R0 (i-1)*n*2+j-1 SLA R1,2 *4 STE F0,TMP(R1) tmp(i,j)=aa(i,j) LA R7,1(R7) j++ ENDDO , enddo j LA R7,1 j=1 A R7,N j=n+1 DO WHILE=(C,R7,LE,N2) do j=n+1 to 2*n LR R1,R6 i BCTR R1,0 -1 MH R1,HN2 *n*2 LR R0,R7 j BCTR R0,0 -1 AR R1,R0 (i-1)*n*2+j-1 SLA R1,2 *4 LE F0,=E'0' 0 STE F0,TMP(R1) tmp(i,j)=0 LA R7,1(R7) j++ ENDDO , enddo j LR R2,R6 i A R2,N i+n LR R1,R6 i BCTR R1,0 -1 MH R1,HN2 *n*2 BCTR R2,0 -1 AR R1,R2 (i+n-1)*n*2+i-1 SLA R1,2 *4 LE F0,=E'1' 1 STE F0,TMP(R1) tmp(i,i+n)=1 LA R6,1(R6) i++ ENDDO , enddo i LA R6,1 i=1 DO WHILE=(C,R6,LE,N) do r=1 to n LR R1,R6 r BCTR R1,0 -1 MH R1,HN2 *n*2 LR R0,R6 r BCTR R0,0 -1 AR R1,R0 (r-1)*n*2+r-1 SLA R1,2 *4 LE F0,TMP(R1) tmp(r,r) STE F0,T t=tmp(r,r) LA R7,1 c=1 DO WHILE=(C,R7,LE,N2) do c=1 to n*2 LR R1,R6 r BCTR R1,0 MH R1,HN2 LR R0,R7 c BCTR R0,0 AR R1,R0 SLA R1,2 LE F0,TMP(R1) LTER F0,F0 BZ *+8 DE F0,T tmp(r,c)/t LR R1,R6 r BCTR R1,0 MH R1,HN2 LR R0,R7 c BCTR R0,0 AR R1,R0 SLA R1,2 STE F0,TMP(R1) tmp(r,c)=tmp(r,c)/t LA R7,1(R7) c++ ENDDO , enddo c LA R8,1 i=1 DO WHILE=(C,R8,LE,N) do i=1 to n IF CR,R8,NE,R6 THEN if i^=r then LR R1,R8 i BCTR R1,0 MH R1,HN2 LR R0,R6 r BCTR R0,0 AR R1,R0 SLA R1,2 LE F0,TMP(R1) tmp(i,r) STE F0,T t=tmp(i,r) LA R7,1 c=1 DO WHILE=(C,R7,LE,N2) do c=1 to n*2 LR R2,R8 i BCTR R2,0 MH R2,HN2 LR R0,R7 c BCTR R0,0 AR R2,R0 SLA R2,2 LE F0,TMP(R2) f0=tmp(i,c) LR R1,R6 r BCTR R1,0 MH R1,HN2 LR R0,R7 c BCTR R0,0 AR R1,R0 SLA R1,2 LE F2,TMP(R1) f2=tmp(r,c) LE F4,T t MER F4,F2 f4=t*tmp(r,c) SER F0,F4 f0=tmp(i,c)-t*tmp(r,c) STE F0,TMP(R2) tmp(i,c)=f0 LA R7,1(R7) c++ ENDDO , enddo c ENDIF , endif LA R8,1(R8) i++ ENDDO , enddo i LA R6,1(R6) r++ ENDDO , enddo r LA R6,1 i=1 DO WHILE=(C,R6,LE,N) do i=1 to n L R7,N n LA R7,1(R7) j=n+1 DO WHILE=(C,R7,LE,N2) do j=n+1 to 2*n LR R2,R7 j S R2,N -n LR R1,R6 i BCTR R1,0 MH R1,HN2 *2*n LR R0,R7 j BCTR R0,0 AR R1,R0 SLA R1,2 LE F0,TMP(R1) tmp(i,j) LR R1,R6 i BCTR R1,0 MH R1,HN *n BCTR R2,0 AR R1,R2 SLA R1,2 STE F0,BB(R1) bb(i,j-n)=tmp(i,j) LA R7,1(R7) j++ ENDDO , enddo j LA R6,1(R6) i++ ENDDO , enddo i LA R6,1 i=1 DO WHILE=(C,R6,LE,N) do i=1 to n LA R3,PG @pg LA R7,1 j=1 DO WHILE=(C,R7,LE,N) do j=1 to n LR R1,R6 i BCTR R1,0 MH R1,HN *n LR R0,R7 j BCTR R0,0 AR R1,R0 SLA R1,2 LE F0,BB(R1) bb(i,j) LA R0,5 number of decimals BAL R14,FORMATF FormatF(bb(i,j)) MVC 0(13,R3),0(R1) edit LA R3,13(R3) @pg+=13 LA R7,1(R7) j++ ENDDO , enddo j XPRNT PG,L'PG print buffer LA R6,1(R6) i++ ENDDO , enddo i L R13,4(0,R13) restore previous savearea pointer RETURN (14,12),RC=0 restore registers from calling save COPY plig\$_FORMATF.MLC format a float
NN EQU 5 N DC A(NN) N2 DC A(NN*2) AA DC E'3.0',E'1.0',E'5.0',E'9.0',E'7.0'
DC E'8.0',E'2.0',E'8.0',E'0.0',E'1.0' DC E'1.0',E'7.0',E'2.0',E'0.0',E'3.0' DC E'0.0',E'1.0',E'7.0',E'0.0',E'9.0' DC E'3.0',E'5.0',E'6.0',E'1.0',E'1.0'
BB DS (NN*NN)E TMP DS (NN*NN*2)E T DS E HN DC AL2(NN) HN2 DC AL2(NN*2) PG DC CL80' ' buffer
REGEQU END GAUSSJOR </lang>
- Output:
0.02863 0.19429 0.13303 -0.05956 -0.25778 -0.00580 -0.03255 0.12109 -0.03803 0.05226 -0.03020 -0.06824 -0.17903 0.06054 0.27187 0.10021 -0.06733 -0.05617 -0.06263 0.09806 0.02413 0.05669 0.12579 0.06824 -0.21726
ALGOL 60
<lang algol60>begin
comment Gauss-Jordan matrix inversion - 22/01/2021; integer n; n:=4; begin procedure rref(m); real array m; begin integer r, c, i; real d, w; for r := 1 step 1 until n do begin d := m[r,r]; if d notequal 0 then for c := 1 step 1 until n*2 do m[r,c] := m[r,c] / d else outstring(1,"inversion impossible!\n"); for i := 1 step 1 until n do if i notequal r then begin w := m[i,r]; for c := 1 step 1 until n*2 do m[i,c] := m[i,c] - w * m[r,c] end end end rref;
procedure inverse(mat,inv); real array mat, inv; begin integer i, j; real array aug[1:n,1:n*2]; for i := 1 step 1 until n do begin for j := 1 step 1 until n do aug[i,j] := mat[i,j]; for j := 1+n step 1 until n*2 do aug[i,j] := 0; aug[i,i+n] := 1 end; rref(aug); for i := 1 step 1 until n do for j := n+1 step 1 until 2*n do inv[i,j-n] := aug[i,j] end inverse;
procedure show(c, m); string c; real array m; begin integer i, j; outstring(1,"matrix "); outstring(1,c); outstring(1,"\n"); for i := 1 step 1 until n do begin for j := 1 step 1 until n do outreal(1,m[i,j]); outstring(1,"\n") end end show;
integer i; real array a[1:n,1:n], b[1:n,1:n], c[1:n,1:n]; i:=1; a[i,1]:= 2; a[i,2]:= 1; a[i,3]:= 1; a[i,4]:= 4; i:=2; a[i,1]:= 0; a[i,2]:=-1; a[i,3]:= 0; a[i,4]:=-1; i:=3; a[i,1]:= 1; a[i,2]:= 0; a[i,3]:=-2; a[i,4]:= 4; i:=4; a[i,1]:= 4; a[i,2]:= 1; a[i,3]:= 2; a[i,4]:= 2; show("a",a); inverse(a,b); show("b",b); inverse(b,c); show("c",c) end
end </lang>
- Output:
matrix a 2 1 1 4 0 -1 0 -1 1 0 -2 4 4 1 2 2 matrix b -0.4 0 0.2 0.4 -0.4 -1.2 0 0.2 0.6 0.4 -0.4 -0.2 0.4 0.2 0 -0.2 matrix c 2 1 1 4 0 -1 0 -1 1 0 -2 4 4 1 2 2
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
C++
<lang cpp>#include <algorithm>
- include <cassert>
- include <iomanip>
- include <iostream>
- include <vector>
template <typename scalar_type> class matrix { public:
matrix(size_t rows, size_t columns) : rows_(rows), columns_(columns), elements_(rows * columns) {}
matrix(size_t rows, size_t columns, scalar_type value) : rows_(rows), columns_(columns), elements_(rows * columns, value) {}
matrix(size_t rows, size_t columns, std::initializer_list<scalar_type> values) : rows_(rows), columns_(columns), elements_(rows * columns) { assert(values.size() <= rows_ * columns_); std::copy(values.begin(), values.end(), elements_.begin()); }
size_t rows() const { return rows_; } size_t columns() const { return columns_; }
scalar_type& operator()(size_t row, size_t column) { assert(row < rows_); assert(column < columns_); return elements_[row * columns_ + column]; } const scalar_type& operator()(size_t row, size_t column) const { assert(row < rows_); assert(column < columns_); return elements_[row * columns_ + column]; }
private:
size_t rows_; size_t columns_; std::vector<scalar_type> elements_;
};
template <typename scalar_type> matrix<scalar_type> product(const matrix<scalar_type>& a,
const matrix<scalar_type>& b) { assert(a.columns() == b.rows()); size_t arows = a.rows(); size_t bcolumns = b.columns(); size_t n = a.columns(); matrix<scalar_type> c(arows, bcolumns); for (size_t i = 0; i < arows; ++i) { for (size_t j = 0; j < n; ++j) { for (size_t k = 0; k < bcolumns; ++k) c(i, k) += a(i, j) * b(j, k); } } return c;
}
template <typename scalar_type> void swap_rows(matrix<scalar_type>& m, size_t i, size_t j) {
size_t columns = m.columns(); for (size_t column = 0; column < columns; ++column) std::swap(m(i, column), m(j, column));
}
// Convert matrix to reduced row echelon form template <typename scalar_type> void rref(matrix<scalar_type>& m) {
size_t rows = m.rows(); size_t columns = m.columns(); for (size_t row = 0, lead = 0; row < rows && lead < columns; ++row, ++lead) { size_t i = row; while (m(i, lead) == 0) { if (++i == rows) { i = row; if (++lead == columns) return; } } swap_rows(m, i, row); if (m(row, lead) != 0) { scalar_type f = m(row, lead); for (size_t column = 0; column < columns; ++column) m(row, column) /= f; } for (size_t j = 0; j < rows; ++j) { if (j == row) continue; scalar_type f = m(j, lead); for (size_t column = 0; column < columns; ++column) m(j, column) -= f * m(row, column); } }
}
template <typename scalar_type> matrix<scalar_type> inverse(const matrix<scalar_type>& m) {
assert(m.rows() == m.columns()); size_t rows = m.rows(); matrix<scalar_type> tmp(rows, 2 * rows, 0); for (size_t row = 0; row < rows; ++row) { for (size_t column = 0; column < rows; ++column) tmp(row, column) = m(row, column); tmp(row, row + rows) = 1; } rref(tmp); matrix<scalar_type> inv(rows, rows); for (size_t row = 0; row < rows; ++row) { for (size_t column = 0; column < rows; ++column) inv(row, column) = tmp(row, column + rows); } return inv;
}
template <typename scalar_type> void print(std::ostream& out, const matrix<scalar_type>& m) {
size_t rows = m.rows(), columns = m.columns(); out << std::fixed << std::setprecision(4); for (size_t row = 0; row < rows; ++row) { for (size_t column = 0; column < columns; ++column) { if (column > 0) out << ' '; out << std::setw(7) << m(row, column); } out << '\n'; }
}
int main() {
matrix<double> m(3, 3, {2, -1, 0, -1, 2, -1, 0, -1, 2}); std::cout << "Matrix:\n"; print(std::cout, m); auto inv(inverse(m)); std::cout << "Inverse:\n"; print(std::cout, inv); std::cout << "Product of matrix and inverse:\n"; print(std::cout, product(m, inv)); std::cout << "Inverse of inverse:\n"; print(std::cout, inverse(inv)); return 0;
}</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
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.
<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"
- 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
<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
<lang go>package main
import (
"fmt" "math"
)
type vector = []float64 type matrix []vector
func (m matrix) print(title string) {
fmt.Println(title) for _, v := range m { fmt.Printf("% f\n", v) }
}
func I(n int) matrix {
b := make(matrix, n) for i := 0; i < n; i++ { b[i] = make(vector, n) b[i][i] = 1 } return b
}
func (m matrix) inverse() matrix {
n := len(m) for _, v := range m { if n != len(v) { panic("Not a square matrix") } } b := I(n) a := make(matrix, n) for i, v := range m { a[i] = make(vector, n) copy(a[i], v) } for k := range a { iMax := 0 max := -1. for i := k; i < n; i++ { row := a[i] // compute scale factor s = max abs in row s := -1. for j := k; j < n; j++ { x := math.Abs(row[j]) if x > s { s = x } } if s == 0 { panic("Irregular matrix") } // scale the abs used to pick the pivot. if abs := math.Abs(row[k]) / s; abs > max { iMax = i max = abs } } if k != iMax { a[k], a[iMax] = a[iMax], a[k] b[k], b[iMax] = b[iMax], b[k] } akk := a[k][k] for j := 0; j < n; j++ { a[k][j] /= akk b[k][j] /= akk } for i := 0; i < n; i++ { if (i != k) { aik := a[i][k] for j := 0; j < n; j++ { a[i][j] -= a[k][j] * aik b[i][j] -= b[k][j] * aik } } } } return b
}
func main() {
a := matrix{{1, 2, 3}, {4, 1, 6}, {7, 8, 9}} a.inverse().print("Inverse of A is:\n")
b := matrix{{2, -1, 0}, {-1, 2, -1}, {0, -1, 2}} b.inverse().print("Inverse of B is:\n")
}</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.printf("%7.3f", elements[row][column]); } System.out.println(); } }
// Returns the inverse of this matrix public Matrix inverse() { assert(rows == columns); // Augment by identity matrix Matrix tmp = new Matrix(rows, columns * 2); for (int row = 0; row < rows; ++row) { System.arraycopy(elements[row], 0, tmp.elements[row], 0, columns); tmp.elements[row][row + columns] = 1; } tmp.toReducedRowEchelonForm(); Matrix inv = new Matrix(rows, columns); for (int row = 0; row < rows; ++row) System.arraycopy(tmp.elements[row], columns, inv.elements[row], 0, columns); return inv; }
// Converts this matrix into reduced row echelon form public void toReducedRowEchelonForm() { for (int row = 0, lead = 0; row < rows && lead < columns; ++row, ++lead) { int i = row; while (elements[i][lead] == 0) { if (++i == rows) { i = row; if (++lead == columns) return; } } swapRows(i, row); if (elements[row][lead] != 0) { double f = elements[row][lead]; for (int column = 0; column < columns; ++column) elements[row][column] /= f; } for (int j = 0; j < rows; ++j) { if (j == row) continue; double f = elements[j][lead]; for (int column = 0; column < columns; ++column) elements[j][column] -= f * elements[row][column]; } } }
// Returns the matrix product of a and b public static Matrix product(Matrix a, Matrix b) { assert(a.columns == b.rows); Matrix result = new Matrix(a.rows, b.columns); for (int i = 0; i < a.rows; ++i) { double[] resultRow = result.elements[i]; double[] aRow = a.elements[i]; for (int j = 0; j < a.columns; ++j) { double[] bRow = b.elements[j]; for (int k = 0; k < b.columns; ++k) resultRow[k] += aRow[j] * bRow[k]; } } return result; }
private void swapRows(int i, int j) { double[] tmp = elements[i]; elements[i] = elements[j]; elements[j] = tmp; }
}</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
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>
Mathematica / Wolfram Language
<lang Mathematica>a = N@{{1, 7, 5, 0}, {5, 8, 6, 9}, {2, 1, 6, 4}, {8, 1, 2, 4}}; elimmat = RowReduce[Join[a, IdentityMatrix[Length[a]], 2]]; inv = elimmat[[All, -Length[a] ;;]]</lang>
- Output:
{{0.0463243, -0.0563948, -0.0367573, 0.163646}, {0.0866062, 0.0684794, -0.133938, -0.020141}, {0.0694864, -0.0845921, 0.194864, -0.00453172}, {-0.149043, 0.137966, 0.00956697, -0.0699899}}
MATLAB
<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
- ———————————————————————————————————————————————————————————————————————————————————————————————————
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
uses ToReducedRowEchelonForm() from Reduced_row_echelon_form#Phix
with javascript_semantics function ToReducedRowEchelonForm(sequence M) integer lead = 1, rowCount = length(M), columnCount = length(M[1]), i for r=1 to rowCount do if lead>=columnCount then exit end if i = r while M[i][lead]=0 do i += 1 if i=rowCount then i = r lead += 1 if lead=columnCount then exit end if end if end while object mr = sq_div(M[i],M[i][lead]) M[i] = M[r] M[r] = mr for j=1 to rowCount do if j!=r then M[j] = sq_sub(M[j],sq_mul(M[j][lead],M[r])) end if end for lead += 1 end for return M end function --</end of copy> function inverse(sequence mat) integer len = length(mat) sequence aug = repeat(repeat(0,2*len),len) for i=1 to len do if length(mat[i])!=len then ?9/0 end if -- "Not a square matrix" for j=1 to len do aug[i][j] = mat[i][j] end for -- augment by identity matrix to right aug[i][i + len] = 1 end for aug = ToReducedRowEchelonForm(aug) sequence inv = repeat(repeat(0,len),len) -- remove identity matrix to left for i=1 to len do for j=len+1 to 2*len do inv[i][j-len] = aug[i][j] end for end for return inv end function constant test = {{ 2, -1, 0}, {-1, 2, -1}, { 0, -1, 2}} pp(inverse(test),{pp_Nest,1})
- Output:
{{0.75,0.5,0.25}, {0.5,1,0.5}, {0.25,0.5,0.75}}
alternate
with javascript_semantics function gauss_jordan_inv(sequence a) integer n = length(a) sequence b = repeat(repeat(0,n),n) for i=1 to n do b[i,i] = 1 end for for k=1 to n do integer lmax = k atom amax = abs(a[k][k]) for l=k+1 to n do atom tmp = abs(a[l][k]) if amax<tmp then {amax,lmax} = {tmp,l} end if end for if k!=lmax then {a[k], a[lmax]} = {a[lmax], a[k]} {b[k], b[lmax]} = {b[lmax], b[k]} end if atom akk = a[k][k] if akk=0 then throw("Irregular matrix") end if for j=1 to n do a[k][j] /= akk b[k][j] /= akk end for for i=1 to n do if i!=k then atom aik = a[i][k] for j=1 to n do a[i][j] -= a[k][j]*aik b[i][j] -= b[k][j]*aik end for end if end for end for return b end function sequence a = {{ 2, -1, 0}, {-1, 2, -1}, { 0, -1, 2}}, inva = gauss_jordan_inv(deep_copy(a)) puts(1,"a = ") pp(a,{pp_Nest,1,pp_Indent,4,pp_IntFmt,"%2d"}) puts(1,"inv(a) = ") pp(inva,{pp_Nest,1,pp_Indent,9,pp_FltFmt,"%4.2f"})
- Output:
a = {{ 2,-1, 0}, {-1, 2,-1}, { 0,-1, 2}} inv(a) = {{0.75,0.50,0.25}, {0.50,1.00,0.50}, {0.25,0.50,0.75}}
PL/I
/* Gauss-Jordan matrix inversion */ G_J: procedure options (main); /* 4 November 2020 */ declare t float; declare (i, j, k, n) fixed binary; open file (sysin) title ('/GAUSSJOR.DAT'); get (n); /* Read in the order of the matrix. */ put skip data (n); begin; declare a(n,n)float, aux(n,n) float; aux = 0; do i = 1 to n; aux(i,i) = 1; end; get (a); /* Read in the matrix. */ put skip list ('The matrix to be inverted is:'); put edit (a) ( skip, (n) F(10,4)); /* Print the matrix. */ do k = 1 to n; /* Divide row k by a(k,k) */ t = a(k,k); a(k,*) = a(k,*) / t; aux(k,*) = aux(k,*) / t; do i = 1 to k-1, k+1 to n; /* Work down the rows. */ t = a(i,k); a(i,*) = a(i,*) - t*a(k,*); aux(i,*) = aux(i,*) - t*aux(k,*); end; end; put skip (2) list ('The inverse is:'); put edit (aux) ( skip, (n) F(10,4)); end; /* of BEGIN block */ end G_J;
Output:
N= 4; The matrix to be inverted is: 8.0000 3.0000 7.0000 5.0000 9.0000 12.0000 10.0000 11.0000 6.0000 2.0000 4.0000 13.0000 14.0000 15.0000 16.0000 17.0000 The inverse is: 0.3590 0.5385 0.0769 -0.5128 -0.0190 0.3419 -0.0199 -0.2004 -0.1833 -0.7009 -0.1425 0.6163 -0.1064 -0.0855 0.0883 0.0779
PowerShell
<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)
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 }
- 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
<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
Rust
<lang rust> fn main() {
let mut a: Vec<Vec<f64>> = vec![vec![1.0, 2.0, 3.0], vec![4.0, 1.0, 6.0], vec![7.0, 8.0, 9.0] ]; let mut b: Vec<Vec<f64>> = vec![vec![2.0, -1.0, 0.0], vec![-1.0, 2.0, -1.0], vec![0.0, -1.0, 2.0] ];
let mut ref_a = &mut a; let rref_a = &mut ref_a; let mut ref_b = &mut b; let rref_b = &mut ref_b;
println!("Matrix A:\n"); print_matrix(rref_a); println!("\nInverse of Matrix A:\n"); print_matrix(&mut matrix_inverse(rref_a)); println!("\n\nMatrix B:\n"); print_matrix(rref_b); println!("\nInverse of Matrix B:\n"); print_matrix(&mut matrix_inverse(rref_b));
}
//Begin Matrix Inversion fn matrix_inverse(matrix: &mut Vec<Vec<f64>>) -> Vec<Vec<f64>>{
let len = matrix.len(); let mut aug = zero_matrix(len, len * 2); for i in 0..len { for j in 0.. len { aug[i][j] = matrix[i][j]; } aug[i][i + len] = 1.0; }
gauss_jordan_general(&mut aug); let mut unaug = zero_matrix(len, len); for i in 0..len { for j in 0..len { unaug[i][j] = aug[i][j+len]; } } unaug
} //End Matrix Inversion
//Begin Generalised Reduced Row Echelon Form fn gauss_jordan_general(matrix: &mut Vec<Vec<f64>>) {
let mut lead = 0; let row_count = matrix.len(); let col_count = matrix[0].len();
for r in 0..row_count { if col_count <= lead { break; } let mut i = r; while matrix[i][lead] == 0.0 { i = i + 1; if row_count == i { i = r; lead = lead + 1; if col_count == lead { break; } } }
let temp = matrix[i].to_owned(); matrix[i] = matrix[r].to_owned(); matrix[r] = temp.to_owned();
if matrix[r][lead] != 0.0 { let div = matrix[r][lead]; for j in 0..col_count { matrix[r][j] = matrix[r][j] / div; } }
for k in 0..row_count { if k != r { let mult = matrix[k][lead]; for j in 0..col_count { matrix[k][j] = matrix[k][j] - matrix[r][j] * mult; } } } lead = lead + 1;
} //matrix.to_owned()
}
fn zero_matrix(rows: usize, cols: usize) -> Vec<Vec<f64>> {
let mut matrix = Vec::with_capacity(cols); for _ in 0..rows { let mut col: Vec<f64> = Vec::with_capacity(rows); for _ in 0..cols { col.push(0.0); } matrix.push(col); } matrix
}
fn print_matrix(mat: &mut Vec<Vec<f64>>) {
for row in 0..mat.len(){ println!("{:?}", mat[row]); }
} </lang>
- Output:
Matrix A: [1.0, 2.0, 3.0] [4.0, 1.0, 6.0] [7.0, 8.0, 9.0] Inverse of Matrix A: [-0.8125, 0.125, 0.1875] [0.12499999999999994, -0.24999999999999997, 0.12499999999999997] [0.5208333333333334, 0.12499999999999999, -0.14583333333333331] Matrix B: [2.0, -1.0, 0.0] [-1.0, 2.0, -1.0] [0.0, -1.0, 2.0] Inverse of Matrix B: [0.75, 0.49999999999999994, 0.24999999999999994] [0.49999999999999994, 0.9999999999999999, 0.4999999999999999] [0.24999999999999997, 0.49999999999999994, 0.7499999999999999]
Sidef
Uses the rref(M) function from Reduced row echelon form.
<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
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
VBScript
To run in console mode with cscript. <lang vb>' Gauss-Jordan matrix inversion - VBScript - 22/01/2021 Option Explicit
Function rref(m)
Dim r, c, i, n, div, wrk n=UBound(m) For r = 1 To n 'row div = m(r,r) If div <> 0 Then For c = 1 To n*2 'col m(r,c) = m(r,c) / div Next 'c Else WScript.Echo "inversion impossible!" WScript.Quit End If For i = 1 To n 'row If i <> r Then wrk = m(i,r) For c = 1 To n*2 m(i,c) = m(i,c) - wrk * m(r,c) Next' c End If Next 'i Next 'r rref = m
End Function 'rref
Function inverse(mat)
Dim i, j, aug, inv, n n = UBound(mat) ReDim inv(n,n), aug(n,2*n) For i = 1 To n For j = 1 To n aug(i,j) = mat(i,j) Next 'j aug(i,i+n) = 1 Next 'i aug = rref(aug) For i = 1 To n For j = n+1 To 2*n inv(i,j-n) = aug(i,j) Next 'j Next 'i inverse = inv
End Function 'inverse
Sub wload(m)
Dim i, j, k k = -1 For i = 1 To n For j = 1 To n k = k + 1 m(i,j) = w(k) Next 'j Next 'i
End Sub 'wload
Sub show(c, m, t)
Dim i, j, buf buf = "Matrix " & c &"=" & vbCrlf & vbCrlf For i = 1 To n For j = 1 To n If t="fixed" Then buf = buf & FormatNumber(m(i,j),6,,,0) & " " Else buf = buf & m(i,j) & " " End If Next 'j buf=buf & vbCrlf Next 'i WScript.Echo buf
End Sub 'show
Dim n, a, b, c, w w = Array( _ 2, 1, 1, 4, _ 0, -1, 0, -1, _ 1, 0, -2, 4, _ 4, 1, 2, 2) n=Sqr(UBound(w)+1) ReDim a(n,n), b(n,n), c(n,n) wload a show "a", a, "simple" b = inverse(a) show "b", b, "fixed" c = inverse(b) show "c", c, "fixed" </lang>
- Output:
Matrix a= 2 1 1 4 0 -1 0 -1 1 0 -2 4 4 1 2 2 Matrix b= -0,400000 0,000000 0,200000 0,400000 -0,400000 -1,200000 0,000000 0,200000 0,600000 0,400000 -0,400000 -0,200000 0,400000 0,200000 -0,000000 -0,200000 Matrix c= 2,000000 1,000000 1,000000 4,000000 0,000000 -1,000000 0,000000 -1,000000 1,000000 0,000000 -2,000000 4,000000 4,000000 1,000000 2,000000 2,000000
Wren
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