ergo
truncation.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 
00041 #ifndef MAT_TRUNCATION
00042 #define MAT_TRUNCATION
00043 #include <limits>
00044 #include <stdexcept>
00045 #include <cmath>
00046 namespace mat { /* Matrix namespace */
00047 
00048   // Stuff for Euclidean norm based truncation
00049   
00050   template<typename Tmatrix, typename Treal>
00051     class EuclTruncationBase {
00052   public:
00053     explicit EuclTruncationBase( Tmatrix & A_ );
00054     Treal run(Treal const threshold);
00055     virtual ~EuclTruncationBase() {}
00056   protected:
00057     virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc, 
00058                                      Treal const threshold ) = 0;
00059     virtual void getFrobSqNorms( std::vector<Treal> & frobsq_norms ) = 0;
00060     virtual void frobThreshLowLevel( Treal const threshold ) = 0;
00061     virtual Interval<Treal> euclIfSmall( Treal const absTol,
00062                                          Treal const threshold ) = 0;
00063     Tmatrix & A; // Matrix to be truncated
00064     Tmatrix E; // Error matrix    
00065   };
00066 
00067   template<typename Tmatrix, typename Treal>
00068     EuclTruncationBase<Tmatrix, Treal>::EuclTruncationBase( Tmatrix & A_ ) 
00069     : A(A_) {
00070     SizesAndBlocks rows;
00071     SizesAndBlocks cols;
00072     A.getRows(rows);
00073     A.getCols(cols);
00074     E.resetSizesAndBlocks(rows, cols);      
00075   }
00076 
00077   template<typename Tmatrix, typename Treal>
00078     Treal EuclTruncationBase<Tmatrix, Treal>::run( Treal const threshold ) {
00079     assert(threshold >= (Treal)0.0);
00080     if (threshold == (Treal)0.0)
00081       return (Treal)0;
00082     std::vector<Treal> frobsq_norms;
00083     this->getFrobSqNorms( frobsq_norms ); /*=======*/
00084     std::sort(frobsq_norms.begin(), frobsq_norms.end());
00085     int low = -1;
00086     int high = frobsq_norms.size() - 1; 
00087     Treal lowFrobTrunc, highFrobTrunc;
00088     this->getFrobTruncBounds( lowFrobTrunc, highFrobTrunc, threshold ); /*=======*/
00089     Treal frobsqSum = 0;
00090     while( low < (int)frobsq_norms.size() - 1 && frobsqSum < lowFrobTrunc ) {
00091       ++low;
00092       frobsqSum += frobsq_norms[low];
00093     }
00094     high = low; /* Removing all tom high is to much */
00095     --low;
00096     while( high < (int)frobsq_norms.size() - 1 && frobsqSum < highFrobTrunc ) {
00097       ++high;
00098       frobsqSum += frobsq_norms[high];
00099     }
00100     // Now we have low and high
00101     int minStep   = int( 0.01 * frobsq_norms.size() ); // Consider elements in chunks of at least 1 percent of all elements at a time to not get too many iterations
00102     minStep = minStep > 0 ? minStep : 1; // step is at least one
00103     int testIndex = high;
00104     int previousTestIndex = high * 2;
00105     // Now, removing everything up to and including testIndex is too much
00106     Interval<Treal> euclEInt(0, threshold * 2);
00107     // We go from above (too many elements in the error matrix) and stop as soon as the error matrix is small enough
00108     while ( euclEInt.upp() > threshold ) {
00109       // Removing everything up to and including testIndex is too much, update high:
00110       high = testIndex;
00111       int stepSize = (int)((high - low) * 0.01); // We can accept that only 99% of elements possible to remove are removed
00112       // stepSize must be at least minStep:
00113       stepSize = stepSize >= minStep ? stepSize : minStep; 
00114       previousTestIndex = testIndex;
00115       testIndex -= stepSize;
00116       // testIndex cannot be smaller than low
00117       testIndex = testIndex > low ? testIndex : low; 
00118       /* Now we have decided the testIndex we would like to
00119          use. However, we must be careful to handle the case when
00120          there are several identical values in the frobsq_norms
00121          list. In that case we use a modified value. */
00122       while(testIndex >= 0 && frobsq_norms[testIndex] == frobsq_norms[testIndex+1])
00123         testIndex--;
00124       /* Note that because of the above while loop, at this point it
00125          is possible that testIndex < low. */
00126       if ( testIndex < 0 ) 
00127         // testIndex == -1, we have to break  
00128         break;
00129       assert( previousTestIndex != testIndex );
00130       Treal currentFrobTrunc = frobsq_norms[testIndex];
00131       frobThreshLowLevel( currentFrobTrunc ); /*=======*/
00132       euclEInt = euclIfSmall( Treal(threshold * 1e-2), threshold ); /*=======*/
00133       // Now we have an interval containing the Euclidean norm of E for the current testIndex
00134     } // end while
00135     Treal euclE; 
00136     if ( testIndex <= -1 ) {
00137       frobThreshLowLevel( (Treal)0.0 ); /*=======*/
00138       euclE = 0;
00139     }
00140     else {
00141       euclE = euclEInt.upp();
00142     }
00143     return euclE;    
00144   } // end run
00145   
00146 
00151   template<typename Tmatrix, typename Treal>
00152     class EuclTruncationSymm : public EuclTruncationBase<Tmatrix, Treal> {
00153   public:
00154     explicit EuclTruncationSymm( Tmatrix & A_ ) 
00155       : EuclTruncationBase<Tmatrix, Treal>(A_) {}
00156   protected:
00157     virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc, 
00158                                      Treal const threshold );
00159     virtual void getFrobSqNorms( std::vector<Treal> & frobsq_norms );
00160     virtual void frobThreshLowLevel( Treal const threshold );
00161     virtual Interval<Treal> euclIfSmall( Treal const absTol,
00162                                          Treal const threshold );
00163   }; // end class EuclTruncationSymm
00164   
00165   template<typename Tmatrix, typename Treal>
00166     void EuclTruncationSymm<Tmatrix, Treal>::getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc, 
00167                                                                  Treal const threshold ) {
00168     /* Divide by 2 because of symmetry */
00169     lowTrunc = (threshold * threshold) / 2;
00170     highTrunc = (threshold * threshold * this->A.get_nrows()) / 2;    
00171   }
00172   
00173   template<typename Tmatrix, typename Treal>
00174     void EuclTruncationSymm<Tmatrix, Treal>::getFrobSqNorms( std::vector<Treal> & frobsq_norms ) {
00175     this->A.getMatrix().getFrobSqLowestLevel(frobsq_norms);
00176   }
00177   
00178   template<typename Tmatrix, typename Treal>
00179     void EuclTruncationSymm<Tmatrix, Treal>::frobThreshLowLevel( Treal const threshold ) {
00180     this->A.getMatrix().frobThreshLowestLevel( threshold, &this->E.getMatrix() );
00181   }
00182   
00183   template<typename Tmatrix, typename Treal>
00184     Interval<Treal> EuclTruncationSymm<Tmatrix, Treal>::euclIfSmall( Treal const absTol,
00185                                                                      Treal const threshold ) {
00186     Treal relTol = std::sqrt(std::sqrt(std::numeric_limits<Treal>::epsilon()));
00187     Interval<Treal> tmpInterval =  mat::euclIfSmall(this->E, absTol, relTol, threshold); 
00188     if ( tmpInterval.length() < 2*absTol )
00189       return Interval<Treal>( tmpInterval.midPoint()-absTol, 
00190                               tmpInterval.midPoint()+absTol );
00191     return tmpInterval;
00192   }
00193 
00200   template<typename Tmatrix, typename TmatrixZ, typename Treal>
00201     class EuclTruncationSymmWithZ : public EuclTruncationSymm<Tmatrix, Treal> {
00202   public:
00203   EuclTruncationSymmWithZ( Tmatrix & A_, TmatrixZ const & Z_ ) 
00204     : EuclTruncationSymm<Tmatrix, Treal>(A_), Z(Z_) {}
00205   protected:
00206     virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc, 
00207                                      Treal const threshold );
00208     // getFrobSqNorms(...)     from EuclTruncationSymm
00209     // frobThreshLowLevel(...) from EuclTruncationSymm
00210     virtual Interval<Treal> euclIfSmall( Treal const absTol,
00211                                          Treal const threshold );
00212     TmatrixZ const & Z;
00213   }; // end class EuclTruncationSymmWithZ
00214   
00215   template<typename Tmatrix, typename TmatrixZ, typename Treal>
00216     void EuclTruncationSymmWithZ<Tmatrix, TmatrixZ, Treal>::getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc, 
00217                                                                                 Treal const threshold ) {
00218     Treal Zfrob = Z.frob();
00219     Treal thresholdTakingZIntoAccount = threshold / (Zfrob * Zfrob);
00220     /* Divide by 2 because of symmetry */
00221     lowTrunc  = thresholdTakingZIntoAccount * thresholdTakingZIntoAccount / 2.0;
00222     highTrunc = std::numeric_limits<Treal>::max();
00223   }
00224   
00225   template<typename Tmatrix, typename TmatrixZ, typename Treal>
00226     Interval<Treal> EuclTruncationSymmWithZ<Tmatrix, TmatrixZ, Treal>::euclIfSmall( Treal const absTol,
00227                                                                                     Treal const threshold ) {
00228     Treal relTol = std::sqrt(std::sqrt(std::numeric_limits<Treal>::epsilon()));
00229     mat::TripleMatrix<Tmatrix, TmatrixZ, Treal> ErrMatTriple( this->E, Z);
00230     Interval<Treal> tmpInterval = mat::euclIfSmall(ErrMatTriple, absTol, relTol, threshold);      
00231     if ( tmpInterval.length() < 2*absTol )
00232       return Interval<Treal>( tmpInterval.midPoint()-absTol, 
00233                               tmpInterval.midPoint()+absTol );
00234     return tmpInterval;
00235   }
00236   
00243   template<typename Tmatrix, typename Treal>
00244     class EuclTruncationSymmElementLevel : public EuclTruncationSymm<Tmatrix, Treal> {
00245   public:
00246   explicit EuclTruncationSymmElementLevel( Tmatrix & A_ ) 
00247     : EuclTruncationSymm<Tmatrix, Treal>(A_) {}
00248   protected:
00249     // getFrobTruncBounds(...) from EuclTruncationSymm
00250     virtual void getFrobSqNorms( std::vector<Treal> & frobsq_norms );
00251     virtual void frobThreshLowLevel( Treal const threshold );
00252     // Interval<Treal> euclIfSmall(...) from EuclTruncationSymm    
00253   }; // end class EuclTruncationSymmElementLevel
00254 
00255   template<typename Tmatrix, typename Treal>
00256     void EuclTruncationSymmElementLevel<Tmatrix, Treal>::getFrobSqNorms( std::vector<Treal> & frobsq_norms ) {
00257     this->A.getMatrix().getFrobSqElementLevel(frobsq_norms);
00258   }
00259 
00260   template<typename Tmatrix, typename Treal>
00261     void EuclTruncationSymmElementLevel<Tmatrix, Treal>::frobThreshLowLevel( Treal const threshold ) {
00262     this->A.getMatrix().frobThreshElementLevel(threshold, &this->E.getMatrix() );
00263   }
00264 
00269   template<typename Tmatrix, typename Treal>
00270     class EuclTruncationGeneral : public EuclTruncationBase<Tmatrix, Treal> {
00271   public:
00272     explicit EuclTruncationGeneral( Tmatrix & A_ ) 
00273       : EuclTruncationBase<Tmatrix, Treal>(A_) {}
00274   protected:
00275     virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc, 
00276                                      Treal const threshold );
00277     virtual void getFrobSqNorms( std::vector<Treal> & frobsq_norms );
00278     virtual void frobThreshLowLevel( Treal const threshold );
00279     virtual Interval<Treal> euclIfSmall( Treal const absTol,
00280                                          Treal const threshold );
00281   }; // end class EuclTruncationGeneral
00282   
00283   template<typename Tmatrix, typename Treal>
00284     void EuclTruncationGeneral<Tmatrix, Treal>::getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc, 
00285                                                                     Treal const threshold ) {
00286     // Try to improve bounds based on the Frobenius norm
00287     /*  ||E||_F^2 <= thres^2  -> 
00288      *  ||E||_F   <= thres    -> 
00289      *  ||E||_2   <= thresh 
00290      */
00291     lowTrunc = (threshold * threshold);
00292     /*  ||E||_F^2 >= thres^2 * n     -> 
00293      *  ||E||_F   >= thres * sqrt(n) -> 
00294      *  ||E||_2   >= thresh 
00295      */
00296     highTrunc = (threshold * threshold * this->A.get_nrows());    
00297   }
00298   
00299   template<typename Tmatrix, typename Treal>
00300     void EuclTruncationGeneral<Tmatrix, Treal>::getFrobSqNorms( std::vector<Treal> & frobsq_norms ) {
00301     this->A.getMatrix().getFrobSqLowestLevel(frobsq_norms);
00302   }
00303   
00304   template<typename Tmatrix, typename Treal>
00305     void EuclTruncationGeneral<Tmatrix, Treal>::frobThreshLowLevel( Treal const threshold ) {
00306     this->A.getMatrix().frobThreshLowestLevel( threshold, &this->E.getMatrix() );
00307   }
00308   
00309   template<typename Tmatrix, typename Treal>
00310     Interval<Treal> EuclTruncationGeneral<Tmatrix, Treal>::euclIfSmall( Treal const absTol,
00311                                                                         Treal const threshold ) {
00312     // FIXME: this should be changed (for all trunc classes) so that
00313     // some relative precision is always requested instead of the input
00314     // absTol which in the current case is not used(!)
00315     mat::ATAMatrix<Tmatrix, Treal> EtE(this->E);    
00316     Treal absTolDummy = std::numeric_limits<Treal>::max(); // Treal(threshold * 1e-2)
00317     Treal relTol = 100 * std::numeric_limits<Treal>::epsilon(); 
00318     Interval<Treal> tmpInterval = mat::euclIfSmall(EtE, absTolDummy, relTol, threshold);
00319     tmpInterval = Interval<Treal>( std::sqrt(tmpInterval.low()), std::sqrt(tmpInterval.upp()) );
00320     if ( tmpInterval.length() < 2*absTol )
00321       return Interval<Treal>( tmpInterval.midPoint()-absTol, 
00322                               tmpInterval.midPoint()+absTol );
00323     return tmpInterval;
00324   }
00325 
00326 
00327 
00328 
00335   template<typename Tmatrix, typename TmatrixB, typename Treal>
00336     class EuclTruncationCongrTransMeasure : public EuclTruncationGeneral<Tmatrix, Treal> {
00337   public:
00338   EuclTruncationCongrTransMeasure( Tmatrix & A_, TmatrixB const & B_ ) 
00339     : EuclTruncationGeneral<Tmatrix, Treal>(A_), B(B_) {}
00340   protected:
00341     virtual void getFrobTruncBounds( Treal & lowTrunc, Treal & highTrunc, 
00342                                      Treal const threshold );
00343     // getFrobSqNorms(...)     from EuclTruncationGeneral
00344     // frobThreshLowLevel(...) from EuclTruncationGeneral
00345     virtual Interval<Treal> euclIfSmall( Treal const absTol,
00346                                          Treal const threshold );
00347     TmatrixB const & B;
00348   }; // end class EuclTruncationCongrTransMeasure
00349   
00350   template<typename Tmatrix, typename TmatrixB, typename Treal>
00351     void EuclTruncationCongrTransMeasure<Tmatrix, TmatrixB, Treal>::getFrobTruncBounds( Treal & lowTrunc, 
00352                                                                                         Treal & highTrunc, 
00353                                                                                         Treal const threshold ) {
00354     Treal Afrob = this->A.frob();
00355     Treal Bfrob = B.frob();
00356     lowTrunc    = -Afrob + std::sqrt( Afrob*Afrob + threshold / Bfrob );
00357     highTrunc   = std::numeric_limits<Treal>::max();
00358   }
00359   
00360   template<typename Tmatrix, typename TmatrixB, typename Treal>
00361     Interval<Treal> EuclTruncationCongrTransMeasure<Tmatrix, TmatrixB, Treal>::euclIfSmall( Treal const absTol,
00362                                                                                             Treal const threshold ) {
00363     Treal relTol = std::sqrt(std::sqrt(std::numeric_limits<Treal>::epsilon()));
00364     mat::CongrTransErrorMatrix<TmatrixB, Tmatrix, Treal> ErrMatTriple( B, this->A, this->E );
00365     Interval<Treal> tmpInterval = mat::euclIfSmall(ErrMatTriple, absTol, relTol, threshold);      
00366     if ( tmpInterval.length() < 2*absTol ) {
00367       return Interval<Treal>( tmpInterval.midPoint()-absTol, 
00368                               tmpInterval.midPoint()+absTol );
00369     }
00370     return tmpInterval;
00371   }
00372 
00373 
00374 } /* end namespace mat */
00375 #endif