Hello,
I want to compute and decompose a Schur complement J.M^{-1}.J^T where J and M are large sparse matrices, and M is symmetric positive definite.
I have not found an efficient way to compute this matrix product.
The inverse of M being not available, I thougt I could use a SparseLDLT decomposition and apply in to J^T:

Code: Select all

SparseLDLT<SparseMatrix<double>,Cholmod> sldlt(M);
sldlt.solveInPlace( J.transpose() );


Unfortunately, solveInPlace works only for dense matrices...
Converting J to a dense matrix is very inefficient.
Computing M^{-1} by solving a series of (dense) unit vectors is also very inefficient.
Does anyone know a good solution ?
Thanks,

well, first of all it seems you are using pretty old, deprecated stuff. I suggest you the 3.1.0-alpha2. It provides a simplified assembly mechanism based on a list of triplet:

http://eigen.tuxfamily.org/dox-devel/Tu ... rseFilling

and many new and more reliable solvers. In your case this becomes:

CholmodDecomposition<SparseMatrix<double>, Lower> solver(M);
A = solver.solve(J.transpose());
A = J * A;

this solution makes use of 'cholmod_spsolve' though, as far as I remember, this function simply uses intermediate dense vectors. The reason is that, in general, it is very unlikely that inv(M) * J is sparse even if M and J are sparse.

Hi,
thanks for the quick answer. I installed the 3.1.0-alpha2. Unfortunately, I have not been able to successfully implement your suggestion. For the following code,

Code: Select all

            typedef Eigen::SparseMatrix<double> Matrix;
            Matrix M,J,A;
            Eigen::CholmodDecomposition<Matrix>  solver(M);
            A = solver.solve(J);


, I get the following error message:

Code: Select all

/ComplianceSolver.cpp:127: erreur : conversion from ‘const Eigen::internal::sparse_solve_retval<Eigen::CholmodDecomposition<Eigen::SparseMatrix<double, 0, int>, 1>, Eigen::SparseMatrix<double, 0, int> >’ to non-scalar type ‘Eigen::SparseMatrix<double, 0, int>’ requested


How can I convert the return type to a sparse matrix ?

config: g++-4.4 (Ubuntu/Linaro 4.4.6-11ubuntu2) 4.4.6

hm... works for me:

#include <Eigen/Sparse>
#include <Eigen/CholmodSupport>

int main()
{
int n = 100;
typedef Eigen::SparseMatrix<double> Matrix;
Matrix M(n,n),J(n,n),A(n,n);
Eigen::CholmodDecomposition<Matrix> solver(M);
A = solver.solve(J);
}

compiled with 3.0.1-alpha2 and gcc 4.4.6:

g++-mp-4.4 cholmod.cpp -I /opt/local/include/ufsparse/ /opt/local/lib/libcholmod.a /opt/local/lib/libcamd.a /opt/local/lib/libamd.a /opt/local/lib/libcolamd.a /opt/local/lib/libmetis.a /opt/local/lib/libccolamd.a ../build-44/blas/libeigen_blas.dylib ../build-44/lapack/libeigen_lapack.dylib -I ../../eigen3.1/

Thanks, I will have another look.
Allow me to raise another point:

ggael wrote:
this solution makes use of 'cholmod_spsolve' though, as far as I remember, this function simply uses intermediate dense vectors. The reason is that, in general, it is very unlikely that inv(M) * J is sparse even if M and J are sparse.

Actually, in a very large class of practical problems including Finite Elements, matrix M is (block-)diagonal and inv(M) * J is sparse. And matrix J^T.inv(M).J is less sparse than J, but it is still very sparse. Typically, each row of J represent a value associated with a cell, M is associated with the nodes, and each row of J^T.inv(M).J has non-null entries only for the nodes of the associated cell or their one-ring neighbors. So a 100% sparse method would probably be worthwile. Is such a sparse method used with SparseLDLT ?

yes I'm aware that in some cases the result is expected to be sparse, and the Schur complement is one typical use case for that. None of the sparse solver libraries we support offers a true solve method with a sparse right hand side, and we also don't offer that yet in Eigen, though this is something we plane to have.

Hi,

I'm also getting the same error. If I try to compile this simple code...

Code: Select all

#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
int main(){
Eigen::SparseMatrix<double> A(10,10);
Eigen::SparseMatrix<double> B(10,10);   
    Eigen::SimplicialLLT<Eigen::SparseMatrix<double> > llt(A);   
    Eigen::SparseMatrix<double> C=llt.solve(B);
    return 0;
}


... I get:


Am I doing something wrong? I'm on Eigen 3.1.1 and I tried both with gcc 4.2 and gcc 4.7.

Thanks!

I would like to add some findings which I can't undertand:
1- The following compiles

Code: Select all

#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
int main(){
Eigen::SparseMatrix<double> A(10,10);   
Eigen::VectorXd B(10);   
    Eigen::SimplicialLLT<Eigen::SparseMatrix<double> > llt(A);   
    Eigen::VectorXd C=llt.solve(B);
    return 0;
}


2- The following also compiles

Code: Select all

#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
int main(){
Eigen::SparseMatrix<double> A(10,10);   
Eigen::MatrixXd B(10,10);   
    Eigen::SimplicialLLT<Eigen::SparseMatrix<double> > llt(A);   
    Eigen::MatrixXd C=llt.solve(B);
    return 0;
}


3- This instead fails

Code: Select all

#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
int main(){
Eigen::SparseMatrix<double> A(10,10);   
Eigen::SparseMatrix<double> B(10,10);   
    Eigen::SimplicialLLT<Eigen::SparseMatrix<double> > llt(A);   
    Eigen::MatrixXd C=llt.solve(B);
    return 0;
}


4- ... but this is compiles ok

Code: Select all

#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
int main(){
Eigen::SparseMatrix<double> A(10,10);   
Eigen::SparseMatrix<double> B(10,10);   
    Eigen::SimplicialLLT<Eigen::SparseMatrix<double> > llt(A);   
    llt.solve(B); // don't assign result
     return 0;
}


So, if B is sparse, what am I supposed to define C as?
Thanks.

The following should do the job:

Code: Select all

Eigen::SparseMatrix<double> C;
C=llt.solve(B);


The problem should be fixed now in 3.1 and devel branches (missing ctor overload)