Hello,

I'm having some troubles with the conjugate gradient iterative solver. I get the following error:

Code: Select all

In file included from /usr/local/include/eigen3/Eigen/Sparse:24:
In file included from /usr/local/include/eigen3/Eigen/IterativeLinearSolvers:32:
/usr/local/include/eigen3/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h:274:21: error: no matching
member function for call to '_solve_sparse'
    dec().derived()._solve_sparse(rhs(),dst);
    ~~~~~~~~~~~~~~~~^~~~~~~~~~~~~
/usr/local/include/eigen3/Eigen/src/misc/SparseSolve.h:46:75: note: in instantiation of function template
specialization
'Eigen::internal::sparse_solve_retval<Eigen::IterativeSolverBase<Eigen::ConjugateGradient<Eigen::SparseMatrix<double,
0, int>, 3, Eigen::DiagonalPreconditioner<double> > >, Eigen::SparseMatrix<double, 0, int>
>::evalTo<Eigen::Matrix<double, -1, 1, 0, -1, 1> >' requested here
    static_cast<const sparse_solve_retval<DecompositionType,Rhs>*>(this)->evalTo(dst);
                                                                          ^
/usr/local/include/eigen3/Eigen/src/Core/ReturnByValue.h:61:42: note: in instantiation of function template
specialization
'Eigen::internal::sparse_solve_retval_base<Eigen::IterativeSolverBase<Eigen::ConjugateGradient<Eigen::SparseMatrix<double,
0, int>, 3, Eigen::DiagonalPreconditioner<double> > >, Eigen::SparseMatrix<double, 0, int>
>::evalTo<Eigen::Matrix<double, -1, 1, 0, -1, 1> >' requested here
    { static_cast<const Derived*>(this)->evalTo(dst); }
                                         ^
/usr/local/include/eigen3/Eigen/src/Core/Matrix.h:311:13: note: in instantiation of function template
specialization
'Eigen::ReturnByValue<Eigen::internal::sparse_solve_retval_base<Eigen::IterativeSolverBase<Eigen::ConjugateGradient<Eigen::SparseMatrix<double,
0, int>, 3, Eigen::DiagonalPreconditioner<double> > >, Eigen::SparseMatrix<double, 0, int> >
>::evalTo<Eigen::Matrix<double, -1, 1, 0, -1, 1> >' requested here
      other.evalTo(*this);
            ^
/Users/Brian/Documents/Thesis/2-Finite_elasticity/C++/mex_lab/mex_HYPER3D.cpp:205:31: note: in instantiation
of function template specialization 'Eigen::Matrix<double, -1, 1, 0, -1,
1>::Matrix<Eigen::internal::sparse_solve_retval_base<Eigen::IterativeSolverBase<Eigen::ConjugateGradient<Eigen::SparseMatrix<double,
0, int>, 3, Eigen::DiagonalPreconditioner<double> > >, Eigen::SparseMatrix<double, 0, int> > >' requested here
    Eigen::VectorXd vDISPTD = solver.solve(FORCE);
                              ^
/usr/local/include/eigen3/Eigen/src/IterativeLinearSolvers/IterativeSolverBase.h:203:8: note: candidate
template ignored: could not match 'SparseMatrix' against 'Matrix'
  void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
       ^
1 error generated.

I really can't figure out what's happening! Here's a sketch of my C++ program:

Code: Select all

//  main.cpp
#include <iostream>
#include <math.h>
#include <ctime>
#include "mex.h"
#include <vector>
#include <eigen3/Eigen/Core>
#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/Sparse>
#include <eigen3/Eigen/IterativeLinearSolvers>

#define t2D(A, i, j, N) (A[ (j)*(N) + (i) ])

typedef Eigen::Triplet<double> T;

