<h3>lac</h3>
<ol>
+ <li> <p> Changed: Solver functions <code class="member">solve</code>
+ are void now. If the solver has not converged within the
+ maximum number of iterations, it throws an exception of type
+ <code class="class">ExcNoConvergence</code>.
+ <br>
+ (GK 2001/03/29)
+ </p>
+
<li> <p>
New: The functions <code class="class">FullMatrix</code>::<code
class="member">mmult</code> and <code
* l2-norm.
*/
template <class MATRIX>
- typename Solver<VECTOR>::ReturnState
+ void
solve (double &value,
const MATRIX &A,
VECTOR &x);
* l2-norm.
*/
template <class MATRIX>
- typename Solver<VECTOR>::ReturnState
+ void
solve (double &value,
const MATRIX &A,
VECTOR &x);
template <class VECTOR>
template <class MATRIX>
-typename Solver<VECTOR>::ReturnState
+void
EigenPower<VECTOR>::solve (double &value,
const MATRIX &A,
VECTOR &x)
deallog.pop();
// Output
- if (conv == SolverControl::failure)
- return exceeded;
- else
- return success;
+ AssertThrow(control().last_check() == SolverControl::success,
+ typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+ control().last_value()));
}
//----------------------------------------------------------------------//
template <class VECTOR>
template <class MATRIX>
-typename Solver<VECTOR>::ReturnState
+void
EigenInverse<VECTOR>::solve (double &value,
const MATRIX &A,
VECTOR &x)
deallog.pop();
// Output
- if (conv == SolverControl::failure)
- return exceeded;
- else
- return success;
+ AssertThrow(control().last_check() == SolverControl::success,
+ typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+ control().last_value()));
}
#endif
* changes and it is still possible in a simple way to give these additional data to
* the @p{SolverSelector} object for each solver which it may use.
*
- * @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann, 1997, 1998, 1999
+ * @author Wolfgang Bangerth, Guido Kanschat, Ralf Hartmann, 1997-2001
*/
template <class Vector = ::Vector<double> >
class Solver
{
public:
- /**
- * Declare possible return values of a
- * solver object.
- */
- enum ReturnState {
- success=0, exceeded, breakdown
- };
-
-
/**
* Constructor. Assign a control
* object which stores the required
* Access to control object.
*/
SolverControl& control() const;
-
+
+ DeclException2(ExcNoConvergence, int, double,
+ << "Iteration did not converge after " << arg1
+ << " steps. Final residual " << arg2);
+
protected:
/**
* Solve primal problem only.
*/
template<class MATRIX, class PRECONDITIONER>
- typename Solver<VECTOR>::ReturnState
+ void
solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
* The iteration loop itself.
*/
template<class MATRIX, class PRECONDITIONER>
- typename Solver<VECTOR>::ReturnState
+ bool
iterate(const MATRIX& A, const PRECONDITIONER& precondition);
};
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState
+bool
SolverBicgstab<VECTOR>::iterate(const MATRIX& A,
const PRECONDITIONER& precondition)
{
//TODO:[GK] Find better breakdown criterion
if (fabs(alpha) > 1.e10)
- return typename Solver<VECTOR>::ReturnState(breakdown);
+ return true;
s.equ(1., r, -alpha, v);
precondition.vmult(z,s);
Vx->add(alpha, y, omega, z);
r.equ(1., s, -omega, t);
- state = control().check(++step, criterion(A, *Vx, *Vb));
+ res = criterion(A, *Vx, *Vb);
+ state = control().check(++step, res);
print_vectors(step, *Vx, r, y);
}
while (state == SolverControl::iterate);
- if (state == SolverControl::success) return success;
- return exceeded;
+ return false;
}
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState
+void
SolverBicgstab<VECTOR>::solve(const MATRIX &A,
VECTOR &x,
const VECTOR &b,
step = 0;
- typename Solver<VECTOR>::ReturnState state = breakdown;
+ bool state;
do
{
if (start(A) == SolverControl::success) break;
state = iterate(A, precondition);
}
- while (state == breakdown);
+ while (state);
memory.free(Vr);
memory.free(Vrbar);
memory.free(Vv);
deallog.pop();
- return state;
+
+ AssertThrow(control().last_check() == SolverControl::success,
+ typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+ control().last_value()));
}
#include <lac/solver.h>
#include <lac/solver_control.h>
+#include <base/exceptions.h>
#include <base/logstream.h>
#include <base/subscriptor.h>
#include <cmath>
* Solver method.
*/
template<class MATRIX, class PRECONDITIONER>
- typename Solver<VECTOR>::ReturnState
+ void
solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState
+void
SolverCG<VECTOR>::solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
memory.free(Vz);
memory.free(VAp);
deallog.pop();
- return success;
+ return;
};
g.scale(-1.);
// Output
deallog.pop();
- if (conv == SolverControl::failure)
- return exceeded;
- else
- return success;
+ AssertThrow(control().last_check() == SolverControl::success,
+ typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+ control().last_value()));
};
virtual State check (const unsigned int step,
const double check_value);
+ /**
+ * Return the result of the last check operation.
+ */
+ State last_check() const;
+
/**
* Return the convergence value of last
* iteration step for which @p{check} was
*/
double tol;
+ /**
+ * Result of last check operation.
+ */
+ State lcheck;
+
/**
* Last value of the convergence criterion.
*/
* Solver method.
*/
template<class MATRIX, class PRECONDITIONER>
- typename Solver<VECTOR>::ReturnState solve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER &precondition);
+ void solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER& precondition);
DeclException1 (ExcTooFewTmpVectors,
int,
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState
-SolverGMRES<VECTOR>::solve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER &precondition)
+void
+SolverGMRES<VECTOR>::solve (const MATRIX& A,
+ VECTOR & x,
+ const VECTOR& b,
+ const PRECONDITIONER& precondition)
{
// this code was written a very
// long time ago by people not
memory.free (tmp_vectors[tmp]);
deallog.pop();
- if (iteration_state)
- return success;
- else
- return exceeded;
+
+ AssertThrow(control().last_check() == SolverControl::success,
+ typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+ control().last_value()));
};
* Solver method.
*/
template<class MATRIX, class PRECONDITIONER>
- typename Solver<VECTOR>::ReturnState
+ void
solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
// Output
deallog.pop ();
- if (conv == SolverControl::failure)
- return exceeded;
-
- return success;
+ AssertThrow(control().last_check() == SolverControl::success,
+ typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+ control().last_value()));
};
* Solver method.
*/
template<class MATRIX, class PRECONDITIONER>
- typename Solver<VECTOR>::ReturnState
+ void
solve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER& precondition);
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER& precondition);
/**
* Interface for derived class.
* The iteration loop itself.
*/
template<class MATRIX, class PRECONDITIONER>
- typename Solver<VECTOR>::ReturnState
+ bool
iterate(const MATRIX& A, const PRECONDITIONER& precondition);
/**
* The current iteration step.
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState
+void
SolverQMRS<VECTOR>::solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
step = 0;
- typename Solver<VECTOR>::ReturnState state = breakdown;
+ bool state;
do
{
deallog << "Restart step " << step << std::endl;
state = iterate(A, precondition);
}
- while (state == breakdown);
+ while (state);
// Deallocate Memory
// Output
deallog.pop();
-
- return state;
+
+ AssertThrow(control().last_check() == SolverControl::success,
+ typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+ control().last_value()));
};
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState
+bool
SolverQMRS<VECTOR>::iterate(const MATRIX& A,
const PRECONDITIONER& precondition)
{
res = v.l2_norm();
if (control().check(step, res) == SolverControl::success)
- return success;
+ return false;
p = v;
//TODO:[GK] Find a really good breakdown criterion. The absolute one detects breakdown instead of convergence
if (fabs(sigma/rho) < additional_data.breakdown)
- return breakdown;
+ return true;
// Step 3
alpha = rho/sigma;
//deallog << "alpha:" << alpha << std::endl;
else
res = sqrt((it+1)*tau);
state = control().check(step,res);
- if (state == SolverControl::success)
- return success;
- else if (state == SolverControl::failure)
- return exceeded;
+ if ((state == SolverControl::success)
+ || (state == SolverControl::failure))
+ return false;
// Step 6
if (fabs(rho) < additional_data.breakdown)
- return breakdown;
+ return true;
// Step 7
rho_old = rho;
precondition.vmult(q,v);
p.sadd(beta,v);
precondition.vmult(q,p);
}
- return exceeded;
+ return false;
}
* Solve $Ax=b$ for $x$.
*/
template<class MATRIX, class PRECONDITIONER>
- typename Solver<VECTOR>::ReturnState solve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER& precondition);
+ void solve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER& precondition);
/**
* Solve $A^Tx=b$ for $x$.
*/
template<class MATRIX, class PRECONDITIONER>
- typename Solver<VECTOR>::ReturnState Tsolve (const MATRIX &A,
- VECTOR &x,
- const VECTOR &b,
- const PRECONDITIONER& precondition);
+ void Tsolve (const MATRIX &A,
+ VECTOR &x,
+ const VECTOR &b,
+ const PRECONDITIONER& precondition);
/**
* Set the damping-coefficient.
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState
+void
SolverRichardson<VECTOR>::solve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
deallog.pop();
// Output
- if (conv == SolverControl::failure)
- return exceeded;
- else
- return success;
+ AssertThrow(control().last_check() == SolverControl::success,
+ typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+ control().last_value()));
}
template<class VECTOR>
template<class MATRIX, class PRECONDITIONER>
-typename Solver<VECTOR>::ReturnState
+void
SolverRichardson<VECTOR>::Tsolve (const MATRIX &A,
VECTOR &x,
const VECTOR &b,
deallog.pop();
// Output
- if (conv == SolverControl::failure)
- return exceeded;
- else
- return success;
+ AssertThrow(control().last_check() == SolverControl::success,
+ typename Solver<VECTOR>::ExcNoConvergence(control().last_step(),
+ control().last_value()));
}
*
*/
template<class Matrix, class Preconditioner>
- typename Solver<Vector>::ReturnState solve(const Matrix &A,
- Vector &x,
- const Vector &b,
- const Preconditioner &precond) const;
+ void solve(const Matrix &A,
+ Vector &x,
+ const Vector &b,
+ const Preconditioner &precond) const;
/**
* Set the additional data. For more
template <class Vector>
template<class Matrix, class Preconditioner>
-typename Solver<Vector>::ReturnState
+void
SolverSelector<Vector>::solve(const Matrix &A,
Vector &x,
const Vector &b,
{
SolverRichardson<Vector> solver(*control,*vector_memory,
richardson_data);
- return solver.solve(A,x,b,precond);
+ solver.solve(A,x,b,precond);
}
else if (solver_name=="cg")
{
SolverCG<Vector> solver(*control,*vector_memory,
cg_data);
- return solver.solve(A,x,b,precond);
+ solver.solve(A,x,b,precond);
}
else if (solver_name=="bicgstab")
{
SolverBicgstab<Vector> solver(*control,*vector_memory,
bicgstab_data);
- return solver.solve(A,x,b,precond);
+ solver.solve(A,x,b,precond);
}
else if (solver_name=="gmres")
{
SolverGMRES<Vector> solver(*control,*vector_memory,
gmres_data);
- return solver.solve(A,x,b,precond);
+ solver.solve(A,x,b,precond);
}
else
Assert(false,ExcSolverDoesNotExist(solver_name));
-
- return Solver<Vector>::breakdown;
};
if (_log_result)
deallog << "Failure step " << step
<< " value " << check_value << std::endl;
+ lcheck = failure;
return failure;
}
if (_log_result)
deallog << "Convergence step " << step
<< " value " << check_value << std::endl;
+ lcheck = success;
return success;
}
-
+ lcheck = iterate;
return iterate;
}
+SolverControl::State
+SolverControl::last_check() const
+{
+ return lcheck;
+}
+
+
double
SolverControl::last_value() const
{
if (_log_result)
deallog << "Convergence step " << step
<< " value " << check_value << std::endl;
+ lcheck = success;
return success;
}
else