ergo
purification_old.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 
00034 #ifndef MAT_PURIFICATION_OLD
00035 #define MAT_PURIFICATION_OLD
00036 #include <iomanip>
00037 #include <math.h>
00038 namespace mat{
00039 
00040 
00041   template<class Treal, class MAT>
00042     void tc2_extra(MAT& D, const int nshells, 
00043                    const int nsteps, const Treal frob_trunc = 0) {
00044     MAT tmp(D);
00045     for (int step = 0; step < nsteps; step++) {
00046       Treal tracediff =  tmp.trace() - nshells;
00047       if (tracediff > 0) {
00048         tmp = (Treal)1.0 * D * D;
00049         D = tmp;
00050         D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
00051       }
00052       else {
00053         tmp = D;
00054         tmp = (Treal)-1.0 * D * D + (Treal)2.0 * tmp;
00055         D = (Treal)1.0 * tmp * tmp;     
00056       }      
00057       D.frob_thresh(frob_trunc);
00058     }      
00059   } 
00060 
00061   template<class Treal, class MAT>
00062     void tc2_extra_auto(MAT& D, const int nshells, 
00063                         int& nsteps, const Treal frob_trunc,
00064                         const int maxiter = 50) {
00065     MAT tmp(D);
00066     Treal froberror;
00067     Treal tracediff =  tmp.trace() - nshells;
00068     nsteps = 0;
00069     do { 
00070       if (tracediff > 0) {
00071         tmp = (Treal)1.0 * D * D;
00072         D = tmp;
00073         D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
00074       }
00075       else {
00076         tmp = D;
00077         tmp = (Treal)-1.0 * D * D + (Treal)2.0 * tmp;
00078         D = (Treal)1.0 * tmp * tmp;     
00079       }
00080       froberror = MAT::frob_diff(D, tmp);
00081       D.frob_thresh(frob_trunc);
00082       tracediff = D.trace() - nshells;
00083 #if 0
00084       std::cout<<"Iteration:"<<nsteps<<"  Froberror: "
00085                <<std::setprecision(8)<<froberror<<
00086         "  Tracediff: "<<tracediff<<std::endl;
00087 #endif
00088       nsteps++;
00089       if (nsteps > maxiter)
00090         throw AcceptableMaxIter("Extra Auto Purification reached maxiter" 
00091                                 " without convergence", maxiter);
00092     } while (froberror > frob_trunc);
00093   } 
00094 
00095   
00096   
00097   template<class Treal, class MAT>
00098     void tc2(MAT& D, int& iterations,
00099              const MAT& H, const int nshells, 
00100              const Treal trace_error_limit = 1e-2,
00101              const Treal frob_trunc = 0,
00102              int* polys= NULL,
00103              const int maxiter = 100) {
00104     MAT tmp(H);
00105     D = H;
00106     iterations = 0;
00107     Treal tracediff =  tmp.trace() - nshells;
00108     Treal tracediff_old = 2.0 * trace_error_limit;
00109     
00110     while (template_blas_fabs(tracediff) > trace_error_limit && 
00111            template_blas_fabs(tracediff_old) > trace_error_limit) {
00112       //      std::cout<<"Iteration:"<<iterations
00113       //               <<"   Tracediff ="<<tracediff<<std::endl;
00114       if (tracediff > 0) {
00115         D = (Treal)1.0 * tmp * tmp;
00116         if (polys)
00117           polys[iterations] = 0;
00118       }
00119       else {
00120         D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
00121         if (polys)
00122           polys[iterations] = 1;
00123       }
00124       D.frob_thresh(frob_trunc);
00125       tmp = D;
00126       tracediff_old = tracediff;
00127       tracediff =  D.trace() - nshells;
00128       iterations++;
00129       if (iterations > maxiter)
00130         throw AcceptableMaxIter("Purification reached maxiter without convergence", maxiter);
00131     }
00132   }
00133 
00134 #if 1
00135   template<class Treal, class MAT>
00136     void tc2_auto(MAT& D, int& iterations,
00137                   const MAT& H, const int nocc, 
00138                   const Treal frob_trunc = 0,
00139                   int* polys= NULL,
00140                   const int maxiter = 100) {
00141     assert(frob_trunc >= 0);
00142     assert(nocc >= 0);
00143     assert(maxiter >= 0);
00144     MAT tmp(H);
00145     D = H;
00146     iterations = 0;
00147     Treal tracediff = 2;
00148     Treal newtracediff = tmp.trace() - nocc;
00149     Treal froberror = 1;
00150     Treal tracechange = 0; /* Initial value has no relevance */
00151     /* Need check for too loose truncation? */
00152     while (2 * froberror > (3 - template_blas_sqrt((Treal)5)) / 2 - frob_trunc || 
00153            template_blas_fabs(tracediff) > 1 || 
00154            template_blas_fabs(tracechange) > (1 - 2 * froberror) * (1 - template_blas_fabs(tracediff))) {
00155       //      std::cout<<"Iteration:"<<iterations
00156       //               <<"   Tracediff ="<<tracediff<<std::endl;
00157       tracediff = newtracediff;
00158       if (tracediff > 0) {
00159         D = (Treal)1.0 * tmp * tmp;
00160         if (polys)
00161           polys[iterations] = 0;
00162       }
00163       else {
00164         D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
00165         if (polys)
00166           polys[iterations] = 1;
00167       }
00168       D.frob_thresh(frob_trunc);
00169       newtracediff = D.trace() - nocc;
00170       tracechange = newtracediff - tracediff;
00171       froberror = MAT::frob_diff(D, tmp);
00172       tmp = D;
00173 #if 0
00174       std::cout<<"Iteration:"<<iterations<<"  epsilon: "
00175                <<std::setprecision(8)<<std::setw(12)<<froberror<<
00176         "  delta: "<<std::setw(12)<<tracediff<<
00177         "  tracechange: "<<std::setw(12)<<tracechange<<std::endl;
00178 #endif
00179       iterations++;
00180       if (iterations > maxiter)
00181         throw AcceptableMaxIter("tc2_auto::Purification reached maxiter"
00182                                 " without convergence", maxiter);
00183     }
00184     
00185     /* Take one step to make sure the eigenvalues stays in */
00186     /* { [ 0 , 2 * epsilon [ , ] 1 - 2 * epsilon , 1] }    */
00187     if (tracediff > 0) {
00188       D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
00189       if (polys)
00190         polys[iterations] = 1;
00191     }
00192     else {
00193       D = (Treal)1.0 * tmp * tmp;
00194       if (polys)
00195         polys[iterations] = 0;
00196     }
00197     iterations++;
00198 
00199     /* Use second order convergence polynomials to converge completely */
00200     D.frob_thresh(frob_trunc);
00201     tracediff = D.trace() - nocc;
00202     do { 
00203       if (tracediff > 0) {
00204         tmp = (Treal)1.0 * D * D;
00205         D = tmp;
00206         D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
00207         if (polys) {
00208           polys[iterations] = 0;
00209           polys[iterations + 1] = 1;
00210         }
00211       }
00212       else {
00213         tmp = D;
00214         tmp = (Treal)-1.0 * D * D + (Treal)2.0 * tmp;
00215         D = (Treal)1.0 * tmp * tmp;
00216         if (polys) {
00217           polys[iterations] = 1;
00218           polys[iterations + 1] = 0;
00219         }
00220       }
00221       froberror = MAT::frob_diff(D, tmp);
00222       D.frob_thresh(frob_trunc);
00223       tracediff = D.trace() - nocc;
00224 #if 0
00225       std::cout<<"Iteration:"<<iterations<<"  Froberror: "
00226                <<std::setprecision(8)<<froberror<<
00227         "  Tracediff: "<<tracediff<<std::endl;
00228 #endif
00229       iterations += 2;
00230       if (iterations > maxiter)
00231         throw AcceptableMaxIter("tc2_auto::Purification reached maxiter" 
00232                                 " in extra part without convergence", maxiter);
00233     } while (froberror > frob_trunc);
00234     
00235 
00236   }
00237 
00238 #endif
00239 
00240 
00241 
00242 #if 0
00243   iterations = 0;
00244   /* D = (F - lmax * I) / (lmin - lmax)                                  */
00245   D = F;
00246   D.add_identity(-lmax);       
00247     D *= (1.0 / (lmin - lmax));
00248     /*****/ clock_t start;
00249     /*****/ float syrktime = 0;
00250     /*****/ float threshtime = 0;
00251     MAT tmp;     
00252     Treal tracediff;
00253     for (int steps = 0; steps < nextra_steps; steps++) {
00254       tmp = D;
00255       tracediff = tmp.trace() - nshells;
00256       while (template_blas_fabs(tracediff) > trace_error_limit) {
00257         iterations++;
00258         /*****/ start = clock();
00259         if (tracediff > 0) 
00260           MAT::syrk('U', false, 1.0, tmp, 0.0, D);
00261         else 
00262           MAT::syrk('U', false, -1.0, tmp, 2.0, D);
00263         /*****/ syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
00264         D.sym_to_nosym();
00265         /*****/ start = clock();
00266         D.frob_thresh(frob_trunc);
00267         /*****/ threshtime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
00268         tmp = D;
00269         tracediff = tmp.trace() - nshells;
00270       } /* end while */
00271       
00272       /* extra steps */
00273       iterations += 2;
00274       if (tracediff > 0) {
00275         /*****/ start = clock();
00276         MAT::syrk('U', false, 1.0, tmp, 0.0, D);
00277         /*****/ syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
00278         D.sym_to_nosym();
00279         tmp = D;
00280         /*****/ start = clock();
00281         MAT::syrk('U', false, -1.0, tmp, 2.0, D);
00282         /*****/ syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
00283       }
00284       else {
00285         /*****/ start = clock();
00286         MAT::syrk('U', false, -1.0, tmp, 2.0, D);
00287         /*****/ syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
00288         D.sym_to_nosym();
00289         tmp = D;
00290         /*****/ start = clock();
00291         MAT::syrk('U', false, 1.0, tmp, 0.0, D);
00292         /*****/ syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
00293       }      
00294       D.sym_to_nosym();
00295         /*****/ start = clock();
00296       D.frob_thresh(frob_trunc);
00297         /*****/ threshtime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
00298     } /* end for */
00300     std::cout<<std::endl<<" total syrktime in tcp ="<<std::setw(5)
00301              <<syrktime<<std::endl;
00302     std::cout<<" total threshtime in tcp ="<<std::setw(5)
00303              <<threshtime<<std::endl;
00304 
00305   }
00306 #endif  
00307 } /* end namespace mat */
00308 
00309 #endif