Skip to content

Commit

Permalink
Merge pull request milc-qcd#18 from milc-qcd/develop
Browse files Browse the repository at this point in the history 
Develop update
  • Loading branch information
@maddyscientist
maddyscientist authored Jan 2, 2020
2 parents 193c3f7 + 4116cdd commit b0afc2e
Show file tree
Hide file tree
Showing 17 changed files with 238 additions and 48 deletions.
3 changes: 3 additions & 0 deletions Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -327,6 +327,9 @@ include ../Make_template_scidac
ifeq ($(strip ${COMPILER}),intel)
INCFFTW = -mkl
LIBFFTW = -mkl
else ifeq ($(strip ${COMPILER}),cray-intel)
INCFFTW = -mkl
LIBFFTW = -mkl
endif

#----------------------------------------------------------------------
Expand Down
2 changes: 1 addition & 1 deletion file_utilities/Make_template
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,7 +133,7 @@ check_gauge::
"DEFINES= -DONLY_GLUON_FILES" \
"EXTRA_OBJECTS= ${G_OBJECTS} check_gauge.o nersc_cksum.o \
gauge_utilities.o \
check_unitarity.o gauge_info_dummy.o d_plaq4.o io_lat4.o"
check_unitarity.o gauge_info_dummy.o d_plaq4.o io_lat4.o io_scidac.o"

