Skip to content

Commit

Permalink
Merge pull request #15 from etmc/quda_work_hmc_clover_massprec
Browse files Browse the repository at this point in the history 
Interface with new 'tm_rho' parameter for the preconditioned twisted …
  • Loading branch information
@Marcogarofalo
Marcogarofalo authored Aug 6, 2021
2 parents 1f9090c + 054916a commit 87276e3
Show file tree
Hide file tree
Showing 2 changed files with 85 additions and 44 deletions.
41 changes: 30 additions & 11 deletions quda_interface.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1127,9 +1127,6 @@ void M_quda(spinor * const P, spinor * const Q) {
_initQuda();

inv_param.kappa = g_kappa;
// IMPORTANT: use opposite TM flavor since gamma5 -> -gamma5 (until LXLYLZT prob. resolved)
inv_param.mu = -g_mu;
inv_param.epsilon = 0.0;

inv_param.twist_flavor = QUDA_TWIST_SINGLET;
inv_param.preserve_source = QUDA_PRESERVE_SOURCE_YES;
Expand All @@ -1143,10 +1140,9 @@ void M_quda(spinor * const P, spinor * const Q) {
1,//even_odd
1e-12,
1000,
0, 0);
//inv_param.dslash_type = QUDA_TWISTED_MASS_DSLASH;
1, 0);

