1

This is related to this question. I have today experimented a bit with Conjugate Gradient, in particular I experimented with max_iterations and tolerance. It is faster but not fast enough. According to the documentation it should be enough to add -fopenmp in the compilation to enable multi-threading.

I have tested using both

`omp_set_num_threads(nbrThreads);
 Eigen::setNbThreads(nbrThreads);`

It makes no difference in time if I use 5 threads or 1 thread, and that I think is a bit strange.

Secondly, another way to speed up the solution using pre-conditoning. When I try to do:

    Eigen::ConjugateGradient<Eigen::SparseMatrix<float>, Eigen::Lower, Eigen::IncompleteCholesky<float>> cg;

I get the error:

    void Eigen::IncompleteCholesky<Scalar, _UpLo, _OrderingType>::_solve_impl(const Rhs&, Dest&) const [with Rhs = Eigen::Matr
ix<float, -1, 1>; Dest = Eigen::Matrix<float, -1, 1>; Scalar = float; int _UpLo = 1; _OrderingType = Eigen::AMDOrdering<int>]: Assertion `m_factorizationIsOk && "factorize() should be called first"' failed.

Given that Eigen::SimplicialLDLT works which is a Cholesky factorization, then the incomplete should also work?

EDIT: Here is how I call cg:

Eigen::ConjugateGradient<Eigen::SparseMatrix<float>, Eigen::Lower, Eigen::IncompleteCholesky<float>> cg;
cg.setTolerance(0.01);
cg.setMaxIterations(50);
cg.analyzePattern(A_tot);
cg.compute(A_tot);
Eigen::VectorXf opt = cg.solveWithGuess(b_tot, rho_current);

Actually when you read about IterativeSolvers here, then IncompleteCholesky is not listed. Although IncompleCholesky is defined here.

Community
  • 1
  • 1
El_Loco
  • 1,716
  • 4
  • 20
  • 35
  • Regarding `IncompleteCholesky`, please add the code calling `cg`. – ggael Feb 09 '17 at 11:32
  • All calls with cg are added now. – El_Loco Feb 09 '17 at 12:17
  • 1
    ILLT probably failed. After the call to compute, check that `cg.preconditioner().info()==Eigen::Success`. (btw compute is a shortcut for analyzePattern and factorize, so no need to call both analyzePattern and compute). If it returns `NumericalIssue`, then you might try to increase the shift with `cg.preconditioner().setInitialShift(a_small_value_but_large_enough_to_get_it_work)`. – ggael Feb 09 '17 at 13:49
  • Thank you, I will try that. – El_Loco Feb 09 '17 at 15:16
  • I tried to change the initialshift, but that did not help. Maybe IncompleteCholesky is not supported? – El_Loco Feb 10 '17 at 08:43

2 Answers2

1

As explained in the documentation, you need to store the full matrix (both upper and lower triangular part), and pass Lower|Upper to ConjugateGradient to get multi-threading with cg.

ggael
  • 28,425
  • 2
  • 65
  • 71
0

The problem seems to be that I had an beta-version of Eigen installed. It runs with pre-conditioner using Eigen 3.3.2

El_Loco
  • 1,716
  • 4
  • 20
  • 35