Gauss-Jordan matrix inversion
Gauss-Jordan matrix inversion
Gauss-Jordan matrix inversion is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
- Task
Invert matrix A using Gauss-Jordan method.
A being an n by n matrix.
C#
<lang csharp> using System; using System.Collections.Generic; using System.Linq; using System.Text; using System.Threading.Tasks;
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 SwapRow(int r1, int r2)
{
if (r1 == r2) return;
double tmp = b[r1];
b[r1] = b[r2];
b[r2] = tmp;
}
internal double norm()
{
throw new NotImplementedException();
}
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]);
}
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 SwapRow(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;
}
}
internal Matrix Inv()
{
if (rows != cols) throw new Exception("rows != cols for Inv");
Matrix M = new Matrix(rows); //unitary
// do partial pivoting
for (int dia = 0; dia < cols; dia++)
{
int max_row = dia;
double max = this[dia, dia];
double tmp;
for (int row = dia + 1; row < rows; row++)
if ((tmp = Math.Abs(this[row, dia])) > max)
{
max_row = row;
max = tmp;
}
SwapRow(dia, max_row);
M.SwapRow(dia, max_row);
}
// change first matrix to diagonal; second matrix follows
for (int row = 0; row < rows; row++)
{
for (int r = 0; r < rows; r++)
if (r != row)
{
double ratio = this[r,row] / this[row,row];
for (int k = 0; k < cols; k++)
{
this[r, k] -= this[row, k] * ratio;
M[r, k] -= M[row, k] * ratio;
}
}
}
//reducing to unit matrix
for (int row = 0; row < rows; row++)
{
double ratio = this[row,row];
for (int col = 0; col < cols; col++)
{
this[row, col] = this[row, col] / ratio;
M[row, col] = M[row, col] / ratio;
}
}
return M;
}
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(3, 3, new double[] { 1, 2, 3, 4, 1, 6, 7, 8, 9 });
M.Inv().print();
}
}
} </lang>
- Output:
-0.8125 0.125 0.1875 0.125 -0.25 0.125 0.520833333333333 0.125 -0.145833333333333
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
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