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.

  • $\begingroup$ Is this just for dense methods, or are you hoping to cover sparse/both? $\endgroup$ – Jesse Chan Sep 1 '13 at 4:15
  • 1
    $\begingroup$ A (slightly) tongue-in-cheek 'answer': yes, this has been done before, it's called Matlab. I suppose that Octave is the nearest open-source equivalent. $\endgroup$ – High Performance Mark Sep 1 '13 at 11:24
  • $\begingroup$ @JLC Preferably sparse. The code I want to use it for is a finite volume gas discharge code. $\endgroup$ – Godric Seer Sep 1 '13 at 13:06

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.

  • $\begingroup$ My idea about the initializer would hopefully avoid element by element access anywhere, but understand what you are saying. More than likely I will end up using one of the two libraries you mentioned. If I limit my general library to just static mesh finite volume solves, it isn't exactly general. $\endgroup$ – Godric Seer Sep 1 '13 at 16:13
  • $\begingroup$ Well, I think what I'm saying is that if your goal is to write a general and widely used finite element library, then being generic is a good approach. But most of us write codes that have a more limited ambition and there genericity is just a waste of time. $\endgroup$ – Wolfgang Bangerth Sep 1 '13 at 18:19
  • $\begingroup$ I am using the word library a bit loosely here. It is more a personal project that decouples the linear algebra from the physics in my own codes. $\endgroup$ – Godric Seer Sep 1 '13 at 18:36

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.

  • $\begingroup$ I would be happy being able to use 2-3 shared memory libraries and 1 good GPU library. If I limit myself to working with finite volume implementations (block diagonal plus known arbitrary sparsity pattern externally) maybe I can simplify it enough to efficiently use several libraries. $\endgroup$ – Godric Seer Sep 1 '13 at 13:28
  • $\begingroup$ From your edit, since I am working with a static mesh, finite volume code, I know my memory layout before I calculate any matrix values, which will definitely help. Also, in regards to you last sentence, I planned on using a config file to tell the code which library to use, so it wouldn't have to be "smart" in any way, other than knowing what libraries were installed. $\endgroup$ – Godric Seer Sep 1 '13 at 13:31
  • $\begingroup$ I see; that does help. Efficient sparse matrix algorithms sometimes utilize knowing the sparsity pattern in advance. $\endgroup$ – Jesse Chan Sep 1 '13 at 21:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.