// some stuff and variables
   
    Eigen::Matrix<double,3,3> I, F, C, L, S;
    Eigen::Matrix<double,3,4> U;
    Eigen::Matrix<double,4,3> D;
    Eigen::Matrix<double,12,1> uE;
    Eigen::Matrix<double,6,12> BD;
    Eigen::Matrix<double,9,12> BG;
    Eigen::Matrix<double,6,6> K;
    Eigen::Matrix<double,9,9> G;
    Eigen::Matrix<double,12,12> GKFe;
    Eigen::Matrix<double,12,1> FORCEe;
    Eigen::Matrix<double,6,1> vecS;
   
    Eigen::SparseMatrix<double> GKF(GDOF,GDOF);
    Eigen::SparseMatrix<double> FORCE(GDOF,1);
    std::vector<T> X;
    std::vector<T> Y;
    X.reserve(12*12*NE);
    Y.reserve(12*NE);
   
    I << 1.0, 0.0, 0.0,
         0.0, 1.0, 0.0,
         0.0, 0.0, 1.0;
   
// some stuff

    GKF.setFromTriplets(X.begin(), X.end());
    FORCE.setFromTriplets(Y.begin(), Y.end());
   
    Eigen::ConjugateGradient<Eigen::SparseMatrix<double>, Eigen::Lower|Eigen::Upper> solver;
    solver.compute(GKF);
   
    if(solver.info()!=Eigen::Success) {
        //mexPrintf(" decomposition failed \n ");
        return;
        }
   
    Eigen::VectorXd vDISPTD = solver.solve(FORCE);  // can't compile when I add this line
   
}

I would appreciate any kind of help, thanks!

Brian.

EDIT: I tried the SimplicialLLT solver and I get some other errors. Is it possible that I didn't install correctly the packages ?

Errors with SimplicialLLT :

Code: Select all

In file included from /usr/local/include/eigen3/Eigen/Sparse:21:
In file included from /usr/local/include/eigen3/Eigen/SparseCholesky:39:
/usr/local/include/eigen3/Eigen/src/SparseCholesky/SimplicialCholesky.h:663:11: error: no matching member
function for call to 'defaultEvalTo'
    this->defaultEvalTo(dst);
    ~~~~~~^~~~~~~~~~~~~
/usr/local/include/eigen3/Eigen/src/misc/SparseSolve.h:46:75: note: in instantiation of function template
specialization
'Eigen::internal::sparse_solve_retval<Eigen::SimplicialCholeskyBase<Eigen::SimplicialLLT<Eigen::SparseMatrix<double,
0, int>, 3, Eigen::AMDOrdering<int> > >, Eigen::SparseMatrix<double, 0, int> >::evalTo<Eigen::Matrix<double,
-1, 1, 0, -1, 1> >' requested here
    static_cast<const sparse_solve_retval<DecompositionType,Rhs>*>(this)->evalTo(dst);
                                                                          ^
/usr/local/include/eigen3/Eigen/src/Core/ReturnByValue.h:61:42: note: in instantiation of function template
specialization
'Eigen::internal::sparse_solve_retval_base<Eigen::SimplicialCholeskyBase<Eigen::SimplicialLLT<Eigen::SparseMatrix<double,
0, int>, 3, Eigen::AMDOrdering<int> > >, Eigen::SparseMatrix<double, 0, int> >::evalTo<Eigen::Matrix<double,
-1, 1, 0, -1, 1> >' requested here
    { static_cast<const Derived*>(this)->evalTo(dst); }
                                         ^
/usr/local/include/eigen3/Eigen/src/Core/Assign.h:529:107: note: in instantiation of function template
specialization
'Eigen::ReturnByValue<Eigen::internal::sparse_solve_retval_base<Eigen::SimplicialCholeskyBase<Eigen::SimplicialLLT<Eigen::SparseMatrix<double,
0, int>, 3, Eigen::AMDOrdering<int> > >, Eigen::SparseMatrix<double, 0, int> >
>::evalTo<Eigen::Matrix<double, -1, 1, 0, -1, 1> >' requested here
  static EIGEN_STRONG_INLINE Derived& evalTo(ActualDerived& dst, const ActualOtherDerived& other) {
  other.evalTo(dst); return dst; }
                                                                                                          ^