diff_colormatrix::
${MAKE} -f ${MAKEFILE} target "MYTARGET= diff_colormatrix" \
Expand Down
2 changes: 1 addition & 1 deletion generic/io_scidac.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ QIO_Writer *open_scidac_output(char *filename, int volfmt,
#ifdef QIO_TRELEASE
QIO_set_trelease(0,QIO_TRELEASE);
#endif
QIO_verbose(QIO_VERB_REG);
QIO_verbose(QIO_VERB_OFF);
outfile = QIO_open_write(xml_write_file, filename, volfmt, layout,
fs, &oflag);
if(outfile == NULL){
Expand Down
21 changes: 12 additions & 9 deletions generic_ks/eigen_stuff_PRIMME.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
#ifdef PRIMME

/* NOTE: The PRIMME release version has clashing definitions for complex functions, so we provide
a modifie version. This must be checked against future releases.
a modified version. This must be checked against future releases.
Also, watch out for the definition of PRIMME_INT = long? */
#include "../include/primme.h"
#include "../include/dslash_ks_redefine.h"
Expand Down Expand Up @@ -220,17 +220,19 @@ int ks_eigensolve_PRIMME(su3_vector **eigVec, double *eigVal,
primme.globalSumReal=par_GlobalSumDouble; /* the wrapper function to do global sums */

primme.matrixMatvec =ks_mxv; /* the matrix on vector product */

// ret = primme_set_method(PRIMME_DEFAULT_MIN_MATVECS, &primme);
ret = primme_set_method(PRIMME_DYNAMIC, &primme);


/* primme.printLevel=3; */
#ifdef PRIMME_PRINTLEVEL
node0_printf("Setting print level to %d\n", PRIMME_PRINTLEVEL);
primme.printLevel = PRIMME_PRINTLEVEL;
#else
node0_printf("Setting print level to default\n");
primme.printLevel = 1;
#endif

// ret = primme_set_method(PRIMME_DEFAULT_MIN_MATVECS, &primme);
ret = primme_set_method(PRIMME_DYNAMIC, &primme);

#ifdef MATVEC_PRECOND
primme.target=primme_largest;
#else
Expand All @@ -242,7 +244,7 @@ int ks_eigensolve_PRIMME(su3_vector **eigVec, double *eigVal,
primme.applyPreconditioner = ks_precond_mxv;
#endif

/* Optimaized Parameter Setting */
/* Optimized Parameter Setting */
primme.correctionParams.robustShifts = 1; // led to faster convergence with 0 for tol=1e-8
primme.locking=1;
#if 1 /* James Osborn's and Xiao-Yong Jin's optimal setting */
Expand All @@ -261,16 +263,17 @@ int ks_eigensolve_PRIMME(su3_vector **eigVec, double *eigVal,

/* Display parameters */
if(this_node==0){
// printf("PRIMME workspace int = %d long int = %ld\n", primme.intWorkSize, primme.realWorkSize); fflush(stdout);
primme_display_params(primme);
printf("PRIMME_INT size is %d\n", sizeof(PRIMME_INT));
// printf("primme_op_datatype size is %d\n", sizeof(primme_op_datatype));
fflush(stdout);
}

/* call the actual EV finder*/
ret = zprimme(evals, (PRIMME_COMPLEX_DOUBLE *)evecs, rnorms, &primme);

if (ret!=0){ /*check return value */
node0_printf("ks_eigensolve_PRIMME: zprimme error\nCall stack:\n");
/** primme_PrintStackTrace(primme); NOT SUPPORTED in 2.0**/
node0_printf("ks_eigensolve_PRIMME: zprimme error %d\n",ret);
fflush(stdout);
exit(1);
}
Expand Down
25 changes: 24 additions & 1 deletion generic_ks/f_meas_current.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,16 @@ thin_source(su3_vector *src, int thinning, int ex, int ey, int ez, int et){
}

/*Write current record for the accumulated average over random sources */
static void
write_tslice_values_begin(char *tag){
node0_printf("BEGIN JTMU%s\n", tag);
}

static void
write_tslice_values_end(char *tag){
node0_printf("END JTMU%s\n", tag);
}

static void
write_tslice_values(char *tag, int jr, Real mass1, Real mass2, Real *j_mu ){
double *jtmu = (double *)malloc(sizeof(double)*4*param.nt);
Expand Down Expand Up @@ -846,6 +856,7 @@ f_meas_current_diff( int n_masses, int nrand, int thinning,
#endif


write_tslice_values_begin("DIFF");
for(int j = 0; j < n_masses; j++){
for(int ir = 0; ir < nr; ir++){
#if 0
Expand All @@ -871,6 +882,7 @@ f_meas_current_diff( int n_masses, int nrand, int thinning,
clear_r_array_field(j_mu[j][ir], NMU);
} /* ir */
} /* j */
write_tslice_values_end("DIFF");

} /* jrand */

Expand Down Expand Up @@ -1406,7 +1418,8 @@ f_meas_current( int n_masses, int nrand, int thinning,
#endif

/* Print the exact low mode contribution to the current density */
if(Nvecs > 0)
if(Nvecs > 0){
write_tslice_values_begin("LOW");
for(int j = 0; j < n_masses; j++){
#ifdef MASS_UDLSC
if((j == 0 && n_masses > 1) || (j == 2 && n_masses > 3))
Expand All @@ -1421,6 +1434,8 @@ f_meas_current( int n_masses, int nrand, int thinning,
#endif
clear_r_array_field(jlow_mu[j], NMU);
} /* j */
write_tslice_values_end("LOW");
}


/* Block solver parameters */
Expand Down Expand Up @@ -1454,6 +1469,7 @@ f_meas_current( int n_masses, int nrand, int thinning,
block_current( n_masses, j_mu, masses, fn_mass, qic, nsrc, nr, gr_even, gr_odd);
#endif

write_tslice_values_begin("HI");
for(int j = 0; j < n_masses; j++){
for(int ir = 0; ir < nr; ir++){
#if 0
Expand Down Expand Up @@ -1481,6 +1497,7 @@ f_meas_current( int n_masses, int nrand, int thinning,
clear_r_array_field(j_mu[j][ir], NMU);
} /* ir */
} /* j */
write_tslice_values_end("HI");
} /* jrand */

for(int j = 0; j < n_masses; j++){
Expand Down Expand Up @@ -1642,6 +1659,7 @@ f_meas_current_diff( int n_masses, int nrand, int thinning,
} /* j */
} /* ex, ey, ez, et */

write_tslice_values_begin("DIFF");
for(j = 0; j < n_masses; j++){
Real mass = ksp[j].mass;
#if 0
Expand All @@ -1656,6 +1674,7 @@ f_meas_current_diff( int n_masses, int nrand, int thinning,
write_tslice_values("DIFF", jrand, ksp[j].mass, 0., j_mu[j]);
clear_r_array_field(j_mu[j], NMU);
} /* j */
write_tslice_values_end("DIFF");

} /* jrand */

Expand Down Expand Up @@ -1764,10 +1783,12 @@ f_meas_current( int n_masses, int nrand, int thinning,
} /* j */
} /* n */

write_tslice_values_begin("LOW");
for(j = 0; j < n_masses; j++){
write_tslice_values("LOW", 0, mass[j], 0., jlow_mu[j]);
clear_r_array_field(jlow_mu[j], NMU);
}
write_tslice_values_end("LOW");

#if 0
for(j = 0; j < n_masses; j++){
Expand All @@ -1785,10 +1806,12 @@ f_meas_current( int n_masses, int nrand, int thinning,
node0_printf("Time for exact low modes %g sec\n", dtime);
#endif

write_tslice_values_begin("");
for(j = 0; j < n_masses; j++){
write_tslice_values("", jrand, mass[j], 0., j_mu[j]);
clear_r_array_field(j_mu[j], NMU);
}
write_tslice_values_end("");

/* Loop over random sources */
for(jrand = 0; jrand < nrand; jrand++){
Expand Down
63 changes: 54 additions & 9 deletions generic_ks/io_helpers_ks_eigen.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -160,13 +160,13 @@ void w_close_ks_eigen(int flag, ks_eigen_file *kseigf){
>1 for seek, read error, or missing data error
*/
int reload_ks_eigen(int flag, char *eigfile, int *Nvecs, double *eigVal,
su3_vector **eigVec, int timing){
su3_vector **eigVec, imp_ferm_links_t *fn, int timing){

register int i, j;
int status = 0;
int serpar;
int qio_status;
int packed;
int packed = 0, file_type = 0;
QIO_Reader *infile;
double dtime = (double)0.0;
char myname[] = "reload_ks_eigen";
Expand All @@ -186,20 +186,36 @@ int reload_ks_eigen(int flag, char *eigfile, int *Nvecs, double *eigVal,
if(flag == RELOAD_SERIAL)serpar = QIO_SERIAL;
else serpar = QIO_PARALLEL;

infile = open_ks_eigen_infile(eigfile, Nvecs, &packed, serpar);
infile = open_ks_eigen_infile(eigfile, Nvecs, &packed, &file_type, serpar);
if(infile == NULL){
node0_printf("ERROR: Can't open %s for reading\n", eigfile);
status = 1;
break;
}
for(int i = 0; i < *Nvecs; i++){
qio_status = read_ks_eigenvector(infile, packed, eigVec[i], &eigVal[i]);
/* Read using MILC format */
if(file_type == 0){
for(int i = 0; i < *Nvecs; i++){
qio_status = read_ks_eigenvector(infile, packed, eigVec[i], &eigVal[i]);
if(qio_status != QIO_SUCCESS){
if(qio_status == QIO_EOF){
node0_printf("WARNING: Premature EOF at %d eigenvectors\n", i);
*Nvecs = i;
} else {
node0_printf("ERROR: Can't read an eigenvector. Error %d\n", qio_status);
status = 1;
}
break;
}
}
} else {
/* Read using QUDA format */
qio_status = read_quda_ks_eigenvectors(infile, eigVec, eigVal, Nvecs, EVEN, fn);
if(qio_status != QIO_SUCCESS){
if(qio_status == QIO_EOF){
node0_printf("WARNING: Premature EOF at %d eigenvectors\n", i);
*Nvecs = i;
node0_printf("WARNING: Premature EOF at eigenvectors\n");
terminate(1);
} else {
node0_printf("ERROR: Can't read an eigenvector. Error %d\n", qio_status);
node0_printf("ERROR: Can't read eigenvectors. Error %d\n", qio_status);
status = 1;
}
break;
Expand Down Expand Up @@ -326,6 +342,33 @@ int save_ks_eigen(int flag, char *savefile, int Nvecs, double *eigVal,

close_ks_eigen_outfile(outfile);
break;

case SAVE_PARTFILE_SCIDAC:
#ifdef PACK_EIGEN
node0_printf("ERROR: Eigenvector packing is incompatible with partfile\n");
node0_printf("Recompile with PACK_EIGEN undefined\n");
status = 1;
break;
#else
outfile = open_ks_eigen_outfile(savefile, Nvecs, QIO_PARTFILE, QIO_SERIAL);
if(outfile == NULL){
node0_printf("ERROR: Can't open %s for writing\n", savefile);
status = 1;
break;
}
for(int i = 0; i < Nvecs; i++){
int status = write_ks_eigenvector(outfile, eigVec[i], eigVal[i], resid[i]);
if(status != QIO_SUCCESS){
node0_printf("ERROR: Can't write eigenvector.\n");
status = 1;
break;
}
}

close_ks_eigen_outfile(outfile);
break;
#endif

default:
node0_printf("%s: Unrecognized save flag.\n", myname);
terminate(1);
Expand Down Expand Up @@ -680,7 +723,7 @@ int ask_starting_ks_eigen(FILE *fp, int prompt, int *flag, char *filename){

static void print_save_options(void){

printf("'forget_ks_eigen', 'save_ascii_ks_eigen', 'save_serial_ks_eigen', or 'save_parallel_ks_eigen'");
printf("'forget_ks_eigen', 'save_ascii_ks_eigen', 'save_serial_ks_eigen', 'save_parallel_ks_eigen', or 'save_partfile_ks_eigen'");
}

/*--------------------------------------------------------------------*/
Expand Down Expand Up @@ -714,6 +757,8 @@ int ask_ending_ks_eigen(FILE *fp, int prompt, int *flag, char *filename){
*flag = SAVE_SERIAL;
else if(strcmp("save_parallel_ks_eigen",savebuf) == 0)
*flag = SAVE_PARALLEL;
else if(strcmp("save_partfile_ks_eigen",savebuf) == 0)
*flag = SAVE_PARTFILE_SCIDAC;
else{
printf("ERROR IN INPUT: ks_eigen outpu command %s is invalid.\n", savebuf);
printf("Choices are ");
Expand Down
Loading

0 comments on commit b0afc2e

Please sign in to comment.