Gauss-Jordan matrix inversion

From Rosetta Code
Revision as of 23:50, 12 March 2018 by rosettacode>Gpapo (Added Julia language)

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

Works with: Julia version 0.6

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

Native implementation: <lang julia>function gaussjordan(A::Matrix)

   size(A, 1) == size(A, 2) || throw(ArgumentError("A must be squared"))
   n = size(A, 1)
   M = [convert(Matrix{float(eltype(A))}, A) I]
   i = 1
   local tmp = Vector{eltype(M)}(2n)
   # forward
   while i ≤ n
       if M[i, i] ≈ 0.0
           local j = i + 1
           while j ≤ n && M[j, i] ≈ 0.0
               j += 1
           end
           if j ≤ n
               tmp     .= M[i, :]
               M[i, :] .= M[j, :]
               M[j, :] .= tmp
           else
               throw(ArgumentError("matrix is singular, cannot compute the inverse"))
           end
       end
       for j in (i + 1):n
           M[j, :] .-= M[j, i] / M[i, i] .* M[i, :]
       end
       i += 1
   end
   i = n
   # backward
   while i ≥ 1
       if M[i, i] ≈ 0.0
           local j = i - 1
           while j ≥ 1 && M[j, i] ≈ 0.0
               j -= 1
           end
           if j ≥ 1
               tmp     .= M[i, :]
               M[i, :] .= M[j, :]
               M[j, :] .= tmp
           else
               throw(ArgumentError("matrix is singular, cannot compute the inverse"))
           end
       end
       for j in (i - 1):-1:1
           M[j, :] .-= M[j, i] / M[i, i] .* M[i, :]
       end
       i -= 1
   end
   M ./= diag(M) # normalize
   return M[:, n+1:2n]

end

@show gaussjordan(A) @assert gaussjordan(A) ≈ inv(A)

A = rand(10, 10) @assert gaussjordan(A) ≈ inv(A)</lang>

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