diff --git a/prism/include/Measures.h b/prism/include/Measures.h index 8bb092eb..1eb923d1 100644 --- a/prism/include/Measures.h +++ b/prism/include/Measures.h @@ -103,4 +103,86 @@ public: }; +/** + * Measure for determining the difference between the upper and lower values in + * an interval iteration. + * In relative mode, the difference is scaled by the mid-point between upper and lower value. + */ +class MeasureSupNormInterval { +private: + // relative mode? + const bool relative; + // the current maximal value + double sup_norm; + +public: + /** Constructor, set relative flag */ + MeasureSupNormInterval(bool relative) : relative(relative) { + reset(); + } + + /** Reset for new measurement */ + void reset() { + sup_norm = 0.0; + } + + /** Relative mode? */ + bool isRelative() const { + return relative; + } + + /** + * Do the measurement for an upper and lower value pair. + * For relative mode, the lower value is used as the divisor. + */ + inline void measure(double lower, double upper) { + double x; + + // special case: one of the values is infinite (this can happen e.g. for non-converging + // iterations when the values grow extremely large and overflow to infinity). + if (isinf(lower) || isinf(upper)) { + x = std::numeric_limits::infinity(); + } else { + // compute difference + // we don't use fabs like for MeasureSupNorm, as we want x to become negative + // in situations where upper < lower (should only happen due to numerical inaccuracies / rounding) + x = upper - lower; + // we clamp to zero for negative values + if (x < 0) + x = 0; + + if (relative && x != 0.0) { + // for relative mode: divide by lower + // taking lower ensures that if the actual value should happen to be + // the lower value, that then the relative precision is satisfied. + // We take the absolute value of the lower to ensure that + // x does not flip signs. + // Note: if lower is 0.0, then x will become +inf, as x!=0 + x /= fabs(lower); + } + } + + // sup_norm = max { x, sup_norm } + if (x > sup_norm) { + sup_norm = x; + } + } + + /** + * Do the measurement for a pair of arrays (lower and upper bounds) of size n. + * For relative mode, the midpoints are used as the divisors. + */ + inline void measure(const double *lower, const double *upper, std::size_t n) { + for (std::size_t i = 0; i < n; i++) { + measure(lower[i], upper[i]); + } + } + + /** Return the measured value */ + double value() const { + return sup_norm; + } + +}; + #endif