Invert a matrix with Crystal via Lapack bindings

inverse_matrix.cr

@[Link(framework: "Accelerate")]
lib LibLapack
  fun sgetrf_(m: Int32*, n: Int32*, a: Void*, lda: Int32*, ipiv: Int32*, info: Int32*)

  fun sgetri_(n: Int32*, a: Void*, lda: Int32*, ipiv: Int32*, work: Void*, lwork: Int32*, info: Int32*)
end

# Matrix A (the one being inverted)
# Matrix [1 2
#         3 4]
# (matrices have to be provided to Lapack as 1D arrays)
# MV: Is there a more elegant way of initializing a StaticArray?
a = StaticArray(Float32, 4).new(0_f32)
a[0] = 1_f32
a[1] = 3_f32
a[2] = 2_f32
a[3] = 4_f32

n = 2
arows = n
acols = n

# Leading dimension of A
lda = 2

p "Matrix A before any calls:"

(0...n).each do |i|
  (0...n).each do |j|
    p "A[#{i}][#{j}]: #{a[n * j + i]}"
  end
end

# Set the dimension of the "workspace" array WORK
# see http://www.netlib.org/lapack-dev/lapack-coding/program-style.html#workspace
LWORK = 100
work = StaticArray(Float32, LWORK).new(0_f32)

# Pivot indices array
ipiv = StaticArray(Int32, 4).new(0)

# Lapack error codes
info = 0

# Calculate LU decomposition of A (and store it in A)
LibLapack.sgetrf_(pointerof(arows), pointerof(acols), pointerof(a) as Void*, pointerof(lda), ipiv, pointerof(info))

p "Matrix A after sgetrf"
(0...2).each do |i|
  (0...2).each do |j|
    p "A[#{i}][#{j}]: #{a[n*j+i]}"
  end
end

p "sgetrf_ returned an error!" if info != 0

lwork = LWORK

LibLapack.sgetri_(pointerof(n), pointerof(a) as Void*, pointerof(lda), ipiv, pointerof(work) as Void*, pointerof(lwork), pointerof(info))

p "Matrix A after sgetri"
(0...2).each do |i|
  (0...2).each do |j|
    p "A[#{i}][#{j}]: #{a[n*j+i]}"
  end
end

p "sgetri_ returned an error!" if info != 0