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