ergo
Interval.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
38#ifndef MAT_INTERVAL
39#define MAT_INTERVAL
40#include <math.h>
41#include <iomanip>
42namespace mat {
43
44 /* FIXME represent interval as midpoint and length to get better precision */
45 template<typename Treal>
46 class Interval {
47 public:
48 explicit Interval(Treal low = 1, Treal upp = -1)
50 }
51 inline bool empty() const {return lowerBound > upperBound;}
52
53 static Interval intersect(Interval const & A, Interval const & B) {
54 if (A.empty())
55 return A;
56 else if (B.empty())
57 return B;
58 else
59 return Interval(A.lowerBound > B.lowerBound ?
60 A.lowerBound : B.lowerBound,
61 A.upperBound < B.upperBound ?
62 A.upperBound : B.upperBound);
63 }
64 inline void intersect(Interval const & other) {
65 if (this->empty())
66 return;
67 else if (other.empty()) {
68 *this = other;
69 return;
70 }
71 else {
72 this->lowerBound = this->lowerBound > other.lowerBound ?
73 this->lowerBound : other.lowerBound;
74 this->upperBound = this->upperBound < other.upperBound ?
75 this->upperBound : other.upperBound;
76 return;
77 }
78 }
79
80 inline void intersect_always_non_empty(Interval const & other) {
81 if (this->empty()) {
82 *this = other;
83 return;
84 }
85 if (other.empty())
86 return;
87
88 Interval intersection = intersect(*this, other);
89 if ( !intersection.empty() ) {
90 *this = intersection;
91 return;
92 }
93 // Ok, the intersection is actually empty
94 Treal tmp_val;
95 if ( this->lowerBound > other.upperBound )
96 tmp_val = (this->lowerBound + other.upperBound) / 2;
97 else if ( this->upperBound < other.lowerBound )
98 tmp_val = (this->upperBound + other.lowerBound) / 2;
99 else
100 assert(0); // This should never happen!
101 this->lowerBound = tmp_val;
102 this->upperBound = tmp_val;
103 return;
104 }
105
109 inline Treal length() const {
110 if (empty())
111 return 0.0;
112 else
113 return upperBound - lowerBound;
114 }
115 inline Treal midPoint() const {
116 assert(!empty());
117 return (upperBound + lowerBound) / 2;
118 }
119 inline bool cover(Treal const value) const {
120 if (empty())
121 return false;
122 else
123 return (value <= upperBound && value >= lowerBound);
124 }
125 inline bool overlap(Interval const & other) const {
126 Interval tmp = intersect(*this, other);
127 return !tmp.empty();
128 }
129
133 inline void increase(Treal const value) {
134 if (empty())
135 throw Failure("Interval<Treal>::increase(Treal const) : "
136 "Attempt to increase empty interval.");
137 lowerBound -= value;
138 upperBound += value;
139 }
140 inline void decrease(Treal const value) {
141 lowerBound += value;
142 upperBound -= value;
143 }
144 inline Treal low() const {return lowerBound;}
145 inline Treal upp() const {return upperBound;}
146
147#if 0
148 inline Interval<Treal> operator*(Interval<Treal> const & other) const {
149 assert(lowerBound >= 0);
150 assert(other.lowerBound >= 0);
152 upperBound * other.upperBound);
153 }
154#endif
155
156 inline Interval<Treal> operator*(Treal const & value) const {
157 if (value < 0)
158 return Interval<Treal>(upperBound * value, lowerBound * value);
159 else
160 return Interval<Treal>(lowerBound * value, upperBound * value);
161 }
162
163 /* Minus operator well defined? */
164 inline Interval<Treal> operator-(Interval<Treal> const & other) const {
165 return *this + (-1.0 * other);
166 // return Interval<Treal>(lowerBound - other.lowerBound,
167 // upperBound - other.upperBound);
168 }
169
170 inline Interval<Treal> operator+(Interval<Treal> const & other) const {
172 upperBound + other.upperBound);
173 }
174 inline Interval<Treal> operator/(Treal const & value) const {
175 if (value < 0)
176 return Interval<Treal>(upperBound / value, lowerBound / value);
177 else
178 return Interval<Treal>(lowerBound / value, upperBound / value);
179 }
180 inline Interval<Treal> operator-(Treal const & value) const {
181 return Interval<Treal>(lowerBound - value, upperBound - value);
182 }
183 inline Interval<Treal> operator+(Treal const & value) const {
184 return Interval<Treal>(lowerBound + value, upperBound + value);
185 }
186
187
188
189 void puriStep(int poly);
190 void invPuriStep(int poly);
191 // The two functions above really not needed now, just special
192 // case of functions below with alpha == 1
193 void puriStep(int poly, Treal alpha);
194 void invPuriStep(int poly, Treal alpha);
195 protected:
198 private:
199 };
200
201#if 0
202 /* Minus operator is not well defined for intervals */
203 template<typename Treal>
204 inline Interval<Treal> operator-(Treal const value,
205 Interval<Treal> const & other) {
206 return Interval<Treal>(value - other.lowerBound,
207 value - other.upperBound);
208 }
209 template<typename Treal>
210 inline Interval<Treal> operator+(Treal const value,
211 Interval<Treal> const & other) {
212 return Interval<Treal>(value + other.lowerBound,
213 value + other.upperBound);
214 }
215#endif
216
217#if 1
218 template<typename Treal>
220 if (other.empty())
221 return other;
222 else {
223 assert(other.low() >= 0);
225 }
226 }
227#endif
228
229 template<typename Treal>
231 if (upperBound > 2.0 || lowerBound < -1.0)
232 throw Failure("Interval<Treal>::puriStep(int) : It is assumed here "
233 "that the interval I is within [-1.0, 2.0]");
234 Interval<Treal> oneInt(-1.0,1.0);
235 Interval<Treal> zeroInt(0.0,2.0);
236 /* Elias note 2010-09-12:
237 Sometimes the polynomial 2*x-x*x makes a non-empty interval in [0,1]
238 become empty because of rounding errors. For some reason this problem
239 showed up when using Intel compiler version 11.1 but not when using
240 version 10.1. Many test cases failed on Kalkyl when using
241 the compiler "icpc (ICC) 11.1 20100414".
242 Anyway, we know that the function f(x) = 2*x-x*x is growing
243 in the interval [0,1] so we know a non-empty interval inside [0,1] should
244 always stay non-empty. To avoid assertion failures in purification,
245 the code below now forces the interval to be non-empty if it became
246 empty due to rounding errors. */
247 bool nonEmptyIntervalInZeroToOne = false;
248 if(upperBound > lowerBound && lowerBound >= 0 && upperBound <= 1)
249 nonEmptyIntervalInZeroToOne = true;
250 if (poly) {
251 /* 2*x - x*x */
252 *this = Interval<Treal>::intersect(*this,oneInt);
253 // assert(upperBound <= 1.0);
254 upperBound = 2 * upperBound - upperBound * upperBound;
255 lowerBound = 2 * lowerBound - lowerBound * lowerBound;
256 if(nonEmptyIntervalInZeroToOne && upperBound < lowerBound) {
257 // Force interval to be non-empty.
258 Treal midPoint = (lowerBound + upperBound) / 2;
259 upperBound = lowerBound = midPoint;
260 }
261 }
262 else { /* x*x */
263 *this = Interval<Treal>::intersect(*this,zeroInt);
264 // assert(lowerBound >= 0.0);
265 upperBound = upperBound * upperBound;
266 lowerBound = lowerBound * lowerBound;
267 }
268 }
269
270 template<typename Treal>
272 if (upperBound > 2.0 || lowerBound < -1.0)
273 throw Failure("Interval<Treal>::puriStep(int) : It is assumed here "
274 "that the interval I is within [-1.0, 1.0]");
275 Interval<Treal> oneInt(-1.0,1.0);
276 Interval<Treal> zeroInt(0.0,2.0);
277 if (poly) {
278 /* 1 - sqrt(1 - x) */
279 *this = Interval<Treal>::intersect(*this,oneInt);
280 // assert(upperBound <= 1.0);
281 upperBound = 1.0 - template_blas_sqrt(1.0 - upperBound);
282 lowerBound = 1.0 - template_blas_sqrt(1.0 - lowerBound);
283 }
284 else { /* sqrt(x) */
285 *this = Interval<Treal>::intersect(*this,zeroInt);
286 // assert(lowerBound >= 0.0);
287 upperBound = template_blas_sqrt(upperBound);
288 lowerBound = template_blas_sqrt(lowerBound);
289 }
290 }
291
292#if 1
293 template<typename Treal>
294 void Interval<Treal>::puriStep(int poly, Treal alpha) {
295 if (upperBound > 2.0 || lowerBound < -1.0)
296 throw Failure("Interval<Treal>::puriStep(int, real) : It is assumed here "
297 "that the interval I is within [-1.0, 2.0]");
298 Interval<Treal> oneInt(-1.0,1.0);
299 Interval<Treal> zeroInt(0.0,2.0);
300 /* Elias note 2010-09-12:
301 Sometimes the polynomial 2*x-x*x makes a non-empty interval in [0,1]
302 become empty because of rounding errors. For some reason this problem
303 showed up when using Intel compiler version 11.1 but not when using
304 version 10.1. Many test cases failed on Kalkyl when using
305 the compiler "icpc (ICC) 11.1 20100414".
306 Anyway, we know that the function f(x) = 2*x-x*x is growing
307 in the interval [0,1] so we know a non-empty interval inside [0,1] should
308 always stay non-empty. To avoid assertion failures in purification,
309 the code below now forces the interval to be non-empty if it became
310 empty due to rounding errors. */
311 bool nonEmptyIntervalInZeroToOne = false;
312 if(upperBound > lowerBound && lowerBound >= 0 && upperBound <= 1)
313 nonEmptyIntervalInZeroToOne = true;
314 if (poly) {
315 /* 2*alpha*x - alpha*alpha*x*x */
316 *this = Interval<Treal>::intersect(*this,oneInt);
317 // assert(upperBound <= 1.0);
318 upperBound = 2 * alpha * upperBound - alpha * alpha * upperBound * upperBound;
319 lowerBound = 2 * alpha * lowerBound - alpha * alpha * lowerBound * lowerBound;
320 if(nonEmptyIntervalInZeroToOne && upperBound < lowerBound) {
321 // Force interval to be non-empty.
322 Treal midPoint = (lowerBound + upperBound) / 2;
323 upperBound = lowerBound = midPoint;
324 }
325 }
326 else { /* ( alpha*x + (1-alpha) )^2 */
327 *this = Interval<Treal>::intersect(*this,zeroInt);
328 // assert(lowerBound >= 0.0);
329 upperBound = ( alpha * upperBound + (1-alpha) ) * ( alpha * upperBound + (1-alpha) );
330 lowerBound = ( alpha * lowerBound + (1-alpha) ) * ( alpha * lowerBound + (1-alpha) );
331 }
332 }
333
334 template<typename Treal>
335 void Interval<Treal>::invPuriStep(int poly, Treal alpha) {
336 if (upperBound > 2.0 || lowerBound < -1.0)
337 throw Failure("Interval<Treal>::invPuriStep(int, real) : It is assumed here "
338 "that the interval I is within [-1.0, 2.0]");
339 Interval<Treal> oneInt(-1.0,1.0);
340 Interval<Treal> zeroInt(0.0,2.0);
341 if (poly) {
342 /* ( 1 - sqrt(1 - x) ) / alpha */
343 *this = Interval<Treal>::intersect(*this,oneInt);
344 // assert(upperBound <= 1.0);
345 upperBound = ( 1.0 - template_blas_sqrt(1.0 - upperBound) ) / alpha;
346 lowerBound = ( 1.0 - template_blas_sqrt(1.0 - lowerBound) ) / alpha;
347 }
348 else { /* ( sqrt(x) - (1-alpha) ) / alpha */
349 *this = Interval<Treal>::intersect(*this,zeroInt);
350 // assert(lowerBound >= 0.0);
351 upperBound = ( template_blas_sqrt(upperBound) - (1-alpha) ) / alpha;
352 lowerBound = ( template_blas_sqrt(lowerBound) - (1-alpha) ) / alpha;
353 }
354 }
355
356#endif
357
358 template<typename Treal>
359 std::ostream& operator<<(std::ostream& s,
360 Interval<Treal> const & in) {
361 if (in.empty())
362 s<<"[empty]";
363 else
364 s<<"["<<in.low()<<", "<<in.upp()<<"]";
365 return s;
366 }
367
368} /* end namespace mat */
369#endif
Definition: Failure.h:57
Definition: Interval.h:46
void puriStep(int poly)
Definition: Interval.h:230
void puriStep(int poly, Treal alpha)
Definition: Interval.h:294
Treal low() const
Definition: Interval.h:144
bool overlap(Interval const &other) const
Definition: Interval.h:125
Interval< Treal > operator-(Treal const &value) const
Definition: Interval.h:180
void increase(Treal const value)
Increases interval with value in both directions.
Definition: Interval.h:133
Interval< Treal > operator/(Treal const &value) const
Definition: Interval.h:174
void intersect_always_non_empty(Interval const &other)
Definition: Interval.h:80
bool cover(Treal const value) const
Definition: Interval.h:119
Treal midPoint() const
Definition: Interval.h:115
Treal length() const
Returns the length of the interval.
Definition: Interval.h:109
static Interval intersect(Interval const &A, Interval const &B)
Definition: Interval.h:53
Treal upperBound
Definition: Interval.h:197
Treal upp() const
Definition: Interval.h:145
bool empty() const
Definition: Interval.h:51
Interval< Treal > operator+(Treal const &value) const
Definition: Interval.h:183
void invPuriStep(int poly)
Definition: Interval.h:271
void decrease(Treal const value)
Definition: Interval.h:140
void intersect(Interval const &other)
Definition: Interval.h:64
Treal lowerBound
Definition: Interval.h:196
Interval< Treal > operator*(Treal const &value) const
Definition: Interval.h:156
void invPuriStep(int poly, Treal alpha)
Definition: Interval.h:335
Interval< Treal > operator+(Interval< Treal > const &other) const
Definition: Interval.h:170
Interval(Treal low=1, Treal upp=-1)
Definition: Interval.h:48
Interval< Treal > operator-(Interval< Treal > const &other) const
Definition: Interval.h:164
#define B
#define A
Definition: allocate.cc:39
XmY< TX, TY > operator-(TX const &AA, TY const &BB)
Substraction of two objects of type TX and TY.
Definition: matrix_proxy.h:277
Interval< Treal > sqrtInt(Interval< Treal > const &other)
Definition: Interval.h:219
XYZpUV< TX, TY, TZ, TU, TV > operator+(XYZ< TX, TY, TZ > const &ABC, XY< TU, TV > const &DE)
Addition of two multiplication proxys XYZ and XY.
Definition: matrix_proxy.h:237
std::ostream & operator<<(std::ostream &s, Interval< Treal > const &in)
Definition: Interval.h:359
Treal template_blas_sqrt(Treal x)