ergo
Perturbation.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_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