Skip to content
Snippets Groups Projects
Commit b4807eba authored by Oliver Sander's avatar Oliver Sander Committed by sander
Browse files

Replace hard-coded 'double' by 'field_type'

[[Imported from SVN: r12436]]
parent ac49582a
Branches
Tags
No related merge requests found
...@@ -53,10 +53,10 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() ...@@ -53,10 +53,10 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve()
std::cout << std::endl; std::cout << std::endl;
} }
double error = std::numeric_limits<double>::max(); field_type error = std::numeric_limits<field_type>::max();
double normOfOldCorrection = 0; field_type normOfOldCorrection = 0;
double totalConvRate = 1; field_type totalConvRate = 1;
this->maxTotalConvRate_ = 0; this->maxTotalConvRate_ = 0;
int convRateCounter = 0; int convRateCounter = 0;
...@@ -67,7 +67,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() ...@@ -67,7 +67,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve()
VectorType q(*x_); // a temporary vector VectorType q(*x_); // a temporary vector
// some local variables // some local variables
double rho,rholast,beta; field_type rho,rholast,beta;
// determine initial search direction // determine initial search direction
p = 0; // clear correction p = 0; // clear correction
...@@ -82,9 +82,9 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() ...@@ -82,9 +82,9 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve()
rholast = p*b; // orthogonalization rholast = p*b; // orthogonalization
double solutionNormSquared = 0; field_type solutionNormSquared = 0;
double solutionTimesCorrection = 0; field_type solutionTimesCorrection = 0;
double correctionNormSquared = rholast; field_type correctionNormSquared = rholast;
// Loop until desired tolerance or maximum number of iterations is reached // Loop until desired tolerance or maximum number of iterations is reached
for (i=0; i<this->maxIterations_ && (error>this->tolerance_ || std::isnan(error)); i++) { for (i=0; i<this->maxIterations_ && (error>this->tolerance_ || std::isnan(error)); i++) {
...@@ -147,7 +147,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() ...@@ -147,7 +147,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve()
// minimize in given search direction p // minimize in given search direction p
matrix_->mv(p,q); // q=Ap 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 // Compute the new iterate. If it is outside of the admissible set
...@@ -159,7 +159,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() ...@@ -159,7 +159,7 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve()
tentativeNewIterate.axpy(alpha,p); // update solution tentativeNewIterate.axpy(alpha,p); // update solution
// Compute length of this tentative new iterate // Compute length of this tentative new iterate
double tentativeLengthSquared = solutionNormSquared field_type tentativeLengthSquared = solutionNormSquared
+ 2*alpha*solutionTimesCorrection + 2*alpha*solutionTimesCorrection
+ alpha*alpha*correctionNormSquared; + alpha*alpha*correctionNormSquared;
...@@ -221,17 +221,17 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve() ...@@ -221,17 +221,17 @@ void TruncatedCGSolver<MatrixType, VectorType>::solve()
// ////////////////////////////////////////////////// // //////////////////////////////////////////////////
// Compute error // Compute error
// ////////////////////////////////////////////////// // //////////////////////////////////////////////////
double oldNorm = errorNorm_->operator()(oldSolution); field_type oldNorm = errorNorm_->operator()(oldSolution);
oldSolution -= *x_; oldSolution -= *x_;
double normOfCorrection = errorNorm_->operator()(oldSolution); field_type normOfCorrection = errorNorm_->operator()(oldSolution);
if (this->useRelativeError_) if (this->useRelativeError_)
error = normOfCorrection / oldNorm; error = normOfCorrection / oldNorm;
else else
error = normOfCorrection; error = normOfCorrection;
double convRate = normOfCorrection / normOfOldCorrection; field_type convRate = normOfCorrection / normOfOldCorrection;
normOfOldCorrection = normOfCorrection; normOfOldCorrection = normOfCorrection;
......
...@@ -86,7 +86,7 @@ public: ...@@ -86,7 +86,7 @@ public:
void setProblem(const MatrixType& matrix, void setProblem(const MatrixType& matrix,
VectorType* x, VectorType* x,
VectorType* rhs, VectorType* rhs,
double trustRegionRadius) field_type trustRegionRadius)
{ {
matrix_ = &matrix; matrix_ = &matrix;
x_ = x; x_ = x;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment