diff --git a/dune/solvers/iterationsteps/amgstep.hh b/dune/solvers/iterationsteps/amgstep.hh index 58f948d54fe327f46a80a514d10a95533a068046..a0f37001b9b4eebdcd53209e1754025c96395c18 100644 --- a/dune/solvers/iterationsteps/amgstep.hh +++ b/dune/solvers/iterationsteps/amgstep.hh @@ -3,7 +3,7 @@ #ifndef ISTL_AMG_STEP_HH #define ISTL_AMG_STEP_HH -/** \file +/** \file \brief A wrapper class for the ISTL AMG implementation */ @@ -28,7 +28,7 @@ class AMGStep /** \brief Use a sequential SSOR for smoothing */ typedef Dune::SeqSSOR<MatrixType,VectorType,VectorType> Smoother; typedef typename Dune::Amg::template SmootherTraits<Smoother>::Arguments SmootherArgs; - + typedef Dune::Amg::AMG<Operator,VectorType,Smoother> AMG; public: @@ -40,8 +40,8 @@ public: \param smootherArgs Arguments for the smoother. See the dune-istl documentation for details \param coarseningCriterion Arguments for the coarsening. See the dune-istl documentation for details */ - AMGStep (const MatrixType* stiffnessMatrix, - VectorType& x, + AMGStep (const MatrixType* stiffnessMatrix, + VectorType& x, VectorType& rhs, const SmootherArgs& smootherArgs, const Criterion& coarseningCriterion) @@ -54,8 +54,8 @@ public: } /** \brief Initialize the iteration step but don't actually build the matrix hierarchy yet */ - AMGStep (const MatrixType* stiffnessMatrix, - VectorType& x, + AMGStep (const MatrixType* stiffnessMatrix, + VectorType& x, VectorType& rhs) { setProblem(*stiffnessMatrix, x, rhs); diff --git a/dune/solvers/iterationsteps/linegsstep.cc b/dune/solvers/iterationsteps/linegsstep.cc index aada568c86d0156375adb1d7c32349f7b39b13a0..3beca0dfc6a643a33f2681b823c2303940a9e906 100755 --- a/dune/solvers/iterationsteps/linegsstep.cc +++ b/dune/solvers/iterationsteps/linegsstep.cc @@ -23,7 +23,7 @@ void LineGSStep<MatrixType, DiscFuncType, BitVectorType >::iterate() const int current_block_size = 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 DiscFuncType permuted_r_i(current_block_size); @@ -42,7 +42,7 @@ void LineGSStep<MatrixType, DiscFuncType, 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); @@ -54,10 +54,10 @@ void LineGSStep<MatrixType, DiscFuncType, BitVectorType >::iterate() // left off-diagonal: if (j>0) tridiagonalMatrix[j][j-1] = 0; - + // diagonal tridiagonalMatrix[j][j] = Dune::ScaledIdentityMatrix<double,BlockSize>(1); - + // right off-diagonal: if (j<current_block_size-1) tridiagonalMatrix[j][j+1] = 0; @@ -67,10 +67,10 @@ void LineGSStep<MatrixType, DiscFuncType, BitVectorType >::iterate() // left off-diagonal: if (j>0) tridiagonalMatrix[j][j-1] = mat[blockStructure_[b_num][j]][blockStructure_[b_num][j-1]]; - + // diagonal tridiagonalMatrix[j][j] = mat[blockStructure_[b_num][j]][blockStructure_[b_num][j]]; - + // right off-diagonal: if (j<current_block_size-1) tridiagonalMatrix[j][j+1] = mat[blockStructure_[b_num][j]][blockStructure_[b_num][j+1]]; diff --git a/dune/solvers/iterationsteps/linegsstep.hh b/dune/solvers/iterationsteps/linegsstep.hh index 1c59163995009e35d3b2b9703d3108de4eca76f7..a3d23bc8a34f306c39f469e5a4d0b0a3e15ce5ae 100755 --- a/dune/solvers/iterationsteps/linegsstep.hh +++ b/dune/solvers/iterationsteps/linegsstep.hh @@ -23,7 +23,7 @@ template<class MatrixType, std::vector<std::vector<unsigned int> > blockStructure_; public: - + //! Default constructor. Doesn't init anything LineGSStep() {} @@ -36,7 +36,7 @@ template<class MatrixType, LineGSStep(const MatrixType& mat, DiscFuncType& x, DiscFuncType& rhs) : Base(mat, x, rhs) {} - + //! Perform one iteration virtual void iterate(); diff --git a/dune/solvers/solvers/iterativesolver.cc b/dune/solvers/solvers/iterativesolver.cc index 419b6a1444f29311579515cf2a5d78890fa9917b..47a31e95139a94130a5a0b56dff802b66691afde 100644 --- a/dune/solvers/solvers/iterativesolver.cc +++ b/dune/solvers/solvers/iterativesolver.cc @@ -26,14 +26,14 @@ void IterativeSolver<VectorType, BitVectorType>::writeIterate(const VectorType& int iterationNumber) const { std::stringstream iSolFilename; - iSolFilename << historyBuffer_ << "/intermediatesolution_" + iSolFilename << historyBuffer_ << "/intermediatesolution_" << std::setw(4) << std::setfill('0') << iterationNumber; - + std::ofstream file(iSolFilename.str().c_str(), std::ios::out|std::ios::binary); if (not(file)) DUNE_THROW(SolverError, "Couldn't open file " << iSolFilename.str() << " for writing"); - + GenericVector::writeBinary(file, iterate); - + file.close(); } diff --git a/dune/solvers/solvers/tcgsolver.cc b/dune/solvers/solvers/tcgsolver.cc index 566cdfe8cbb2a8269c0810fdc4c201dee0158bf1..e6b63438452e3d86e5004a2777758aad93b49532 100644 --- a/dune/solvers/solvers/tcgsolver.cc +++ b/dune/solvers/solvers/tcgsolver.cc @@ -65,13 +65,13 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() VectorType p(*x_); // the search direction VectorType q(*x_); // a temporary vector - + // some local variables field_type rho,rholast,beta; - + // determine initial search direction p = 0; // clear correction - + // apply preconditioner if (preconditioner_) { preconditioner_->setProblem(*matrix_,p,b); @@ -85,13 +85,13 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() field_type solutionNormSquared = 0; field_type solutionTimesCorrection = 0; field_type correctionNormSquared = rholast; - + // Loop until desired tolerance or maximum number of iterations is reached for (i=0; i<this->maxIterations_ && (error>this->tolerance_ || std::isnan(error)); i++) { - + // Backup of the current solution for the error computation later on VectorType oldSolution = *x_; - + // /////////////////////////////////////////////////////////////////////// // Perform one iteration step. The current search direction is in p // /////////////////////////////////////////////////////////////////////// @@ -120,7 +120,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() } if (correctionNormSquared <= 0) { - + if (this->historyBuffer_!="") this->writeIterate(*x_, i); return; @@ -131,7 +131,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() if (energyNormSquared <= 0) { field_type tau = positiveRoot(correctionNormSquared, - 2*solutionTimesCorrection, + 2*solutionTimesCorrection, solutionNormSquared - trustRegionRadius_*trustRegionRadius_); // add scaled correction to current iterate @@ -160,7 +160,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() // Compute length of this tentative new iterate field_type tentativeLengthSquared = solutionNormSquared - + 2*alpha*solutionTimesCorrection + + 2*alpha*solutionTimesCorrection + alpha*alpha*correctionNormSquared; // std::cout << "tentative radius^2: " << trustRegionScalarProduct(tentativeNewIterate,tentativeNewIterate) @@ -169,12 +169,12 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() if (tentativeLengthSquared >= trustRegionRadius_*trustRegionRadius_) { - field_type tau = positiveRoot(correctionNormSquared, - 2*solutionTimesCorrection, + field_type tau = positiveRoot(correctionNormSquared, + 2*solutionTimesCorrection, solutionNormSquared - trustRegionRadius_*trustRegionRadius_); - - - + + + //*x_ = eta_j; // return eta + tau*delta_j x_->axpy(tau, p); @@ -192,7 +192,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() // determine new search direction q = 0; // clear correction - + // apply preconditioner if (preconditioner_) { preconditioner_->setProblem(*matrix_,q,b); @@ -200,7 +200,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() q = preconditioner_->getSol(); } else q = b; - + rho = q*b; // orthogonalization beta = rho/rholast; // scaling factor p *= beta; // scale old search direction @@ -215,11 +215,11 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() // ////////////////////////////////////////////////// // write iteration to file, if requested // ////////////////////////////////////////////////// - if (this->historyBuffer_!="") + if (this->historyBuffer_!="") this->writeIterate(*x_, i); // ////////////////////////////////////////////////// - // Compute error + // Compute error // ////////////////////////////////////////////////// field_type oldNorm = (errorNorm_) ? errorNorm_->operator()(oldSolution) : oldSolution.two_norm(); oldSolution -= *x_; diff --git a/dune/solvers/solvers/tcgsolver.hh b/dune/solvers/solvers/tcgsolver.hh index 703d5db78122417840fafffab0793f9a30265a44..82f649a4df0a370dddf8669dfcf3799292df3ef3 100644 --- a/dune/solvers/solvers/tcgsolver.hh +++ b/dune/solvers/solvers/tcgsolver.hh @@ -11,7 +11,7 @@ #include <dune/solvers/iterationsteps/lineariterationstep.hh> #include <dune/solvers/norms/norm.hh> -/** \brief A truncated conjugate gradient solver +/** \brief A truncated conjugate gradient solver * */ template <class MatrixType, class VectorType> @@ -43,12 +43,12 @@ class TruncatedCGSolver : public IterativeSolver<VectorType> if (trustRegionNormMatrix_) return Dune::MatrixVector::Axy(*trustRegionNormMatrix_, b, a); - + return a*b; } -public: +public: /** \brief Constructor taking all relevant data */ TruncatedCGSolver(const MatrixType* matrix, @@ -64,7 +64,7 @@ public: bool useRelativeError=true) : IterativeSolver<VectorType>(tolerance,maxIterations,verbosity,useRelativeError), matrix_(matrix), x_(x), rhs_(rhs), - preconditioner_(preconditioner), + preconditioner_(preconditioner), errorNorm_(errorNorm), trustRegionRadius_(trustRegionRadius), trustRegionNormMatrix_(trustRegionNormMatrix) @@ -80,12 +80,12 @@ public: bool useRelativeError=true) : IterativeSolver<VectorType>(tolerance,maxIterations,verbosity,useRelativeError), matrix_(NULL), x_(NULL), rhs_(NULL), - preconditioner_(preconditioner), + preconditioner_(preconditioner), errorNorm_(errorNorm), trustRegionRadius_(0), trustRegionNormMatrix_(trustRegionNormMatrix) {} - + void setProblem(const MatrixType& matrix, VectorType* x, VectorType* rhs, @@ -96,11 +96,11 @@ public: rhs_ = rhs; trustRegionRadius_ = trustRegionRadius; } - + /** \brief Loop, call the iteration procedure * and monitor convergence */ virtual void solve(); - + /** \brief Checks whether all relevant member variables are set * \exception SolverError if the iteration step is not set up properly */ @@ -111,10 +111,10 @@ public: VectorType* x_; VectorType* rhs_; - + //! The iteration step used by the algorithm LinearIterationStep<MatrixType,VectorType>* preconditioner_; - + //! The norm used to measure convergence Norm<VectorType>* errorNorm_;