Recursive linear algebra algorithms

[Patrick Alken <alken@colorado.edu>]

Hi all,


   I have recently learned of a project called ReLAPACK 
(https://github.com/HPAC/ReLAPACK, paper here: 
https://arxiv.org/abs/1602.06763) which implements a number of LAPACK 
algorithms (such as LU, Cholesky, Sylvester equations) using recursive 
methods which can use Level 3 BLAS calls. The paper shows that most of 
these algorithms out-perform the block Level 3 algorithms in LAPACK. The 
main advantage is that LAPACK block algorithms require fixing the block 
size ahead of time, which may not be optimal for a given architecture, 
while the recursive methods don't require a block size parameter.

The recursive methods do however require a "base case" - i.e. at what 
size matrix should it switch to the Level 2 BLAS based algorithms. 
ReLAPACK fixes this currently at N=24.

Anyway, the recursive Cholesky variant is quite straightforward to 
implement, and I have already coded it for GSL (both the decomposition 
and inversion). I did some tests for N=10,000 with ATLAS BLAS and found 
that it runs faster than DPOTRF from LAPACK. This fast Cholesky 
decomposition will improve the performance also for the generalized 
symmetric definite eigensolvers, and the least squares modules (linear 
and nonlinear).

So GSL now has a competitive Cholesky solver, which I think should make 
many GSL users happy :)

Work is currently underway to implement the recursive pivoted LU 
decomposition in GSL.

Unfortunately the ReLAPACK authors state that the QR algorithm is not 
amenable to recursive methods, so the block QR seems to still be the 
best choice. It would be nice to implement this for GSL, in case any 
volunteers are looking for a project ;)

Patrick