Santiago Palladino

Santiago Palladino

Linear algebra in Crystal from LAPACK

2 min
Oct 30 2015
crystal, math
2 min
Oct 30 2015

This week Martín started a project for implementing linear algebra routines in Crystal by binding to LAPACK, a de-facto standard library for such operations. You can check out the project in github at mverzilli/crystalla, as in Crystal Linear Algebra.

The project makes use of Crystal’s simple C bindings to leverage LAPACK power. The first operation implemented for testing was computing the inverse of a matrix, which first computes de LU decomposition of the matrix and then calculates the inverse from there.

A simple benchmark shows that the implementation works quite fast with a 10000 x 10000 float random matrix, taking less than 40 seconds:

require "benchmark"
require "./src/crystalla"

cols = [] of Array(Float64)
range = 10000

r = Random.new
(1..range).each do |i|
  row = [] of Float64
  (1..range).each do |j|
    row.push r.next_float
  end
  cols.push row
end

m = Crystalla::Matrix.columns cols
m.invert!

We decided to contrast this against numpy, a widely-used python package for numeric operations, which also leverages LAPACK for linear algebra.

import numpy

SIZE = 10000
m = numpy.random.rand(SIZE, SIZE)
numpy.linalg.inv(m)

A time on this script yields one order of magnitude higher than Crystal:

real  3m38.331s
user  3m36.148s
sys   0m1.817s

We expect to continue adding basic functionality to crystalla over the following weeks. Feel free to contribute with your own additions!