Skip to content
Snippets Groups Projects
Commit 90236d45 authored by Christian Engwer's avatar Christian Engwer
Browse files

Merge branch 'const_correctness_in_iterative_solvers' into 'master' 

Const correctness in iterative solvers

See merge request !421
parents 13b13352 1c8084b7
No related branches found
No related tags found
1 merge request!421Const correctness in iterative solvers
Pipeline #54591 failed
Stage: test
Hide whitespace changes
Inline Side-by-side
CHANGELOG.md
+ 4
− 0
View file @ 90236d45
# Master (will become release 2.9)
- Add `const` qualifier to `LinearOperator` and `ScalarProduct` in
`IterativeSolver`. In particular, the constructors of iterative solvers have
changed.
- Solvers are more robust if used with multiple right-hand sides and one lane starts with the exact solution.
- Added a function to write nested matrices as SVG objects: `writeSVGMatrix(...)`
......
dune/istl/solver.hh
+ 8
− 8
View file @ 90236d45
......@@ -225,7 +225,7 @@ namespace Dune
<li> 2 : print line for each iteration </li>
</ul>
*/
IterativeSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int maxit, int verbose) :
IterativeSolver (const LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int maxit, int verbose) :
_op(stackobject_to_shared_ptr(op)),
_prec(stackobject_to_shared_ptr(prec)),
_sp(new SeqScalarProduct<X>),
......@@ -257,7 +257,7 @@ namespace Dune
<li> 2 : print line for each iteration </li>
</ul>
*/
IterativeSolver (LinearOperator<X,Y>& op, ScalarProduct<X>& sp, Preconditioner<X,Y>& prec,
IterativeSolver (const LinearOperator<X,Y>& op, const ScalarProduct<X>& sp, Preconditioner<X,Y>& prec,
scalar_real_type reduction, int maxit, int verbose) :
_op(stackobject_to_shared_ptr(op)),
_prec(stackobject_to_shared_ptr(prec)),
......@@ -285,7 +285,7 @@ namespace Dune
See \ref ISTL_Factory for the ParameterTree layout and examples.
*/
IterativeSolver (std::shared_ptr<LinearOperator<X,Y> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
IterativeSolver (std::shared_ptr<const LinearOperator<X,Y> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
IterativeSolver(op,std::make_shared<SeqScalarProduct<X>>(),prec,
configuration.get<real_type>("reduction"),
configuration.get<int>("maxit"),
......@@ -308,7 +308,7 @@ namespace Dune
See \ref ISTL_Factory for the ParameterTree layout and examples.
*/
IterativeSolver (std::shared_ptr<LinearOperator<X,Y> > op, std::shared_ptr<ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
IterativeSolver (std::shared_ptr<const LinearOperator<X,Y> > op, std::shared_ptr<const ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
IterativeSolver(op,sp,prec,
configuration.get<scalar_real_type>("reduction"),
configuration.get<int>("maxit"),
......@@ -335,8 +335,8 @@ namespace Dune
<li> 2 : print line for each iteration </li>
</ul>
*/
IterativeSolver (std::shared_ptr<LinearOperator<X,Y>> op,
std::shared_ptr<ScalarProduct<X>> sp,
IterativeSolver (std::shared_ptr<const LinearOperator<X,Y>> op,
std::shared_ptr<const ScalarProduct<X>> sp,
std::shared_ptr<Preconditioner<X,Y>> prec,
scalar_real_type reduction, int maxit, int verbose) :
_op(op),
......@@ -499,9 +499,9 @@ namespace Dune
};
protected:
std::shared_ptr<LinearOperator<X,Y>> _op;
std::shared_ptr<const LinearOperator<X,Y>> _op;
std::shared_ptr<Preconditioner<X,Y>> _prec;
std::shared_ptr<ScalarProduct<X>> _sp;
std::shared_ptr<const ScalarProduct<X>> _sp;
scalar_real_type _reduction;
int _maxit;
int _verbose;
......
dune/istl/solvers.hh
+ 37
− 37
View file @ 90236d45
......@@ -212,12 +212,12 @@ namespace Dune {
/*!
\brief Constructor to initialize a CG solver.
\copydetails IterativeSolver::IterativeSolver(LinearOperator<X,Y>&, Preconditioner<X,Y>&, real_type, int, int)
\copydetails IterativeSolver::IterativeSolver(const LinearOperator<X,Y>&, Preconditioner<X,Y>&, real_type, int, int)
\param condition_estimate Whether to calculate an estimate of the condition number.
The estimate is given in the InverseOperatorResult returned by apply().
This is only supported for float and double field types.
*/
CGSolver (LinearOperator<X,X>& op, Preconditioner<X,X>& prec,
CGSolver (const LinearOperator<X,X>& op, Preconditioner<X,X>& prec,
scalar_real_type reduction, int maxit, int verbose, bool condition_estimate) : IterativeSolver<X,X>(op, prec, reduction, maxit, verbose),
condition_estimate_(condition_estimate)
{
......@@ -229,12 +229,12 @@ namespace Dune {
/*!
\brief Constructor to initialize a CG solver.
\copydetails IterativeSolver::IterativeSolver(LinearOperator<X,Y>&, ScalarProduct<X>&, Preconditioner<X,Y>&, real_type, int, int)
\copydetails IterativeSolver::IterativeSolver(const LinearOperator<X,Y>&, const ScalarProduct<X>&, Preconditioner<X,Y>&, real_type, int, int)
\param condition_estimate Whether to calculate an estimate of the condition number.
The estimate is given in the InverseOperatorResult returned by apply().
This is only supported for float and double field types.
*/
CGSolver (LinearOperator<X,X>& op, ScalarProduct<X>& sp, Preconditioner<X,X>& prec,
CGSolver (const LinearOperator<X,X>& op, const ScalarProduct<X>& sp, Preconditioner<X,X>& prec,
scalar_real_type reduction, int maxit, int verbose, bool condition_estimate) : IterativeSolver<X,X>(op, sp, prec, reduction, maxit, verbose),
condition_estimate_(condition_estimate)
{
......@@ -246,12 +246,12 @@ namespace Dune {
/*!
\brief Constructor to initialize a CG solver.
\copydetails IterativeSolver::IterativeSolver(std::shared_ptr<LinearOperator<X,Y>>, std::shared_ptr<ScalarProduct<X>>, std::shared_ptr<Preconditioner<X,Y>>, real_type, int, int)
\copydetails IterativeSolver::IterativeSolver(std::shared_ptr<const LinearOperator<X,Y>>, std::shared_ptr<ScalarProduct<X>>, std::shared_ptr<Preconditioner<X,Y>>, real_type, int, int)
\param condition_estimate Whether to calculate an estimate of the condition number.
The estimate is given in the InverseOperatorResult returned by apply().
This is only supported for float and double field types.
*/
CGSolver (std::shared_ptr<LinearOperator<X,X>> op, std::shared_ptr<ScalarProduct<X>> sp,
CGSolver (std::shared_ptr<const LinearOperator<X,X>> op, std::shared_ptr<ScalarProduct<X>> sp,
std::shared_ptr<Preconditioner<X,X>> prec,
scalar_real_type reduction, int maxit, int verbose, bool condition_estimate)
: IterativeSolver<X,X>(op, sp, prec, reduction, maxit, verbose),
......@@ -842,10 +842,10 @@ namespace Dune {
/*!
\brief Set up RestartedGMResSolver solver.
\copydoc LoopSolver::LoopSolver(L&,P&,double,int,int)
\copydoc LoopSolver::LoopSolver(const L&,P&,double,int,int)
\param restart number of GMRes cycles before restart
*/
RestartedGMResSolver (LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int restart, int maxit, int verbose) :
RestartedGMResSolver (const LinearOperator<X,Y>& op, Preconditioner<X,Y>& prec, scalar_real_type reduction, int restart, int maxit, int verbose) :
IterativeSolver<X,Y>::IterativeSolver(op,prec,reduction,maxit,verbose),
_restart(restart)
{}
......@@ -853,10 +853,10 @@ namespace Dune {
/*!
\brief Set up RestartedGMResSolver solver.
\copydoc LoopSolver::LoopSolver(L&,S&,P&,double,int,int)
\copydoc LoopSolver::LoopSolver(const L&, const S&,P&,double,int,int)
\param restart number of GMRes cycles before restart
*/
RestartedGMResSolver (LinearOperator<X,Y>& op, ScalarProduct<X>& sp, Preconditioner<X,Y>& prec, scalar_real_type reduction, int restart, int maxit, int verbose) :
RestartedGMResSolver (const LinearOperator<X,Y>& op, const ScalarProduct<X>& sp, Preconditioner<X,Y>& prec, scalar_real_type reduction, int restart, int maxit, int verbose) :
IterativeSolver<X,Y>::IterativeSolver(op,sp,prec,reduction,maxit,verbose),
_restart(restart)
{}
......@@ -864,7 +864,7 @@ namespace Dune {
/*!
\brief Constructor.
\copydoc IterativeSolver::IterativeSolver(L&,S&,P&,const ParameterTree&)
\copydoc IterativeSolver::IterativeSolver(const L&, const S&,P&,const ParameterTree&)
Additional parameter:
ParameterTree Key | Meaning
......@@ -873,12 +873,12 @@ namespace Dune {
See \ref ISTL_Factory for the ParameterTree layout and examples.
*/
RestartedGMResSolver (std::shared_ptr<LinearOperator<X,Y> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
RestartedGMResSolver (std::shared_ptr<const LinearOperator<X,Y> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
IterativeSolver<X,Y>::IterativeSolver(op,prec,configuration),
_restart(configuration.get<int>("restart"))
{}
RestartedGMResSolver (std::shared_ptr<LinearOperator<X,Y> > op, std::shared_ptr<ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
RestartedGMResSolver (std::shared_ptr<const LinearOperator<X,Y> > op, std::shared_ptr<const ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
IterativeSolver<X,Y>::IterativeSolver(op,sp,prec,configuration),
_restart(configuration.get<int>("restart"))
{}
......@@ -886,11 +886,11 @@ namespace Dune {
/*!
\brief Set up RestartedGMResSolver solver.
\copydoc LoopSolver::LoopSolver(std::shared_ptr<L>,std::shared_ptr<S>,std::shared_ptr<P>,double,int,int)
\copydoc LoopSolver::LoopSolver(std::shared_ptr<const L>,std::shared_ptr<const S>,std::shared_ptr<P>,double,int,int)
\param restart number of GMRes cycles before restart
*/
RestartedGMResSolver (std::shared_ptr<LinearOperator<X,Y>> op,
std::shared_ptr<ScalarProduct<X>> sp,
RestartedGMResSolver (std::shared_ptr<const LinearOperator<X,Y>> op,
std::shared_ptr<const ScalarProduct<X>> sp,
std::shared_ptr<Preconditioner<X,Y>> prec,
scalar_real_type reduction, int restart, int maxit, int verbose) :
IterativeSolver<X,Y>::IterativeSolver(op,sp,prec,reduction,maxit,verbose),
......@@ -1323,10 +1323,10 @@ private:
/*!
\brief Set up nonlinear preconditioned conjugate gradient solver.
\copydoc LoopSolver::LoopSolver(L&,P&,double,int,int)
\copydoc LoopSolver::LoopSolver(const L&,P&,double,int,int)
\param restart number of GMRes cycles before restart
*/
GeneralizedPCGSolver (LinearOperator<X,X>& op, Preconditioner<X,X>& prec, scalar_real_type reduction, int maxit, int verbose, int restart = 10) :
GeneralizedPCGSolver (const LinearOperator<X,X>& op, Preconditioner<X,X>& prec, scalar_real_type reduction, int maxit, int verbose, int restart = 10) :
IterativeSolver<X,X>::IterativeSolver(op,prec,reduction,maxit,verbose),
_restart(restart)
{}
......@@ -1334,11 +1334,11 @@ private:
/*!
\brief Set up nonlinear preconditioned conjugate gradient solver.
\copydoc LoopSolver::LoopSolver(L&,S&,P&,double,int,int)
\copydoc LoopSolver::LoopSolver(const L&, const S&,P&,double,int,int)
\param restart When to restart the construction of
the Krylov search space.
*/
GeneralizedPCGSolver (LinearOperator<X,X>& op, ScalarProduct<X>& sp, Preconditioner<X,X>& prec, scalar_real_type reduction, int maxit, int verbose, int restart = 10) :
GeneralizedPCGSolver (const LinearOperator<X,X>& op, const ScalarProduct<X>& sp, Preconditioner<X,X>& prec, scalar_real_type reduction, int maxit, int verbose, int restart = 10) :
IterativeSolver<X,X>::IterativeSolver(op,sp,prec,reduction,maxit,verbose),
_restart(restart)
{}
......@@ -1347,7 +1347,7 @@ private:
/*!
\brief Constructor.
\copydoc IterativeSolver::IterativeSolver(L&,S&,P&,const ParameterTree&)
\copydoc IterativeSolver::IterativeSolver(const L&,const S&,P&,const ParameterTree&)
Additional parameter:
ParameterTree Key | Meaning
......@@ -1356,24 +1356,24 @@ private:
See \ref ISTL_Factory for the ParameterTree layout and examples.
*/
GeneralizedPCGSolver (std::shared_ptr<LinearOperator<X,X> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
GeneralizedPCGSolver (std::shared_ptr<const LinearOperator<X,X> > op, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
IterativeSolver<X,X>::IterativeSolver(op,prec,configuration),
_restart(configuration.get<int>("restart"))
{}
GeneralizedPCGSolver (std::shared_ptr<LinearOperator<X,X> > op, std::shared_ptr<ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
GeneralizedPCGSolver (std::shared_ptr<const LinearOperator<X,X> > op, std::shared_ptr<const ScalarProduct<X> > sp, std::shared_ptr<Preconditioner<X,X> > prec, const ParameterTree& configuration) :
IterativeSolver<X,X>::IterativeSolver(op,sp,prec,configuration),
_restart(configuration.get<int>("restart"))
{}
/*!
\brief Set up nonlinear preconditioned conjugate gradient solver.
\copydoc LoopSolver::LoopSolver(std::shared_ptr<L>,std::shared_ptr<S>,std::shared_ptr<P>,double,int,int)
\copydoc LoopSolver::LoopSolver(std::shared_ptr<const L>,std::shared_ptr<const S>,std::shared_ptr<P>,double,int,int)
\param restart When to restart the construction of
the Krylov search space.
*/
GeneralizedPCGSolver (std::shared_ptr<LinearOperator<X,X>> op,
std::shared_ptr<ScalarProduct<X>> sp,
GeneralizedPCGSolver (std::shared_ptr<const LinearOperator<X,X>> op,
std::shared_ptr<const ScalarProduct<X>> sp,
std::shared_ptr<Preconditioner<X,X>> prec,
scalar_real_type reduction, int maxit, int verbose,
int restart = 10) :
......@@ -1511,31 +1511,31 @@ private:
using IterativeSolver<X,X>::apply;
/*!
\brief Constructor to initialize a RestartedFCG solver.
\copydetails IterativeSolver::IterativeSolver(LinearOperator<X,Y>&, Preconditioner<X,Y>&, real_type, int, int, int)
\copydetails IterativeSolver::IterativeSolver(const LinearOperator<X,Y>&, Preconditioner<X,Y>&, real_type, int, int, int)
\param mmax is the maximal number of previous vectors which are orthogonalized against the new search direction.
*/
RestartedFCGSolver (LinearOperator<X,X>& op, Preconditioner<X,X>& prec,
RestartedFCGSolver (const LinearOperator<X,X>& op, Preconditioner<X,X>& prec,
scalar_real_type reduction, int maxit, int verbose, int mmax = 10) : IterativeSolver<X,X>(op, prec, reduction, maxit, verbose), _mmax(mmax)
{
}
/*!
\brief Constructor to initialize a RestartedFCG solver.
\copydetails IterativeSolver::IterativeSolver(LinearOperator<X,Y>&, ScalarProduct<X>&, Preconditioner<X,Y>&, real_type, int, int,int)
\copydetails IterativeSolver::IterativeSolver(const LinearOperator<X,Y>&, const ScalarProduct<X>&, Preconditioner<X,Y>&, real_type, int, int,int)
\param mmax is the maximal number of previous vectors which are orthogonalized against the new search direction.
*/
RestartedFCGSolver (LinearOperator<X,X>& op, ScalarProduct<X>& sp, Preconditioner<X,X>& prec,
RestartedFCGSolver (const LinearOperator<X,X>& op, const ScalarProduct<X>& sp, Preconditioner<X,X>& prec,
scalar_real_type reduction, int maxit, int verbose, int mmax = 10) : IterativeSolver<X,X>(op, sp, prec, reduction, maxit, verbose), _mmax(mmax)
{
}
/*!
\brief Constructor to initialize a RestartedFCG solver.
\copydetails IterativeSolver::IterativeSolver(std::shared_ptr<LinearOperator<X,Y>>, std::shared_ptr<ScalarProduct<X>>, std::shared_ptr<Preconditioner<X,Y>>, real_type, int, int,int)
\copydetails IterativeSolver::IterativeSolver(std::shared_ptr<const LinearOperator<X,Y>>, std::shared_ptr<const ScalarProduct<X>>, std::shared_ptr<Preconditioner<X,Y>>, real_type, int, int,int)
\param mmax is the maximal number of previous vectors which are orthogonalized against the new search direction.
*/
RestartedFCGSolver (std::shared_ptr<LinearOperator<X,X>> op,
std::shared_ptr<ScalarProduct<X>> sp,
RestartedFCGSolver (std::shared_ptr<const LinearOperator<X,X>> op,
std::shared_ptr<const ScalarProduct<X>> sp,
std::shared_ptr<Preconditioner<X,X>> prec,
scalar_real_type reduction, int maxit, int verbose,
int mmax = 10)
......@@ -1545,7 +1545,7 @@ private:
/*!
\brief Constructor.
\copydoc IterativeSolver::IterativeSolver(L&,S&,P&,const ParameterTree&)
\copydoc IterativeSolver::IterativeSolver(const L&, const S&,P&,const ParameterTree&)
Additional parameter:
ParameterTree Key | Meaning
......@@ -1554,14 +1554,14 @@ private:
See \ref ISTL_Factory for the ParameterTree layout and examples.
*/
RestartedFCGSolver (std::shared_ptr<LinearOperator<X,X>> op,
RestartedFCGSolver (std::shared_ptr<const LinearOperator<X,X>> op,
std::shared_ptr<Preconditioner<X,X>> prec,
const ParameterTree& config)
: IterativeSolver<X,X>(op, prec, config), _mmax(config.get("mmax", 10))
{}
RestartedFCGSolver (std::shared_ptr<LinearOperator<X,X>> op,
std::shared_ptr<ScalarProduct<X>> sp,
RestartedFCGSolver (std::shared_ptr<const LinearOperator<X,X>> op,
std::shared_ptr<const ScalarProduct<X>> sp,
std::shared_ptr<Preconditioner<X,X>> prec,
const ParameterTree& config)
: IterativeSolver<X,X>(op, sp, prec, config), _mmax(config.get("mmax", 10))
......
dune/istl/test/solvertest.cc
+ 1
− 1
View file @ 90236d45
......@@ -68,7 +68,7 @@ int main(int argc, char** argv)
typedef Dune::MatrixAdapter<BCRSMat,BVector,BVector> Operator;
BCRSMat mat;
Operator fop(mat);
const Operator fop(mat);
BVector b(N*N), x(N*N);
setupLaplacian(mat,N);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment