Commit be9468dc authored by Jakub Klinkovský's avatar Jakub Klinkovský
Browse files

Fixed bug in GMRES solver 

The call to the nextIteration() method must occur exactly once per
iteration, otherwise incrementing the total iterations counter in
wrong places might cause serious bugs. For example if minIterations is
equal to 1 and the initial residue is smaller than convergenceResidue,
the loop on line 154 in tnlGMRESSolver_impl.h will be entered, but the
m-loop on line 179 will be skipped and on line 263 update() will be
called with empty vectors, which will lead to division by zero.
parent 11be86b5

src/solvers/linear/krylov/tnlGMRESSolver_impl.h

+2 −5
Original line number Diff line number Diff line
@@ -151,7 +151,7 @@ bool tnlGMRESSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector& 
   this->setResidue( beta / normb );



   tnlSharedVector< RealType, DeviceType, IndexType > vi, vk;

   while( this->nextIteration() )

   while( this->checkNextIteration() )

   {

      const IndexType m = restarting;

      for( IndexType i = 0; i < m + 1; i ++ )
@@ -251,10 +251,7 @@ bool tnlGMRESSolver< Matrix, Preconditioner > :: solve( const Vector& b, Vector& 
                             sn[ i ] );



         this->setResidue( fabs( s[ i + 1 ] ) / normb );

         // The nextIteration() method should not be called here, because it increments

         // the iteration counter. The condition below is slightly different than in

         // nextIteration(), because it first increments and then compares.

         if( this->getIterations() >= this->getMinIterations() && this->getResidue() < this->getConvergenceResidue() ) {

         if( ! this->checkNextIteration() ) {

            update( i, m, _H, _s, _v, x );

            this->refreshSolverMonitor( true );

            return this->checkConvergence();

src/solvers/tnlIterativeSolver.h

+2 −0
Original line number Diff line number Diff line
@@ -65,6 +65,8 @@ class tnlIterativeSolver 


   bool nextIteration();



   bool checkNextIteration();



   bool checkConvergence();



   void refreshSolverMonitor( bool force = false );

src/solvers/tnlIterativeSolver_impl.h

+11 −3
Original line number Diff line number Diff line
@@ -89,10 +89,18 @@ void tnlIterativeSolver< Real, Index> :: resetIterations() 


template< typename Real, typename Index >

bool tnlIterativeSolver< Real, Index> :: nextIteration()

{

   // this->checkNextIteration() must be called before the iteration counter is incremented

   bool result = this->checkNextIteration();

   this->currentIteration++;

   return result;

}



template< typename Real, typename Index >

bool tnlIterativeSolver< Real, Index> :: checkNextIteration()

{

   // TODO: fix

   //tnlAssert( solverMonitor, );

   this->currentIteration++;

   if( this->solverMonitor )

   {

      solverMonitor->setIterations( this->currentIteration );
@@ -103,8 +111,8 @@ bool tnlIterativeSolver< Real, Index> :: nextIteration() 


   if( std::isnan( this->getResidue() ) || 

       this->getIterations() > this->getMaxIterations()  ||

       ( this->getResidue() > this->getDivergenceResidue() && this->getIterations() > this->getMinIterations() ) ||

       ( this->getResidue() < this->getConvergenceResidue() && this->getIterations() > this->minIterations ) ) 

       ( this->getResidue() > this->getDivergenceResidue() && this->getIterations() >= this->getMinIterations() ) ||

       ( this->getResidue() < this->getConvergenceResidue() && this->getIterations() >= this->getMinIterations() ) )

      return false;

   return true;

}