CODEWITHSUNDEEP
c++linear-algebranumericarmadillodecomposition
alberto.cuoci
alberto.cuoci

Reputation: 87

Armadillo C++ LU decomposition

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?

Upvotes: 2

Views: 2825

Answers (1)

Stefano M
Stefano M

Reputation: 4809

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.

  1. 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);
  2. 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.

Upvotes: 3

Related Questions