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 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