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

Use the field_type instead of hard-coded 'double'

[[Imported from SVN: r11761]]
parent e0b40e9f
Branches
Tags
No related merge requests found
......@@ -7,7 +7,7 @@ DiscFuncType TrustRegionGSStep<OperatorType, DiscFuncType>::getSol()
template<class OperatorType, class DiscFuncType>
inline
double TrustRegionGSStep<OperatorType, DiscFuncType>::
typename DiscFuncType::field_type TrustRegionGSStep<OperatorType, DiscFuncType>::
residual(unsigned int index, int blockIndex) const
{
const OperatorType& mat = *this->mat_;
......@@ -18,7 +18,7 @@ residual(unsigned int index, int blockIndex) const
/* The following loop computes
* \f[ sum_i = \sum_{j \ne i} A_{ij}x_j \f]
*/
double r = (*this->rhs_)[index][blockIndex];
field_type r = (*this->rhs_)[index][blockIndex];
ColumnIterator cIt = mat[index].template begin();
ColumnIterator cEndIt = mat[index].template end();
......@@ -62,23 +62,23 @@ void TrustRegionGSStep<OperatorType, DiscFuncType>::iterate()
continue;
// Compute residual
double r = residual(i, j);
const double& diag = mat[i][i][j][j];
field_type r = residual(i, j);
const field_type& diag = mat[i][i][j][j];
// Find local minimum
double& x = (*this->x_)[i][j];
field_type& x = (*this->x_)[i][j];
if (diag > 0) {
//printf("%d %d konvex\n", i, j);
//double oldEnergy = computeEnergy(mat, *this->x_, *this->rhs_);
// 1d problem is convex
#ifndef NDEBUG
double oldX = x;
field_type oldX = x;
#endif
x = r / diag;
#ifndef NDEBUG
double energyX = 0.5*diag*x*x - r*x;
double energyOldX = 0.5*diag*oldX*oldX - r*oldX;
field_type energyX = 0.5*diag*x*x - r*x;
field_type energyOldX = 0.5*diag*oldX*oldX - r*oldX;
assert(energyX <= energyOldX + 1e-6);
#endif
if ( x < obstacles[i].lower(j)) {
......@@ -104,11 +104,11 @@ void TrustRegionGSStep<OperatorType, DiscFuncType>::iterate()
//printf("%d %d konkav\n", i, j);
// 1d problem is concave or linear.
// Minimum is attained at one of the boundaries
double lBound = obstacles[i].lower(j);
double uBound = obstacles[i].upper(j);
field_type lBound = obstacles[i].lower(j);
field_type uBound = obstacles[i].upper(j);
double lValue = 0.5*diag*lBound*lBound - r*lBound;
double uValue = 0.5*diag*uBound*uBound - r*uBound;
field_type lValue = 0.5*diag*lBound*lBound - r*lBound;
field_type uValue = 0.5*diag*uBound*uBound - r*uBound;
#if 0
double max = 0.5*diag*(r/diag)*(r/diag) - r*r/diag;
......
......@@ -37,7 +37,7 @@
virtual DiscFuncType getSol();
double residual(unsigned int index, int blockIndex) const;
field_type residual(unsigned int index, int blockIndex) const;
};
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment