36 #ifndef PUR_PURIFICATION_SCALED
37 #define PUR_PURIFICATION_SCALED
45 template<
typename Tmatrix>
59 real
const toleratedEigenvalError,
62 real
const toleratedSubspaceError,
81 void mInfo(std::ostream & file)
const;
83 void mTime(std::ostream & file)
const;
88 static int get_poly( real
const homo, real
const lumo ) {
91 return homo + lumo < 1;
96 real
const toleratedEigenvalError,
107 int const steps_left,
108 real
const acc_error,
109 real
const toleratedSubspaceError );
120 bool & homo_was_computed,
121 bool & lumo_was_computed );
163 template<
typename Tmatrix>
169 real const toleratedEigenvalError,
170 real const toleratedSubspaceError,
173 bool const use_scaling)
175 puri_steps(max_steps+1),
177 tolEigenvalError(toleratedEigenvalError),
178 tolSubspaceError(toleratedSubspaceError),
180 normTruncation(normForTruncation),
181 use_scaling(use_scaling),
183 lumo_is_computed(0) {
189 X.add_identity(-lmax);
190 X *= ((
real)1.0 / (lmin - lmax));
198 throw "homo interval empty in Purification_scaled constructor due to incorrect input.";
200 throw "lumo interval empty in Purification_scaled constructor due to incorrect input.";
201 homo = (homo - lmax) / (lmin - lmax);
202 lumo = (lumo - lmax) / (lmin - lmax);
210 steps_left, acc_error,
214 double wall_sec_thresh = t_thr.toc();
225 real const ONE = 1.0;
228 double wall_sec_square = t_sqr.
toc();
229 double wall_sec_total = t_tot.
toc();
232 X.trace(),
X2.trace(),
242 template<
typename Tmatrix>
255 alpha = 2 / ( 1+homo.
upp() );
258 alpha = 2 / ( 2-lumo.
low() );
263 X2 *= (
real)-(alpha*alpha);
268 X2 *= (
real)(alpha*alpha);
269 X2 += ((
real)2.0 * alpha * (1-alpha)) *
X;
270 X2.add_identity( (1-alpha)*(1-alpha) );
285 assert(acc_error >= 0);
287 steps_left, acc_error,
291 double wall_sec_thresh = t_thr.toc();
298 real const ONE = 1.0;
301 double wall_sec_square = t_sqr.
toc();
302 bool homo_was_computed =
false;
303 bool lumo_was_computed =
false;
308 double wall_sec_XmX2norm = t_norm.
toc();
310 double wall_sec_total = t_tot.
toc();
313 X.trace(),
X2.trace(),
321 if ( homo_was_computed ) {
325 if ( lumo_was_computed ) {
332 template<
typename Tmatrix>
336 real const toleratedEigenvalError,
345 while ((1 - homo.
low() > toleratedEigenvalError &&
346 lumo.
upp() > toleratedEigenvalError) && steps < max_steps) {
352 alpha = 2 / ( 1+homo.
upp() );
355 alpha = 2 / ( 2-lumo.
low() );
364 template<
typename Tmatrix>
367 int const steps_left,
368 real const acc_error,
369 real const toleratedSubspaceError ) {
370 real acceptable_error = toleratedSubspaceError - acc_error;
381 real smallest_factor = 1.01;
382 real steps_left_pessimistic = steps_left > smallest_factor ? steps_left : smallest_factor;
383 assert(acceptable_error >= 0);
384 return ( gap * acceptable_error/steps_left_pessimistic ) /
385 ( 1 + acceptable_error/steps_left_pessimistic );
401 template<
typename Tmatrix>
405 bool & homo_was_computed,
406 bool & lumo_was_computed ) {
407 homo_was_computed =
false;
408 lumo_was_computed =
false;
409 bool computeAccurately =
false;
410 bool homo_may_be_computed =
false;
411 bool lumo_may_be_computed =
false;
412 if ( homo.
low() > 0.5 && homo.
low() + lumo.
low() < 1 ) {
413 homo_may_be_computed =
true;
417 computeAccurately =
true;
419 if ( lumo.
upp() < 0.5 && homo.
upp() + lumo.
upp() > 1 ) {
420 lumo_may_be_computed =
true;
424 computeAccurately =
true;
430 real const XmX2ENIsSmallValue = 0.2475;
433 if ( computeAccurately )
434 normXmX2 = Tmatrix::diffIfSmall(
X,
X2,
438 normXmX2 = Tmatrix::diffIfSmall(
X,
X2,
444 real norm_upper_bound = Tmatrix::mixed_diff(
X,
X2, diffAcc);
445 norm_upper_bound = norm_upper_bound < normXmX2.
upp() ? norm_upper_bound : normXmX2.
upp();
450 real lowerBound = normXmX2.
low();
465 if (homo_may_be_computed && !lumo.
overlap(lumoTmp) ) {
467 if ( computeAccurately && normXmX2.
length() <= 2.1*diffAcc )
468 homo_was_computed =
true;
470 if (lumo_may_be_computed && !homo.
overlap(homoTmp)) {
472 if ( computeAccurately && normXmX2.
length() <= 2.1*diffAcc )
473 lumo_was_computed =
true;
478 template<
typename Tmatrix>
482 bool & homo_was_computed,
483 bool & lumo_was_computed ) {
484 homo_was_computed =
false;
485 lumo_was_computed =
false;
486 bool computeAccurately =
false;
487 bool homoIsComputable =
false;
488 bool lumoIsComputable =
false;
489 if ( homo.
upp() + lumo.
upp() < 1 ) {
492 if (homo.
low() > 0.5) {
493 homoIsComputable =
true;
495 computeAccurately =
true;
498 if ( homo.
low() + lumo.
low() > 1 ){
501 if (lumo.
upp() < 0.5) {
502 lumoIsComputable =
true;
504 computeAccurately =
true;
512 real const XmX2ENIsSmallValue = 0.2475;
515 if ( computeAccurately )
516 normXmX2 = Tmatrix::diffIfSmall(
X,
X2,
520 normXmX2 = Tmatrix::diffIfSmall(
X,
X2,
525 real lowerBound = normXmX2.
low();
531 if (homoIsComputable) {
535 if ( computeAccurately && normXmX2.
length() <= 2.1*diffAcc )
536 homo_was_computed =
true;
538 if (lumoIsComputable) {
543 if ( computeAccurately && normXmX2.
length() <= 2.1*diffAcc )
544 lumo_was_computed =
true;
549 template<
typename Tmatrix>
551 real total_error = 0;
552 for (
unsigned int step = 0; step <=
current_step; step++) {
555 total_error += normE / (gap - normE);
560 template<
typename Tmatrix>
562 for (
unsigned int step =
current_step; step >= 1; step-- )
565 template<
typename Tmatrix>
567 for (
unsigned int step =
current_step; step >= 1; step-- )
571 template<
typename Tmatrix>
579 template<
typename Tmatrix>
584 throw "Attempt to get homo and lumo eigenvalue intervals when purification has not converged.";
591 hoF = homo * (lmin - lmax) + lmax;
592 luF = lumo * (lmin - lmax) + lmax;
595 template<
typename Tmatrix>
598 s<<
"% PURIFICATION INFO IN MATLAB/OCTAVE FILE FORMAT"<<std::endl;
601 s<<
"n = "<<
n<<
";"<<std::endl;
606 s<<std::setprecision(15);
607 s<<
"% traceX and traceX2 \n";
609 for (
unsigned int ind = 0; ind <=
current_step; ++ind) {
616 for (
unsigned int ind = 0; ind <=
current_step; ++ind) {
620 s<<
"% chosenThresh actualThresh \n";
622 for (
unsigned int ind = 0; ind <=
current_step; ++ind) {
624 <<
puri_steps[ind].actual_thresh << std::endl;
627 s<<
"% nnzX nnzX2 \n";
629 for (
unsigned int ind = 0; ind <=
current_step; ++ind) {
634 s<<
"% lumo_low lumo_upp homo_low homo_upp \n";
635 s<<
"homo_lumo_eigs = [";
636 for (
unsigned int ind = 0; ind <=
current_step; ++ind) {
637 s<<
puri_steps[ind].eig_lumo_orig.low() <<
" "
640 <<
puri_steps[ind].eig_homo_orig.upp() << std::endl;
645 <<
"plot(1:nIter,(homo_lumo_eigs(:,3)+homo_lumo_eigs(:,4))/2,'-b',"
646 <<
"1:nIter,(homo_lumo_eigs(:,1)+homo_lumo_eigs(:,2))/2,'-r')\n"
647 <<
"legend('HOMO','LUMO'),xlabel('Iteration')\n"
649 <<
"for ind = 1:nIter \n"
650 <<
" plot([ind ind],[homo_lumo_eigs(ind,3) homo_lumo_eigs(ind,4)], '-b')\n"
651 <<
" plot([ind ind],[homo_lumo_eigs(ind,1) homo_lumo_eigs(ind,2)], '-r')\n"
653 <<
"axis([0 nIter 0 1])\n";
657 <<
"plot(1:nIter,100*nnzVals(:,1)/(n*n),'o-r',...\n"
658 <<
"1:nIter, 100*nnzVals(:,2)/(n*n),'x-b')\n"
659 <<
"legend('nnz(X)','nnz(X^2)'),xlabel('Iteration') \n"
660 <<
"ylabel('Percentage')\n"
661 <<
"axis([0 nIter 0 100])\n";
663 <<
"semilogy(1:nIter,thresh(:,1),'x-r',1:nIter,thresh(:,2),'o-b')\n"
664 <<
"xlabel('Iteration'), ylabel('Threshold')\n"
665 <<
"legend('Chosen threshold', 'Actual threshold')\n"
666 <<
"axis([0 nIter min(thresh(:,2))/10 max(thresh(:,1))*10])\n";
670 template<
typename Tmatrix>
673 s<<
"% PURIFICATION TIMINGS IN MATLAB/OCTAVE FILE FORMAT"<<std::endl;
676 s<<
"n = "<<
n<<
";"<<std::endl;
679 s<<
"wall_sec_thresh = [";
680 for (
unsigned int ind = 0; ind <=
current_step; ++ind) {
685 s<<
"wall_sec_square = [";
686 for (
unsigned int ind = 0; ind <=
current_step; ++ind) {
691 s<<
"wall_sec_XmX2norm = [";
692 for (
unsigned int ind = 0; ind <=
current_step; ++ind) {
697 s<<
"wall_sec_total = [";
698 for (
unsigned int ind = 0; ind <=
current_step; ++ind) {
704 s<<
"puriTime = [wall_sec_square wall_sec_thresh wall_sec_XmX2norm];\n";
705 s<<
"figure; bar(puriTime(:,1:3),'stacked')"<<std::endl<<
706 "legend('Matrix Square', 'Truncation', '||X-X^2||'),"<<
707 " xlabel('Iteration'), ylabel('Time (seconds)')"<<std::endl;
708 s<<
"figure; plot(wall_sec_total,'-x')"<<std::endl;
711 template<
typename Tmatrix>
721 if(gettimeofday(&tv, NULL) != 0)
722 throw std::runtime_error(
"Error in get_wall_seconds(), in gettimeofday().");
723 double seconds = tv.tv_sec + (double)tv.tv_usec / 1000000;
std::string getNormTypeString(normType nType)
Definition: matInclude.cc:55
bool overlap(Interval const &other) const
Definition: Interval.h:123
Treal upp() const
Definition: Interval.h:143
Definition: Purification_scaled.h:46
Treal length() const
Returns the length of the interval.
Definition: Interval.h:107
bool converged()
Definition: Purification_scaled.h:572
real const tolEigenvalError
Tolerated error in eigenvalues at convergence.
Definition: Purification_scaled.h:152
static Interval intersect(Interval const &A, Interval const &B)
Definition: Interval.h:51
static double get_wall_seconds()
Definition: bench_gemm_only.cc:47
double ticTime
Definition: Purification_scaled.h:718
Proxy structs used by the matrix API.
void mInfo(std::ostream &file) const
Definition: Purification_scaled.h:597
std::vector< Step< real > > puri_steps
Definition: Purification_scaled.h:150
double toc()
Definition: Purification_scaled.h:716
Tmatrix X2
Definition: Purification_scaled.h:148
Time()
Definition: Purification_scaled.h:714
void increase(Treal const value)
Increases interval with value in both directions.
Definition: Interval.h:131
bool empty() const
Definition: Interval.h:49
Definition: allocate.cc:30
bool use_scaling
Definition: Purification_scaled.h:158
static int get_poly(real const homo, real const lumo)
Definition: Purification_scaled.h:88
mat::Interval< real > const & eigFInt
Definition: Purification_scaled.h:151
void purify()
Definition: Purification_scaled.h:243
bool homo_is_computed
Definition: Purification_scaled.h:159
void puriStep(int poly)
Definition: Interval.h:228
Definition: matInclude.h:135
static real get_threshold(real const gap, int const steps_left, real const acc_error, real const toleratedSubspaceError)
Compute acceptable error due to truncation in an iteration based on the current band gap...
Definition: Purification_scaled.h:366
mat::normType normTruncation
Definition: Purification_scaled.h:157
unsigned int current_step
Definition: Purification_scaled.h:156
void propagate_homo_information()
Go through step vector and improve homo information in each step possible.
Definition: Purification_scaled.h:561
Treal low() const
Definition: Interval.h:142
void propagate_lumo_information()
Go through step vector and improve lumo information in each step possible.
Definition: Purification_scaled.h:566
static int estimated_steps_left(mat::Interval< real > eig_homo, mat::Interval< real > eig_lumo, real const toleratedEigenvalError, int const max_steps, bool const use_scaling)
Estimates number of remaining iterations.
Definition: Purification_scaled.h:334
Interval< Treal > sqrtInt(Interval< Treal > const &other)
Definition: Interval.h:217
Purification_scaled(Tmatrix &F_and_D, mat::Interval< real > const &eigFInt, mat::Interval< real > const &hoF, mat::Interval< real > const &luF, real const toleratedEigenvalError, real const toleratedSubspaceError, int const max_steps, mat::normType normForTruncation, bool const use_scaling)
Definition: Purification_scaled.h:165
void get_homo_lumo_intervals(mat::Interval< real > &hoF, mat::Interval< real > &luF)
Computed eigenvalues of F.
Definition: Purification_scaled.h:581
void mTime(std::ostream &file) const
Definition: Purification_scaled.h:672
Definition: Purification_scaled.h:712
ergo_real real
Definition: cubature_rules.h:33
real accumulated_error()
Total accumulated error due to removal of small matrix elements up to current iteration.
Definition: Purification_scaled.h:550
Tmatrix & X
Definition: Purification_scaled.h:147
static double get_wall_seconds()
Definition: Purification_scaled.h:719
void tic()
Definition: Purification_scaled.h:715
bool lumo_is_computed
Definition: Purification_scaled.h:160
Copyright(c) Emanuel Rubensson 2006.
Definition: matInclude.h:135
real const tolSubspaceError
Tolerated subspace error.
Definition: Purification_scaled.h:154
void intersect_always_non_empty(Interval const &other)
Definition: Interval.h:78
Tmatrix::real real
Definition: Purification_scaled.h:48
int n
Definition: Purification_scaled.h:149
Definition: Purification_scaled.h:44
normType
Definition: matInclude.h:135
void improve_homo_lumo_based_on_normXmX2(mat::Interval< real > &homo, mat::Interval< real > &lumo, bool &homo_was_computed, bool &lumo_was_computed)
Computes interval containing spectral norm of X-X^2 matrix.
Definition: Purification_scaled.h:403