Skip to content

Commit

Permalink
WIP: interface for calling two-flavour QPhiX CG and BiCGstab from inv…
Browse files Browse the repository at this point in the history 
…ert_doublet_eo in place, solver does not converge yet
  • Loading branch information
@kostrzewa
kostrzewa committed Jun 30, 2017
1 parent 9941c2c commit a3c43a7
Show file tree
Hide file tree
Showing 10 changed files with 253 additions and 146 deletions.
8 changes: 4 additions & 4 deletions invert_clover_eo.c
View file
Expand Up @@ -66,7 +66,7 @@ int invert_clover_eo(spinor * const Even_new, spinor * const Odd_new,
const int solver_flag, const int rel_prec, const int even_odd_flag,
solver_params_t solver_params,
su3 *** gf, matrix_mult Qsq, matrix_mult Qm,
const ExternalInverter inverter, const SloppyPrecision sloppy, const CompressionType compression) {
const ExternalInverter external_inverter, const SloppyPrecision sloppy, const CompressionType compression) {
int iter;

if(even_odd_flag) {
Expand All @@ -75,7 +75,7 @@ int invert_clover_eo(spinor * const Even_new, spinor * const Odd_new,
}

#ifdef QUDA
if( inverter==QUDA_INVERTER ) {
if( external_inverter==QUDA_INVERTER ) {
return invert_eo_quda(Even_new, Odd_new, Even, Odd,
precision, max_iter,
solver_flag, rel_prec,
Expand Down Expand Up @@ -109,10 +109,10 @@ int invert_clover_eo(spinor * const Even_new, spinor * const Odd_new,

/* Here we invert the hermitean operator squared */
#ifdef TM_USE_QPHIX
if( inverter==QPHIX_INVERTER ) {
if( external_inverter==QPHIX_INVERTER ) {
// QPhiX inverts M(mu)M(mu)^dag or M(mu), no gamma_5 multiplication required
iter = invert_eo_qphix_oneflavour(Odd_new, g_spinor_field[DUM_DERI],
precision, max_iter,
max_iter, precision,
solver_flag, rel_prec,
solver_params,
sloppy,
Expand Down
2 changes: 1 addition & 1 deletion invert_clover_eo.h
View file
Expand Up @@ -13,6 +13,6 @@ int invert_clover_eo(spinor * const Even_new, spinor * const Odd_new,
const int solver_flag, const int rel_prec,
const int even_odd_flag, solver_params_t solver_params,
su3 *** gf, matrix_mult Qsq, matrix_mult Qm,
const ExternalInverter inverter, const SloppyPrecision sloppy, const CompressionType compression);
const ExternalInverter external_inverter, const SloppyPrecision sloppy, const CompressionType compression);

#endif
87 changes: 48 additions & 39 deletions invert_doublet_eo.c
View file
Expand Up @@ -50,7 +50,9 @@
#ifdef QUDA
# include "quda_interface.h"
#endif

#ifdef TM_USE_QPHIX
#include "qphix_interface.h"
#endif

#ifdef HAVE_GPU
# include"GPU/cudadefs.h"
Expand All @@ -63,20 +65,19 @@ extern su3* g_trafo;
# endif
#endif


int invert_doublet_eo(spinor * const Even_new_s, spinor * const Odd_new_s,
spinor * const Even_new_c, spinor * const Odd_new_c,
spinor * const Even_s, spinor * const Odd_s,
spinor * const Even_c, spinor * const Odd_c,
const double precision, const int max_iter,
const int solver_flag, const int rel_prec,
solver_params_t solver_params, const ExternalInverter inverter,
solver_params_t solver_params, const ExternalInverter external_inverter,
const SloppyPrecision sloppy, const CompressionType compression) {

int iter = 0;

#ifdef QUDA
if( inverter==QUDA_INVERTER ) {
if( external_inverter==QUDA_INVERTER ) {
return invert_doublet_eo_quda( Even_new_s, Odd_new_s, Even_new_c, Odd_new_c,
Even_s, Odd_s, Even_c, Odd_c,
precision, max_iter,
Expand Down Expand Up @@ -118,43 +119,51 @@ int invert_doublet_eo(spinor * const Even_new_s, spinor * const Odd_new_s,
printf("# Using CG for TMWILSON flavour doublet!\n");
fflush(stdout);
}
gamma5(g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI], VOLUME/2);
gamma5(g_spinor_field[DUM_DERI+1], g_spinor_field[DUM_DERI+1], VOLUME/2);

if ( external_inverter == NO_EXT_INV ){
gamma5(g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI], VOLUME/2);
gamma5(g_spinor_field[DUM_DERI+1], g_spinor_field[DUM_DERI+1], VOLUME/2);

#ifdef HAVE_GPU
if (usegpu_flag) { // GPU, mixed precision solver
# if ( defined TM_USE_MPI && defined PARALLELT )
iter = mixedsolve_eo_nd(Odd_new_s, Odd_new_c, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1],
max_iter, precision, rel_prec);
# elif ( !defined TM_USE_MPI && !defined PARALLELT )
iter = mixedsolve_eo_nd(Odd_new_s, Odd_new_c, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1],
max_iter, precision, rel_prec);
# else
printf("MPI and/or PARALLELT are not appropriately set for the GPU implementation. Aborting...\n");
exit(-1);
# endif
}
else { // CPU, conjugate gradient
iter = cg_her_nd(Odd_new_s, Odd_new_c, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1],
max_iter, precision, rel_prec,
VOLUME/2, &Qtm_pm_ndpsi);
}
if (usegpu_flag) { // GPU, mixed precision solver
# if ( defined TM_USE_MPI && defined PARALLELT )
iter = mixedsolve_eo_nd(Odd_new_s, Odd_new_c, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1],
max_iter, precision, rel_prec);
# elif ( !defined TM_USE_MPI && !defined PARALLELT )
iter = mixedsolve_eo_nd(Odd_new_s, Odd_new_c, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1],
max_iter, precision, rel_prec);
# else
printf("MPI and/or PARALLELT are not appropriately set for the GPU implementation. Aborting...\n");
exit(-1);
# endif
}
else { // CPU, conjugate gradient
iter = cg_her_nd(Odd_new_s, Odd_new_c, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1],
max_iter, precision, rel_prec,
VOLUME/2, &Qtm_pm_ndpsi);
}
#else // CPU, conjugate gradient
if(solver_flag == RGMIXEDCG){
iter = rg_mixed_cg_her_nd(Odd_new_s, Odd_new_c, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1],
solver_params, max_iter, precision, rel_prec, VOLUME/2,
&Qtm_pm_ndpsi, &Qtm_pm_ndpsi_32);
}
else {
iter = cg_her_nd(Odd_new_s, Odd_new_c, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1],
max_iter, precision, rel_prec, VOLUME/2, &Qtm_pm_ndpsi);
}
if(solver_flag == RGMIXEDCG){
iter = rg_mixed_cg_her_nd(Odd_new_s, Odd_new_c, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1],
solver_params, max_iter, precision, rel_prec, VOLUME/2,
&Qtm_pm_ndpsi, &Qtm_pm_ndpsi_32);
}
else {
iter = cg_her_nd(Odd_new_s, Odd_new_c, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1],
max_iter, precision, rel_prec, VOLUME/2, &Qtm_pm_ndpsi);
}
#endif


Qtm_dagger_ndpsi(Odd_new_s, Odd_new_c,
Odd_new_s, Odd_new_c);
Qtm_dagger_ndpsi(Odd_new_s, Odd_new_c,
Odd_new_s, Odd_new_c);
} // if(NO_EXT_INV)
#ifdef TM_USE_QPHIX
else if (external_inverter == QPHIX_INVERTER ) {
// using QPhiX, we invert M M^dagger y = b, so we don't need gamma_5 multiplications
iter = invert_eo_qphix_twoflavour(Odd_new_s, Odd_new_c, g_spinor_field[DUM_DERI], g_spinor_field[DUM_DERI+1],
max_iter, precision, solver_flag, rel_prec,
solver_params, sloppy, compression);
// and it multiplies y internally by M^dagger, returning M^{-1} b as required
}
#endif // TM_USE_QPHIX

