Hello,

I use BiCGSTAB solver. When I use solve method, it works fine. But when I use solveWithGuess it cannot solve the problem (comes out with 100 iterations and not Success). It is somehow strange because providing an initial guess should help the solver. Here is my code:

Code: Select all

void projection::solve_p(domain *d, psolver *solver)
{
    double *rhs = get_pressure_rhs(d);

    Eigen::VectorXd b(d->n), x(d->n), x0(d->n);
    for (int i = 0; i < d->n; ++i)
    {
        b[i] = rhs[i];
        x0[i] = d->p[i];
    }
    delete[] rhs;

    x = solver->solveWithGuess(b, x0);
    // x = solver->solve(b);

    for (int i = 0; i < d->n; ++i) d->p[i] = x[i];
}


What should I do?

which Eigen version? what about the residual with x0 and the resulting x?

I am using Eigen 3. I am using Ubuntu package.
I did not check residuals etc., because in the first place I expect a decrease in the number of iterations.

In Eigen3.2 we have a restart strategy such that if your initial guess is problematic, we restart with another guess. Checking the residuals will also help to understand what is wrong with your guess. It's easy to compute: (A*x-b).norm()/b.norm()

The norm of x0 is something between 0.0055 and 0.0065 in different time steps. Meanwhile, the norm of produced x is between 5e-17 to 3e-15.

melmi wrote:the norm of produced x is between 5e-17 to 3e-15.

it cannot be better with double, so it should report Success. Could you try with Eigen 3.2 because I see there has been a few fixes in BiCGSTAB.

Thank you ggael. I installed version 3.2 and now everything works fine and in a logical manner.

I've been having a similar problem with BiCGSTAB.
Here is an example where I solve a system and use the solution as a subsequent guess.
On the first solve I get a very reasonable answer with a very low residual (5.3084e-011) in a short time
On solvewithguess however, I get a completely unreasonable answer with a huge residual (8.29119e+016) in about 10x the time as the original solve.
I would have expected solvewithguess to return (about) the same result after the first iteration. Why doesn't this happen?
This problem persists in Eigen 3.2.1 and 3.2.2

Code: Select all

   VectorXd PSolver::Solve(SpMat& mat, VectorXd& b)
   {
      Eigen::BiCGSTAB<SpMat, Eigen::IncompleteLUT<SpMat::Scalar>> solver;
      
      try{
         solver.compute(mat);
      }
      catch (std::exception e)
      {
         std::cout << e.what() << std::endl;
         throw e;
      }

      if (solver.info() != Success)
      {
         std::cout << solver.info();
         throw std::exception("Error computing matrix");
      }

      VectorXd res = solver.solve(b);
      std::cout << (mat*res - b).norm() / b.norm() << std::endl;
      res = solver.solveWithGuess(b, res);
      std::cout << (mat*res - b).norm() / b.norm() << std::endl;

      if (solver.info() != Success)
      {
         std::cout << solver.info();
         throw std::exception("Error solving matrix");
      }
      return res;
   }

We cannot help without the matrix or a way to reproduce the issue. See this example (http://eigen.tuxfamily.org/dox-devel/gr ... tml#title3) to export 'mat' and 'b' to files.

nevermind:

Code: Select all

https://bitbucket.org/eigen/eigen/commits/51c2effe19bd/
Summary:     Do not apply the preconditioner before starting the iterations as this might destroy a very good initial guess.


Sorry, I missed this response earlier. Thank you! Bug is fixed!

I am facing a similar issue.

@ggael, What is your point about "Do not apply the preconditioner before starting the iterations as this might destroy a very good initial guess." I am not able to figure out where tvt173 is applying it before starting the iterations. I am seeing only compute being called, which computes the preconditioner.

@tvt173, how did you fix the bug?

This is the summary of the commit that fixed the bug within Eigen! Download the latest version (https://bitbucket.org/eigen/eigen/get/3.2.tar.gz) to get the fix.

Thanks for your reply. I am using the latest version of Eigen and still getting the error. Here is a toy example:

Code: Select all

    #include <Eigen/IterativeLinearSolvers>
    #include <iostream>
    #include <cstdlib>
   
    int main()
    {
     int n = 7;
   
     Eigen::SparseMatrix<double> A;
     A.resize(n,n);
   
     Eigen::VectorXd x(n), xWithGuess(n), b(n), x0(n);
   
     for (int i = 0; i != n; ++i) {
      for (int j = 0; j != n; ++j) {
    //   if (i == j) {
        A.insert(i,j) = rand()%100 + 1;
    //   }
      }
      b[i] = rand()%100 + 1;
      x0[i] = 0;
     }
   
     Eigen::BiCGSTAB<Eigen::SparseMatrix<double> > cg;
     Eigen::BiCGSTAB<Eigen::SparseMatrix<double> > cgWithGuess;
   
     cg.setMaxIterations(1);
     cg.setTolerance(0.000001);
     
     cg.compute(A);
     x = cg.solve(b);
   
     cgWithGuess.setMaxIterations(1);
     cgWithGuess.setTolerance(0.000001);
   
     cgWithGuess.compute(A);
     xWithGuess = cgWithGuess.solveWithGuess(b, x0);
   
     for (int i = 0; i != n; ++i) {
      if (x[i] != xWithGuess[i]) {
       std::cout<<"problem case "<<i<<" "<<x[i]<<" "<<xWithGuess[i]<<std::endl;
      }
     }
   
     return 0;
    }

Here is a typical output:
problem case 0 83.66 0.22452
problem case 1 21.7196 0.170796
problem case 2 -7.39731 0.065385
problem case 3 12.3621 -0.0844874
problem case 4 -32.8539 0.0880718
problem case 5 -28.2271 0.0705384
problem case 6 -14.0711 0.272096

How is your example relevant? It does only one iteration starting from a different initial position (b versus 0).