Best C++ Matrix Library for sparse unitary matrices

Question

I am looking for a good (in the best case actively maintained) C++ matrix library. Thereby it should be templated, because I want to use a complex of rationals as numerical type. The matrices what I am dealing with are mainly sparse and unitary.

Can you please suggest libraries and also give a small explaination why to use them, because I know how to find them, but I cannot really decide what is suitable for me because I am missing the experience with them.

EDIT:

The main operations I am dealing with are matrix multiplication, scalar multiplication with a vector and kronecker product. The size of the matrices is exponential and I wanna at least be able to deal with matrices up to 1024x1024 entries.

Solution 1

Many people doing "serious" matrix stuff, rely on BLAS, adding LAPACK / ATLAS (normal matrices) or UMFPACK (sparse matrices) for more advanced math. The reason is that this code is well-tested, stable, reliable, and quite fast. Furthermore, you can buy them directly from a vendor (e.g. Intel MKL) tuned towards your architecture, but also get them for free. uBLAS mentioned in Manuel's answer is probably the standard C++ BLAS implementation. And if you need something like LAPACK later on, there are bindings to do so.

However, none of these standard libraries (BLAS / LAPACK / ATLAS or uBLAS + bindings + LAPACK / ATLAS) ticks your box for being templated and easy to use (unless uBLAS is all you'll ever need). Actually, I must admit, that I tend to call the C / Fortran interface directly when I use a BLAS / LAPACK implementation, since I often don't see much additional advantage in the uBLAS + bindings combination.

If I a need a simple-to-use, general-purpose C++ matrix library, I tend to use Eigen (I used to use NewMat in the past). Advantages:

• quite fast on Intel architecture, probably the fastest for smaller matrices
• nice interface
• almost everything you expect from a matrix library
• you can easily add new types

Disadvantages (IMO):

• single-processor [Edit: partly fixed in Eigen 3.0]
• slower for larger matrices and some advanced math than ATLAS or Intel MKL (e.g. LU decomposition) [Edit: also improved in Eigen 3.0]
• only experimental support for sparse matrices [Edit: improved in upcoming version 3.1].

Edit: The upcoming Eigen 3.1 allows some functions to use the Intel MKL (or any other BLAS / LAPACK implementation).

Solution 2

Boost uBLAS, because it's passed the Boost filter.

There are a few template libs that support sparse matrices, so it's really hard to come up with a better rationale if you're not more specific about your needs.

Solution 3

You should also try MLT and HASEM Matrix C++ Library. The last one is very well documented.

Comments

[Manuel]

One obvious advantage of a C++ wrapper is the performance gain provided by expression templates

[stephan]

@Manuel: agree. Avoiding temps is the main promise of matrix template libraries like oonumerics.org/blitz, osl.iu.edu/research/mtl, etc. But you can do most of this in C code, too, since BLAS routines do not copy matrices but work in-place where possible.

[Manuel]

you can spare some temporaries by doing everything in place but there are other things that expression templates optimize. For example, in a C library if you do add(A, add(B, add(C, D))) then the intermediate results must be computed, and each step involves a loop, so in total you have 3 loops. With expression templates the result is only evaluated once: 1 loop.

[Sard]

Eigen seems to use right to left multiplication which is a bit strange.

[rcollyer]

uBLAS isn't a wrapper, it's wholly implemented in c++.

[stephan]

@rcollyer: thanks for your input. This was indeed wrongly worded. I am surprised that you are the first person to correct me on this error. The only apology I can give is that I have always used uBLAS in combination with the bindings (hence my "wrapper" perception), which (as mentioned) wasn't my favorite approach. But you are right -- and I should have remembered, because one of the reasons we opted against uBLAS was that uBLAS prod was substantially slower for larger dense matrices than ATLAS gemm (IIRC, Intel MKL overloads prod with gemm calls for this reason).

[rcollyer]

I've never used the wrapper functionality myself, as I rejected its use because I need access to the internal representation, nor could I figure out how to bolt on the functionality I needed. So, I went with codesourcery's VSIPL++ implementation.