Hi:

I am using the iterative solver for dense matrix. I have a matrix A with type MatrixXd, and a vector b with type VectorXd. I simply type

x = ConjugateGradient<MatrixXd>(A).solve(b);

It just returns a vector of all entires -1.#IND. When I use the direct solver, I would get the correct result. Could anyone tell what is the problem here? Thank you!

Is A symmetric positive definite?

Yes. I also tried BICGSTAB, has the same issue. Hence I don't the matrix A cause the problem.

Which eigen version? What about A.diagonal().cwiseAbs().minCoeff() ?

I believe I am using the latest version Eigen 3.2.5, and A.diagonal().cwiseAbs.minCoeff() works well.

What does is returns ? something really positive or close to zero??

Could you try with the devel branch? (http://bitbucket.org/eigen/eigen/get/default.zip)

It returns 0.

ok, that's the problem. I see that a fix of the devel branch did not go to the 3.2 branch. So it should definitely work with the devel branch. I'll backport the fix for the next 3.2.6 release.

Thank you!

Backport done: https://bitbucket.org/eigen/eigen/commits/b32d157299b3/
If you don't feel too adventurous better use the 3.2 branch: https://bitbucket.org/eigen/eigen/get/3.2.zip

I tried this version, still got the invalid return. What is the problem here? Did I do something wrong here, since I didn't see other people ask this question. Thank you.

could you share your matrix and right hand side vector?

I am calculating the boundary integral equation in a circle, for discretization, I am using the corrected trapezoidal rule. Now I restrict my number of nodes to 10, and here is my matrix:

A =
0 -0.287157 -0.245947 1.29213 -1.36476 0.757617 -1.36476 1.29213
-0.245947 -0.287157
-0.287157 0 -0.287157 -0.245947 1.29213 -1.36476 0.757617 -1.36476
1.29213 -0.245947
-0.245947 -0.287157 0 -0.287157 -0.245947 1.29213 -1.36476 0.757617
-1.36476 1.29213
1.29213 -0.245947 -0.287157 0 -0.287157 -0.245947 1.29213 -1.36476
0.757617 -1.36476
-1.36476 1.29213 -0.245947 -0.287157 0 -0.287157 -0.245947 1.29213
-1.36476 0.757617
0.757617 -1.36476 1.29213 -0.245947 -0.287157 0 -0.287157 -0.245947
1.29213 -1.36476
-1.36476 0.757617 -1.36476 1.29213 -0.245947 -0.287157 0 -0.287157
-0.245947 1.29213
1.29213 -1.36476 0.757617 -1.36476 1.29213 -0.245947 -0.287157 0
-0.287157 -0.245947
-0.245947 1.29213 -1.36476 0.757617 -1.36476 1.29213 -0.245947 -0.287157
0 -0.287157
-0.287157 -0.245947 1.29213 -1.36476 0.757617 -1.36476 1.29213 -0.245947
-0.287157 0

b =
-0.330006
-0.357292
-0.119652
0.0582414
0.163355
0.213277
0.216202
0.172588
0.0753447
-0.0920579

I can also attach code if it helps. Thank you.

Works for me (see below). Make sure you are using the latest 3.2 branch or, since the diagonal is zero, you can also use the IdentityPreconditioner.

btw, what is the motivation for using a conjugate gradient on such problem? It is small and dense, so better use LDLT that will be much faster.

Code: Select all

#include <Eigen/Sparse>
#include <Eigen/Dense>
#include <iostream>

using namespace std;
using namespace Eigen;

int main() {
    int N = 10;
    Eigen::MatrixXd A(N,N);
    VectorXd b(N);
    A <<
      0,-0.287157,-0.245947,1.29213,-1.36476,0.757617,-1.36476,1.29213,-0.245947,-0.287157,
      -0.287157,0,-0.287157,-0.245947,1.29213,-1.36476,0.757617,-1.36476,1.29213,-0.245947,
      -0.245947,-0.287157,0,-0.287157,-0.245947,1.29213,-1.36476,0.757617,-1.36476,1.29213,
      1.29213,-0.245947,-0.287157,0,-0.287157,-0.245947,1.29213,-1.36476,0.757617,-1.36476,
      -1.36476,1.29213,-0.245947,-0.287157,0,-0.287157,-0.245947,1.29213,-1.36476,0.757617,
      0.757617,-1.36476,1.29213,-0.245947,-0.287157,0,-0.287157,-0.245947,1.29213,-1.36476,
      -1.36476,0.757617,-1.36476,1.29213,-0.245947,-0.287157,0,-0.287157,-0.245947,1.29213,
      1.29213,-1.36476,0.757617,-1.36476,1.29213,-0.245947,-0.287157,0,-0.287157,-0.245947,
      -0.245947,1.29213,-1.36476,0.757617,-1.36476,1.29213,-0.245947,-0.287157,0,-0.287157,
      -0.287157,-0.245947,1.29213,-1.36476,0.757617,-1.36476,1.29213,-0.245947,-0.287157,0;

    b << -0.330006,-0.357292,-0.119652,0.0582414,0.163355,0.213277,0.216202,0.172588,0.0753447,-0.0920579;
   
    VectorXd x = VectorXd::Zero(N);
   
    std::cout << "Eigenvalues: " << SelfAdjointEigenSolver<MatrixXd>(A).eigenvalues().transpose() << "\n";

    {
      ConjugateGradient< MatrixXd, Upper|Lower > cg;
      cg.compute(A);
      x = cg.solve(b);
      std::cout << x.transpose() << " -> " << (A*x-b).norm() / b.norm() << " in " << cg.iterations() << " iterations" << std::endl;
    }
   
    {
      ConjugateGradient< MatrixXd, Upper|Lower,  IdentityPreconditioner> cg;
      cg.compute(A);
      x = cg.solve(b);
      std::cout << x.transpose() << " -> " << (A*x-b).norm() / b.norm() << " in " << cg.iterations() << " iterations" << std::endl;
    }   
}


I use the Eigen 3.2.5. I copy your code, and I find out that if I add IdentityPreconditioner, then it works. But if I don't add IdentityPreconditioner, I still get the invalid answer.

For the later question, the real problem can be huge, I set the N to 10 is just for debugging.

Thank you!