diff --git a/dune/solvers/iterationsteps/multigridstep.cc b/dune/solvers/iterationsteps/multigridstep.cc
index df031e7e54781008316108bdaf3482fd2fb906ae..5a73ce1460a1011e9b0ae1ae088f64dcb1760d79 100644
--- a/dune/solvers/iterationsteps/multigridstep.cc
+++ b/dune/solvers/iterationsteps/multigridstep.cc
@@ -9,20 +9,20 @@
#include <dune/solvers/solvers/quadraticipopt.hh>
#endif
-template <class MatrixType, class VectorType, class BitVectorType>
-VectorType MultigridStep<MatrixType, VectorType, BitVectorType>::
-getSol()
-{
- return this->x_[numLevels_-1];
-}
-
-template<class MatrixType, class VectorType, class BitVectorType>
-inline
-const MatrixType* MultigridStep<MatrixType, VectorType, BitVectorType>::
-getMatrix()
-{
- return this->mat_.back().get();
-}
+//template <class MatrixType, class VectorType, class BitVectorType>
+//VectorType MultigridStep<MatrixType, VectorType, BitVectorType>::
+//getSol()
+//{
+// return *(this->x_);
+//}
+
+//template<class MatrixType, class VectorType, class BitVectorType>
+//inline
+//const MatrixType* MultigridStep<MatrixType, VectorType, BitVectorType>::
+//getMatrix()
+//{
+// return this->mat_;
+//}
template<class MatrixType, class VectorType, class BitVectorType>
inline
@@ -34,6 +34,23 @@ setMGType(int mu, int nu1, int nu2)
nu2_ = nu2;
}
+/** \brief Sets up the MultigridStep for computation
+ *
+ * Resizes all hierarchy containers to appropriate size
+ * Sets up:
+ * \li matrix hierarchy
+ * \li ignore nodes hierarchy
+ * \li the base solver
+ *
+ * \b NOTE: All data, i.e.
+ * \li transfer operators
+ * \li fine level matrix, solution vector and right hand side
+ * \li smoothers
+ * have to be set prior to the call to preprocess().
+ * If any is changed subsequently, preprocess() has to be called again before iterate().
+ * There is one exception. The right hand side may be changed by setRhs() without a subsequent call to preprocess().
+ * This is in order to be able to compute solutions to multiple rhs w/o having to reassemble the same matrix hierarchy time and again.
+ */
template<class MatrixType, class VectorType, class BitVectorType>
void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
{
@@ -41,36 +58,76 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
std::cout << "Warning: Preprocess has already been called without calling iterate afterwards!" << std::endl;
preprocessCalled = true;
+ numLevels_ = mgTransfer_.size() + 1;
+
this->level_ = this->numLevels_-1;
-
+
// //////////////////////////////////////////////////////////
// Check that transfer operators have been supplied
// //////////////////////////////////////////////////////////
- if (this->mgTransfer_.size()!=this->numLevels_-1)
- DUNE_THROW(SolverError, "You requested " << this->numLevels_ << " levels "
- << "but you supplied " << this->mgTransfer_.size() << " transfer operators!");
-
for (size_t i=0; i<this->mgTransfer_.size(); i++)
if (this->mgTransfer_[i] == NULL)
DUNE_THROW(SolverError, "You have not supplied a multigrid restriction operator "
"for level " << i);
+
+ // //////////////////////////////////////////////////////////
+ // Resize hierarchy containers to current numLevels_
+ // and set fine level data where needed.
+ // //////////////////////////////////////////////////////////
+ for (int i=0; i<int(ignoreNodesHierarchy_.size())-1; i++)
+ if (ignoreNodesHierarchy_[i])
+ delete(ignoreNodesHierarchy_[i]);
+
+ xHierarchy_.resize(numLevels_);
+ rhsHierarchy_.resize(numLevels_);
+ matrixHierarchy_.resize(numLevels_);
+ ignoreNodesHierarchy_.resize(numLevels_);
+
+ for (int i=0; i<int(matrixHierarchy_.size())-1; i++)
+ {
+ matrixHierarchy_[i].reset();
+ xHierarchy_[i].reset();
+ ignoreNodesHierarchy_[i] = NULL;
+ }
+
+ matrixHierarchy_.back() = this->mat_;
+ xHierarchy_.back() = Dune::stackobject_to_shared_ptr(*(this->x_));
+ rhsHierarchy_.back() = *(this->rhs_);
+
+ presmoother_.resize(numLevels_);
+ postsmoother_.resize(numLevels_);
+
+ typename SmootherCache::iterator end=levelWiseSmoothers_.end();
+ for (size_t i=0; i<numLevels_; ++i)
+ {
+ typename SmootherCache::iterator levelSmoother = levelWiseSmoothers_.find(i);
+ if (levelSmoother != end)
+ presmoother_[i] = postsmoother_[i] = levelSmoother->second;
+ else
+ {
+ presmoother_[i] = presmootherDefault_;
+ postsmoother_[i] = postsmootherDefault_;
+ }
+ }
+
// /////////////////////////////////////////////////////
// Assemble the complete hierarchy of matrices
// /////////////////////////////////////////////////////
for (int i=this->numLevels_-2; i>=0; i--)
{
- this->mat_[i] = Dune::shared_ptr<MatrixType>(new MatrixType);
+ this->matrixHierarchy_[i] = Dune::shared_ptr<MatrixType>(new MatrixType);
+ this->xHierarchy_[i] = Dune::shared_ptr<VectorType>(new VectorType);
// Compute which entries are present in the (sparse) coarse grid stiffness
// matrices.
- this->mgTransfer_[i]->galerkinRestrictSetOccupation(*this->mat_[i+1], *std::const_pointer_cast<MatrixType>(this->mat_[i]));
+ this->mgTransfer_[i]->galerkinRestrictSetOccupation(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
// setup matrix
- this->mgTransfer_[i]->galerkinRestrict(*this->mat_[i+1], *std::const_pointer_cast<MatrixType>(this->mat_[i]));
+ this->mgTransfer_[i]->galerkinRestrict(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
// Set solution vector sizes for the lower levels
- GenericVector::resize(this->x_[i],*this->mat_[i]);
+ GenericVector::resize(*(this->xHierarchy_[i]),*this->matrixHierarchy_[i]);
}
// /////////////////////////////////////////////////////
@@ -78,7 +135,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
// /////////////////////////////////////////////////////
for (size_t i=0; i<ignoreNodesHierarchy_.size()-1; i++)
delete(ignoreNodesHierarchy_[i]);
-
+
if (this->ignoreNodes_!=0)
{
ignoreNodesHierarchy_[this->numLevels_-1] = const_cast<BitVectorType*>(this->ignoreNodes_);
@@ -87,7 +144,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
ignoreNodesHierarchy_[i] = new BitVectorType();
this->mgTransfer_[i]->restrictToFathers(*ignoreNodesHierarchy_[i+1], *ignoreNodesHierarchy_[i]);
}
- } else
+ } else
DUNE_THROW(SolverError, "We need a set of nodes to ignore");
// /////////////////////////////////////////////
@@ -110,7 +167,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
typedef LinearIterationStep<MatrixType, VectorType> SmootherType;
assert(dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_));
- dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->setProblem(*(this->mat_[0]), this->x_[0], this->rhs_[0]);
+ dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]);
dynamic_cast<SmootherType*>(loopBaseSolver->iterationStep_)->ignoreNodes_ = ignoreNodesHierarchy_[0];
}
@@ -119,7 +176,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
QuadraticIPOptSolver<MatrixType,VectorType>* ipoptBaseSolver = dynamic_cast<QuadraticIPOptSolver<MatrixType,VectorType>*> (this->basesolver_);
- ipoptBaseSolver->setProblem(*(this->mat_[0]), this->x_[0], this->rhs_[0]);
+ ipoptBaseSolver->setProblem(*(this->matrixHierarchy_[0]), *this->xHierarchy_[0], this->rhsHierarchy_[0]);
}
#endif
else {
@@ -129,7 +186,6 @@ else {
}
-
template<class MatrixType, class VectorType, class BitVectorType>
void MultigridStep<MatrixType, VectorType, BitVectorType>::nestedIteration()
{
@@ -137,7 +193,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::nestedIteration()
{
iterate();
- mgTransfer_[level_]->prolong(x_[level_], x_[level_+1]);
+ mgTransfer_[level_]->prolong(*xHierarchy_[level_], *xHierarchy_[level_+1]);
}
iterate();
@@ -158,9 +214,9 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate()
int& level = this->level_;
// Define references just for ease of notation
- std::vector<Dune::shared_ptr<const MatrixType> > const &mat = this->mat_;
- std::vector<VectorType>& x = this->x_;
- std::vector<VectorType>& rhs = this->rhs_;
+ std::vector<Dune::shared_ptr<const MatrixType> > const &mat = this->matrixHierarchy_;
+ std::vector<Dune::shared_ptr<VectorType> >& x = this->xHierarchy_;
+ std::vector<VectorType>& rhs = this->rhsHierarchy_;
// Solve directly if we're looking at the coarse problem
if (level == 0) {
@@ -170,7 +226,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate()
// Presmoothing
- presmoother_[level]->setProblem(*(this->mat_[level]), x[level], rhs[level]);
+ presmoother_[level]->setProblem(*(this->matrixHierarchy_[level]), *x[level], rhs[level]);
presmoother_[level]->ignoreNodes_ = ignoreNodesHierarchy_[level];
for (int i=0; i<this->nu1_; i++)
@@ -182,7 +238,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate()
// compute residual
// fineResidual = rhs[level] - mat[level] * x[level];
VectorType fineResidual = rhs[level];
- mat[level]->mmv(x[level], fineResidual);
+ mat[level]->mmv(*x[level], fineResidual);
// restrict residual
this->mgTransfer_[level-1]->restrict(fineResidual, rhs[level-1]);
@@ -192,7 +248,7 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate()
GenericVector::truncate(rhs[level-1], *ignoreNodesHierarchy_[level-1]);
// Choose all zeros as the initial correction
- x[level-1] = 0;
+ *x[level-1] = 0;
// ///////////////////////////////////////
// Recursively solve the coarser system
@@ -206,14 +262,16 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::iterate()
// add correction to the presmoothed solution
VectorType tmp;
- this->mgTransfer_[level-1]->prolong(x[level-1], tmp);
- x[level] += tmp;
+ this->mgTransfer_[level-1]->prolong(*x[level-1], tmp);
+ *x[level] += tmp;
// Postsmoothing
- postsmoother_[level]->setProblem(*(mat[level]), x[level], rhs[level]);
+ postsmoother_[level]->setProblem(*(mat[level]), *x[level], rhs[level]);
postsmoother_[level]->ignoreNodes_ = ignoreNodesHierarchy_[level];
for (int i=0; i<this->nu2_; i++)
postsmoother_[level]->iterate();
+// *(this->x_) = xHierarchy_.back();
+
}
diff --git a/dune/solvers/iterationsteps/multigridstep.hh b/dune/solvers/iterationsteps/multigridstep.hh
index 71ccd26fc6b39af8f4acfce779388d040ee95e0a..d315df2c797e0387d16f6a94a0c2951a49e645c4 100644
--- a/dune/solvers/iterationsteps/multigridstep.hh
+++ b/dune/solvers/iterationsteps/multigridstep.hh
@@ -2,6 +2,7 @@
#define MULTIGRID_STEP_HH
#include <vector>
+#include <map>
#include <dune/common/bitsetvector.hh>
#include <dune/solvers/transferoperators/multigridtransfer.hh>
@@ -37,22 +38,41 @@
LinearIterationStep<MatrixType, VectorType>* preSmoother,
LinearIterationStep<MatrixType, VectorType>* postSmoother,
Solver* baseSolver,
- const BitVectorType* ignoreNodes) :
+ const BitVectorType* ignoreNodes)
+ DUNE_DEPRECATED_MSG("The number of levels is no longer set explicitely, but instead inferred from the number of transfer operators.\n Just erase the number of levels from the argument list in your constructor call.") :
LinearIterationStep<MatrixType, VectorType>(mat, x, rhs),
presmoother_(numLevels), postsmoother_(numLevels),
- mat_(numLevels),
+ matrixHierarchy_(numLevels),
ignoreNodesHierarchy_(numLevels),
- x_(numLevels),
- rhs_(numLevels),
+ xHierarchy_(numLevels),
+ rhsHierarchy_(numLevels),
preprocessCalled(false)
{
numLevels_ = numLevels;
level_ = numLevels-1;
- mat_[level_] = &mat;
- x_[level_] = x;
- rhs_[level_] = rhs;
+ mu_ = mu;
+ nu1_ = nu1;
+ nu2_ = nu2;
+
+ setSmoother(preSmoother,postSmoother);
+
+ basesolver_ = baseSolver;
+
+ this->ignoreNodes_ = ignoreNodes;
+ }
+ MultigridStep(const MatrixType& mat,
+ VectorType& x,
+ VectorType& rhs,
+ int mu, int nu1, int nu2,
+ LinearIterationStep<MatrixType, VectorType>* preSmoother,
+ LinearIterationStep<MatrixType, VectorType>* postSmoother,
+ Solver* baseSolver,
+ const BitVectorType* ignoreNodes) :
+ LinearIterationStep<MatrixType, VectorType>(mat, x, rhs),
+ preprocessCalled(false)
+ {
mu_ = mu;
nu1_ = nu1;
nu2_ = nu2;
@@ -64,29 +84,35 @@
this->ignoreNodes_ = ignoreNodes;
}
+
MultigridStep(const MatrixType& mat,
VectorType& x,
- VectorType& rhs, int numLevels) :
+ VectorType& rhs, int numLevels)
+ DUNE_DEPRECATED_MSG("The number of levels is no longer set explicitely, but instead inferred from the number of transfer operators.\n Just erase the number of levels from the argument list in your constructor call.") :
LinearIterationStep<MatrixType, VectorType>(mat, x, rhs),
presmoother_(numLevels),postsmoother_(numLevels),
basesolver_(0),
- mat_(numLevels),
+ matrixHierarchy_(numLevels),
ignoreNodesHierarchy_(numLevels),
- x_(numLevels),
- rhs_(numLevels),
+ xHierarchy_(numLevels),
+ rhsHierarchy_(numLevels),
preprocessCalled(false)
{
numLevels_ = numLevels;
level_ = numLevels-1;
-
- mat_[level_] = Dune::stackobject_to_shared_ptr(mat);
- x_[level_] = x;
- rhs_[level_] = rhs;
}
+ MultigridStep(const MatrixType& mat,
+ VectorType& x,
+ VectorType& rhs) :
+ LinearIterationStep<MatrixType, VectorType>(mat, x, rhs),
+ basesolver_(0),
+ preprocessCalled(false)
+ {}
- virtual ~MultigridStep() {
- for (size_t i=0; i<mat_.size()-1; i++)
+ virtual ~MultigridStep()
+ {
+ for (size_t i=0; i<matrixHierarchy_.size()-1; i++)
delete(ignoreNodesHierarchy_[i]);
}
@@ -94,12 +120,13 @@
VectorType& x,
VectorType& rhs,
int numLevels)
+ DUNE_DEPRECATED_MSG("Use setProblem(const MatrixType& mat,VectorType& x,const VectorType& rhs).\n The number of levels is no longer set explicitely, but instead inferred from the number of transfer operators.")
{
- setNumberOfLevels(numLevels);
setProblem(mat, x, rhs);
}
virtual void setNumberOfLevels(int numLevels)
+ DUNE_DEPRECATED_MSG("The number of levels is no longer set explicitely, but instead inferred from the number of transfer operators.\n The functionality of this function is now performed in preprocess().")
{
numLevels_ = numLevels;
@@ -107,31 +134,28 @@
if (ignoreNodesHierarchy_[i])
delete(ignoreNodesHierarchy_[i]);
- mat_.resize(numLevels);
+ matrixHierarchy_.resize(numLevels);
ignoreNodesHierarchy_.resize(numLevels);
- for (int i=0; i<int(mat_.size())-1; i++)
+ for (int i=0; i<int(matrixHierarchy_.size())-1; i++)
{
- mat_[i].reset();
+ matrixHierarchy_[i].reset();
ignoreNodesHierarchy_[i] = NULL;
}
presmoother_.resize(numLevels);
postsmoother_.resize(numLevels);
- x_.resize(numLevels);
- rhs_.resize(numLevels);
+ xHierarchy_.resize(numLevels);
+ rhsHierarchy_.resize(numLevels);
}
- virtual void setProblem(const MatrixType& mat,
- VectorType& x,
+ virtual void setProblem(const MatrixType& mat,
+ VectorType& x,
const VectorType& rhs)
{
level_ = numLevels_-1;
-
- mat_[level_] = Dune::stackobject_to_shared_ptr(mat);
- x_[level_] = x;
- rhs_[level_] = rhs;
+ LinearIterationStep<MatrixType, VectorType, BitVectorType>::setProblem(mat,x,rhs);
// Preprocess must be called again, to create the new matrix hierarchy
preprocessCalled = false;
@@ -139,7 +163,9 @@
void setRhs(const VectorType& rhs)
{
- rhs_[numLevels_-1] = rhs;
+ level_ = numLevels_-1;
+ this->rhs_ = &rhs;
+ rhsHierarchy_.back() = rhs;
}
template <class DerivedTransferHierarchy>
@@ -158,9 +184,9 @@
virtual void postprocess();
- virtual VectorType getSol();
+// virtual VectorType getSol();
- virtual const MatrixType* getMatrix();
+// virtual const MatrixType* getMatrix();
virtual int level() const {return level_;}
@@ -173,36 +199,41 @@
virtual void setMGType(int mu, int nu1, int nu2);
/** \brief Set the smoother iteration step */
- virtual void setSmoother(LinearIterationStep<MatrixType, VectorType>* smoother) {
- for (size_t i=0; i<presmoother_.size(); i++)
- presmoother_[i] = postsmoother_[i] = Dune::stackobject_to_shared_ptr(*smoother);
+ virtual void setSmoother(LinearIterationStep<MatrixType, VectorType>* smoother)
+ {
+ presmootherDefault_ = postsmootherDefault_ = Dune::stackobject_to_shared_ptr(*smoother);
+
+ levelWiseSmoothers_.clear();
}
/** \brief Set the smoother iteration step from a smart pointer*/
- virtual void setSmoother(Dune::shared_ptr<LinearIterationStep<MatrixType, VectorType> > smoother) {
- for (size_t i=0; i<presmoother_.size(); i++)
- presmoother_[i] = postsmoother_[i] = smoother;
+ virtual void setSmoother(Dune::shared_ptr<LinearIterationStep<MatrixType, VectorType> > smoother)
+ {
+ presmootherDefault_ = postsmootherDefault_ = smoother;
+
+ levelWiseSmoothers_.clear();
}
/** \brief Set pre- and post smoothers individually */
virtual void setSmoother(LinearIterationStep<MatrixType, VectorType>* preSmoother,
- LinearIterationStep<MatrixType, VectorType>* postSmoother) {
- for (size_t i=0; i<presmoother_.size(); i++) {
- presmoother_[i] = Dune::stackobject_to_shared_ptr(*preSmoother);
- postsmoother_[i] = Dune::stackobject_to_shared_ptr(*postSmoother);
- }
+ LinearIterationStep<MatrixType, VectorType>* postSmoother)
+ {
+ presmootherDefault_ = Dune::stackobject_to_shared_ptr(*preSmoother);
+ postsmootherDefault_ = Dune::stackobject_to_shared_ptr(*postSmoother);
+
+ levelWiseSmoothers_.clear();
}
/** \brief Set the smoother iteration step for a particular level */
- virtual void setSmoother(LinearIterationStep<MatrixType, VectorType>* smoother,
- int level) {
- presmoother_[level] = postsmoother_[level] = Dune::stackobject_to_shared_ptr(*smoother);
+ virtual void setSmoother(LinearIterationStep<MatrixType, VectorType>* smoother, std::size_t level)
+ {
+ levelWiseSmoothers_[level] = Dune::stackobject_to_shared_ptr(*smoother);
}
/** \brief Set the smoother iteration step for a particular level, from a smart pointer */
- virtual void setSmoother(Dune::shared_ptr<LinearIterationStep<MatrixType, VectorType> > smoother,
- int level) {
- presmoother_[level] = postsmoother_[level] = smoother;
+ virtual void setSmoother(Dune::shared_ptr<LinearIterationStep<MatrixType, VectorType> > smoother, std::size_t level)
+ {
+ levelWiseSmoothers_[level] = smoother;
}
protected:
@@ -232,23 +263,29 @@
//! The total number of levels of the underlying grid
int numLevels_;
- public:
//! The linear operators on each level
- std::vector<Dune::shared_ptr<const MatrixType> > mat_;
+ std::vector<Dune::shared_ptr<const MatrixType> > matrixHierarchy_;
protected:
//! Flags specifying the dirichlet nodes on each level
std::vector<BitVectorType*> ignoreNodesHierarchy_;
public:
- std::vector<VectorType> x_;
+ std::vector<Dune::shared_ptr<VectorType> > xHierarchy_;
- std::vector<VectorType> rhs_;
+ std::vector<VectorType> rhsHierarchy_;
//protected:
std::vector<MultigridTransfer<VectorType, BitVectorType, MatrixType>* > mgTransfer_;
protected:
+
+ Dune::shared_ptr<LinearIterationStep<MatrixType, VectorType> > presmootherDefault_;
+ Dune::shared_ptr<LinearIterationStep<MatrixType, VectorType> > postsmootherDefault_;
+ typedef std::map<std::size_t, Dune::shared_ptr<LinearIterationStep<MatrixType, VectorType> > > SmootherCache;
+ SmootherCache levelWiseSmoothers_;
+
+
bool preprocessCalled;
};