diff --git a/dune/solvers/solvers/tcgsolver.cc b/dune/solvers/solvers/tcgsolver.cc index c46f0b6c3e829d2b39392216d49c932e7c868ef6..787dfc0b661bb344bf0c26303ef1d8e5219a78fc 100644 --- a/dune/solvers/solvers/tcgsolver.cc +++ b/dune/solvers/solvers/tcgsolver.cc @@ -53,10 +53,10 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() std::cout << std::endl; } - double error = std::numeric_limits<double>::max(); + field_type error = std::numeric_limits<field_type>::max(); - double normOfOldCorrection = 0; - double totalConvRate = 1; + field_type normOfOldCorrection = 0; + field_type totalConvRate = 1; this->maxTotalConvRate_ = 0; int convRateCounter = 0; @@ -67,7 +67,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() VectorType q(*x_); // a temporary vector // some local variables - double rho,rholast,beta; + field_type rho,rholast,beta; // determine initial search direction p = 0; // clear correction @@ -82,9 +82,9 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() rholast = p*b; // orthogonalization - double solutionNormSquared = 0; - double solutionTimesCorrection = 0; - double correctionNormSquared = rholast; + 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++) { @@ -147,7 +147,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() // minimize in given search direction p matrix_->mv(p,q); // q=Ap - double alpha = rholast/(p*q); // scalar product + field_type alpha = rholast/(p*q); // scalar product // //////////////////////////////////////////////////////////////////////// // Compute the new iterate. If it is outside of the admissible set @@ -159,7 +159,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() tentativeNewIterate.axpy(alpha,p); // update solution // Compute length of this tentative new iterate - double tentativeLengthSquared = solutionNormSquared + field_type tentativeLengthSquared = solutionNormSquared + 2*alpha*solutionTimesCorrection + alpha*alpha*correctionNormSquared; @@ -221,17 +221,17 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() // ////////////////////////////////////////////////// // Compute error // ////////////////////////////////////////////////// - double oldNorm = errorNorm_->operator()(oldSolution); + field_type oldNorm = errorNorm_->operator()(oldSolution); oldSolution -= *x_; - double normOfCorrection = errorNorm_->operator()(oldSolution); + field_type normOfCorrection = errorNorm_->operator()(oldSolution); if (this->useRelativeError_) error = normOfCorrection / oldNorm; else error = normOfCorrection; - double convRate = normOfCorrection / normOfOldCorrection; + field_type convRate = normOfCorrection / normOfOldCorrection; normOfOldCorrection = normOfCorrection; diff --git a/dune/solvers/solvers/tcgsolver.hh b/dune/solvers/solvers/tcgsolver.hh index a050eb785039e0ae9bf9b52fd7fcae2b843e2033..da48e8e87f958a3f010701eea523bad24e6ca8b3 100644 --- a/dune/solvers/solvers/tcgsolver.hh +++ b/dune/solvers/solvers/tcgsolver.hh @@ -86,7 +86,7 @@ public: void setProblem(const MatrixType& matrix, VectorType* x, VectorType* rhs, - double trustRegionRadius) + field_type trustRegionRadius) { matrix_ = &matrix; x_ = x;