Armadillo C++ LU decomposition

Question

I am using the Armadillo C++ library for solving linear systems of medium/large dimensions (1000-5000 equations).

Since I have to solve different linear systems

AX=b

in which A is always the same and B changes, I would like to LU factorize A only once and reuse the LU factorization with different b. Unfortunately I do not know how to perform this kind of operations in Armadillo.

What I did was just the LU factorization of the A matrix:

arma::mat A;
// ... fill the A matrix ...
arma::mat P,L,U;
arma::lu(L, U, P, A);

But now I would like to use the matrices P, L and U to solve several linear systems with different b vectors.

Could you help me please?

Maybe you could tell us what you have tried already? Perhaps you can edit your question to include a [SSCCE](http://sscce.org/) of your attempt? — Some programmer dude, Jul 30 '13 at 07:27

score 3 · Answer 1 · edited Apr 13 '17 at 12:53

Since A = P.t()*L*U (where equality is only approximate due to rounding errors), solving for x in P.t()*L*U*x = b requires to permute rows of B and performing forward and back substitution:

x = solve(trimatu(U), solve(trimatl(L), P*b) );

Due to the lack of a true triangular solver in armadillo, and a fast way to perform row permutation, this procedure will not be very efficient, with respect to a direct call to the relevant computational LAPACK subroutines.

General advice is to avoid explicit LU decomposition in higher level libraries, like armadillo.

if all different b's are known at the same time, store them as columns in a rectangular matrix B and X = solve(A,B);
if the different b's are known one at a time, then precomputing AINV = A.i(); and x = AINV*b; will be more efficient if the number of different r.h.s. vectors is big enough. See this answer to a similar question.

Armadillo C++ LU decomposition

1 Answers1