Merge pull request #17 from daschaich/rescale

Rescale

4d_Q16/susy/Make_template

@@ -144,13 +144,13 @@ susy_eig::
 
 susy_cheb::
 	${MAKE} -f ${MAKEFILE} target "MYTARGET= $@" \
-	"DEFINES = ${DEFINES} -DPHI_ALGORITHM -DCHEB " \
+	"DEFINES = ${DEFINES} -DPHI_ALGORITHM -DCHEB -DRESCALE " \
 	"LAPACK = -llapack -lblas " \
 	"EXTRA_OBJECTS = control_cheb.o chebyshev.o z2source.o "
 
 susy_mode::
 	${MAKE} -f ${MAKEFILE} target "MYTARGET= $@" \
-	"DEFINES = ${DEFINES} -DPHI_ALGORITHM -DMODE " \
+	"DEFINES = ${DEFINES} -DPHI_ALGORITHM -DMODE -DRESCALE " \
 	"LAPACK = -llapack -lblas " \
 	"EXTRA_OBJECTS = control_mode.o modenumber.o mode_coeffs.o z2source.o "
 

4d_Q16/susy/README

@@ -30,6 +30,7 @@ Eigenvalue measurements (susy_eig) require PRIMME (https://github.com/primme/pri
 -DEIG switches on the PRIMME eigenvalue calculation
 -DCHEB switches on Chebyshev spectral density computations
 -DMODE switches on stochastic eigenmode number computations
+-DRESCALE switches to more symmetric fermion operator
 -DPL_CORR switches on the Polyakov loop correlator calculation (NOT CURRENTLY IN USE)
 -DPUREGAUGE switches off the fermions (FOR TESTING)
 

4d_Q16/susy/congrad_multi.c

@@ -40,7 +40,6 @@ int congrad_multi(Twist_Fermion *src, Twist_Fermion **psim,
   double *beta_im1 = malloc(sizeof *beta_im1 * Norder);
   double *alpha    = malloc(sizeof *alpha * Norder);
   double rsqj;
-  complex ctmp;
   Twist_Fermion **pm = malloc(sizeof(Twist_Fermion*) * Norder);
   for (i = 1; i < Norder; i++)    // !!!
     pm[i] = malloc(sizeof(Twist_Fermion) * sites_on_node);
@@ -86,10 +85,8 @@ int congrad_multi(Twist_Fermion *src, Twist_Fermion **psim,
 
     // beta_i[0] = -(r, r) / (pm, Mpm)
     cd = 0;
-    FORALLSITES(i, s) {
-      ctmp = TF_dot(&(pm0[i]), &(mpm[i]));
-      cd += ctmp.real;
-    }
+    FORALLSITES(i, s)
+      TF_rdot_sum(&(pm0[i]), &(mpm[i]), &cd);
     g_doublesum(&cd);
 
     beta_i[0] = -rsq / cd;
@@ -155,7 +152,7 @@ int congrad_multi(Twist_Fermion *src, Twist_Fermion **psim,
       floatvark[j] = (Real)alpha[j];
     }
     FORALLSITES(i, s) {
-      scalar_mult_TF(&(rm[i]),floatvar, &(mpm[i]));
+      scalar_mult_TF(&(rm[i]), floatvar, &(mpm[i]));
       scalar_mult_add_TF(&(mpm[i]), &(pm0[i]), floatvar2, &(pm0[i]));
       for (j = 1; j < Norder; j++) {
         if (converged[j] == 0) {
@@ -215,7 +212,7 @@ int congrad_multi(Twist_Fermion *src, Twist_Fermion **psim,
     source_norm = 0;                // Re-using for convenience
     FORALLSITES(i, s) {             // Add shift.psi and subtract src
       scalar_mult_sum_TF(&(psim[j][i]), shift[j],  &(mpm[i]));
-      scalar_mult_sum_TF(&(src[i]), -1.0, &(mpm[i]));
+      dif_TF(&(src[i]), &(mpm[i]));
       source_norm += (double)magsq_TF(&(mpm[i]));
     }
     g_doublesum(&source_norm);

4d_Q16/susy/control_cheb.c

@@ -10,7 +10,7 @@
 #include "susy_includes.h"
 
 #ifndef CHEB
-#error "Don't use control_mode unless compiling with -DMODE!"
+#error "Don't use control_cheb unless compiling with -DCHEB!"
 #endif
 #ifdef MODE
 #error "-DCHEB and -DMODE will clobber each other!"
@@ -93,7 +93,6 @@ int main(int argc, char *argv[]) {
   node0_printf("RUNNING COMPLETED\n");
   dtime += dclock();
   node0_printf("\nTime = %.4g seconds\n", dtime);
-  node0_printf("total_iters = %d\n", total_iters);
   fflush(stdout);
 
   normal_exit(0);         // Needed by at least some clusters

4d_Q16/susy/defines.h

@@ -27,6 +27,7 @@
 #define SV                // Site/vector terms in action
 #define VP                // Vector/plaquette terms in action
 #define QCLOSED           // Q-closed terms in action
+//#define RESCALE           // Rescale for more symmetric fermion operator
 //#define DIMREDUCE         // Consider dimensionally reduced systems
 //#define DEBUG_CHECK       // Print lambdas, offsets, etc.
 

4d_Q16/susy/setup.c

@@ -124,6 +124,9 @@ void make_fields() {
   node0_printf("Single-trace scalar potential\n");
 #else
   node0_printf("Double-trace scalar potential\n");
+#endif
+#ifdef RESCALE
+  node0_printf("Rescaled fermion operator\n");
 #endif
   Real size = (Real)(2.0 * sizeof(complex));
   FIELD_ALLOC(tr_eta, complex);

4d_Q16/susy/update_h.c

@@ -268,6 +268,7 @@ double gauge_force(Real eps) {
 // All called by assemble_fermion_force below
 // First Q-closed piece: chi_ab D_c chi_de epsilon_{abcde}
 // Note factor of -1/2
+#ifdef QCLOSED
 void F1Q(matrix *plaq_sol[NPLAQ], matrix *plaq_psol[NPLAQ]) {
   register int i, opp_a, opp_b;
   register site *s;
@@ -359,13 +360,15 @@ void F1Q(matrix *plaq_sol[NPLAQ], matrix *plaq_psol[NPLAQ]) {
     flip = (flip + 1) % 2;
   }
 }
+#endif
 // -----------------------------------------------------------------
 
 
 
 // -----------------------------------------------------------------
 // Second Q-closed piece
 // Note factor of -1/2
+#ifdef QCLOSED
 void F2Q(matrix *plaq_sol[NPLAQ], matrix *plaq_psol[NPLAQ]) {
   register int i, opp_a, opp_b;
   register site *s;
@@ -457,6 +460,7 @@ void F2Q(matrix *plaq_sol[NPLAQ], matrix *plaq_psol[NPLAQ]) {
     flip = (flip + 1) % 2;
   }
 }
+#endif
 // -----------------------------------------------------------------
 
 
@@ -469,6 +473,7 @@ void F2Q(matrix *plaq_sol[NPLAQ], matrix *plaq_psol[NPLAQ]) {
 // Assume compute_plaqdet() has already been run
 // Appropriate adjoints set up in assemble_fermion_force
 // A bit more code reuse may be possible
+#ifdef SV
 void detF(matrix *eta, matrix *psi[NUMLINK], int sign) {
   register int i;
   register site *s;
@@ -644,6 +649,7 @@ void detF(matrix *eta, matrix *psi[NUMLINK], int sign) {
   free(inv_term);
   free(adj_term);
 }
+#endif
 // -----------------------------------------------------------------
 
 
@@ -686,6 +692,13 @@ void assemble_fermion_force(Twist_Fermion *sol, Twist_Fermion *psol) {
       mat_copy(&(sol[i].Fplaq[mu]), &(plaq_src[mu][i]));
       adjoint(&(psol[i].Fplaq[mu]), &(plaq_dest[mu][i]));
     }
+
+    // Optionally rescale to make fermion operator more symmetric
+    // psol[i].Fsite already rescaled in fermion_op
+    // !!!TODO: Not yet conserving...
+#ifdef RESCALE
+    scalar_mult_matrix(&(site_src[i]), 2.0, &(site_src[i]));
+#endif
   }
 
 #ifdef SV
@@ -730,6 +743,11 @@ void assemble_fermion_force(Twist_Fermion *sol, Twist_Fermion *psol) {
     }
     cleanup_gather(mtag[mu]);
   }
+#else
+  FORALLDIR(mu) {                         // Zero force collectors
+    FORALLSITES(i, s)
+      clear_mat(&(s->f_U[mu]));
+  }
 #endif
 #ifdef VP
   // Now calculate DU on chi_{munu} D_mu(U) psi_nu
@@ -875,6 +893,7 @@ void assemble_fermion_force(Twist_Fermion *sol, Twist_Fermion *psol) {
   }
 #endif
 
+#ifdef SV
   // Plaquette determinant contributions if G is non-zero
   if (doG) {
     // First connect link_src with site_dest[DIMF - 1]^dag (LtoS)
@@ -883,6 +902,7 @@ void assemble_fermion_force(Twist_Fermion *sol, Twist_Fermion *psol) {
     // Second connect site_src[DIMF - 1] with link_dest^dag (StoL)
     detF(site_src, link_dest, MINUS);
   }
+#endif
 
 #ifdef QCLOSED
   if (NUMLINK != 5) {

4d_Q16/susy/utilities.c

@@ -772,13 +772,25 @@ void fermion_op(Twist_Fermion *src, Twist_Fermion *dest, int sign) {
     node0_printf("Error: incorrect sign in fermion_op: %d\n", sign);
     terminate(1);
   }
-  FORALLSITES(i, s)
+  FORALLSITES(i, s) {
+    // Optionally rescale to make fermion operator more symmetric
+#ifdef RESCALE
+    scalar_mult_matrix(&(site_src[i]), 2.0, &(site_src[i]));
+#endif
     tr_eta[i] = trace(&(site_src[i]));
+  }
 
   // Assemble separate routines for each term in the fermion operator
 #ifdef VP
   Dplus(link_src, plaq_dest);             // Overwrites plaq_dest
   Dminus(plaq_src, link_dest);            // Overwrites link_dest
+#else
+  FORALLSITES(i, s) {                     // Zero link_dest and plaq_dest
+    FORALLDIR(mu)
+      clear_mat(&(link_dest[mu][i]));
+    for (mu = 0; mu < NPLAQ; mu++)
+      clear_mat(&(plaq_dest[mu][i]));
+  }
 #endif
 
 #ifdef SV
@@ -794,13 +806,22 @@ void fermion_op(Twist_Fermion *src, Twist_Fermion *dest, int sign) {
   // Link-to-site plaquette determinant contribution if G is non-zero
   if (doG)
     detLtoS(link_src, site_dest);         // Adds to site_dest
+#else
+  FORALLSITES(i, s)                       // Zero site_dest
+    clear_mat(&(site_dest[i]));
 #endif
 
 #ifdef QCLOSED
   DbminusPtoP(plaq_src, plaq_dest);       // Adds to plaq_dest
   DbplusPtoP(plaq_src, plaq_dest);        // Adds to plaq_dest
 #endif
 
+  // Optionally rescale to make fermion operator more symmetric
+#ifdef RESCALE
+  FORALLSITES(i, s)
+    scalar_mult_matrix(&(site_dest[i]), 2.0, &(site_dest[i]));
+#endif
+
   // Copy local plaquette, link and site fermions into dest TwistFermion
   if (sign == 1) {
     FORALLSITES(i, s) {
@@ -835,7 +856,7 @@ void DSq(Twist_Fermion *src, Twist_Fermion *dest) {
 
   fermion_op(src, tempTF, PLUS);
   fermion_op(tempTF, dest, MINUS);
-  if (fmass > IMAG_TOL) {
+  if (fmass > IMAG_TOL) {           // Assume fmass non-negative
     Real fmass2 = fmass * fmass;
     FORALLSITES(i, s)
       scalar_mult_sum_TF(&(src[i]), fmass2, &(dest[i]));