Skip to content
Snippets Groups Projects
Commit 1e0a32bb authored by Jonathan Youett's avatar Jonathan Youett
Browse files

Remove trailing whitespaces

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