Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Invert a matrix with Crystal via Lapack bindings
@[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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment