Inverse of a Matrix using Gauss Jordan Elimination on Julia

Question

I want to make an algorithm that uses the Gauss-Jordan method to compute the inverse of a Matrix. I am more familiar with Matlab so I wrote this Matlab algorithm first. It works fine but I cannot make it to work on Julia

function [Inv,P]=GJinv(A)
  format rat
  n=length(A);
  P=eye(n); 
  Inv=[A,eye(n)]; 
  for k=1:n
        
        [~,r] = max(abs(A(k:end,k)));
        r = n-(n-k+1)+r;    
        #αλλάζω γραμμές#
        Inv([k r],:) = Inv([r k],:);
        P([k r],:) = P([r k],:);
       
        Inv(k,:) = Inv(k,:) / Inv(k,k);
        for i=1:n
          if i~=k
            Inv(i,:)=Inv(i,:)-(Inv(i,k)*Inv(k,:)); 
          endif
        endfor
  endfor
  Inv=Inv(1:n,n+1:2*n);
endfunction

My work so far:

using LinearAlgebra
function GJinv(A)
    (n,b)=size(A)
    Inv=[A,I(n)] 
    for k=1:n
          
          local ma=maximum(abs.(A[k:n,k:k])) 
          local pos = findfirst( x -> x == ma, abs.(A[k:n,k:n]))
          local r=pos[1] ; local r= n-(n-k+1)+r 
          #αλλάζω γραμμές#
          Inv[[k, r],:] = Inv[[r, k],:]
          
          Inv[k,:] = Inv[k,:] / Inv[k,k]
          for i=1:n
            if i!=k
              Inv[i,:]=Inv[i,:]-(Inv[i,k]*Inv[k,:]) 
            end
          end
    end
    Inv=Inv[1:n,n+1:2*n]
    print(Inv)
end
#TEST
A=[2 2 3; 4 5 6; 1 2 4]
GJinv(A)    

I am getting the following error:

ERROR: LoadError: InexactError: Int64(1.25)
Stacktrace:
  [1] Int64
    @ .\float.jl:912 [inlined]
  [2] convert
    @ .\number.jl:7 [inlined]
  [3] setindex!
    @ .\array.jl:1024 [inlined]
  [4] macro expansion
    @ .\multidimensional.jl:960 [inlined]
  [5] macro expansion
    @ .\cartesian.jl:64 [inlined]
  [6] _unsafe_setindex!(::IndexLinear, ::Matrix{Int64}, ::Vector{Float64}, ::Int64, ::Base.Slice{Base.OneTo{Int64}})
    @ Base .\multidimensional.jl:955
  [7] _setindex!
    @ .\multidimensional.jl:944 [inlined]
  [8] setindex!
    @ .\abstractarray.jl:1396 [inlined]
  [9] GJinv(A::Matrix{Int64})
    @\GJinv.jl:13
 [10] top-level scope
    @\GJinv.jl:24
in expression starting at @\GJinv.jl:24

Any help will be appreciated.

Answer 1

You were bitten by two problems of the syntax of MATLAB not being quite the same as Julia's syntax.

First, if you declare a matrix with all integer values in Julia the matrix will have an element type of Int, which will give you problems later when you divide, since in Julia the return type of

/(x::Int, y::Int)::Float64

is floating point

x / y

which you cannot then assign back to an Int matrix element. Second, as mentioned above by S Zakharov,

[A, I(n)]

will give you a Vector of Matrix in Julia, but you want the horizontal (hcat) matrix formed from two matrices.

So consider this correction (I fixed the syntax but did not change your algebra afaik):

using LinearAlgebra

function GJinv(squaremat)
    n, b = size(squaremat)
    @assert n == b # number of rows must == number of cols
    A = float.(squaremat)
    Inv = hcat(A, I(n))
    for k = 1:n
        ma = maximum(abs.(A[k:n,k:k])) 
        r = k + findfirst( x -> x == ma, abs.(A[k:n,k:n]))[1] - 1
        #αλλάζω γραμμές#
        Inv[[k, r],:] = Inv[[r, k],:]
        Inv[k,:] = Inv[k,:] / Inv[k,k]
        for i=1:n
            if i!=k
                Inv[i,:]=Inv[i,:]-(Inv[i,k]*Inv[k,:]) 
            end
        end
    end
    Inv = round.(Inv[1:n,n+1:2*n], digits = 6) # rounding makes it print better
    print(Inv)
end

 #TEST
 A=[2 2 3; 4 5 6; 1 2 4]

 GJinv(A)  # [1.6 -0.4 -0.6; -2.0 1.0 0.0; 0.6 -0.4 0.4]