00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
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
00113
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;
00151
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
00156
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
00186
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
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
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 }
00271
00272
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 }
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 }
00308
00309 #endif