34 #ifndef MAT_PURIFICATION_OLD
35 #define MAT_PURIFICATION_OLD
41 template<
class Treal,
class MAT>
43 const int nsteps,
const Treal frob_trunc = 0) {
45 for (
int step = 0; step < nsteps; step++) {
46 Treal tracediff = tmp.trace() - nshells;
48 tmp = (Treal)1.0 * D * D;
50 D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
54 tmp = (Treal)-1.0 * D * D + (Treal)2.0 * tmp;
55 D = (Treal)1.0 * tmp * tmp;
57 D.frob_thresh(frob_trunc);
61 template<
class Treal,
class MAT>
63 int& nsteps,
const Treal frob_trunc,
64 const int maxiter = 50) {
67 Treal tracediff = tmp.trace() - nshells;
71 tmp = (Treal)1.0 * D * D;
73 D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
77 tmp = (Treal)-1.0 * D * D + (Treal)2.0 * tmp;
78 D = (Treal)1.0 * tmp * tmp;
80 froberror = MAT::frob_diff(D, tmp);
81 D.frob_thresh(frob_trunc);
82 tracediff = D.trace() - nshells;
84 std::cout<<
"Iteration:"<<nsteps<<
" Froberror: "
85 <<std::setprecision(8)<<froberror<<
86 " Tracediff: "<<tracediff<<std::endl;
91 " without convergence", maxiter);
92 }
while (froberror > frob_trunc);
97 template<
class Treal,
class MAT>
98 void tc2(MAT& D,
int& iterations,
99 const MAT& H,
const int nshells,
100 const Treal trace_error_limit = 1e-2,
101 const Treal frob_trunc = 0,
103 const int maxiter = 100) {
107 Treal tracediff = tmp.trace() - nshells;
108 Treal tracediff_old = 2.0 * trace_error_limit;
115 D = (Treal)1.0 * tmp * tmp;
117 polys[iterations] = 0;
120 D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
122 polys[iterations] = 1;
124 D.frob_thresh(frob_trunc);
126 tracediff_old = tracediff;
127 tracediff = D.trace() - nshells;
129 if (iterations > maxiter)
130 throw AcceptableMaxIter(
"Purification reached maxiter without convergence", maxiter);
135 template<
class Treal,
class MAT>
137 const MAT& H,
const int nocc,
138 const Treal frob_trunc = 0,
140 const int maxiter = 100) {
141 assert(frob_trunc >= 0);
143 assert(maxiter >= 0);
148 Treal newtracediff = tmp.trace() - nocc;
150 Treal tracechange = 0;
157 tracediff = newtracediff;
159 D = (Treal)1.0 * tmp * tmp;
161 polys[iterations] = 0;
164 D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
166 polys[iterations] = 1;
168 D.frob_thresh(frob_trunc);
169 newtracediff = D.trace() - nocc;
170 tracechange = newtracediff - tracediff;
171 froberror = MAT::frob_diff(D, tmp);
174 std::cout<<
"Iteration:"<<iterations<<
" epsilon: "
175 <<std::setprecision(8)<<std::setw(12)<<froberror<<
176 " delta: "<<std::setw(12)<<tracediff<<
177 " tracechange: "<<std::setw(12)<<tracechange<<std::endl;
180 if (iterations > maxiter)
182 " without convergence", maxiter);
188 D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
190 polys[iterations] = 1;
193 D = (Treal)1.0 * tmp * tmp;
195 polys[iterations] = 0;
200 D.frob_thresh(frob_trunc);
201 tracediff = D.trace() - nocc;
204 tmp = (Treal)1.0 * D * D;
206 D = (Treal)-1.0 * tmp * tmp + (Treal)2.0 * D;
208 polys[iterations] = 0;
209 polys[iterations + 1] = 1;
214 tmp = (Treal)-1.0 * D * D + (Treal)2.0 * tmp;
215 D = (Treal)1.0 * tmp * tmp;
217 polys[iterations] = 1;
218 polys[iterations + 1] = 0;
221 froberror = MAT::frob_diff(D, tmp);
222 D.frob_thresh(frob_trunc);
223 tracediff = D.trace() - nocc;
225 std::cout<<
"Iteration:"<<iterations<<
" Froberror: "
226 <<std::setprecision(8)<<froberror<<
227 " Tracediff: "<<tracediff<<std::endl;
230 if (iterations > maxiter)
232 " in extra part without convergence", maxiter);
233 }
while (froberror > frob_trunc);
246 D.add_identity(-lmax);
247 D *= (1.0 / (lmin - lmax));
250 float threshtime = 0;
253 for (
int steps = 0; steps < nextra_steps; steps++) {
255 tracediff = tmp.trace() - nshells;
262 MAT::syrk(
'U',
false, -1.0, tmp, 2.0, D);
263 syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
266 D.frob_thresh(frob_trunc);
267 threshtime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
269 tracediff = tmp.trace() - nshells;
277 syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
281 MAT::syrk(
'U',
false, -1.0, tmp, 2.0, D);
282 syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
286 MAT::syrk(
'U',
false, -1.0, tmp, 2.0, D);
287 syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
292 syrktime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
296 D.frob_thresh(frob_trunc);
297 threshtime += ((float)(clock()-start))/(CLOCKS_PER_SEC);
300 std::cout<<std::endl<<
" total syrktime in tcp ="<<std::setw(5)
301 <<syrktime<<std::endl;
302 std::cout<<
" total threshtime in tcp ="<<std::setw(5)
303 <<threshtime<<std::endl;
void tc2_extra(MAT &D, const int nshells, const int nsteps, const Treal frob_trunc=0)
Definition: purification_old.h:42
static void syrk(const char *uplo, const char *trans, const int *n, const int *k, const T *alpha, const T *A, const int *lda, const T *beta, T *C, const int *ldc)
Definition: gblas.h:332
void tc2(MAT &D, int &iterations, const MAT &H, const int nshells, const Treal trace_error_limit=1e-2, const Treal frob_trunc=0, int *polys=NULL, const int maxiter=100)
Definition: purification_old.h:98
Treal template_blas_fabs(Treal x)
void tc2_extra_auto(MAT &D, const int nshells, int &nsteps, const Treal frob_trunc, const int maxiter=50)
Definition: purification_old.h:62
void tc2_auto(MAT &D, int &iterations, const MAT &H, const int nocc, const Treal frob_trunc=0, int *polys=NULL, const int maxiter=100)
Definition: purification_old.h:136
Treal template_blas_sqrt(Treal x)