ergo
|
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek. 00004 * 00005 * This program is free software: you can redistribute it and/or modify 00006 * it under the terms of the GNU General Public License as published by 00007 * the Free Software Foundation, either version 3 of the License, or 00008 * (at your option) any later version. 00009 * 00010 * This program is distributed in the hope that it will be useful, 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 * GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License 00016 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * Primary academic reference: 00019 * KohnâSham Density Functional Theory Electronic Structure Calculations 00020 * with Linearly Scaling Computational Time and Memory Usage, 00021 * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek, 00022 * J. Chem. Theory Comput. 7, 340 (2011), 00023 * <http://dx.doi.org/10.1021/ct100611z> 00024 * 00025 * For further information about Ergo, see <http://www.ergoscf.org>. 00026 */ 00027 00036 #ifndef MAT_INTERVAL 00037 #define MAT_INTERVAL 00038 #include <math.h> 00039 #include <iomanip> 00040 namespace mat { 00041 00042 /* FIXME represent interval as midpoint and length to get better precision */ 00043 template<typename Treal> 00044 class Interval { 00045 public: 00046 explicit Interval(Treal low = 1, Treal upp = -1) 00047 : lowerBound(low), upperBound(upp) { 00048 } 00049 inline bool empty() const {return lowerBound > upperBound;} 00050 00051 static Interval intersect(Interval const & A, Interval const & B) { 00052 if (A.empty()) 00053 return A; 00054 else if (B.empty()) 00055 return B; 00056 else 00057 return Interval(A.lowerBound > B.lowerBound ? 00058 A.lowerBound : B.lowerBound, 00059 A.upperBound < B.upperBound ? 00060 A.upperBound : B.upperBound); 00061 } 00062 inline void intersect(Interval const & other) { 00063 if (this->empty()) 00064 return; 00065 else if (other.empty()) { 00066 *this = other; 00067 return; 00068 } 00069 else { 00070 this->lowerBound = this->lowerBound > other.lowerBound ? 00071 this->lowerBound : other.lowerBound; 00072 this->upperBound = this->upperBound < other.upperBound ? 00073 this->upperBound : other.upperBound; 00074 return; 00075 } 00076 } 00077 00078 inline void intersect_always_non_empty(Interval const & other) { 00079 if (this->empty()) { 00080 *this = other; 00081 return; 00082 } 00083 if (other.empty()) 00084 return; 00085 00086 Interval intersection = intersect(*this, other); 00087 if ( !intersection.empty() ) { 00088 *this = intersection; 00089 return; 00090 } 00091 // Ok, the intersection is actually empty 00092 Treal tmp_val; 00093 if ( this->lowerBound > other.upperBound ) 00094 tmp_val = (this->lowerBound + other.upperBound) / 2; 00095 else if ( this->upperBound < other.lowerBound ) 00096 tmp_val = (this->upperBound + other.lowerBound) / 2; 00097 else 00098 assert(0); // This should never happen! 00099 this->lowerBound = tmp_val; 00100 this->upperBound = tmp_val; 00101 return; 00102 } 00103 00107 inline Treal length() const { 00108 if (empty()) 00109 return 0.0; 00110 else 00111 return upperBound - lowerBound; 00112 } 00113 inline Treal midPoint() const { 00114 assert(!empty()); 00115 return (upperBound + lowerBound) / 2; 00116 } 00117 inline bool cover(Treal const value) const { 00118 if (empty()) 00119 return false; 00120 else 00121 return (value <= upperBound && value >= lowerBound); 00122 } 00123 inline bool overlap(Interval const & other) const { 00124 Interval tmp = intersect(*this, other); 00125 return !tmp.empty(); 00126 } 00127 00131 inline void increase(Treal const value) { 00132 if (empty()) 00133 throw Failure("Interval<Treal>::increase(Treal const) : " 00134 "Attempt to increase empty interval."); 00135 lowerBound -= value; 00136 upperBound += value; 00137 } 00138 inline void decrease(Treal const value) { 00139 lowerBound += value; 00140 upperBound -= value; 00141 } 00142 inline Treal low() const {return lowerBound;} 00143 inline Treal upp() const {return upperBound;} 00144 00145 #if 0 00146 inline Interval<Treal> operator*(Interval<Treal> const & other) const { 00147 assert(lowerBound >= 0); 00148 assert(other.lowerBound >= 0); 00149 return Interval<Treal>(lowerBound * other.lowerBound, 00150 upperBound * other.upperBound); 00151 } 00152 #endif 00153 00154 inline Interval<Treal> operator*(Treal const & value) const { 00155 if (value < 0) 00156 return Interval<Treal>(upperBound * value, lowerBound * value); 00157 else 00158 return Interval<Treal>(lowerBound * value, upperBound * value); 00159 } 00160 00161 /* Minus operator well defined? */ 00162 inline Interval<Treal> operator-(Interval<Treal> const & other) const { 00163 return *this + (-1.0 * other); 00164 // return Interval<Treal>(lowerBound - other.lowerBound, 00165 // upperBound - other.upperBound); 00166 } 00167 00168 inline Interval<Treal> operator+(Interval<Treal> const & other) const { 00169 return Interval<Treal>(lowerBound + other.lowerBound, 00170 upperBound + other.upperBound); 00171 } 00172 inline Interval<Treal> operator/(Treal const & value) const { 00173 if (value < 0) 00174 return Interval<Treal>(upperBound / value, lowerBound / value); 00175 else 00176 return Interval<Treal>(lowerBound / value, upperBound / value); 00177 } 00178 inline Interval<Treal> operator-(Treal const & value) const { 00179 return Interval<Treal>(lowerBound - value, upperBound - value); 00180 } 00181 inline Interval<Treal> operator+(Treal const & value) const { 00182 return Interval<Treal>(lowerBound + value, upperBound + value); 00183 } 00184 00185 00186 00187 void puriStep(int poly); 00188 void invPuriStep(int poly); 00189 // The two functions above really not needed now, just special 00190 // case of functions below with alpha == 1 00191 void puriStep(int poly, Treal alpha); 00192 void invPuriStep(int poly, Treal alpha); 00193 protected: 00194 Treal lowerBound; 00195 Treal upperBound; 00196 private: 00197 }; 00198 00199 #if 0 00200 /* Minus operator is not well defined for intervals */ 00201 template<typename Treal> 00202 inline Interval<Treal> operator-(Treal const value, 00203 Interval<Treal> const & other) { 00204 return Interval<Treal>(value - other.lowerBound, 00205 value - other.upperBound); 00206 } 00207 template<typename Treal> 00208 inline Interval<Treal> operator+(Treal const value, 00209 Interval<Treal> const & other) { 00210 return Interval<Treal>(value + other.lowerBound, 00211 value + other.upperBound); 00212 } 00213 #endif 00214 00215 #if 1 00216 template<typename Treal> 00217 Interval<Treal> sqrtInt(Interval<Treal> const & other) { 00218 if (other.empty()) 00219 return other; 00220 else { 00221 assert(other.low() >= 0); 00222 return Interval<Treal>(template_blas_sqrt(other.low()), template_blas_sqrt(other.upp())); 00223 } 00224 } 00225 #endif 00226 00227 template<typename Treal> 00228 void Interval<Treal>::puriStep(int poly) { 00229 if (upperBound > 2.0 || lowerBound < -1.0) 00230 throw Failure("Interval<Treal>::puriStep(int) : It is assumed here " 00231 "that the interval I is within [-1.0, 2.0]"); 00232 Interval<Treal> oneInt(-1.0,1.0); 00233 Interval<Treal> zeroInt(0.0,2.0); 00234 /* Elias note 2010-09-12: 00235 Sometimes the polynomial 2*x-x*x makes a non-empty interval in [0,1] 00236 become empty because of rounding errors. For some reason this problem 00237 showed up when using Intel compiler version 11.1 but not when using 00238 version 10.1. Many test cases failed on Kalkyl when using 00239 the compiler "icpc (ICC) 11.1 20100414". 00240 Anyway, we know that the function f(x) = 2*x-x*x is growing 00241 in the interval [0,1] so we know a non-empty interval inside [0,1] should 00242 always stay non-empty. To avoid assertion failures in purification, 00243 the code below now forces the interval to be non-empty if it became 00244 empty due to rounding errors. */ 00245 bool nonEmptyIntervalInZeroToOne = false; 00246 if(upperBound > lowerBound && lowerBound >= 0 && upperBound <= 1) 00247 nonEmptyIntervalInZeroToOne = true; 00248 if (poly) { 00249 /* 2*x - x*x */ 00250 *this = Interval<Treal>::intersect(*this,oneInt); 00251 // assert(upperBound <= 1.0); 00252 upperBound = 2 * upperBound - upperBound * upperBound; 00253 lowerBound = 2 * lowerBound - lowerBound * lowerBound; 00254 if(nonEmptyIntervalInZeroToOne && upperBound < lowerBound) { 00255 // Force interval to be non-empty. 00256 Treal midPoint = (lowerBound + upperBound) / 2; 00257 upperBound = lowerBound = midPoint; 00258 } 00259 } 00260 else { /* x*x */ 00261 *this = Interval<Treal>::intersect(*this,zeroInt); 00262 // assert(lowerBound >= 0.0); 00263 upperBound = upperBound * upperBound; 00264 lowerBound = lowerBound * lowerBound; 00265 } 00266 } 00267 00268 template<typename Treal> 00269 void Interval<Treal>::invPuriStep(int poly) { 00270 if (upperBound > 2.0 || lowerBound < -1.0) 00271 throw Failure("Interval<Treal>::puriStep(int) : It is assumed here " 00272 "that the interval I is within [-1.0, 1.0]"); 00273 Interval<Treal> oneInt(-1.0,1.0); 00274 Interval<Treal> zeroInt(0.0,2.0); 00275 if (poly) { 00276 /* 1 - sqrt(1 - x) */ 00277 *this = Interval<Treal>::intersect(*this,oneInt); 00278 // assert(upperBound <= 1.0); 00279 upperBound = 1.0 - template_blas_sqrt(1.0 - upperBound); 00280 lowerBound = 1.0 - template_blas_sqrt(1.0 - lowerBound); 00281 } 00282 else { /* sqrt(x) */ 00283 *this = Interval<Treal>::intersect(*this,zeroInt); 00284 // assert(lowerBound >= 0.0); 00285 upperBound = template_blas_sqrt(upperBound); 00286 lowerBound = template_blas_sqrt(lowerBound); 00287 } 00288 } 00289 00290 #if 1 00291 template<typename Treal> 00292 void Interval<Treal>::puriStep(int poly, Treal alpha) { 00293 if (upperBound > 2.0 || lowerBound < -1.0) 00294 throw Failure("Interval<Treal>::puriStep(int, real) : It is assumed here " 00295 "that the interval I is within [-1.0, 2.0]"); 00296 Interval<Treal> oneInt(-1.0,1.0); 00297 Interval<Treal> zeroInt(0.0,2.0); 00298 /* Elias note 2010-09-12: 00299 Sometimes the polynomial 2*x-x*x makes a non-empty interval in [0,1] 00300 become empty because of rounding errors. For some reason this problem 00301 showed up when using Intel compiler version 11.1 but not when using 00302 version 10.1. Many test cases failed on Kalkyl when using 00303 the compiler "icpc (ICC) 11.1 20100414". 00304 Anyway, we know that the function f(x) = 2*x-x*x is growing 00305 in the interval [0,1] so we know a non-empty interval inside [0,1] should 00306 always stay non-empty. To avoid assertion failures in purification, 00307 the code below now forces the interval to be non-empty if it became 00308 empty due to rounding errors. */ 00309 bool nonEmptyIntervalInZeroToOne = false; 00310 if(upperBound > lowerBound && lowerBound >= 0 && upperBound <= 1) 00311 nonEmptyIntervalInZeroToOne = true; 00312 if (poly) { 00313 /* 2*alpha*x - alpha*alpha*x*x */ 00314 *this = Interval<Treal>::intersect(*this,oneInt); 00315 // assert(upperBound <= 1.0); 00316 upperBound = 2 * alpha * upperBound - alpha * alpha * upperBound * upperBound; 00317 lowerBound = 2 * alpha * lowerBound - alpha * alpha * lowerBound * lowerBound; 00318 if(nonEmptyIntervalInZeroToOne && upperBound < lowerBound) { 00319 // Force interval to be non-empty. 00320 Treal midPoint = (lowerBound + upperBound) / 2; 00321 upperBound = lowerBound = midPoint; 00322 } 00323 } 00324 else { /* ( alpha*x + (1-alpha) )^2 */ 00325 *this = Interval<Treal>::intersect(*this,zeroInt); 00326 // assert(lowerBound >= 0.0); 00327 upperBound = ( alpha * upperBound + (1-alpha) ) * ( alpha * upperBound + (1-alpha) ); 00328 lowerBound = ( alpha * lowerBound + (1-alpha) ) * ( alpha * lowerBound + (1-alpha) ); 00329 } 00330 } 00331 00332 template<typename Treal> 00333 void Interval<Treal>::invPuriStep(int poly, Treal alpha) { 00334 if (upperBound > 2.0 || lowerBound < -1.0) 00335 throw Failure("Interval<Treal>::invPuriStep(int, real) : It is assumed here " 00336 "that the interval I is within [-1.0, 2.0]"); 00337 Interval<Treal> oneInt(-1.0,1.0); 00338 Interval<Treal> zeroInt(0.0,2.0); 00339 if (poly) { 00340 /* ( 1 - sqrt(1 - x) ) / alpha */ 00341 *this = Interval<Treal>::intersect(*this,oneInt); 00342 // assert(upperBound <= 1.0); 00343 upperBound = ( 1.0 - template_blas_sqrt(1.0 - upperBound) ) / alpha; 00344 lowerBound = ( 1.0 - template_blas_sqrt(1.0 - lowerBound) ) / alpha; 00345 } 00346 else { /* ( sqrt(x) - (1-alpha) ) / alpha */ 00347 *this = Interval<Treal>::intersect(*this,zeroInt); 00348 // assert(lowerBound >= 0.0); 00349 upperBound = ( template_blas_sqrt(upperBound) - (1-alpha) ) / alpha; 00350 lowerBound = ( template_blas_sqrt(lowerBound) - (1-alpha) ) / alpha; 00351 } 00352 } 00353 00354 #endif 00355 00356 template<typename Treal> 00357 std::ostream& operator<<(std::ostream& s, 00358 Interval<Treal> const & in) { 00359 if (in.empty()) 00360 s<<"[empty]"; 00361 else 00362 s<<"["<<in.low()<<", "<<in.upp()<<"]"; 00363 return s; 00364 } 00365 00366 } /* end namespace mat */ 00367 #endif