ergo
Interval.h
Go to the documentation of this file.
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