Hi,

I have a system of equations that I would like to solve. The system comes from some FEM formulation with Dirichlet boundary conditions imposed by Lagrange multipliers. Numerically this translates in solving the following block problem:

Code: Select all

[Kqq Kpq^T]*[Xq]=[Fq]
[Kpq 0]     [Lp] [Fp]


where Kqq is symmetric and positive-definite. I need to solve this system many times with different Fq, Fp and Kpq, for a fixed Kqq. This is how I solve the problem if I use dense matrices:

Code: Select all

Eigen::MatrixXd Kqq, Kpq;
Eigen::MatrixXd Fp, Fq;
...
Eigen::LLT<Eigen::MatrixXd> llt(Kqq); // STORE THIS ONCE
Eigen::PartialPivLU<Eigen::MatrixXd> pplu(Kpq*llt.solve(Kpq.transpose()));
Eigen::VectorXd Lp = pplu.solve(Kpq*llt.solve(Fq)-Fp);
Eigen::VectorXd Xq = llt.solve(Fq-Kpq.transpose()*Lp);


Now, my question is: how do I do the same with sparse matrices Kqq and Kpq?

Code: Select all

Eigen::SparseLLT<Eigen::SparseMatrix<double> > llt(Kqq); // OK
llt.solve(Kpq.transpose()); // FAILS


Apart from other minor issues the problem here seems to be the fact that llt.solve() takes vector-rhs only as opposed to matrix-rhs. Is my understanding correct? If so, is there any efficient workaround that I could use to overcome this issue?

Thank you for your help!

hm... you should rather use the new SimplicialCholesky:

SimplicialCholesky<SparseMAtrix<double>, Upper or Lower> llt(Kqq);
res = llt.Solve(Kpq.transpose());

if that still fails, try to evaluate the transposed expression:

res = llt.Solve(Kpq.transpose().evaluate());

Thanks, I didn't notice the class SimplicialCholesky.

However, apart from the .transpose() problem I now find something strange. It seems that if I do:

Code: Select all

Eigen::SparseMatrix<double> Kqq(50,50);
Eigen::SparseMatrix<double> Kqp(50,20);   
Eigen::SimplicialCholesky<Eigen::SparseMatrix<double> > llt(Kqq);
llt.solve(Kqp) // FAILS HERE;


then I get:

error: no matching function for call to 'Eigen::SimplicialCholesky<Eigen::SparseMatrix<double, 0, int>, 1>::solve(Eigen::SparseMatrix<double, 0, int>&)'

However if I do:

Code: Select all

Eigen::SparseMatrix<double> Kqq(50,50);
Eigen::MatrixXd Kqp(50,20);
Eigen::SimplicialCholesky<Eigen::SparseMatrix<double> > llt(Kqq);
llt.solve(Kqp) // OK;


my code compiles fine.

Am I doing something wrong? I'm on Eigen 3.0.2. Is this dev branch stuff?

Thanks agin!

Indeed, solving with a sparse rhs is not implemented yet because of laziness/lack of time, however everything is there to implement it quickly.

Is Kqp really sparse in your case? because there is no built-in sparse LU solver yet in Eigen. You will have to either rely on our wrapper to UmfPack or SuperLU, or use the BiCGSTAB iterative solver that I'll commit very soon.

I see. Yes, my Kqp is definitely sparse. So far I've been using another library for this but I was hoping to get Eigen to do the job (the rest of the code uses Eigen). Ok, I'll just wait for this to be available then. What is the target release for the sparse module? 3.1?

Thanks.