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
00028
00029
00039 #ifndef HEADER_PURIFICATION_SP2
00040 #define HEADER_PURIFICATION_SP2
00041
00042 #include "purification_general.h"
00043
00044
00045
00050 template<typename MatrixType>
00051 class Purification_sp2 : public PurificationGeneral<MatrixType>
00052 {
00053 public:
00054
00055 typedef typename PurificationGeneral<MatrixType>::real real;
00056 typedef typename PurificationGeneral<MatrixType>::IntervalType IntervalType;
00057 typedef typename PurificationGeneral<MatrixType>::NormType NormType;
00058
00059 typedef typename PurificationGeneral<MatrixType>::VectorTypeInt VectorTypeInt;
00060 typedef typename PurificationGeneral<MatrixType>::VectorTypeReal VectorTypeReal;
00061
00062 typedef generalVector VectorType;
00063
00064 Purification_sp2() : PurificationGeneral<MatrixType>() {}
00065
00066 void set_init_params()
00067 {
00068 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen SP2 purification method");
00069 this->info.method = 1;
00070
00071 this->gammaStopEstim = (3 - template_blas_sqrt((real)5)) / 2.0;
00072 }
00073
00074 void get_poly(const int it, int& poly);
00075 void set_poly(const int it);
00076
00077 void truncate_matrix(real& threshold, int it);
00078
00079 void estimate_number_of_iterations(int& numit);
00080 void purify_X(const int it);
00081 void purify_bounds(const int it);
00082 void save_other_iter_info(IterationInfo& iter_info, int it);
00083 void apply_inverse_poly_vector(const int it, VectorTypeReal& bounds_from_it);
00084
00085 void return_constant_C(const int it, real& Cval);
00086
00087
00088 real apply_poly(const int it, real x);
00089 void apply_poly_to_eigs(const int it, real& homo, real& lumo);
00090 real compute_derivative(const int it, real x, real& DDf);
00091 };
00092
00093 template<typename MatrixType>
00094 void Purification_sp2<MatrixType>::set_poly(const int it)
00095 {
00096 assert((int)this->VecPoly.size() > it);
00097
00098
00099 if (this->VecPoly[it] == -1)
00100 {
00101 real Xtrace = this->X.trace();
00102 real Xsqtrace = this->Xsq.trace();
00103
00104 if ((template_blas_fabs(Xsqtrace - this->nocc) <
00105 template_blas_fabs(2 * Xtrace - Xsqtrace - this->nocc))
00106 ||
00107 (it % 2
00108 &&
00109 (template_blas_fabs(Xsqtrace - this->nocc) ==
00110 template_blas_fabs(2 * Xtrace - Xsqtrace - this->nocc))
00111 ))
00112 {
00113 this->VecPoly[it] = 1;
00114 }
00115 else
00116 {
00117 this->VecPoly[it] = 0;
00118 }
00119 }
00120 }
00121
00122
00123 template<typename MatrixType>
00124 void Purification_sp2<MatrixType>::get_poly(const int it, int& poly)
00125 {
00126 assert((int)this->VecPoly.size() > it);
00127 assert(this->VecPoly[it] != -1);
00128 poly = this->VecPoly[it];
00129 }
00130
00131
00132 template<typename MatrixType>
00133 void Purification_sp2<MatrixType>::purify_X(const int it)
00134 {
00135 #ifdef DEBUG_OUTPUT
00136 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Purify X...");
00137 #endif
00138 int poly;
00139
00140 set_poly(it);
00141 get_poly(it, poly);
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157 if (poly == 0)
00158 {
00159 this->Xsq *= ((real) - 1.0);
00160 this->X *= (real)2.0;
00161 this->Xsq += this->X;
00162 }
00163
00164 this->Xsq.transfer(this->X);
00165 }
00166
00167
00168 template<typename MatrixType>
00169 void Purification_sp2<MatrixType>::purify_bounds(const int it)
00170 {
00171 #ifdef DEBUG_OUTPUT
00172 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Change homo and lumo bounds according to the chosen polynomial VecPoly = %d", this->VecPoly[it]);
00173 #endif
00174 real homo_low, homo_upp, lumo_upp, lumo_low;
00175 int poly;
00176
00177 get_poly(it, poly);
00178
00179 if (poly == 1)
00180 {
00181
00182 homo_low = 2 * this->homo_bounds.low() - this->homo_bounds.low() * this->homo_bounds.low();
00183 homo_upp = 2 * this->homo_bounds.upp() - this->homo_bounds.upp() * this->homo_bounds.upp();
00184 lumo_low = this->lumo_bounds.low() * this->lumo_bounds.low();
00185 lumo_upp = this->lumo_bounds.upp() * this->lumo_bounds.upp();
00186 this->homo_bounds = IntervalType(homo_low, homo_upp);
00187 this->lumo_bounds = IntervalType(lumo_low, lumo_upp);
00188 }
00189 else
00190 {
00191
00192 lumo_low = 2 * this->lumo_bounds.low() - this->lumo_bounds.low() * this->lumo_bounds.low();
00193 lumo_upp = 2 * this->lumo_bounds.upp() - this->lumo_bounds.upp() * this->lumo_bounds.upp();
00194 homo_low = this->homo_bounds.low() * this->homo_bounds.low();
00195 homo_upp = this->homo_bounds.upp() * this->homo_bounds.upp();
00196 this->homo_bounds = IntervalType(homo_low, homo_upp);
00197 this->lumo_bounds = IntervalType(lumo_low, lumo_upp);
00198 }
00199
00200 IntervalType zero_one(0, 1);
00201 this->homo_bounds.intersect(zero_one);
00202 this->lumo_bounds.intersect(zero_one);
00203
00204 #ifdef DEBUG_OUTPUT
00205 if (this->homo_bounds.empty())
00206 {
00207 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval homo_bounds is empty.");
00208 }
00209 if (this->lumo_bounds.empty())
00210 {
00211 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval lumo_bounds is empty.");
00212 }
00213
00214 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "1-homo: [ %lf , %lf ],", this->homo_bounds.low(), this->homo_bounds.upp());
00215 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "lumo: [ %lf , %lf ].", this->lumo_bounds.low(), this->lumo_bounds.upp());
00216 #endif
00217 }
00218
00219
00220
00221
00222 template<typename MatrixType>
00223 void Purification_sp2<MatrixType>::return_constant_C(const int it, real& Cval)
00224 {
00225 Cval = C_SP2;
00226 }
00227
00228
00229
00230
00231 template<typename MatrixType>
00232 void Purification_sp2<MatrixType>::truncate_matrix(real& threshold, int it)
00233 {
00234 real allowed_error = this->error_per_it;
00235
00236 threshold = 0;
00237 if (allowed_error == 0)
00238 {
00239 return;
00240 }
00241
00242 assert((int)this->VecGap.size() > it);
00243 #ifdef DEBUG_OUTPUT
00244 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Truncate matrix X: ");
00245 #endif
00246
00247 real tau;
00248
00249 if (this->VecGap[it] > 0)
00250 {
00251 #ifdef DEBUG_OUTPUT
00252 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "VecGap[ %d ] = %e , ", it, this->VecGap[it]);
00253 #endif
00254 tau = (allowed_error * this->VecGap[it]) / (1 + allowed_error);
00255 }
00256 else
00257 {
00258 tau = allowed_error * 0.01;
00259 }
00260
00261
00262 #ifdef USE_CHUNKS_AND_TASKS
00263 threshold = (this->X).thresh_frob(tau);
00264 #else
00265 threshold = (this->X).thresh(tau, this->normPuriTrunc);
00266 #endif
00267
00268 #ifdef DEBUG_OUTPUT
00269 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "tau = %e", tau);
00270 #endif
00271 }
00272
00273
00274
00275
00276
00277 template<typename MatrixType>
00278 void Purification_sp2<MatrixType>::estimate_number_of_iterations(int& estim_num_iter)
00279 {
00280 int it = 1;
00281 int maxit_tmp = this->maxit + this->additional_iterations;
00282 real x, y;
00283 real epsilon = this->get_epsilon();
00284
00285
00286 this->check_stopping_criterion_iter = 2;
00287
00288 int max_size = maxit_tmp + 1 + this->additional_iterations + 2;
00289
00290 this->VecPoly.clear();
00291 this->VecPoly.resize(max_size, -1);
00292
00293 this->VecGap.clear();
00294 this->VecGap.resize(max_size, -1);
00295
00296
00297 x = this->lumo_bounds.upp();
00298 y = this->homo_bounds.upp();
00299
00300
00301
00302 if (1 - x - y <= 0)
00303 {
00304 #ifdef DEBUG_OUTPUT
00305 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "VecGap cannot be computed. Set estimated number of iteration to the maxit.");
00306 #endif
00307 estim_num_iter = this->maxit;
00308 return;
00309 }
00310
00311 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "INIT LUMO: [ %.12lf , %.12lf ]", (double)this->lumo_bounds.low(), (double)this->lumo_bounds.upp());
00312 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "INIT HOMO: [ %.12lf , %.12lf ]", (double)this->homo_bounds.low(), (double)this->homo_bounds.upp());
00313
00314
00315
00316 this->VecPoly[0] = -1;
00317 this->VecGap[0] = 1 - x - y;
00318
00319 estim_num_iter = -1;
00320
00321 while (it <= maxit_tmp)
00322 {
00323
00324 if ((x > y) || (it % 2 && (x == y)))
00325 {
00326 x *= x;
00327 y = 2 * y - y * y;
00328 this->VecPoly[it] = 1;
00329 }
00330 else
00331 {
00332 x = 2 * x - x * x;
00333 y *= y;
00334 this->VecPoly[it] = 0;
00335 }
00336
00337 this->VecGap[it] = 1 - x - y;
00338
00339
00340 if ((estim_num_iter == -1) &&
00341 (x - x * x < epsilon) && (y - y * y < epsilon) &&
00342 (this->VecPoly[it] != this->VecPoly[it - 1]))
00343 {
00344 estim_num_iter = it;
00345 maxit_tmp = it + this->additional_iterations;
00346
00347 if (this->normPuriStopCrit == mat::frobNorm)
00348 {
00349 estim_num_iter += 2;
00350 maxit_tmp += 2;
00351 }
00352 }
00353
00354 ++it;
00355 }
00356
00357
00358
00359 if ((estim_num_iter == -1) && (it == maxit_tmp + 1))
00360 {
00361 #ifdef DEBUG_OUTPUT
00362 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "maxit = %d number of iterations is reached in estimate_number_of_iteration()", this->maxit);
00363 #endif
00364 estim_num_iter = this->maxit;
00365 }
00366
00367
00368 this->VecPoly.resize(maxit_tmp + 1);
00369 this->VecGap.resize(maxit_tmp + 1);
00370
00371 #ifdef DEBUG_OUTPUT
00372 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Sequence of polynomials VecPoly: ");
00373 for (int i = 0; i < (int)this->VecPoly.size(); ++i)
00374 {
00375 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "%d ", this->VecPoly[i]);
00376 }
00377 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "");
00378 #endif
00379 }
00380
00381
00382
00383
00384 template<typename MatrixType>
00385 void Purification_sp2<MatrixType>::save_other_iter_info(IterationInfo& iter_info, int it)
00386 {
00387 assert((int)this->VecPoly.size() > it);
00388 assert((int)this->VecGap.size() > it);
00389
00390 iter_info.poly = this->VecPoly[it];
00391 iter_info.gap = this->VecGap[it];
00392 }
00393
00394
00395
00396
00397
00398 template<typename MatrixType>
00399 void Purification_sp2<MatrixType>::apply_inverse_poly_vector(const int it, VectorTypeReal& bounds_from_it)
00400 {
00401 int poly;
00402 for (int i = it; i >= 1; i--)
00403 {
00404 get_poly(i, poly);
00405
00406 if (poly == 1)
00407 {
00408 bounds_from_it[0] = template_blas_sqrt(bounds_from_it[0]);
00409 bounds_from_it[1] = template_blas_sqrt(bounds_from_it[1]);
00410
00411 bounds_from_it[2] = bounds_from_it[2] / (1 + template_blas_sqrt(1 - bounds_from_it[2]));
00412 bounds_from_it[3] = bounds_from_it[3] / (1 + template_blas_sqrt(1 - bounds_from_it[3]));
00413 }
00414 else
00415 {
00416 bounds_from_it[0] = bounds_from_it[0] / (1 + template_blas_sqrt(1 - bounds_from_it[0]));
00417 bounds_from_it[1] = bounds_from_it[1] / (1 + template_blas_sqrt(1 - bounds_from_it[1]));
00418
00419 bounds_from_it[2] = template_blas_sqrt(bounds_from_it[2]);
00420 bounds_from_it[3] = template_blas_sqrt(bounds_from_it[3]);
00421 }
00422 }
00423 }
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451 template<typename MatrixType>
00452 typename Purification_sp2<MatrixType>::real
00453 Purification_sp2<MatrixType>::apply_poly(const int it, real x)
00454 {
00455 assert(it >= 0);
00456 if (it == 0)
00457 {
00458 return x;
00459 }
00460
00461 real fx;
00462 int poly;
00463 get_poly(it, poly);
00464
00465 if (poly == 1)
00466 {
00467 fx = x * x;
00468 }
00469 else
00470 {
00471 fx = 2 * x - x * x;
00472 }
00473
00474 return fx;
00475 }
00476
00477
00478 template<typename MatrixType>
00479 void Purification_sp2<MatrixType>::apply_poly_to_eigs(const int it, real& homo, real& lumo)
00480 {
00481 assert(it >= 0);
00482 if (it == 0)
00483 {
00484 return;
00485 }
00486
00487 int poly;
00488 get_poly(it, poly);
00489
00490 if (poly == 1)
00491 {
00492 homo = 2 * homo - homo * homo;
00493 lumo = lumo * lumo;
00494 }
00495 else
00496 {
00497 homo = homo * homo;
00498 lumo = 2 * lumo - lumo * lumo;
00499 }
00500 }
00501
00502
00503 template<typename MatrixType>
00504 typename Purification_sp2<MatrixType>::real
00505 Purification_sp2<MatrixType>::compute_derivative(const int it, real x, real& DDf)
00506 {
00507 assert(it > 0);
00508
00509 real Df;
00510 real a;
00511 int poly;
00512
00513 a = x;
00514 Df = 1;
00515 DDf = 1;
00516
00517 for (int i = 1; i <= it; i++)
00518 {
00519 get_poly(i, poly);
00520
00521 if (poly == 1)
00522 {
00523 DDf = 2 * Df * Df + (2 * a) * DDf;
00524 Df *= 2 * a;
00525 a = a * a;
00526 }
00527 else
00528 {
00529 DDf = -2 * Df * Df + (2 - 2 * a) * DDf;
00530 Df *= 2 - 2 * a;
00531 a = 2 * a - a * a;
00532 }
00533 }
00534
00535
00536
00537 return Df;
00538 }
00539
00540
00541 #endif //HEADER_PURIFICATION_SP2