/usr/local/include/eigen3/Eigen/src/Core/Assign.h:585:65: note: in instantiation of function template
specialization 'Eigen::internal::assign_selector<Eigen::Matrix<double, -1, 1, 0, -1, 1>,
Eigen::internal::sparse_solve_retval_base<Eigen::SimplicialCholeskyBase<Eigen::SimplicialLLT<Eigen::SparseMatrix<double,
0, int>, 3, Eigen::AMDOrdering<int> > >, Eigen::SparseMatrix<double, 0, int> >, false,
false>::evalTo<Eigen::Matrix<double, -1, 1, 0, -1, 1>,
Eigen::ReturnByValue<Eigen::internal::sparse_solve_retval_base<Eigen::SimplicialCholeskyBase<Eigen::SimplicialLLT<Eigen::SparseMatrix<double,
0, int>, 3, Eigen::AMDOrdering<int> > >, Eigen::SparseMatrix<double, 0, int> > > >' requested here
  return internal::assign_selector<Derived,OtherDerived,false>::evalTo(derived(), other.derived());
                                                                ^
/usr/local/include/eigen3/Eigen/src/Core/PlainObjectBase.h:421:20: note: in instantiation of function
template specialization 'Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>
>::operator=<Eigen::internal::sparse_solve_retval_base<Eigen::SimplicialCholeskyBase<Eigen::SimplicialLLT<Eigen::SparseMatrix<double,
0, int>, 3, Eigen::AMDOrdering<int> > >, Eigen::SparseMatrix<double, 0, int> > >' requested here
      return Base::operator=(func);
                   ^
/usr/local/include/eigen3/Eigen/src/Core/Matrix.h:190:20: note: in instantiation of function template
specialization 'Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>
>::operator=<Eigen::internal::sparse_solve_retval_base<Eigen::SimplicialCholeskyBase<Eigen::SimplicialLLT<Eigen::SparseMatrix<double,
0, int>, 3, Eigen::AMDOrdering<int> > >, Eigen::SparseMatrix<double, 0, int> > >' requested here
      return Base::operator=(func);
                   ^
/Users/Brian/Documents/Thesis/2-Finite_elasticity/C++/mex_lab/mex_HYPER3D.cpp:208:13: note: in instantiation
of function template specialization 'Eigen::Matrix<double, -1, 1, 0, -1,
1>::operator=<Eigen::internal::sparse_solve_retval_base<Eigen::SimplicialCholeskyBase<Eigen::SimplicialLLT<Eigen::SparseMatrix<double,
0, int>, 3, Eigen::AMDOrdering<int> > >, Eigen::SparseMatrix<double, 0, int> > >' requested here
    vDISPTD = solver.solve(FORCE);
            ^
/usr/local/include/eigen3/Eigen/src/misc/SparseSolve.h:51:17: note: candidate template ignored: could not
match 'SparseMatrix' against 'Matrix'
    inline void defaultEvalTo(SparseMatrix<DestScalar,DestOptions,DestIndex>& dst) const
                ^
1 error generated.

EDIT2: I think I found where the errors come from. Looks like there is a problem with me trying to construct a sparse vector this way:

Code: Select all

Eigen::SparseMatrix<double> FORCE(GDOF,1);

I replaced it by

Code: Select all

 Eigen::SparseVector<double> FORCE(GDOF);

and filled it with a for loop (i) and FORCE.coeffRef(i). However it still does not work, is it because the linear solvers do not accept sparse vectors b in Ax = b ? Thanks again.

This works in 3.3, but not in 3.2 because the rhs is sparse whereas the result is dense. Better use a dense vector for the RHS anyway.

Hello, ok thanks for your help.