General Linear Algebra Wrapper Library

Question

I am currently mulling over the idea of taking a code I currently work with and rebuilding it from the ground up to allow for the use of more efficient programming and numerical techniques.

In the early stages of the design I decided I would like to write a wrapper class for matrices and vectors. This would allow me to (hopefully) switch out lin. alg. libraries and solvers without actually touching the code. I would simply write a class that has virtual functions for BLAS functionality along with some methods to internalize the initialization (allowing the main code to not worry about initializing in parallel). Then, when I want to make my code work with a new library, I simply write a new matrix class that implements all those functions using the efficient library functions.

I would be working with a finite-volume code, so I would be using sparse matrices with a preknown block memory layout. I am also looking at approximately a million unknowns although I would like this to not matter (so long as enough memory was available). The matrix wrapper would also get instructions from a config file to know which library, memory layout, and solvers to use, so it would not have to do any optimization on its own. It would only check if the requested library was installed.

To me this seems like a fairly straight forward concept, even if the implementation won't be easy. As such, my question comes in a couple parts.

First off, has this been tried/implemented before, and if so where? I would much rather be able to grab an open source library that start from nothing. Any reference or papers would be greatly appreciated.

And secondly, if I implement this sort of system, what sort of issues should I expect? I can already foresee that making it work transparently with multiple memory models (shared memory, GPU, or distributed memory) will be difficult, but I don't know what other pitfalls to look out for in designing a fairly general linear algebra interface.

Answer 1

deal.II has such wrappers for our own linear algebra implementation, for PETSc and for Trilinos. It's a lot of work to get all of this right because every library provides slightly different abstractions. As far as performance is concerned, deal.II does not use virtual functions since the overhead of calling a virtual function just to read a single element of a matrix or vector is likely non-negligible.

My recommendation would be to not go down this road. It's simply too much work. Rather, my suggestion would be to carefully choose which library to use and then stick with it. I don't know what you plan on doing with your code, but in all likelihood you're going to want to produce results, not spend forever on implementing interfaces. I say that because among the big two (PETSc and Trilinos), there isn't so much to be gained by making things generic -- they're both very good, and the performance differences are much smaller than what you can gain by changing your algorithms. It seems to me that you can only lose by spending a lot of time making everything generic.

Answer 2

Edit: Most of this answer had dealt with dense matrices, I've edited it more for sparse matrices, which is tougher. As far as I can tell, a lot of the issue with sparsity is taking advantage of caching and memory layout the way that you can with dense matrices. Storage and retrieval of sparse matrices entries is a lot more irregular.

I think there are a couple examples out there. There's Matlab and all its open source equivalents (as HighPerformanceMark mentioned), and Eigen is another library you might look at. I know Trilinos has some examples of this as well - Epetra takes an abstract sparse matrix for iterative and direct solvers. They deal with storage/retrieval issues in memory by having you specify an estimated number of nonzeros per row.

A thought on your note on multiple memory models - there seems to be something like BLAS/LAPACK for most architectures (Eigen, LAPACK/BLAS for serial, ViennaCL, CuBLAS, MAGMA, etc for GPU and OpenMP, Trilinos distributed, etc). It might be possible (though tricky?) to write a matrix class/solver which calls different packages depending on your architecture.