Skip to content
Snippets Groups Projects
Commit 1212aa08 authored by Oliver Sander's avatar Oliver Sander
Browse files

Remove old debugging code

parent 565ab669
Branches
No related tags found
No related merge requests found
...@@ -37,9 +37,7 @@ residual(unsigned int index, int blockIndex) const ...@@ -37,9 +37,7 @@ residual(unsigned int index, int blockIndex) const
} }
//std::cout << "+++ index " << index << " rhs: " << (*this->rhs_)[index] << std::endl;
return r; return r;
} }
template<class MatrixType, class VectorType> template<class MatrixType, class VectorType>
...@@ -49,10 +47,6 @@ void TrustRegionGSStep<MatrixType, VectorType>::iterate() ...@@ -49,10 +47,6 @@ void TrustRegionGSStep<MatrixType, VectorType>::iterate()
const MatrixType& mat = *this->mat_; const MatrixType& mat = *this->mat_;
const std::vector<BoxConstraint<field_type,blocksize> >& obstacles = *this->obstacles_; const std::vector<BoxConstraint<field_type,blocksize> >& obstacles = *this->obstacles_;
// Debug: check for energy decrease
//std::cout << "----------------------------------\n";
//double oldEnergy = computeEnergy(mat, *this->x_, *this->rhs_);
for (size_t i=0; i<this->x_->size(); i++) { for (size_t i=0; i<this->x_->size(); i++) {
// scalar gauss-seidel // scalar gauss-seidel
...@@ -69,39 +63,17 @@ void TrustRegionGSStep<MatrixType, VectorType>::iterate() ...@@ -69,39 +63,17 @@ void TrustRegionGSStep<MatrixType, VectorType>::iterate()
field_type& x = (*this->x_)[i][j]; field_type& x = (*this->x_)[i][j];
if (diag > 0) { if (diag > 0) {
//printf("%d %d konvex\n", i, j);
//double oldEnergy = computeEnergy(mat, *this->x_, *this->rhs_);
// 1d problem is convex // 1d problem is convex
#ifndef NDEBUG
field_type oldX = x;
#endif
x = r / diag; x = r / diag;
#ifndef NDEBUG
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)) { if ( x < obstacles[i].lower(j)) {
// assert that projection increases functional value
// double trueMin = 0.5*diag*x*x - r*x;
// double proMin = 0.5*diag*obstacles[i].lower(j)*obstacles[i].lower(j) - r*obstacles[i].lower(j);
//assert(trueMin < proMin);
x = obstacles[i].lower(j); x = obstacles[i].lower(j);
} }
else if ( x > obstacles[i].upper(j)) { else if ( x > obstacles[i].upper(j)) {
// assert that projection increases functional value
// double trueMin = 0.5*diag*x*x - r*x;
// double proMin = 0.5*diag*obstacles[i].upper(j)*obstacles[i].upper(j) - r*obstacles[i].upper(j);
//assert(trueMin < proMin);
x = obstacles[i].upper(j); x = obstacles[i].upper(j);
} }
// Debug: check for energy decrease
//double energy = computeEnergy(mat, *this->x_, *this->rhs_);
//assert(oldEnergy+1e-6>=energy);
} else { } else {
//printf("%d %d konkav\n", i, j);
// 1d problem is concave or linear. // 1d problem is concave or linear.
// Minimum is attained at one of the boundaries // Minimum is attained at one of the boundaries
field_type lBound = obstacles[i].lower(j); field_type lBound = obstacles[i].lower(j);
...@@ -110,25 +82,12 @@ void TrustRegionGSStep<MatrixType, VectorType>::iterate() ...@@ -110,25 +82,12 @@ void TrustRegionGSStep<MatrixType, VectorType>::iterate()
field_type lValue = 0.5*diag*lBound*lBound - r*lBound; field_type lValue = 0.5*diag*lBound*lBound - r*lBound;
field_type uValue = 0.5*diag*uBound*uBound - r*uBound; field_type uValue = 0.5*diag*uBound*uBound - r*uBound;
#if 0
double max = 0.5*diag*(r/diag)*(r/diag) - r*r/diag;
assert(max >= lValue);
assert(max >= uValue);
#endif
x = (lValue < uValue) ? lBound : uBound; x = (lValue < uValue) ? lBound : uBound;
// Debug: check for energy decrease
//double energy = computeEnergy(mat, *this->x_, *this->rhs_);
//assert(oldEnergy+1e-6>=energy);
} }
} }
} }
// Debug: check for energy decrease
//double energy = computeEnergy(mat, *this->x_, *this->rhs_);
//assert(oldEnergy+1e-6>=energy);
} }
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment