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_SP2ACC
00040 #define HEADER_PURIFICATION_SP2ACC
00041
00042 #include "purification_general.h"
00043
00044
00045
00050 template<typename MatrixType>
00051 class Purification_sp2acc : 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_sp2acc() : PurificationGeneral<MatrixType>() {}
00065
00066 virtual void set_init_params()
00067 {
00068 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen SP2 ACCELERATED purification method");
00069 this->info.method = 2;
00070
00071 this->gammaStopEstim = 6 - 4 * template_blas_sqrt((real)2);
00072
00073 this->check_stopping_criterion_iter = -1;
00074 }
00075
00076 virtual void get_poly(const int it, int& poly, real& alpha);
00077 virtual void set_poly(const int it);
00078
00079 virtual void truncate_matrix(real& threshold, int it);
00080
00081 virtual void estimate_number_of_iterations(int& numit);
00082 virtual void purify_X(const int it);
00083 virtual void purify_bounds(const int it);
00084 virtual void save_other_iter_info(IterationInfo& iter_info, int it);
00085 virtual void apply_inverse_poly_vector(const int it, VectorTypeReal& bounds_from_it);
00086
00087 virtual void return_constant_C(const int it, real& Cval);
00088
00089 virtual real apply_poly(const int it, real x);
00090 virtual void apply_poly_to_eigs(const int it, real& homo, real& lumo);
00091 virtual real compute_derivative(const int it, real x, real& DDf);
00092
00093
00094
00095
00096 VectorTypeReal VecAlpha;
00097
00098
00099 static const real deltaTurnOffAcc;
00100 };
00101
00102 template<typename MatrixType>
00103 const typename Purification_sp2acc<MatrixType>::real
00104 Purification_sp2acc<MatrixType>::deltaTurnOffAcc = 0.01;
00105
00106
00107
00108 template<typename MatrixType>
00109 void Purification_sp2acc<MatrixType>::set_poly(const int it)
00110 {
00111 assert((int)this->VecPoly.size() > it);
00112
00113
00114 if (this->VecPoly[it] == -1)
00115 {
00116 real Xtrace = this->X.trace();
00117 real Xsqtrace = this->Xsq.trace();
00118
00119 real delta = deltaTurnOffAcc;
00120
00121
00122 if ((this->check_stopping_criterion_iter == -1) && (this->lumo_bounds.low() < delta) && (this->homo_bounds.low() < delta))
00123 {
00124 #ifdef DEBUG_OUTPUT
00125 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Outer bounds of homo and lumo are less then %e: ", delta);
00126 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "lumo_out = %e, homo_out = %e ", this->lumo_bounds.low(), this->homo_bounds.low());
00127 #endif
00128
00129 this->lumo_bounds = IntervalType(0, this->lumo_bounds.upp());
00130 this->homo_bounds = IntervalType(0, this->homo_bounds.upp());
00131
00132
00133 if (it == 1)
00134 {
00135 this->check_stopping_criterion_iter = it + 1;
00136 }
00137 else
00138 {
00139 this->check_stopping_criterion_iter = it + 2;
00140 }
00141 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Start to check stopping criterion on iteration %d", this->check_stopping_criterion_iter);
00142 }
00143
00144 if ((template_blas_fabs(Xsqtrace - this->nocc) <
00145 template_blas_fabs(2 * Xtrace - Xsqtrace - this->nocc))
00146 ||
00147 (it % 2
00148 &&
00149 (template_blas_fabs(Xsqtrace - this->nocc) ==
00150 template_blas_fabs(2 * Xtrace - Xsqtrace - this->nocc))
00151 ))
00152 {
00153 this->VecPoly[it] = 1;
00154 VecAlpha[it] = 2 / (2 - this->lumo_bounds.low());
00155 }
00156 else
00157 {
00158 this->VecPoly[it] = 0;
00159 VecAlpha[it] = 2 / (2 - this->homo_bounds.low());
00160 }
00161 #ifdef DEBUG_OUTPUT
00162 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Acceleration parameter: alpha = %lf", VecAlpha[it]);
00163 #endif
00164 }
00165 }
00166
00167
00168 template<typename MatrixType>
00169 void Purification_sp2acc<MatrixType>::get_poly(const int it, int& poly, real& alpha)
00170 {
00171 assert((int)this->VecPoly.size() > it);
00172 assert(this->VecPoly[it] != -1);
00173
00174
00175 assert(this->VecAlpha[it] != -1);
00176
00177 poly = this->VecPoly[it];
00178 alpha = VecAlpha[it];
00179 }
00180
00181
00182 template<typename MatrixType>
00183 void Purification_sp2acc<MatrixType>::purify_X(const int it)
00184 {
00185 #ifdef DEBUG_OUTPUT
00186 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Purify X...");
00187 #endif
00188 real alpha_tmp;
00189 int poly;
00190
00191 set_poly(it);
00192
00193 get_poly(it, poly, alpha_tmp);
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208 if (poly == 1)
00209 {
00210 if (alpha_tmp != 1)
00211 {
00212
00213
00214
00215
00216
00217
00218 this->X *= ((real)2.0 * (1 - alpha_tmp) * alpha_tmp);
00219 this->X.add_identity((real)(1 - alpha_tmp) * (1 - alpha_tmp));
00220 this->Xsq *= ((real)alpha_tmp * alpha_tmp);
00221 this->Xsq += this->X;
00222
00223
00224 }
00225 else
00226 {
00227
00228 }
00229 }
00230 else
00231 {
00232 if (alpha_tmp != 1)
00233 {
00234 this->X *= ((real) - 2.0 * alpha_tmp);
00235 this->Xsq *= ((real) - alpha_tmp * alpha_tmp);
00236 this->Xsq -= this->X;
00237 }
00238 else
00239 {
00240 this->Xsq *= ((real) - 1.0);
00241 this->X *= (real)2.0;
00242 this->Xsq += this->X;
00243
00244
00245 }
00246 }
00247
00248 this->Xsq.transfer(this->X);
00249 }
00250
00251
00252 template<typename MatrixType>
00253 void Purification_sp2acc<MatrixType>::purify_bounds(const int it)
00254 {
00255 #ifdef DEBUG_OUTPUT
00256 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Change homo and lumo bounds according to the chosen polynomial VecPoly = %d", this->VecPoly[it]);
00257 #endif
00258
00259 real homo_low, homo_upp, lumo_upp, lumo_low;
00260 real alpha_tmp;
00261 int poly;
00262
00263 get_poly(it, poly, alpha_tmp);
00264
00265 if (poly == 1)
00266 {
00267
00268 homo_low = 2 * alpha_tmp * this->homo_bounds.low() - alpha_tmp * alpha_tmp * this->homo_bounds.low() * this->homo_bounds.low();
00269 homo_upp = 2 * alpha_tmp * this->homo_bounds.upp() - alpha_tmp * alpha_tmp * this->homo_bounds.upp() * this->homo_bounds.upp();
00270 lumo_low = (1 - alpha_tmp + alpha_tmp * this->lumo_bounds.low());
00271 lumo_low *= lumo_low;
00272 lumo_upp = (1 - alpha_tmp + alpha_tmp * this->lumo_bounds.upp());
00273 lumo_upp *= lumo_upp;
00274
00275 this->homo_bounds = IntervalType(homo_low, homo_upp);
00276 this->lumo_bounds = IntervalType(lumo_low, lumo_upp);
00277 }
00278 else
00279 {
00280
00281 lumo_low = 2 * alpha_tmp * this->lumo_bounds.low() - alpha_tmp * alpha_tmp * this->lumo_bounds.low() * this->lumo_bounds.low();
00282 lumo_upp = 2 * alpha_tmp * this->lumo_bounds.upp() - alpha_tmp * alpha_tmp * this->lumo_bounds.upp() * this->lumo_bounds.upp();
00283 homo_low = (1 - alpha_tmp + alpha_tmp * this->homo_bounds.low());
00284 homo_low *= homo_low;
00285 homo_upp = (1 - alpha_tmp + alpha_tmp * this->homo_bounds.upp());
00286 homo_upp *= homo_upp;
00287
00288 this->homo_bounds = IntervalType(homo_low, homo_upp);
00289 this->lumo_bounds = IntervalType(lumo_low, lumo_upp);
00290 }
00291
00292 IntervalType zero_one(0, 1);
00293 this->homo_bounds.intersect(zero_one);
00294 this->lumo_bounds.intersect(zero_one);
00295
00296 #ifdef DEBUG_OUTPUT
00297 if (this->homo_bounds.empty())
00298 {
00299 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval homo_bounds is empty.");
00300 }
00301 if (this->lumo_bounds.empty())
00302 {
00303 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Interval lumo_bounds is empty.");
00304 }
00305
00306
00307 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "1-homo: [ %g , %g ],", this->homo_bounds.low(), this->homo_bounds.upp());
00308 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "lumo: [ %g , %g ].", this->lumo_bounds.low(), this->lumo_bounds.upp());
00309 #endif
00310 }
00311
00312
00313
00314
00315 template<typename MatrixType>
00316 void Purification_sp2acc<MatrixType>::return_constant_C(const int it, real& Cval)
00317 {
00318 assert(it >= 1);
00319
00320 real alpha1 = VecAlpha[it - 1];
00321 real alpha2 = VecAlpha[it];
00322
00323 Cval = -1;
00324
00325 #ifdef DEBUG_OUTPUT
00326 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "alpha1 = %.4e , alpha2 = %.4e", alpha1, alpha2);
00327 #endif
00328
00329 if (it < 2)
00330 {
00331 return;
00332 }
00333
00334 if (((alpha1 == 1) && (alpha2 == 1)))
00335 {
00336 #ifdef DEBUG_OUTPUT
00337 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Check SP2 stopping criterion.");
00338 #endif
00339 Cval = C_SP2;
00340 return;
00341 }
00342 #ifdef DEBUG_OUTPUT
00343 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Do not check stopping criterion.");
00344 #endif
00345
00346
00347
00348
00349 #if 0
00350
00351 real homo_low = this->info.Iterations[it - 2].homo_bound_low;
00352 real homo_upp = this->info.Iterations[it - 2].homo_bound_upp;
00353 real lumo_low = this->info.Iterations[it - 2].lumo_bound_low;
00354 real lumo_upp = this->info.Iterations[it - 2].lumo_bound_upp;
00355 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "lumo [%.16e, %e], homo [%.16e, %e]", lumo_low, lumo_upp, homo_low, homo_upp);
00356
00357
00358 if ((homo_upp > 0.5) || (lumo_upp > 0.5))
00359 {
00360 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Inner bounds interval do not contain 0.5. Skip iteration. ");
00361 Cval = -1;
00362 return;
00363 }
00364
00365 a = std::max(lumo_low, homo_low);
00366
00367 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Chosen a = %g", a);
00368
00369 if (a <= 0)
00370 {
00371
00372 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Cannot compute constant C since a = %g when alpha1 = %g"
00373 " and alpha2 = %g", a, alpha1, alpha2);
00374 Cval = -1;
00375 return;
00376 }
00377
00378 real C1;
00379 C1 = -7.88 + 11.6 * alpha1 + 0.71 * alpha2;
00380 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Local maximum of g1: %g", C1);
00381
00382 Cval = C1;
00383
00384 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "*********** C = %g ************", Cval);
00385 #endif
00386 }
00387
00388
00389
00390
00391 template<typename MatrixType>
00392 void Purification_sp2acc<MatrixType>::truncate_matrix(real& threshold, int it)
00393 {
00394 real allowed_error = this->error_per_it;
00395
00396 threshold = 0;
00397 if (allowed_error == 0)
00398 {
00399 return;
00400 }
00401
00402 assert((int)this->VecGap.size() > it);
00403 #ifdef DEBUG_OUTPUT
00404 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Truncate matrix X: ");
00405 #endif
00406 real tau;
00407
00408 if (this->VecGap[it] > 0)
00409 {
00410 #ifdef DEBUG_OUTPUT
00411 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "VecGap[ %d ] = %e , ", it, this->VecGap[it]);
00412 #endif
00413 tau = (allowed_error * this->VecGap[it]) / (1 + allowed_error);
00414 }
00415 else
00416 {
00417 tau = allowed_error * 0.01;
00418 }
00419
00420
00421
00422 #ifdef USE_CHUNKS_AND_TASKS
00423 threshold = (this->X).thresh_frob(tau);
00424 #else
00425 threshold = (this->X).thresh(tau, this->normPuriTrunc);
00426 #endif
00427
00428 #ifdef DEBUG_OUTPUT
00429 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "tau = %e", tau);
00430 #endif
00431 }
00432
00433
00434
00435
00436
00437
00438 template<typename MatrixType>
00439 void Purification_sp2acc<MatrixType>::estimate_number_of_iterations(int& estim_num_iter)
00440 {
00441 int it = 1;
00442 int maxit_tmp = this->maxit + this->additional_iterations;
00443 real x, y, x_out, y_out;
00444 real alpha_tmp;
00445 real epsilon = this->get_epsilon();
00446
00447 int max_size = maxit_tmp + 1 + this->additional_iterations + 2;
00448
00449 this->VecPoly.clear();
00450 this->VecPoly.resize(max_size, -1);
00451
00452 this->VecGap.clear();
00453 this->VecGap.resize(max_size, -1);
00454
00455 VecAlpha.clear();
00456 VecAlpha.resize(max_size, -1);
00457
00458
00459 x = this->lumo_bounds.upp();
00460 y = this->homo_bounds.upp();
00461
00462
00463 if (1 - x - y <= 0)
00464 {
00465 #ifdef DEBUG_OUTPUT
00466 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "VecGap cannot be computed. Set estimated number of iteration to the maxit.");
00467 #endif
00468 estim_num_iter = this->maxit;
00469 return;
00470 }
00471
00472 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "INIT LUMO: [ %.12lf , %.12lf ]", (double)this->lumo_bounds.low(), (double)this->lumo_bounds.upp());
00473 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "INIT HOMO: [ %.12lf , %.12lf ]", (double)this->homo_bounds.low(), (double)this->homo_bounds.upp());
00474
00475
00476 x_out = this->lumo_bounds.low();
00477 y_out = this->homo_bounds.low();
00478
00479
00480
00481 real delta = deltaTurnOffAcc;
00482
00483 this->VecPoly[0] = -1;
00484 this->VecGap[0] = 1 - x - y;
00485
00486 estim_num_iter = -1;
00487
00488 while (it <= maxit_tmp)
00489 {
00490
00491 if ((x > y) || (it % 2 && (x == y)))
00492 {
00493 alpha_tmp = 2 / (2 - x_out);
00494
00495 x = (1 - alpha_tmp + alpha_tmp * x);
00496 x *= x;
00497 y = 2 * alpha_tmp * y - alpha_tmp * alpha_tmp * y * y;
00498
00499 x_out = (1 - alpha_tmp + alpha_tmp * x_out);
00500 x_out *= x_out;
00501 y_out = 2 * alpha_tmp * y_out - alpha_tmp * alpha_tmp * y_out * y_out;
00502
00503 this->VecPoly[it] = 1;
00504 }
00505 else
00506 {
00507 alpha_tmp = 2 / (2 - y_out);
00508
00509 x = 2 * alpha_tmp * x - alpha_tmp * alpha_tmp * x * x;
00510 y = (1 - alpha_tmp + alpha_tmp * y);
00511 y *= y;
00512
00513 x_out = 2 * alpha_tmp * x_out - alpha_tmp * alpha_tmp * x_out * x_out;
00514 y_out = (1 - alpha_tmp + alpha_tmp * y_out);
00515 y_out *= y_out;
00516
00517 this->VecPoly[it] = 0;
00518 }
00519
00520 VecAlpha[it] = alpha_tmp;
00521 this->VecGap[it] = 1 - x - y;
00522
00523
00524 if ((x_out < delta) && (y_out < delta) && (this->check_stopping_criterion_iter == -1))
00525 {
00526
00527 if (it == 1)
00528 {
00529 this->check_stopping_criterion_iter = it + 1;
00530 }
00531 else
00532 {
00533 this->check_stopping_criterion_iter = it + 2;
00534 }
00535 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "Start to check stopping criterion on iteration %d", this->check_stopping_criterion_iter);
00536 x_out = 0;
00537 y_out = 0;
00538 }
00539
00540
00541
00542 if ((estim_num_iter == -1) && (it >= this->check_stopping_criterion_iter) &&
00543 (x - x * x < epsilon) && (y - y * y < epsilon) &&
00544 (this->VecPoly[it] != this->VecPoly[it - 1]))
00545 {
00546 estim_num_iter = it;
00547 maxit_tmp = it + this->additional_iterations;
00548
00549 if (this->normPuriStopCrit == mat::frobNorm)
00550 {
00551 estim_num_iter += 2;
00552 maxit_tmp += 2;
00553 }
00554 }
00555
00556 ++it;
00557 }
00558
00559 if ((estim_num_iter == -1) && (it == maxit_tmp + 1))
00560 {
00561 #ifdef DEBUG_OUTPUT
00562 do_output(LOG_CAT_INFO, LOG_AREA_DENSFROMF, "maxit = %d number of iterations is reached in estimate_number_of_iteration()", this->maxit);
00563 #endif
00564 estim_num_iter = this->maxit;
00565 }
00566
00567 this->VecPoly.resize(maxit_tmp + 1);
00568 this->VecGap.resize(maxit_tmp + 1);
00569 this->VecAlpha.resize(maxit_tmp + 1);
00570 }
00571
00572
00573
00574
00575 template<typename MatrixType>
00576 void Purification_sp2acc<MatrixType>::save_other_iter_info(IterationInfo& iter_info, int it)
00577 {
00578 assert((int)this->VecPoly.size() > it);
00579 assert((int)this->VecGap.size() > it);
00580 assert((int)this->VecAlpha.size() > it);
00581
00582 iter_info.poly = this->VecPoly[it];
00583 iter_info.gap = this->VecGap[it];
00584 iter_info.alpha = VecAlpha[it];
00585 }
00586
00587
00588
00589
00590 template<typename MatrixType>
00591 void Purification_sp2acc<MatrixType>::apply_inverse_poly_vector(const int it, VectorTypeReal& bounds_from_it)
00592 {
00593 real tau;
00594 real alpha_tmp;
00595 int poly;
00596
00597 for (int i = it; i >= 1; i--)
00598 {
00599 tau = 0;
00600
00601 get_poly(i, poly, alpha_tmp);
00602
00603 if (poly == 1)
00604 {
00605 bounds_from_it[0] = template_blas_sqrt(bounds_from_it[0]);
00606 bounds_from_it[0] = (bounds_from_it[0] - 1 + alpha_tmp) / alpha_tmp - tau;
00607 bounds_from_it[1] = template_blas_sqrt(bounds_from_it[1]);
00608 bounds_from_it[1] = (bounds_from_it[1] - 1 + alpha_tmp) / alpha_tmp + tau;
00609
00610 bounds_from_it[2] = bounds_from_it[2] / (1 + template_blas_sqrt(1 - bounds_from_it[2]));
00611 bounds_from_it[2] = bounds_from_it[2] / alpha_tmp - tau;
00612 bounds_from_it[3] = bounds_from_it[3] / (1 + template_blas_sqrt(1 - bounds_from_it[3]));
00613 bounds_from_it[3] = bounds_from_it[3] / alpha_tmp + tau;
00614 }
00615 else
00616 {
00617 bounds_from_it[0] = bounds_from_it[0] / (1 + template_blas_sqrt(1 - bounds_from_it[0]));
00618 bounds_from_it[0] = bounds_from_it[0] / alpha_tmp - tau;
00619 bounds_from_it[1] = bounds_from_it[1] / (1 + template_blas_sqrt(1 - bounds_from_it[1]));
00620 bounds_from_it[1] = bounds_from_it[1] / alpha_tmp + tau;
00621
00622 bounds_from_it[2] = template_blas_sqrt(bounds_from_it[2]);
00623 bounds_from_it[2] = (bounds_from_it[2] - 1 + alpha_tmp) / alpha_tmp - tau;
00624 bounds_from_it[3] = template_blas_sqrt(bounds_from_it[3]);
00625 bounds_from_it[3] = (bounds_from_it[3] - 1 + alpha_tmp) / alpha_tmp + tau;
00626 }
00627 }
00628 }
00629
00630
00631 template<typename MatrixType>
00632 typename Purification_sp2acc<MatrixType>::real
00633 Purification_sp2acc<MatrixType>::apply_poly(const int it, real x)
00634 {
00635 assert(it >= 0);
00636 if (it == 0)
00637 {
00638 return x;
00639 }
00640
00641 real fx;
00642 int poly;
00643 real alpha_tmp;
00644 get_poly(it, poly, alpha_tmp);
00645
00646 if (poly == 1)
00647 {
00648 fx = (1 - alpha_tmp + alpha_tmp * x) * (1 - alpha_tmp + alpha_tmp * x);
00649 }
00650 else
00651 {
00652 fx = 2 * alpha_tmp * x - alpha_tmp * alpha_tmp * x * x;
00653 }
00654
00655 return fx;
00656 }
00657
00658
00659 template<typename MatrixType>
00660 void Purification_sp2acc<MatrixType>::apply_poly_to_eigs(const int it, real& homo, real& lumo)
00661 {
00662 assert(it >= 0);
00663 if (it == 0)
00664 {
00665 return;
00666 }
00667
00668 int poly;
00669 real alpha_tmp;
00670 get_poly(it, poly, alpha_tmp);
00671
00672 if (poly == 1)
00673 {
00674 homo = 2 * alpha_tmp * homo - alpha_tmp * alpha_tmp * homo * homo;
00675 lumo = (1 - alpha_tmp + alpha_tmp * lumo) * (1 - alpha_tmp + alpha_tmp * lumo);
00676 }
00677 else
00678 {
00679 homo = (1 - alpha_tmp + alpha_tmp * homo) * (1 - alpha_tmp + alpha_tmp * homo);
00680 lumo = 2 * alpha_tmp * lumo - alpha_tmp * alpha_tmp * lumo * lumo;
00681 }
00682 }
00683
00684
00685 template<typename MatrixType>
00686 typename Purification_sp2acc<MatrixType>::real
00687 Purification_sp2acc<MatrixType>::compute_derivative(const int it, real x, real& DDf)
00688 {
00689 assert(it > 0);
00690
00691 real Df;
00692 real temp, a, b;
00693 int poly;
00694 real alpha;
00695
00696 a = x;
00697 Df = 1;
00698 DDf = -1;
00699
00700 for (int i = 1; i <= it; i++)
00701 {
00702 temp = a;
00703
00704 get_poly(i, poly, alpha);
00705
00706 if (poly == 1)
00707 {
00708 a = ((1 - alpha) + alpha * temp) * ((1 - alpha) + alpha * temp);
00709 b = 2 * alpha * ((1 - alpha) + alpha * temp);
00710 }
00711 else
00712 {
00713 a = 2 * alpha * temp - alpha * alpha * temp * temp;
00714 b = 2 * alpha - 2 * alpha * alpha * temp;
00715 }
00716 Df *= b;
00717 }
00718
00719 return Df;
00720 }
00721
00722
00723 #endif //HEADER_PURIFICATION_SP2ACC