inv_param.solution_type = QUDA_MATPCDAG_MATPC_SOLUTION;
inv_param.solution_type = QUDA_MATPC_SOLUTION;
inv_param.solve_type = QUDA_NORMOP_PC_SOLVE;
inv_param.matpc_type = QUDA_MATPC_ODD_ODD_ASYMMETRIC;
inv_param.dagger = QUDA_DAG_NO;
Expand All @@ -1159,8 +1155,6 @@ void M_quda(spinor * const P, spinor * const Q) {
reorder_spinor_eo_toQuda( (double*)spinorIn, inv_param.cpu_prec, 0 ,1);

// multiply


MatQuda(spinorOut, spinorIn, &inv_param);

// reorder spinor
Expand Down Expand Up @@ -1189,10 +1183,16 @@ void _setOneFlavourSolverParam(const double kappa, const double c_sw, const doub
// IMPORTANT: use opposite TM flavor since gamma5 -> -gamma5 (until LXLYLZT prob. resolved)
inv_param.mu = -mu/2./kappa;
inv_param.clover_coeff = c_sw*kappa;
inv_param.clover_rho = 0.0;
inv_param.tm_rho = -g_mu3/2./kappa;
inv_param.compute_clover_inverse = 1;
inv_param.compute_clover = 1;
}
else if( fabs(mu) > 0.0 ) {
inv_param.clover_coeff = 0.0;
inv_param.clover_rho = 0.0;
inv_param.tm_rho = 0.0;

inv_param.twist_flavor = QUDA_TWIST_SINGLET;
inv_param.dslash_type = QUDA_TWISTED_MASS_DSLASH;
inv_param.matpc_type = QUDA_MATPC_EVEN_EVEN_ASYMMETRIC;
Expand All @@ -1201,16 +1201,35 @@ void _setOneFlavourSolverParam(const double kappa, const double c_sw, const doub
inv_param.mu = -mu/2./kappa;
}
else if( c_sw > 0.0 ) {
inv_param.twist_flavor = QUDA_TWIST_NO;
inv_param.dslash_type = QUDA_CLOVER_WILSON_DSLASH;
// when we are dealing with the 'rho' mass preconditioning parameter
// the way this is implemented in QUDA is not consistent with the
// way it is implemented in tmLQCD
// To get agreement, we use the twisted clover operator also
// for Wilson clover fermions in this case.
if( fabs(g_mu3) > 0.0 ){
inv_param.twist_flavor = QUDA_TWIST_SINGLET;
inv_param.dslash_type = QUDA_TWISTED_CLOVER_DSLASH;
inv_param.mu = 0.0;
inv_param.clover_rho = 0.0;
inv_param.tm_rho = -g_mu3/2./kappa;
} else {
inv_param.twist_flavor = QUDA_TWIST_NO;
inv_param.dslash_type = QUDA_CLOVER_WILSON_DSLASH;
}
inv_param.matpc_type = QUDA_MATPC_EVEN_EVEN;
inv_param.solution_type = QUDA_MAT_SOLUTION;
inv_param.clover_order = QUDA_PACKED_CLOVER_ORDER;
inv_param.clover_coeff = c_sw*kappa;
inv_param.clover_rho = 0.0;
inv_param.compute_clover_inverse = 1;
inv_param.compute_clover = 1;
}
else {
inv_param.mu = 0.0;
inv_param.clover_coeff = 0.0;
inv_param.clover_rho = 0.0;
inv_param.tm_rho = 0.0;

inv_param.twist_flavor = QUDA_TWIST_NO;
inv_param.dslash_type = QUDA_WILSON_DSLASH;
if( single_parity_solve ){
Expand All @@ -1220,7 +1239,6 @@ void _setOneFlavourSolverParam(const double kappa, const double c_sw, const doub
// twisted mass
inv_param.dslash_type = QUDA_TWISTED_MASS_DSLASH;
inv_param.twist_flavor = QUDA_TWIST_SINGLET;
inv_param.mu = 0.0;
}
inv_param.matpc_type = QUDA_MATPC_EVEN_EVEN;
inv_param.solution_type = QUDA_MAT_SOLUTION;
Expand Down Expand Up @@ -1745,6 +1763,7 @@ int invert_eo_MMd_quda(spinor * const out,

// now we invert \hat{M}^{-} to get the inverse of \hat{Q}^{-} in the end
inv_param.mu = -inv_param.mu;
inv_param.tm_rho = -inv_param.tm_rho;

if(solver_flag == MG){
// flip the sign of the coarse operator and update the setup
Expand Down
88 changes: 55 additions & 33 deletions solver/monomial_solve.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,7 @@ int solve_degenerate(spinor * const P, spinor * const Q, solver_params_t solver_
if ( solver_params.external_inverter == QUDA_INVERTER){

int QmQp = (f == Qsw_pm_psi || f == Qtm_pm_psi);

gamma5(temp[0], Q, VOLUME/2);
iteration_count = invert_eo_MMd_quda(P, //spinor * const Odd_new,
temp[0],
Expand All @@ -140,67 +140,89 @@ int solve_degenerate(spinor * const P, spinor * const Q, solver_params_t solver_
//// init_solver_field(&tempE, VOLUMEPLUSRAND/2, 2);
//// //point like source only if mpi=1
//// for(int x =0; x < (VOLUMEPLUSRAND/2);x++){
//// tempE[0][x].s0.c0=((double)rand())/RAND_MAX; //0.0;
//// tempE[0][x].s0.c1=((double)rand())/RAND_MAX; //0.0;
//// tempE[0][x].s0.c2=((double)rand())/RAND_MAX; //0.0;
//// tempE[0][x].s1.c0=((double)rand())/RAND_MAX; //0.0;
//// tempE[0][x].s1.c1=((double)rand())/RAND_MAX; //0.0;
//// tempE[0][x].s1.c2=((double)rand())/RAND_MAX; //0.0;
//// tempE[0][x].s2.c0=((double)rand())/RAND_MAX; //0.0;
//// tempE[0][x].s2.c1=((double)rand())/RAND_MAX; //0.0;
//// tempE[0][x].s2.c2=((double)rand())/RAND_MAX; //0.0;
//// tempE[0][x].s3.c0=((double)rand())/RAND_MAX; //0.0;
//// tempE[0][x].s3.c1=((double)rand())/RAND_MAX; //0.0;
//// tempE[0][x].s3.c2=((double)rand())/RAND_MAX; //0.0;
//// // tempE[0][x].s0.c0=0.0;
//// // tempE[0][x].s0.c1=0.0;
//// // tempE[0][x].s0.c2=0.0;
//// // tempE[0][x].s1.c0=0.0;
//// // tempE[0][x].s1.c1=0.0;
//// // tempE[0][x].s1.c2=0.0;
//// // tempE[0][x].s2.c0=0.0;
//// // tempE[0][x].s2.c1=0.0;
//// // tempE[0][x].s2.c2=0.0;
//// // tempE[0][x].s3.c0=0.0;
//// // tempE[0][x].s3.c1=0.0;
//// // tempE[0][x].s3.c2=0.0;
////
//// // random
//// tempE[0][x].s0.c0=((double)rand())/RAND_MAX;
//// tempE[0][x].s0.c1=((double)rand())/RAND_MAX;
//// tempE[0][x].s0.c2=((double)rand())/RAND_MAX;
//// tempE[0][x].s1.c0=((double)rand())/RAND_MAX;
//// tempE[0][x].s1.c1=((double)rand())/RAND_MAX;
//// tempE[0][x].s1.c2=((double)rand())/RAND_MAX;
//// tempE[0][x].s2.c0=((double)rand())/RAND_MAX;
//// tempE[0][x].s2.c1=((double)rand())/RAND_MAX;
//// tempE[0][x].s2.c2=((double)rand())/RAND_MAX;
//// tempE[0][x].s3.c0=((double)rand())/RAND_MAX;
//// tempE[0][x].s3.c1=((double)rand())/RAND_MAX;
//// tempE[0][x].s3.c2=((double)rand())/RAND_MAX;
//// }
//// // set something other than component (0,0) to 1.0
//// // tempE[0][0].s0.c0=1.0;
//// // tempE[0][0].s1.c0=1.0;
//// // tempE[0][0].s2.c0=1.0;
//// // tempE[0][0].s3.c0=1.0;
//// // tempE[0][0].s0.c1=1.0;
//// // tempE[0][0].s0.c1=1.0;
//// // tempE[0][0].s0.c2=1.0;
//// // tempE[0][0].s1.c1=1.0;
//// // tempE[0][0].s2.c1=1.0;
//// // tempE[0][0].s3.c1=1.0;

//// // just in case: copy the source
//// assign(tempE[1], tempE[0], VOLUMEPLUSRAND/2);

//// // Qhat_oo = gamma_5 * M^hat_oo * in
//// // calling MatQuda in M_quda should apply M^hat_oo
//// M_quda(P, tempE[0]); // quda changes the source
//// mul_gamma5(P, VOLUME/2);

//// // for now, use Qtm_plus_psi explicitly, this can be generalised later and
//// // placed into a proper driver for testing the various operators between tmLQCD
//// // and QUDA
//// // If I were to use f() here, it would be applying Q^hat_plus Q^hat_minus, which complicates matters
//// // as it involves more hops
//// Qtm_plus_psi(temp[0], tempE[1]);
//// //Qsw_pm_psi(temp[0], tempE[0]);
//// if( f == Qtm_pm_psi ){
//// Mtm_plus_psi(temp[0], tempE[1]);
//// } else if ( f == Qsw_pm_psi ){
//// Msw_plus_psi(temp[0], tempE[1]);
//// }

//// // almost certainly we need to account for the gamma basis
//// for (int ix=0; ix < (VOLUME/2); ix++){
//// spinor *hp=((spinor*)temp[0]) + ix;
//// spinor *dp=((spinor*)P) + ix;
//// double r=creal((hp)->s0.c0)-creal((dp)->s0.c0);
//// printf("ix=%d, r=%.3e\n"
//// "tmLQCD=(%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e)\n"
//// " quda=(%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e)\n",
//// "re tmLQCD=(%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e)\n"
//// "re quda=(%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e)\n"
//// "im tmLQCD=(%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e)\n"
//// "im quda=(%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e), (%.3e,%.3e,%.3e)\n",
//// ix,r,
//// creal((hp)->s0.c0), creal((hp)->s0.c1), creal((hp)->s0.c2),
//// creal((hp)->s1.c0), creal((hp)->s1.c1), creal((hp)->s1.c2),
//// creal((hp)->s2.c0), creal((hp)->s2.c1), creal((hp)->s2.c2),
//// creal((hp)->s3.c0), creal((hp)->s3.c1), creal((hp)->s3.c2),

//// creal((dp)->s0.c0), creal((dp)->s0.c1), creal((dp)->s0.c2),
//// creal((dp)->s1.c0), creal((dp)->s1.c1), creal((dp)->s1.c2),
//// creal((dp)->s2.c0), creal((dp)->s2.c1), creal((dp)->s2.c2),
//// creal((dp)->s3.c0), creal((dp)->s3.c1), creal((dp)->s3.c2)
//// creal((dp)->s3.c0), creal((dp)->s3.c1), creal((dp)->s3.c2),
////
//// cimag((hp)->s0.c0), cimag((hp)->s0.c1), cimag((hp)->s0.c2),
//// cimag((hp)->s1.c0), cimag((hp)->s1.c1), cimag((hp)->s1.c2),
//// cimag((hp)->s2.c0), cimag((hp)->s2.c1), cimag((hp)->s2.c2),
//// cimag((hp)->s3.c0), cimag((hp)->s3.c1), cimag((hp)->s3.c2),

//// cimag((dp)->s0.c0), cimag((dp)->s0.c1), cimag((dp)->s0.c2),
//// cimag((dp)->s1.c0), cimag((dp)->s1.c1), cimag((dp)->s1.c2),
//// cimag((dp)->s2.c0), cimag((dp)->s2.c1), cimag((dp)->s2.c2),
//// cimag((dp)->s3.c0), cimag((dp)->s3.c1), cimag((dp)->s3.c2)
//// );
//// }
//// printf("\n\n\n");
//// print_spinor_similar_components(temp[0], P, VOLUME/2, 1e-4);
//// finalize_solver(tempE,2);
//// exit(1);
//// //////////////////////////////////////////////////////////// end of the test to be removed

//// finalize_solver(tempE,2);

} else
#endif
#ifdef TM_USE_QPHIX
Expand Down

0 comments on commit 87276e3

Please sign in to comment.