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 { 37 typename std::remove_const<T>::type sum; 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 121