stablesum.h

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the class library                   */
/*       SoPlex --- the Sequential object-oriented simPlex.                  */
/*                                                                           */
/*    Copyright (C) 1996-2022 Konrad-Zuse-Zentrum                            */
/*                            fuer Informationstechnik Berlin                */
/*                                                                           */
/*  SoPlex is distributed under the terms of the ZIB Academic Licence.       */
/*                                                                           */
/*  You should have received a copy of the ZIB Academic License              */
/*  along with SoPlex; see the file COPYING. If not email to soplex@zib.de.  */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

#ifndef _SOPLEX_STABLE_SUM_H_
#define _SOPLEX_STABLE_SUM_H_
#include <type_traits>

// #define CHECK_STABLESUM  // double check the stable sum computation

namespace soplex
{

template <typename T>
class StableSum
{
   typename std::remove_const<T>::type sum;

public:
   StableSum() : sum(0) {}
   StableSum(const T& init) : sum(init) {}

   void operator+=(const T& input)
   {
      sum += input;
   }

   void operator-=(const T& input)
   {
      sum -= input;
   }

   operator typename std::remove_const<T>::type() const
   {
      return sum;
   }

};

template <>
class StableSum<double>
{
   double sum = 0;
   double c = 0;
#ifdef CHECK_STABLESUM
   double checksum = 0;
#endif

public:
   StableSum() = default;
   StableSum(double init) : sum(init), c(0) {}

   void operator+=(double input)
   {
#if defined(_MSC_VER) || defined(__INTEL_COMPILER)
#pragma float_control( precise, on )
#endif
#ifdef CHECK_STABLESUM
      checksum += input;
#endif
      double t = sum + input;
      double z = t - sum;
      double y = (sum - (t - z)) + (input - z);
      c += y;

      sum = t;
   }

   void operator-=(double input)
   {
      (*this) += -input;
   }

   operator double() const
   {
#ifdef CHECK_STABLESUM

      if(spxAbs(checksum - (sum + c)) >= 1e-6)
         printf("stablesum viol: %g, rel: %g, checksum: %g\n", spxAbs(checksum - (sum + c)),
                spxAbs(checksum - (sum + c)) / MAXIMUM(1.0, MAXIMUM(spxAbs(checksum), spxAbs(sum + c))), checksum);

      assert(spxAbs(checksum - (sum + c)) < 1e-6);
#endif
      return sum + c;
   }
};

/// Output operator.
template < class T >
std::ostream& operator<<(std::ostream& s, const StableSum<T>& sum)
{
   s << static_cast<T>(sum);

   return s;
}

}

#endif