102 template <
typename T >
109 const Real& zero_cutoff,
112 int n_min,
int n_max,
118 const Real& delta_cycle,
119 const Real& gamma_factor)
123 int N5 = psi_all.size1();
131 const Subset&
s =
A.subset();
136 int N_eig_index = N_eig - 1;
137 int N_eig_minus_one = N_eig_index;
140 for(
n=0;
n <
N5;
n++) {
141 psi[
n][
s] = psi_all[N_eig_index][
n];
144 QDPIO::cout <<
"ritz: a" << std::endl;
147 if( N_eig_minus_one > 0 ) {
155 for(
n=0;
n <
N5;
n++) {
159 QDPIO::cout <<
"ritz: b" << std::endl;
168 if( N_eig_minus_one > 0 ) {
174 QDPIO::cout <<
"ritz: c" << std::endl;
178 for(
n=1;
n <
N5;
n++) {
179 mu += innerProductReal(
psi[
n], Apsi[
n],
s);
187 for(
n=0;
n <
N5;
n++) {
192 QDPIO::cout <<
"ritz: d" << std::endl;
202 QDPIO::cout <<
"ritz: e" << std::endl;
210 QDPIO::cout <<
"Starting Ritz: N_eig=" << N_eig <<
", mu = " <<
mu
211 <<
", g2_0 = " << g2_0 << std::endl;
216 bool minItersDoneP = n_min <= 0;
217 bool maxItersDoneP = n_max <= 0;
221 bool deltaCycleConvP;
227 rsd_a_sq = Rsd_a * Rsd_a;
228 rsd_r_sq = Rsd_r * Rsd_r *
mu *
mu;
234 CGConvP = toBool( g2 < rsd_r_sq ) || toBool(g2 < rsd_a_sq );
240 zeroFoundP = toBool( fabs(
mu) < zero_cutoff );
244 Double delta_cycle_err = delta_cycle;
246 if( Kalk_Sim ==
false )
251 convP = minItersDoneP && ( CGConvP || zeroFoundP );
271 convP = minItersDoneP && ( KSConvP || CGConvP || zeroFoundP );
298 if( N_eig_minus_one > 0 ) {
305 dd =
Double(innerProductReal(
p[0], Ap[0],
s));
306 for(
n=1;
n <
N5;
n++) {
307 dd +=
Double(innerProductReal(
p[
n], Ap[
n],
s));
321 if( toBool (
a >= s3) )
339 if( toBool( s2 > 0) )
341 s2 = half * (
one+s2);
343 s2 = -half * s3 / dd;
347 s2 = half * (
one-s2);
349 dd = -half * s3 / s2;
366 for(
n=0;
n <
N5;
n++)
368 psi[
n][
s] *= Real(ct);
370 Apsi[
n][
s] *= Real(ct);
371 Apsi[
n][
s] += Ap[
n] * Real(st);
375 for(
n=0;
n <
N5;
n++) {
376 Ap[
n][
s] = Apsi[
n] -
psi[
n]*lambda;
386 minItersDoneP = (
k >= n_min );
387 maxItersDoneP = (
k >= n_max );
390 rsd_a_sq = Rsd_a*Rsd_a;
391 rsd_r_sq = Rsd_r*Rsd_r*
mu*
mu;
397 CGConvP = toBool( g2 < rsd_a_sq ) || toBool( g2 < rsd_r_sq);
403 zeroFoundP = toBool( fabs(
mu) < zero_cutoff );
405 if( Kalk_Sim ==
false ) {
410 convP = minItersDoneP && ( maxItersDoneP || CGConvP || zeroFoundP );
436 deltaCycleConvP =
false;
437 Double g2_g0_ratio = g2 / g2_0;
438 KSConvP = deltaCycleConvP || toBool( g2_g0_ratio <=
Double(gamma_factor)) ;
439 convP = minItersDoneP && ( maxItersDoneP || CGConvP || KSConvP || zeroFoundP );
449 if( N_eig_minus_one > 0 ) {
454 dd = sqrt(norm2(
psi,
s));
456 for(
n=0;
n <
N5;
n++) {
457 psi[
n][
s] /= Real(dd);
465 QDPIO::cout <<
"Converged at iter=" <<
k <<
", lambda = " << lambda
466 <<
", rsd | mu | = " << rsd <<
", || g || = "
467 << sqrt(g2) <<
" || x || - 1 = " <<
d << std::endl;
471 QDPIO::cout <<
"KS: gamma = "<< gamma_factor <<
", || g ||^2/|| g_0 ||^2="
473 <<
", delta_cycle_err=" << delta_cycle_err << std::endl;
474 QDPIO::cout <<
"KS: CGConvP = " << CGConvP <<
", KSConvP = " << KSConvP <<
" deltaCycleConvP = " << deltaCycleConvP << std::endl;
478 QDPIO::cout <<
"ritz: some convergence" << std::endl;
484 mu = innerProductReal(
psi[0], Apsi[0],
s);
485 for(
n=1;
n <
N5;
n++) {
486 mu += innerProductReal(
psi[
n], Apsi[
n],
s);
491 QDPIO::cout <<
"Mu-s at convergence: old " << s1 <<
" vs " <<
mu << std::endl;
495 for(
n=0;
n <
N5;
n++) {
496 psi_all[N_eig_index][
n][
s] =
psi[
n];
500 for(
n=0;
n <
N5;
n++) {
501 Ap[
n][
s] = Apsi[
n] -
psi[
n]*lambda;
504 dd = sqrt(norm2(Ap,
s));
505 final_grad = Real(dd);
522 if( toBool( ct >
Double(1)) ) {
525 QDPIO::cout <<
"Restart at iter " <<
k <<
" since beta = " <<
b << std::endl;
526 for(
n=0;
n <
N5;
n++) {
533 for(
n=1;
n <
N5;
n++) {
537 xp = cmplx(Real(real(dxp)), Real(imag(dxp)));
540 for(
n=0;
n <
N5;
n++) {
549 if(
k % n_renorm == 0 ) {
552 if( N_eig_minus_one > 0 ) {
557 dd = sqrt(norm2(
psi,
s));
559 for(
n=0;
n <
N5;
n++) {
560 psi[
n][
s] /= Real(dd);
574 mu = innerProductReal(
psi[0], Apsi[0],
s);
575 for(
n=1;
n <
N5;
n++) {
576 mu += innerProductReal(
psi[
n], Apsi[
n],
s);
581 for(
n=0;
n <
N5;
n++) {
583 Ap[
n][
s] -=
psi[
n] * lambda;
590 if( N_eig_minus_one > 0 ) {
603 for(
n=1;
n <
N5;
n++) {
607 dxp = ct*(cmplx(g2,
Double(0)) - dxp);
608 xp = cmplx(Real(real(dxp)), Real(imag(dxp)));
610 for(
n=0;
n <
N5;
n++) {
611 p[
n][
s] += Ap[
n] * xp;
614 else if (ProjApsiP && N_eig_minus_one > 0)
617 if( N_eig_minus_one > 0 )
629 for(
n=0;
n <
N5;
n++) {
630 psi_all[N_eig_index][
n][
s] =
psi[
n];
634 final_grad = sqrt(g2);
635 QDPIO::cerr <<
"too many CG/Ritz iterations: n_count=" <<
n_count
636 <<
", rsd_r =" << sqrt(rsd_r_sq) <<
" rsd_a=" << Rsd_a <<
", ||g||=" << sqrt(g2) <<
", p2=" << p2
637 <<
", lambda" << lambda << std::endl;
647 multi2d<LatticeFermion>& psi_all,
651 const Real& zero_cutoff,
654 int n_min,
int n_max,
660 const Real& delta_cycle,
661 const Real& gamma_factor)
663 RitzArray_t(
A, lambda, psi_all, N_eig, Rsd_r, Rsd_a, zero_cutoff,
664 n_renorm, n_min, n_max,
MaxCG,
665 ProjApsiP,
n_count, final_grad,
666 Kalk_Sim, delta_cycle, gamma_factor);
Primary include file for CHROMA library code.
Linear Operator to arrays.
Gramm-Schmidt orthogonolization.
void GramSchmArray(multi2d< LatticeFermion > &psi, const int Npsi, const multi2d< LatticeFermion > &vec, const int Nvec, const Subset &sub)
Gramm-Schmidt orthogonolization.
void RitzArray_t(const LinearOperatorArray< T > &A, Real &lambda, multi2d< T > &psi_all, int N_eig, const Real &Rsd_r, const Real &Rsd_a, const Real &zero_cutoff, int n_renorm, int n_min, int n_max, int MaxCG, bool ProjApsiP, int &n_count, Real &final_grad, bool Kalk_Sim, const Real &delta_cycle, const Real &gamma_factor)
Minimizes the Ritz functional with a CG based algorithm.
BinaryReturn< C1, C2, FnInnerProduct >::Type_t innerProduct(const QDPSubType< T1, C1 > &s1, const QDPType< T2, C2 > &s2)
Asqtad Staggered-Dirac operator.
const WilsonTypeFermAct< multi1d< LatticeFermion > > Handle< const ConnectState > const multi1d< Real > enum InvType invType const multi1d< Real > int MaxCG
void Ritz(const LinearOperator< LatticeFermion > &A, Real &lambda, multi1d< LatticeFermion > &psi_all, int N_eig, const Real &Rsd_r, const Real &Rsd_a, const Real &zero_cutoff, int n_renorm, int n_min, int n_max, int MaxCG, bool ProjApsiP, int &n_count, Real &final_grad, bool Kalk_Sim, const Real &delta_cycle, const Real &gamma_factor)
multi1d< LatticeFermion > s(Ncb)
FloatingPoint< double > Double