How to re-use the coefficient matrix decomposition result when solving linear systems by Eigen C++

Question

My problem needs to solve dense symmetric linear systems something like:

A x = b, A y = x, A z = y+x,...in sequence.

In Eigen C++, if I take advantage of the symmetry of A by using:

x=A.ldlt().solve(b);
y=A.ldlt().solve(x);
z=A.ldlt().solve(y+x);

The same matrix A has to be Cholesky factorized many times.

Similar issue may exists for householderQR and LU algorithms.

How can I reuse the coefficient matrix decomposition result for other systems?

It seems this is equivalent to ask whether there is a simple back-substitution solver for upper triangular linear systems.

Answer 1

You need to declare a LDLT object like this:

LDLT<MatrixXd> ldlt(A);
x = ldlt.solve(b);
y = ldlt.solve(x);
...

Answer 2

Given that the accepted answer has typos and does not give full code, I decided to post mine.

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
    MatrixXd A;
    A.resize(3,3);
    VectorXd b;
    b.resize(3);

    A <<
      13, 5, 7 ,
      5 , 9, 3 ,
      7 , 3, 11;
    b << 3, 3, 4;

    cout << "Here is the matrix A:\n" << A << endl;
    cout << "Here is the vector b:\n" << b << endl;
    LDLT<MatrixXd> L = A.ldlt();
    VectorXd x = L.solve(b);
    cout << "The solution is:\n" << x << endl;
    cout << "reconstructed b is:\n"  << A*x << endl; 

}