/* Reconstruct the even sites */
Hopping_Matrix(EO, g_spinor_field[DUM_DERI], Odd_new_s);
Expand Down Expand Up @@ -188,12 +197,12 @@ int invert_cloverdoublet_eo(spinor * const Even_new_s, spinor * const Odd_new_s,
spinor * const Even_c, spinor * const Odd_c,
const double precision, const int max_iter,
const int solver_flag, const int rel_prec, solver_params_t solver_params,
const ExternalInverter inverter, const SloppyPrecision sloppy, const CompressionType compression) {
const ExternalInverter external_inverter, const SloppyPrecision sloppy, const CompressionType compression) {

int iter = 0;

#ifdef QUDA
if( inverter==QUDA_INVERTER ) {
if( external_inverter==QUDA_INVERTER ) {
return invert_doublet_eo_quda( Even_new_s, Odd_new_s, Even_new_c, Odd_new_c,
Even_s, Odd_s, Even_c, Odd_c,
precision, max_iter,
Expand Down
4 changes: 2 additions & 2 deletions invert_doublet_eo.h
View file
Expand Up @@ -40,7 +40,7 @@ int invert_doublet_eo(spinor * const Even_new_s, spinor * const Odd_new_s,
spinor * const Even_c, spinor * const Odd_c,
const double precision, const int max_iter,
const int solver_flag, const int rel_prec, solver_params_t solver_params,
const ExternalInverter inverter, const SloppyPrecision sloppy, const CompressionType compression);
const ExternalInverter extenral_inverter, const SloppyPrecision sloppy, const CompressionType compression);


/* This is the full matrix multiplication */
Expand All @@ -57,5 +57,5 @@ int invert_cloverdoublet_eo(spinor * const Even_new_s, spinor * const Odd_new_s,
spinor * const Even_c, spinor * const Odd_c,
const double precision, const int max_iter,
const int solver_flag, const int rel_prec, solver_params_t solver_params,
const ExternalInverter inverter, const SloppyPrecision sloppy, const CompressionType compression);
const ExternalInverter external_inverter, const SloppyPrecision sloppy, const CompressionType compression);
#endif
8 changes: 4 additions & 4 deletions invert_eo.c
View file
Expand Up @@ -86,12 +86,12 @@ int invert_eo(spinor * const Even_new, spinor * const Odd_new,
const int solver_flag, const int rel_prec,
const int sub_evs_flag, const int even_odd_flag,
const int no_extra_masses, double * const extra_masses, solver_params_t solver_params, const int id,
const ExternalInverter inverter, const SloppyPrecision sloppy, const CompressionType compression ) {
const ExternalInverter external_inverter, const SloppyPrecision sloppy, const CompressionType compression ) {

int iter = 0;

#ifdef QUDA
if( inverter==QUDA_INVERTER ) {
if( external_inverter==QUDA_INVERTER ) {
return invert_eo_quda(Even_new, Odd_new, Even, Odd,
precision, max_iter,
solver_flag, rel_prec,
Expand Down Expand Up @@ -159,10 +159,10 @@ int invert_eo(spinor * const Even_new, spinor * const Odd_new,
/* matrix to get the odd sites */

#ifdef TM_USE_QPHIX
if( inverter==QPHIX_INVERTER ) {
if( external_inverter==QPHIX_INVERTER ) {
// QPhiX inverts M(mu)M(mu)^dag or M(mu), no gamma_5 source multiplication required
iter = invert_eo_qphix_oneflavour(Odd_new, g_spinor_field[DUM_DERI],
precision, max_iter,
max_iter, precision,
solver_flag, rel_prec,
solver_params,
sloppy,
Expand Down
2 changes: 1 addition & 1 deletion invert_eo.h
View file
Expand Up @@ -36,6 +36,6 @@ int invert_eo(spinor * const Even_new, spinor * const Odd_new,
const int solver_flag, const int rel_prec,
const int sub_evs_flag, const int even_odd_flag,
const int no_extra_masses, double * const extra_masses, solver_params_t solver_params, const int id,
const ExternalInverter inverter, const SloppyPrecision sloppy, const CompressionType compression );
const ExternalInverter external_inverter, const SloppyPrecision sloppy, const CompressionType compression );

#endif

0 comments on commit a3c43a7

Please sign in to comment.