From 52fc747dd49744ab1e02a9ae56ac1f8de41984d3 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 6 Apr 2009 14:36:12 +0000
Subject: [PATCH] moved to new module dune-solvers
[[Imported from SVN: r2351]]
---
dune-solvers/solvers/cgsolver.cc | 165 +++++++++++++++++++++++++++++++
dune-solvers/solvers/cgsolver.hh | 61 ++++++++++++
2 files changed, 226 insertions(+)
create mode 100644 dune-solvers/solvers/cgsolver.cc
create mode 100644 dune-solvers/solvers/cgsolver.hh
diff --git a/dune-solvers/solvers/cgsolver.cc b/dune-solvers/solvers/cgsolver.cc
new file mode 100644
index 00000000..7a5565b9
--- /dev/null
+++ b/dune-solvers/solvers/cgsolver.cc
@@ -0,0 +1,165 @@
+#include <cmath>
+#include <limits>
+#include <iomanip>
+
+template <class MatrixType, class VectorType>
+void CGSolver<MatrixType, VectorType>::check() const
+{
+ if (!preconditioner_)
+ DUNE_THROW(SolverError, "You need to supply a preconditioner to a CG solver!");
+
+ preconditioner_->check();
+
+ if (!errorNorm_)
+ DUNE_THROW(SolverError, "You need to supply a norm to a CG solver!");
+
+}
+
+template <class MatrixType, class VectorType>
+void CGSolver<MatrixType, VectorType>::solve()
+{
+
+ int i;
+
+ VectorType& b = *rhs_;
+
+ // Check whether the solver is set up properly
+ check();
+
+ // set up preconditioner
+ preconditioner_->preprocess();
+
+ if (this->verbosity_ != NumProc::QUIET)
+ std::cout << "--- CGSolver ---\n";
+
+ if (this->verbosity_ == NumProc::FULL)
+ {
+ std::cout << " iter";
+ std::cout << " error";
+ std::cout << " rate";
+ std::string header = preconditioner_->getOutput();
+ std::cout << header;
+ std::cout << std::endl;
+ std::cout << "-----";
+ std::cout << "---------------";
+ std::cout << "---------";
+ for(int i=0; i<header.size(); ++i)
+ std::cout << "-";
+ std::cout << std::endl;
+ }
+
+ double error = std::numeric_limits<double>::max();
+
+ double normOfOldCorrection = 0;
+ double totalConvRate = 1;
+ maxTotalConvRate_ = 0;
+ int convRateCounter = 0;
+
+ // overwrite b with defect
+ matrix_->mmv(*x_,b);
+
+ VectorType p(*x_); // the search direction
+ VectorType q(*x_); // a temporary vector
+
+ // some local variables
+ double rho,rholast,lambda,alpha,beta;
+
+ // determine initial search direction
+ p = 0; // clear correction
+
+ // apply preconditioner
+ preconditioner_->setProblem(*matrix_,p,b);
+ preconditioner_->iterate();
+ p = preconditioner_->getSol();
+ rholast = p*b; // orthogonalization
+
+ // 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
+ // minimize in given search direction p
+ matrix_->mv(p,q); // q=Ap
+ alpha = p*q; // scalar product
+ lambda = rholast/alpha; // minimization
+ x_->axpy(lambda,p); // update solution
+ rhs_->axpy(-lambda,q); // update defect
+
+ // determine new search direction
+ q = 0; // clear correction
+
+ // apply preconditioner
+ preconditioner_->setProblem(*matrix_,q,b);
+ preconditioner_->iterate();
+ q = preconditioner_->getSol();
+
+ rho = q*b; // orthogonalization
+ beta = rho/rholast; // scaling factor
+ p *= beta; // scale old search direction
+ p += q; // orthogonalization with correction
+ rholast = rho; // remember rho for recurrence
+
+ // write iteration to file, if requested
+ if (this->historyBuffer_!="") {
+ VectorType intermediateSol = preconditioner_->getSol();
+
+ std::stringstream iSolFilename;
+ iSolFilename << this->historyBuffer_ << "/intermediatesolution_" << std::setw(4) << std::setfill('0') << i;
+
+ std::ofstream file(iSolFilename.str().c_str(), std::ios::out|std::ios::binary);
+ if (not(file))
+ DUNE_THROW(SolverError, "Couldn't open file " << iSolFilename << " for writing");
+
+ GenericVector::writeBinary(file, intermediateSol);
+
+ file.close();
+ }
+
+ // Compute error
+ double oldNorm = errorNorm_->operator()(oldSolution);
+ oldSolution -= *x_;
+
+ double normOfCorrection = errorNorm_->operator()(oldSolution);
+
+ if (this->useRelativeError_)
+ error = normOfCorrection / oldNorm;
+ else
+ error = normOfCorrection;
+
+ double convRate = normOfCorrection / normOfOldCorrection;
+
+ normOfOldCorrection = normOfCorrection;
+
+ if (!isinf(convRate) && !isnan(convRate)) {
+ totalConvRate *= convRate;
+ maxTotalConvRate_ = std::max(maxTotalConvRate_, std::pow(totalConvRate, 1/((double)convRateCounter+1)));
+ convRateCounter++;
+ }
+
+ // Output
+ if (this->verbosity_ == NumProc::FULL) {
+ std::cout << std::setw(5) << i;
+
+ std::cout << std::setiosflags(std::ios::scientific);
+ std::cout << std::setw(15) << std::setprecision(7) << error;
+ std::cout << std::resetiosflags(std::ios::scientific);
+
+ std::cout << std::setiosflags(std::ios::fixed);
+ std::cout << std::setw(9) << std::setprecision(5) << convRate;
+ std::cout << std::resetiosflags(std::ios::fixed);
+
+ std::cout << std::endl;
+ }
+
+ }
+
+
+ if (this->verbosity_ != NumProc::QUIET) {
+ std::cout << "maxTotalConvRate: " << maxTotalConvRate_ << ", "
+ << i << " iterations performed\n";
+ std::cout << "--------------------\n";
+ }
+
+}
diff --git a/dune-solvers/solvers/cgsolver.hh b/dune-solvers/solvers/cgsolver.hh
new file mode 100644
index 00000000..096a5920
--- /dev/null
+++ b/dune-solvers/solvers/cgsolver.hh
@@ -0,0 +1,61 @@
+#ifndef CG_SOLVER_HH
+#define CG_SOLVER_HH
+
+#include <dune/ag-common/iterativesolver.hh>
+#include <dune/ag-common/lineariterationstep.hh>
+#include <dune/ag-common/norms/norm.hh>
+
+/** \brief A conjugate gradient solver
+ *
+ */
+template <class MatrixType, class VectorType>
+class CGSolver : public IterativeSolver<VectorType>
+{
+ static const int blocksize = VectorType::block_type::size;
+
+public:
+
+ /** \brief Constructor taking all relevant data */
+
+ CGSolver(const MatrixType* matrix,
+ VectorType* x,
+ VectorType* rhs,
+ LinearIterationStep<MatrixType,VectorType>* preconditioner,
+ int maxIterations,
+ double tolerance,
+ Norm<VectorType>* errorNorm,
+ NumProc::VerbosityMode verbosity,
+ bool useRelativeError=true)
+ : IterativeSolver<VectorType>(tolerance,maxIterations,verbosity,useRelativeError),
+ matrix_(matrix), x_(x), rhs_(rhs),
+ preconditioner_(preconditioner), errorNorm_(errorNorm)
+ {}
+
+ /** \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
+ */
+ virtual void check() const;
+
+ const MatrixType* matrix_;
+
+ VectorType* x_;
+
+ VectorType* rhs_;
+
+ //! The iteration step used by the algorithm
+ LinearIterationStep<MatrixType,VectorType>* preconditioner_;
+
+ //! The norm used to measure convergence
+ Norm<VectorType>* errorNorm_;
+
+ double maxTotalConvRate_;
+
+};
+
+#include "cgsolver.cc"
+
+#endif
--
GitLab