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 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
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.
<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>
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 <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
<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)
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/19REXX
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
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