diff --git a/dune/solvers/iterationsteps/minimalpolynomialextrapolationstep.hh b/dune/solvers/iterationsteps/minimalpolynomialextrapolationstep.hh
index bc3af39da29e80a03b696c96d07693bb056478fa..d05acc6757e18cb36e7922cc07c8fb5625e3dc40 100644
--- a/dune/solvers/iterationsteps/minimalpolynomialextrapolationstep.hh
+++ b/dune/solvers/iterationsteps/minimalpolynomialextrapolationstep.hh
@@ -14,23 +14,23 @@
#include <dune/solvers/iterationsteps/iterationstep.hh>
/** \brief A single step of the Minimal Polynomial Extrapolation method
- *
+ *
* Minimal Polynomial Extrapolation (MPE) is a general method for sequence acceleration.
* It takes a converging sequence (actually, it even works for some nonconverging ones)
* and interpolates the iterates such as to obtain a second sequence which converges
* faster.
- *
+ *
* In this implementation, the original sequence is produced by an IterationStep object.
*/
template<class VectorType,
class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> >
- class MinimalPolynomialExtrapolationStep
+ class MinimalPolynomialExtrapolationStep
: public IterationStep<VectorType, BitVectorType>
{
public:
-
- /** \brief Constructor
+
+ /** \brief Constructor
* \param iterationStep The iteration step object that creates the sequence being accelerated
* \param restart Restart the method after this number of iterations
*/
@@ -46,19 +46,19 @@ public:
/** \brief Retrieve the solution */
virtual VectorType getSol();
- /** \brief To be called before starting to iterate
+ /** \brief To be called before starting to iterate
\note This calls the preprocess method for the dependant iteration step class, too!
*/
virtual void preprocess();
-
+
IterationStep<VectorType,BitVectorType>* iterationStep_;
/** \brief Restart the method after this number of iterations */
int restart_;
-
+
/** \brief The history of iterates of the original sequence */
std::vector<VectorType> xHistory_;
-
+
/** \brief The history of differences between the iterates of the original sequence */
std::vector<VectorType> U_;
@@ -76,10 +76,10 @@ inline
void MinimalPolynomialExtrapolationStep<VectorType, BitVectorType>::preprocess()
{
iterationStep_->preprocess();
-
+
xHistory_.resize(1);
xHistory_[0] = *this->x_;
-
+
U_.resize(0);
}
@@ -89,13 +89,13 @@ void MinimalPolynomialExtrapolationStep<VectorType, BitVectorType>::iterate()
{
// append a copy of the last element to the history, for the starting value
xHistory_.push_back(xHistory_.back());
-
+
// iterate once
iterationStep_->x_ = &xHistory_.back();
iterationStep_->iterate();
//std::cout << "x[" << xHistory_.size()-1 << "]:\n" << xHistory_.back() << std::endl;
-
+
// update the history of sequence differences
VectorType diff = xHistory_.back();
diff -= xHistory_[xHistory_.size()-2];
@@ -103,30 +103,30 @@ void MinimalPolynomialExtrapolationStep<VectorType, BitVectorType>::iterate()
// current number of iterates used for interpolation
const size_t k= U_.size()-1;
-
+
/** Compute weights c by solving the system
* \f[ U^T_{k-1} U_{k-1} c = -U_{k-1} u_k \f]
*/
typedef Dune::Matrix<Dune::FieldMatrix<double,1,1> > ScalarMatrixType;
typedef Dune::BlockVector<Dune::FieldVector<double,1> > ScalarVectorType;
-
+
// set up the matrix
ScalarMatrixType UTU(k,k);
for (size_t i=0; i<k; i++)
for (size_t j=0; j<k; j++)
UTU[i][j] = U_[i] * U_[j];
-
+
// set up the right hand side
ScalarVectorType rhs(k);
for (size_t i=0; i<k; i++)
rhs[i] = U_[i] * U_[k];
rhs *= -1;
-
+
// solve the system
ScalarVectorType c(k);
- c = 0;
-
+ c = 0;
+
// Technicality: turn the matrix into a linear operator
Dune::MatrixAdapter<ScalarMatrixType,ScalarVectorType,ScalarVectorType> op(UTU);
@@ -147,32 +147,32 @@ void MinimalPolynomialExtrapolationStep<VectorType, BitVectorType>::iterate()
// The last coefficient of the c array is '1'
c.resize(c.size()+1, true);
c[c.size()-1] = 1;
-
+
//std::cout << "c:\n" << c << std::endl;
-
+
/////////////////////////////////////////////////////////////
// Compute the new accelerated iterate
/////////////////////////////////////////////////////////////
VectorType& newIterate = *this->x_;
newIterate = 0;
double cSum = std::accumulate(c.begin(), c.end(), double(0));
-
+
for (size_t i=1; i<=k+1; i++)
newIterate.axpy(c[i-1]/cSum, xHistory_[i]);
-
+
//std::cout << "y:\n" << newIterate << std::endl;
-
+
// Acceleration methods for nonlinear problems should be restarted
// every now and then
if (k==restart_) {
-
+
std::cout << "Restarting MPE..." << std::endl;
U_.resize(0);
xHistory_.resize(1);
xHistory_[0] = *this->x_;
-
+
}
-
+
}
#endif
diff --git a/dune/solvers/iterationsteps/projectedlinegsstep.cc b/dune/solvers/iterationsteps/projectedlinegsstep.cc
index f0fe5b099ac4d2fd5340eb2adf4b6d0bf6ec533b..b852c1afdd17ac7fcb365aaccffffe7182be1ef2 100755
--- a/dune/solvers/iterationsteps/projectedlinegsstep.cc
+++ b/dune/solvers/iterationsteps/projectedlinegsstep.cc
@@ -22,7 +22,7 @@ solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
// /////////////////////////////
x = 0;
-
+
// //////////////////////////////////////
// The TNNMG loop
// //////////////////////////////////////
@@ -35,12 +35,12 @@ solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
for (int iteration=0; iteration<10; iteration++) {
VectorType oldX = x;
-
+
// ///////////////////////////
// Nonlinear presmoothing
// ///////////////////////////
VectorType rhsCopy = rhs; // need a non-const version of rhs
-
+
ProjectedBlockGSStep<LocalMatrixType, VectorType> nonlinearSmootherStep(matrix, x, rhsCopy);
nonlinearSmootherStep.ignoreNodes_ = &ignoreNodes;
@@ -49,11 +49,11 @@ solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
for (int i=0; i<3; i++)
nonlinearSmootherStep.iterate();
-
+
// Compute residual
VectorType residual = rhs;
matrix.mmv(x,residual);
-
+
// ///////////////////////////
// Truncation
// ///////////////////////////
@@ -70,10 +70,10 @@ solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
// left off-diagonal:
if (i>0)
truncatedMatrix[i][i-1] = 0;
-
+
// diagonal
truncatedMatrix[i][i] = Dune::ScaledIdentityMatrix<double,BlockSize>(1);
-
+
// right off-diagonal:
if (i<truncatedMatrix.N()-1)
truncatedMatrix[i][i+1] = 0;
@@ -84,11 +84,11 @@ solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
}
}
-
+
// ///////////////////////////
// Linear correction
// ///////////////////////////
-
+
VectorType linearCorrection(residual.size());
truncatedMatrix.solve(linearCorrection, residual); // the direct linear solution algorithm
@@ -96,53 +96,53 @@ solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
// Line search in the direction of the projected coarse
// grid correction to ensure monotonicity.
// //////////////////////////////////////////////////////////
-
+
// L2-projection of the correction onto the defect obstacle
for (size_t i=0; i<linearCorrection.size(); i++) {
-
+
for (int j=0; j<BlockSize; j++) {
-
+
linearCorrection[i][j] = std::max(linearCorrection[i][j], localObstacles[i].lower(j) - x[i][j]);
linearCorrection[i][j] = std::min(linearCorrection[i][j], localObstacles[i].upper(j) - x[i][j]);
-
+
}
-
+
}
-
+
// Construct obstacles in the direction of the projected correction
BoxConstraint<field_type,1> lineSearchObs;
for (size_t i=0; i<linearCorrection.size(); i++) {
-
+
for (int j=0; j<BlockSize; j++) {
-
+
if (linearCorrection[i][j] > 0) {
-
+
// This division can cause nan on some platforms...
if (!isnan( (localObstacles[i].lower(j)-x[i][j]) / linearCorrection[i][j]) )
lineSearchObs.lower(0) = std::max(lineSearchObs.lower(0),
(localObstacles[i].lower(j)-x[i][j]) / linearCorrection[i][j]);
-
+
if (!isnan( (localObstacles[i].upper(j)-x[i][j]) / linearCorrection[i][j]) )
- lineSearchObs.upper(0) = std::min(lineSearchObs.upper(0),
+ lineSearchObs.upper(0) = std::min(lineSearchObs.upper(0),
(localObstacles[i].upper(j)-x[i][j]) / linearCorrection[i][j]);
}
-
+
if (linearCorrection[i][j] < 0) {
if (!isnan( (localObstacles[i].upper(j)-x[i][j]) / linearCorrection[i][j]) )
- lineSearchObs.lower(0) = std::max(lineSearchObs.lower(0),
+ lineSearchObs.lower(0) = std::max(lineSearchObs.lower(0),
(localObstacles[i].upper(j)-x[i][j]) / linearCorrection[i][j]);
if (!isnan( (localObstacles[i].lower(j)-x[i][j]) / linearCorrection[i][j]) )
- lineSearchObs.upper(0) = std::min(lineSearchObs.upper(0),
+ lineSearchObs.upper(0) = std::min(lineSearchObs.upper(0),
(localObstacles[i].lower(j)-x[i][j]) / linearCorrection[i][j]);
}
-
+
}
-
+
}
// abort when the linear correction has a small norm
field_type correctionNormSquared = EnergyNorm<LocalMatrixType,VectorType>::normSquared(linearCorrection,truncatedMatrix);
-
+
// Line search
field_type alpha = (residual*linearCorrection) / correctionNormSquared;
alpha = std::max(std::min(alpha, lineSearchObs.upper(0)), lineSearchObs.lower(0));
@@ -187,7 +187,7 @@ void ProjectedLineGSStep<MatrixType, VectorType, BitVectorType>::iterate()
const int current_block_size = this->blockStructure_[b_num].size();
- //! compute and save the residuals for the curent block:
+ //! compute and save the residuals for the curent block:
// determine the (permuted) residuals r[p(i)],..., r[p(i+current_block_size-1)]
// this means that we determine the residuals for the current block
VectorType permuted_r_i(current_block_size);
@@ -206,7 +206,7 @@ void ProjectedLineGSStep<MatrixType, VectorType, BitVectorType>::iterate()
// (Later: x_now[p(i+k)] = x_now[p(i+k)] + v[p(i+k)] )
// //////////////////////////////////////////////////////////////////////////////////////////////////
- // Copy the linear system for the current line/block into a tridiagonal matrix
+ // Copy the linear system for the current line/block into a tridiagonal matrix
// //////////////////////////////////////////////////////////////////////////////////////////////////
Dune::BTDMatrix<typename MatrixType::block_type> tridiagonalMatrix(current_block_size);
@@ -218,10 +218,10 @@ void ProjectedLineGSStep<MatrixType, VectorType, BitVectorType>::iterate()
// left off-diagonal:
if (j>0)
tridiagonalMatrix[j][j-1] = 0;
-
+
// diagonal
tridiagonalMatrix[j][j] = Dune::ScaledIdentityMatrix<field_type,BlockSize>(1);
-
+
// left off-diagonal:
if (j<current_block_size-1)
tridiagonalMatrix[j][j+1] = 0;
@@ -231,10 +231,10 @@ void ProjectedLineGSStep<MatrixType, VectorType, BitVectorType>::iterate()
// left off-diagonal:
if (j>0)
tridiagonalMatrix[j][j-1] = matrix[this->blockStructure_[b_num][j]][this->blockStructure_[b_num][j-1]];
-
+
// diagonal
tridiagonalMatrix[j][j] = matrix[this->blockStructure_[b_num][j]][this->blockStructure_[b_num][j]];
-
+
// left off-diagonal:
if (j<current_block_size-1)
tridiagonalMatrix[j][j+1] = matrix[this->blockStructure_[b_num][j]][this->blockStructure_[b_num][j+1]];
diff --git a/dune/solvers/iterationsteps/projectedlinegsstep.hh b/dune/solvers/iterationsteps/projectedlinegsstep.hh
index a03aae099d14315184d6d0bf06febf98a23ded8f..e02ab0df30fee5bdc902ef19e52c64969bac22af 100755
--- a/dune/solvers/iterationsteps/projectedlinegsstep.hh
+++ b/dune/solvers/iterationsteps/projectedlinegsstep.hh
@@ -15,34 +15,34 @@ template<class MatrixType,
class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> >
class ProjectedLineGSStep : public LineGSStep<MatrixType, VectorType, BitVectorType>
{
-
+
typedef typename VectorType::block_type VectorBlock;
-
+
typedef typename VectorType::field_type field_type;
-
+
enum {BlockSize = VectorBlock::dimension};
-
+
public:
-
+
//! Default constructor. Doesn't init anything
ProjectedLineGSStep() {}
-
+
//! Constructor for usage of Permutation Manager
ProjectedLineGSStep( const std::vector<std::vector<unsigned int> >& blockStructure)
: LineGSStep<MatrixType, VectorType, BitVectorType>( blockStructure )
{}
-
+
/** \brief Solve one tridiagonal system */
void solveLocalSystem(const Dune::BTDMatrix<typename MatrixType::block_type>& matrix,
VectorType& x,
const VectorType& rhs,
const std::vector<BoxConstraint<field_type,BlockSize> >& localObstacles) const;
-
+
//! Perform one iteration
virtual void iterate();
-
+
const std::vector<BoxConstraint<field_type,BlockSize> >* obstacles_;
-
+
};