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 1 1
1 1 0
1 0 0]
@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.0 0.0 1.0; 0.0 1.0 -1.0; 1.0 -1.0 0.0] inv(A) = [0.0 0.0 1.0; 0.0 1.0 -1.0; 1.0 -1.0 0.0] gaussjordan(A) = [0.0 0.0 1.0; -0.0 1.0 -1.0; 1.0 -1.0 -0.0]