Most efficient way to solve a system of linear equations

Question

I have an (n x n) symmetric matrix A and a (n x 1) vector B. Basically, I just need to solve Ax = b for x. The issue is that A is going to potentially be massive. So I'm looking for the most efficient algorithm for solving linear equations in C++. I looked over the Eigen library. Apparently it has an SVD method, but I've been told it's slow. Solving x=inverse(A)*b also seems like it would be suboptimal. Is uBLAS faster? Are there any more efficient methods? Thanks.

Edit: matrix A is positive definite and not sparse.

Solution

The best way to solve a system of linear equations of the form Ax = b is to do the following.

• decompose A into the format A = M1 * M2 (where M1 and M2 are triangular)
• Solve M1 * y = b for y using back substitution
• Solve M2 * x = y for x using back substitution

For square matrices, step 1 would use LU Decomposition.

For non square matrices, step 1 would use QR Decomposition.

If matrix A is positive definite and not sparse you'd use Cholesky Decomposition for the first step.

If you want to use eigen, you will have to first decompose it and then triangular solve it.

If this is still slow, thankfully, there are numerous linear algebra libraries available that can help improve the performance. The routine you should be looking for is dpotrs. Some libraries that have this implemented are as follows:

• Netlib's lapack: Documentation and Download (free)
• Intel's MKL: Documentation and Download (Free for non commercial use).
• AMD's ACML: Download (free)
• PLASMA: Download (free, multi core optimized)
• MAGMA : Download (free, Implemented in CUDA, OpenCL)
• CULA: Download (freemium, Implemented in CUDA).

If you are using eigen in the overall project, you can interface the LAPACK routine you need as described here.