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.