Scippy

SoPlex

Sequential object-oriented simPlex

stablesum.h
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (c) 1996-2023 Zuse Institute Berlin (ZIB) */
7 /* */
8 /* Licensed under the Apache License, Version 2.0 (the "License"); */
9 /* you may not use this file except in compliance with the License. */
10 /* You may obtain a copy of the License at */
11 /* */
12 /* http://www.apache.org/licenses/LICENSE-2.0 */
13 /* */
14 /* Unless required by applicable law or agreed to in writing, software */
15 /* distributed under the License is distributed on an "AS IS" BASIS, */
16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
17 /* See the License for the specific language governing permissions and */
18 /* limitations under the License. */
19 /* */
20 /* You should have received a copy of the Apache-2.0 license */
21 /* along with SoPlex; see the file LICENSE. If not email to soplex@zib.de. */
22 /* */
23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
24 
25 #ifndef _SOPLEX_STABLE_SUM_H_
26 #define _SOPLEX_STABLE_SUM_H_
27 #include <type_traits>
28 
29 // #define SOPLEX_CHECK_STABLESUM // double check the stable sum computation
30 
31 namespace soplex
32 {
33 
34 template <typename T>
35 class StableSum
36 {
38 
39 public:
40  StableSum() : sum(0) {}
41  StableSum(const T& init) : sum(init) {}
42 
43  void operator+=(const T& input)
44  {
45  sum += input;
46  }
47 
48  void operator-=(const T& input)
49  {
50  sum -= input;
51  }
52 
53  operator typename std::remove_const<T>::type() const
54  {
55  return sum;
56  }
57 
58 };
59 
60 template <>
61 class StableSum<double>
62 {
63  double sum = 0;
64  double c = 0;
65 #ifdef SOPLEX_CHECK_STABLESUM
66  double checksum = 0;
67 #endif
68 
69 public:
70  StableSum() = default;
71  StableSum(double init) : sum(init), c(0) {}
72 
73  void operator+=(double input)
74  {
75 #if defined(_MSC_VER) || defined(__INTEL_COMPILER)
76 #pragma float_control( precise, on )
77 #endif
78 #ifdef SOPLEX_CHECK_STABLESUM
79  checksum += input;
80 #endif
81  double t = sum + input;
82  double z = t - sum;
83  double y = (sum - (t - z)) + (input - z);
84  c += y;
85 
86  sum = t;
87  }
88 
89  void operator-=(double input)
90  {
91  (*this) += -input;
92  }
93 
94  operator double() const
95  {
96 #ifdef SOPLEX_CHECK_STABLESUM
97 
98  if(spxAbs(checksum - (sum + c)) >= 1e-6)
99  printf("stablesum viol: %g, rel: %g, checksum: %g\n", spxAbs(checksum - (sum + c)),
100  spxAbs(checksum - (sum + c)) / SOPLEX_MAX(1.0, SOPLEX_MAX(spxAbs(checksum), spxAbs(sum + c))),
101  checksum);
102 
103  assert(spxAbs(checksum - (sum + c)) < 1e-6);
104 #endif
105  return sum + c;
106  }
107 };
108 
109 /// Output operator.
110 template < class T >
111 std::ostream& operator<<(std::ostream& s, const StableSum<T>& sum)
112 {
113  s << static_cast<T>(sum);
114 
115  return s;
116 }
117 
118 }
119 
120 #endif
void operator+=(double input)
Definition: stablesum.h:73
void operator-=(double input)
Definition: stablesum.h:89
void operator+=(const T &input)
Definition: stablesum.h:43
int printf(const S &format_str, const Args &... args)
Definition: printf.h:666
type
Definition: core.h:687
Everything should be within this namespace.
std::remove_const< T >::type sum
Definition: stablesum.h:37
R spxAbs(R a)
Definition: spxdefines.h:393
void operator-=(const T &input)
Definition: stablesum.h:48
#define SOPLEX_MAX(x, y)
Definition: spxdefines.h:297
StableSum(const T &init)
Definition: stablesum.h:41