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_PERTURBATION 00037 #define MAT_PERTURBATION 00038 namespace per { 00039 template<typename Treal, typename Tmatrix, typename Tvector> 00040 class Perturbation { 00041 public: 00042 Perturbation(std::vector<Tmatrix *> const & F, 00044 std::vector<Tmatrix *> & D, 00046 mat::Interval<Treal> const & gap, 00047 mat::Interval<Treal> const & allEigs, 00053 Treal const deltaMax, 00054 Treal const errorTol, 00055 mat::normType const norm, 00056 Tvector & vect 00057 ); 00058 void perturb() { 00059 dryRun(); 00060 run(); 00061 } 00062 00063 void checkIdempotencies(std::vector<Treal> & idemErrors); 00064 template<typename TmatNoSymm> 00065 void checkCommutators(std::vector<Treal> & commErrors, 00066 TmatNoSymm const & dummyMat); 00067 void checkMaxSubspaceError(Treal & subsError); 00068 00069 protected: 00070 /* This is input from the beginning */ 00071 std::vector<Tmatrix *> const & F; 00072 std::vector<Tmatrix *> & X; 00073 mat::Interval<Treal> gap; 00074 mat::Interval<Treal> const & allEigs; 00075 Treal deltaMax; 00076 Treal errorTol; 00077 mat::normType const norm; 00078 Tvector & vect; 00079 00080 /* These variables are set in the dry run. */ 00081 int nIter; 00082 std::vector<Treal> threshVal; 00083 std::vector<Treal> sigma; 00084 00095 void dryRun(); 00096 void run(); 00097 private: 00098 00099 }; 00100 00101 template<typename Treal, typename Tmatrix, typename Tvector> 00102 Perturbation<Treal, Tmatrix, Tvector>:: 00103 Perturbation(std::vector<Tmatrix *> const & F_in, 00104 std::vector<Tmatrix *> & X_in, 00105 mat::Interval<Treal> const & gap_in, 00106 mat::Interval<Treal> const & allEigs_in, 00107 Treal const deltaMax_in, 00108 Treal const errorTol_in, 00109 mat::normType const norm_in, 00110 Tvector & vect_in) 00111 : F(F_in), X(X_in), gap(gap_in), allEigs(allEigs_in), 00112 deltaMax(deltaMax_in), errorTol(errorTol_in), norm(norm_in), 00113 vect(vect_in) { 00114 if (!X.empty()) 00115 throw "Perturbation constructor: D vector is expected to be empty (size==0)"; 00116 for (unsigned int iMat = 0; iMat < F.size(); ++iMat) 00117 X.push_back(new Tmatrix(*F[iMat])); 00118 00119 Treal lmin = allEigs.low(); 00120 Treal lmax = allEigs.upp(); 00121 00122 /***** Initial linear transformation of matrix sequence. */ 00123 typename std::vector<Tmatrix *>::iterator matIt = X.begin(); 00124 /* Scale to [0, 1] interval and negate */ 00125 (*matIt)->add_identity(-lmax); 00126 *(*matIt) *= ((Treal)1.0 / (lmin - lmax)); 00127 matIt++; 00128 /* ...and derivatives: */ 00129 for ( ; matIt != X.end(); matIt++ ) 00130 *(*matIt) *= ((Treal)-1.0 / (lmin - lmax)); 00131 /* Compute transformed gap */ 00132 gap = (gap - lmax) / (lmin - lmax); 00133 } 00134 00135 template<typename Treal, typename Tmatrix, typename Tvector> 00136 void Perturbation<Treal, Tmatrix, Tvector>::dryRun() { 00137 Treal errorTolPerIter; 00138 int nIterGuess = 0; 00139 nIter = 1; 00140 Treal lumo; 00141 Treal homo; 00142 Treal m; 00143 Treal g; 00144 while (nIterGuess < nIter) { 00145 nIterGuess++; 00146 errorTolPerIter = 0.5 * errorTol /nIterGuess; 00147 nIter = 0; 00148 mat::Interval<Treal> gapTmp(gap); 00149 sigma.resize(0); 00150 threshVal.resize(0); 00151 while (gapTmp.low() > 0.5 * errorTol || gapTmp.upp() < 0.5 * errorTol) { 00152 lumo = gapTmp.low(); 00153 homo = gapTmp.upp(); 00154 m = gapTmp.midPoint(); 00155 g = gapTmp.length(); 00156 if (m > 0.5) { 00157 lumo = lumo*lumo; 00158 homo = homo*homo; 00159 sigma.push_back(-1); 00160 } 00161 else { 00162 lumo = 2*lumo - lumo*lumo; 00163 homo = 2*homo - homo*homo; 00164 sigma.push_back(1); 00165 } 00166 /* Compute threshold value necessary to converge. */ 00167 Treal forceConvThresh = template_blas_fabs(1-2*m) * g / 10; 00168 /* We divide by 10 > 2 so that this loop converges at some point. */ 00169 /* Compute threshold value necessary to maintain accuracy in subspace.*/ 00170 Treal subspaceThresh = errorTolPerIter * (homo-lumo) / (1+errorTolPerIter); 00171 /* Choose the most restrictive threshold of the two. */ 00172 threshVal.push_back(forceConvThresh < subspaceThresh ? 00173 forceConvThresh : subspaceThresh); 00174 homo -= threshVal.back(); 00175 lumo += threshVal.back(); 00176 gapTmp = mat::Interval<Treal>(lumo, homo); 00177 if (gapTmp.empty()) 00178 throw "Perturbation<Treal, Tmatrix, Tvector>::dryRun() : Perturbation iterations will fail to converge; Gap is too small or desired accuracy too high."; 00179 nIter++; 00180 } /* end 2nd while */ 00181 } /* end 1st while */ 00182 /* Now, we have nIter, threshVal, and sigma. */ 00183 } 00184 00185 template<typename Treal, typename Tmatrix, typename Tvector> 00186 void Perturbation<Treal, Tmatrix, Tvector>::run() { 00187 Treal const ONE = 1.0; 00188 mat::SizesAndBlocks rowsCopy; 00189 X.front()->getRows(rowsCopy); 00190 mat::SizesAndBlocks colsCopy; 00191 X.front()->getCols(colsCopy); 00192 Tmatrix tmpMat; 00193 // tmpMat.resetSizesAndBlocks(rowsCopy, colsCopy); 00194 int nMatrices; 00195 Treal threshValPerOrder; 00196 Treal chosenThresh; 00197 for (int iter = 0; iter < nIter; iter++) { 00198 std::cout<<"\n\nInside outer loop iter = "<<iter 00199 <<" nIter = "<<nIter 00200 <<" sigma = "<<sigma[iter]<<std::endl; 00201 /* Number of matrices increases by 1 in each iteration: */ 00202 X.push_back(new Tmatrix); 00203 nMatrices = X.size(); 00204 X[nMatrices-1]->resetSizesAndBlocks(rowsCopy, colsCopy); 00205 /* Compute threshold value for each order. */ 00206 threshValPerOrder = threshVal[iter] / nMatrices; 00207 /* Loop through all nonzero orders. */ 00208 std::cout<<"Entering inner loop nMatrices = "<<nMatrices<<std::endl; 00209 for (int j = nMatrices-1 ; j >= 0 ; --j) { 00210 std::cout<<"Inside inner loop j = "<<j<<std::endl; 00211 std::cout<<"X[j]->eucl() (before compute) = "<<X[j]->eucl(vect,1e-7)<<std::endl; 00212 std::cout<<"X[j]->frob() (before compute) = "<<X[j]->frob()<<std::endl; 00213 tmpMat = Treal(Treal(1.0)+sigma[iter]) * (*X[j]); 00214 std::cout<<"tmpMat.eucl() (before for) = "<<tmpMat.eucl(vect,1e-7)<<std::endl; 00215 std::cout<<"tmpMat.frob() (before for) = "<<tmpMat.frob()<<std::endl; 00216 for (int k = 0; k <= j; k++) { 00217 /* X[j] = X[j] - sigma * X[k] * X[j-k] */ 00218 Tmatrix::ssmmUpperTriangleOnly(-sigma[iter], *X[k], *X[j-k], 00219 ONE, tmpMat); 00220 } /* End 3rd for */ 00221 std::cout<<"tmpMat.eucl() (after for) = "<<tmpMat.eucl(vect,1e-7)<<std::endl; 00222 *X[j] = tmpMat; 00223 00224 /* Truncate tmpMat, remove if it becomes zero. */ 00225 chosenThresh = threshValPerOrder / pow(deltaMax, Treal(j)); 00226 std::cout<<"X[j]->eucl() (before thresh) = "<<X[j]->eucl(vect,1e-7)<<std::endl; 00227 std::cout<<"Chosen thresh: "<<chosenThresh<<std::endl; 00228 Treal actualThresh = X[j]->thresh(chosenThresh, vect, norm); 00229 std::cout<<"X[j]->eucl() (after thresh) = "<<X[j]->eucl(vect,1e-7)<<std::endl; 00230 #if 1 00231 /* If the current matrix is zero AND 00232 * it is the last matrix 00233 */ 00234 if (*X[j] == 0 && int(X.size()) == j+1) { 00235 std::cout<<"DELETION: j = "<<j<<" X.size() = "<<X.size() 00236 <<" X[j] = "<<X[j]<< " X[j]->frob() = "<<X[j]->frob() 00237 <<std::endl; 00238 delete X[j]; 00239 X.pop_back(); 00240 } 00241 else 00242 std::cout<<"NO DELETION: j = "<<j<<" X.size() = "<<X.size() 00243 <<" X[j] = "<<X[j]<< " X[j]->frob() = "<<X[j]->frob() 00244 <<std::endl; 00245 #endif 00246 00247 } /* End 2nd for (Loop through orders) */ 00248 } /* End 1st for (Loop through iterations) */ 00249 } /* End run() */ 00250 00251 00252 template<typename Treal, typename Tmatrix, typename Tvector> 00253 void Perturbation<Treal, Tmatrix, Tvector>:: 00254 checkIdempotencies(std::vector<Treal> & idemErrors) { 00255 Tmatrix tmpMat; 00256 Treal const ONE = 1.0; 00257 unsigned int j; 00258 for (unsigned int m = 0; m < X.size(); ++m) { 00259 tmpMat = (-ONE) * (*X[m]); 00260 for (unsigned int i = 0; i <= m; ++i) { 00261 j = m - i; 00262 /* TMP = TMP + X[i] * X[j] */ 00263 Tmatrix::ssmmUpperTriangleOnly(ONE, *X[i], *X[j], ONE, tmpMat); 00264 } 00265 /* TMP is symmetric! */ 00266 idemErrors.push_back(tmpMat.eucl(vect,1e-10)); 00267 } 00268 } 00269 00270 template<typename Treal, typename Tmatrix, typename Tvector> 00271 template<typename TmatNoSymm> 00272 void Perturbation<Treal, Tmatrix, Tvector>:: 00273 checkCommutators(std::vector<Treal> & commErrors, 00274 TmatNoSymm const & dummyMat) { 00275 mat::SizesAndBlocks rowsCopy; 00276 X.front()->getRows(rowsCopy); 00277 mat::SizesAndBlocks colsCopy; 00278 X.front()->getCols(colsCopy); 00279 TmatNoSymm tmpMat; 00280 tmpMat.resetSizesAndBlocks(rowsCopy, colsCopy); 00281 Treal const ONE = 1.0; 00282 unsigned int j; 00283 for (unsigned int m = 0; m < X.size(); ++m) { 00284 tmpMat = 0; 00285 std::cout<<"New loop\n"; 00286 for (unsigned int i = 0; i <= m && i < F.size(); ++i) { 00287 j = m - i; 00288 std::cout<<i<<", "<<j<<std::endl; 00289 /* TMP = TMP + F[i] * X[j] - X[j] * F[i] */ 00290 tmpMat += ONE * (*F[i]) * (*X[j]); 00291 tmpMat += -ONE * (*X[j]) * (*F[i]); 00292 } 00293 /* TMP is not symmetric! */ 00294 commErrors.push_back(tmpMat.frob()); 00295 } 00296 } 00297 00298 00299 template<typename Treal, typename Tmatrix, typename Tvector> 00300 void Perturbation<Treal, Tmatrix, Tvector>:: 00301 checkMaxSubspaceError(Treal & subsError) { 00302 Treal const ONE = 1.0; 00303 Tmatrix XdeltaMax(*F.front()); 00304 for (unsigned int ind = 1; ind < F.size(); ++ind) 00305 XdeltaMax += pow(deltaMax, Treal(ind)) * (*F[ind]); 00306 /***** Initial linear transformation of matrix. */ 00307 Treal lmin = allEigs.low(); 00308 Treal lmax = allEigs.upp(); 00309 /* Scale to [0, 1] interval and negate */ 00310 XdeltaMax.add_identity(-lmax); 00311 XdeltaMax *= ((Treal)1.0 / (lmin - lmax)); 00312 00313 Tmatrix X2; 00314 for (int iter = 0; iter < nIter; iter++) { 00315 X2 = ONE * XdeltaMax * XdeltaMax; 00316 if (sigma[iter] == Treal(1.0)) { 00317 XdeltaMax *= 2.0; 00318 XdeltaMax -= X2; 00319 } 00320 else { 00321 XdeltaMax = X2; 00322 } 00323 } /* End of for (Loop through iterations) */ 00324 00325 Tmatrix DdeltaMax(*X.front()); 00326 for (unsigned int ind = 1; ind < X.size(); ++ind) 00327 DdeltaMax += pow(deltaMax, Treal(ind)) * (*X[ind]); 00328 subsError = Tmatrix::eucl_diff(XdeltaMax,DdeltaMax, 00329 vect, errorTol *1e-2); 00330 } 00331 00332 00333 00334 } /* end namespace mat */ 00335 #endif 00336