433-milc Report

**433-milc**
Exec Time (%)	Codelet Name	Error (%)
99.8	control_main_32	99.8
1.1	gauge_stuff_g_measure_345	1.7
22.2	path_product_path_product_49	7.8
1.3	com_vanilla_start_gather_from_temp_581	11.5
96.1	update_update_40	6.2
24.9	gauge_stuff_imp_gauge_force_230	13.4
1.2	path_product_path_product_156	36.4
22.2	path_product_path_product_49	7.8
1.3	com_vanilla_start_gather_from_temp_581	11.5
46.6	quark_stuff_eo_fermion_force_946	6.6
10.7	quark_stuff_u_shift_fermion_1452	7.8
2.0	quark_stuff_u_shift_fermion_1466	4.6
10.1	quark_stuff_u_shift_fermion_1460	19.8
19.9	quark_stuff_add_force_to_mom_1523	8.6
1.3	com_vanilla_start_gather_from_temp_581	11.5
22.2	path_product_path_product_49	7.8
1.3	com_vanilla_start_gather_from_temp_581	11.5
17.8	quark_stuff_load_fatlinks_566	1.7
3.5	quark_stuff_compute_gen_staple_1368	4.9
6.9	quark_stuff_compute_gen_staple_1396	100.0
4.2	quark_stuff_compute_gen_staple_1375	8.3
1.3	quark_stuff_compute_gen_staple_1406	100.0
1.3	com_vanilla_start_gather_from_temp_581	11.5
1.5	update_u_update_u_40	1.5
22.2	path_product_path_product_49	7.8
1.3	com_vanilla_start_gather_from_temp_581	11.5
17.8	quark_stuff_load_fatlinks_566	1.7
3.5	quark_stuff_compute_gen_staple_1368	4.9
6.9	quark_stuff_compute_gen_staple_1396	100.0
4.2	quark_stuff_compute_gen_staple_1375	8.3
1.3	quark_stuff_compute_gen_staple_1406	100.0
1.3	com_vanilla_start_gather_from_temp_581	11.5

Exec Time : 1.330863e+10 cycles

Matching : 97.2%

                /****** quark_stuff.c  -- ******************/
/* MIMD version 6 */
/* quark action stuff for improved action
* D.T. 1/28/98, starting from gauge_stuff.c
* K.O. 3/99 Added optimized fattening for Asq actions
* D.T. 4/99 Combine force calculations for both mass quarks
* K.O. 4/99 Optimized force for Asq action
* S.G. 7/01, modified to use t_longlink and t_fatlink
* C.D. 10/02, consolidated quark_stuff.c and quark_stuff_tmp.c
*
* This code combines quark_stuff.c and quark_stuff_tmp.c
* with DSLASH_TMP_LINKS defined, puts links in field-major variables
* t_longlink and t_fatlink. Otherwise, puts them in the site structure.
*
* In this directory, assume all paths connect even to odd sites, etc.
* Tabulate "backwards" paths (e.g. "XDOWN" is backward path to "XUP")
* as separate parity transforms of the fundamental paths.  They will
* generally need a negative sign in Dslash.  See bottom for a long
* comment on sign conventions.
*/

/**#define FFTIME**/
/**#define LLTIME**/

#include "generic_ks_includes.h"	/* definitions files and prototypes */


#define NULL_FP -1 /* NULL field_offset to be used in the optimized version *
                    * of the load_fatlinks subroutine */

    /* Specify paths in orientation in which they appear in the
       forward part of the x component of dslash().  Rotations and
       reflections will be automatically included. Be careful
       about signs of coefficients.  See long comment at bottom. */
#include <quark_action.h>
    /* Include file specifies the basic paths */

#define GOES_FORWARDS(dir) (dir<=TUP)
#define GOES_BACKWARDS(dir) (dir>TUP)

void printpath( int *path, int length );
void path_transport( field_offset src, field_offset dest, int parity,
    int *dir, int length );
void path_transport_hwv( field_offset src, field_offset dest, int parity,
    int *dir, int length );
#ifdef  ASQ_OPTIMIZED_FATTENING
void compute_gen_staple(field_offset staple, int mu, int nu,
                        field_offset link, double coef ) ;
#endif
#ifdef  ASQ_OPTIMIZED_FORCE
#ifndef FN
BOMB THE COMPILE
#endif
void u_shift_fermion(su3_vector *src, su3_vector *dest, int dir ) ;
void add_force_to_mom(su3_vector *back, su3_vector *forw, int dir, double coef);
void side_link_force(int mu, int nu, double coeff, su3_vector *Path,
		     su3_vector *Path_nu, su3_vector *Path_mu, 
		     su3_vector *Path_numu) ;

void u_shift_hw_fermion(half_wilson_vector *src, 
			half_wilson_vector *dest, int dir ) ;
void add_3f_force_to_mom(half_wilson_vector *back,
			 half_wilson_vector *forw, int dir, double coeff[2]) ;
void side_link_3f_force(int mu, int nu, double coeff[2], 
			half_wilson_vector *Path   , 
			half_wilson_vector *Path_nu, 
			half_wilson_vector *Path_mu, 
			half_wilson_vector *Path_numu) ;
#endif

int path_num[MAX_BASIC_PATHS];	/* number of rotations/reflections for each 
					kind */
static double act_path_coeff[MAX_BASIC_PATHS]; /* actual path coefficient     *
                                               * it is equal to path_coeff   *
                                               * if not tadpole improvement  *
                                               * is specified                *
                                               * or path_coeff*u_0^(L-1) when*
                                               * tadpole improvement is      *
                                               * specified                   */

/* Array of structures, for each rotation and reflection of each kind of
	path.  */
struct {
    int dir[MAX_LENGTH];	/* directions in path */
    int length;		/* length of path */
    double coeff;	/* coefficient, including minus sign if backwards */
    double forwback;	/* +1 if in forward Dslash, -1 if in backward */
}   q_paths[MAX_NUM];
int num_q_paths;	/* number of paths in dslash */
int num_basic_paths;	/* number of paths before rotation/reflection */

int is_path_equal( int *path1, int* path2, int length );
int add_basic_path( int *vec, int length, double coeff );

/* Make table of paths in action */
void make_path_table() {

    int i,j;
#ifdef TADPOLE_IMPROVE
    int k;
#endif

    /* table of directions, 1 for each kind of path */
    /**int path_ind[MAX_BASIC_PATHS][MAX_LENGTH];**/
    /* table of coefficients in action, for each path */

    node0_printf("%s\n",quark_action_description);
    num_q_paths = 0;
    num_basic_paths = 0;


    /* add rots. and reflects to table, print out the path coefficients */
    node0_printf("path coefficients: npath  path_coeff  multiplicity\n");
    for(j=0;j<quark_action_npaths;j++) {
	double this_coeff;
	this_coeff = path_coeff[j];
#ifdef TADPOLE_IMPROVE
	for(k=1;k< path_length_in[j];k++)this_coeff /= u0;
#endif
	act_path_coeff[j] = this_coeff ;
	i = add_basic_path( path_ind[j], path_length_in[j],
	    this_coeff );
	node0_printf("                    %d      %e     %d\n",
	    j,this_coeff,i);
    }
}

/* add rotations and reflections of a path to the table.  Return
   multiplicity of paths added */
int add_basic_path( int *basic_vec, int length, double coeff ) {

    int perm[8],pp[8],ir[4];
    int j,path_num;
    int vec[MAX_LENGTH];
    int flag;

    path_num = 0;
    /* now fill the long table with all rotations and reflections
	of the fundamental path.  The path presented to us is for
        the positive x component of dslash, so if the x coordinate
        is reflected it will appear with a negative sign. */
      /* permutations */
      for(perm[0]=0;perm[0]<4;perm[0]++)
      for(perm[1]=0;perm[1]<4;perm[1]++)
      for(perm[2]=0;perm[2]<4;perm[2]++)
      for(perm[3]=0;perm[3]<4;perm[3]++){
	if(perm[0] != perm[1] && perm[0] != perm[2] 
	  && perm[0] != perm[3] && perm[1] != perm[2]
	  && perm[1] != perm[3] && perm[2] != perm[3] ) {
	  /* reflections*/
	  for(ir[0]=0;ir[0]<2;ir[0]++)
	  for(ir[1]=0;ir[1]<2;ir[1]++)
	  for(ir[2]=0;ir[2]<2;ir[2]++)
	  for(ir[3]=0;ir[3]<2;ir[3]++){
	    for(j=0;j<4;j++){
	      pp[j]=perm[j];

	      if(ir[j] == 1) pp[j]=OPP_DIR(pp[j]);
	      pp[OPP_DIR(j)]=OPP_DIR(pp[j]);
	    }
	    /* create new vector*/
	    for(j=0;j<length;j++) vec[j]=pp[basic_vec[j]];
	    for(j=length;j<MAX_LENGTH;j++) vec[j]=NODIR;

            flag=0;
	    /* check if it's a new set: */
	    for(j=0;j<num_q_paths;j++){
	      flag = is_path_equal( vec, q_paths[j].dir, MAX_LENGTH );
	      if(flag==1)break;
	    }
	    if(flag == 0 ){
	      if(num_q_paths>=MAX_NUM){
		node0_printf("OOPS: MAX_NUM too small\n");
		exit(0);
	      }
	      q_paths[num_q_paths].length=length;
	      for(j=0;j<MAX_LENGTH;j++) q_paths[num_q_paths].dir[j]=vec[j];
		/* remember to copy NODIR's, or comparison will get confused */
	      if(ir[0]==0){
		q_paths[num_q_paths].coeff =  coeff;
		q_paths[num_q_paths].forwback =  +1;
	      }
	      else{
		q_paths[num_q_paths].coeff = -coeff;
		q_paths[num_q_paths].forwback = -1;
	      }
	      num_q_paths++;
	      path_num++;
	      /**node0_printf("ADD PATH %d:  rx=%d ",num_q_paths-1,ir[0]);
		 printpath( vec, length );**/
	    }

	  } /* end reflection*/
        } /* end permutation if block */
      } /* end permutation */
    num_basic_paths++;
    return(path_num);
} /* add_basic_path */

/* parallel transport a vector to the current site along a path.
   For example, if the path is "XUP", bring in the vector from
   the +X direction to the current site. If KS phases are in lattice,
   this transport will automatically include them. 
   OK for src and dest to be the same.  OK for length=0.  */
/* NOT OPTIMIZED at the moment - do lots of extra copying.  Use temp.
   vectors rather than stuff in lattice.h */
void path_transport( field_offset src, field_offset dest, int parity,
    int *dir, int length ){
    register int i;
    register site *s;
    msg_tag *mtag0;
    int j;
    su3_vector *tmp_src,*tmp_dest,*tmp_work; /*source, dest and workspace*/
    su3_vector *tmp_pt; /* scratch */
    int tmp_parity, tmp_otherparity; /* parity for this step */

  if( length > 0 ){
    tmp_src = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));
    tmp_dest = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));
    tmp_work = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));

    for( j=length-1; j>=0; j-- ){
	/* figure out parities for this step */
	if( j%2==0 ){
	    tmp_parity = parity;
	    switch(tmp_parity){
		case EVEN: tmp_otherparity=ODD; break;
		case ODD: tmp_otherparity=EVEN; break;
		case EVENANDODD: tmp_otherparity=EVENANDODD; break;
	    }
	}
	else { /* odd # step */
	    tmp_otherparity = parity;
	    switch(tmp_otherparity){
		case EVEN: tmp_parity=ODD; break;
		case ODD: tmp_parity=EVEN; break;
		case EVENANDODD: tmp_parity=EVENANDODD; break;
	    }
	}

	if( j==length-1 ){
	    FORSOMEPARITY(i,s,tmp_otherparity){
	        tmp_src[i] = *(su3_vector *)F_PT(s,src);
	    }
	}

	if( GOES_FORWARDS(dir[j]) ) {
	    mtag0 = start_gather_from_temp( tmp_src, sizeof(su3_vector),
	        dir[j], tmp_parity, gen_pt[0] );
	    wait_gather(mtag0);
	    FORSOMEPARITY(i,s,tmp_parity){
		mult_su3_mat_vec( &(s->link[dir[j]]),
		    (su3_vector *)(gen_pt[0][i]),
		    &(tmp_dest[i]) );
	    }
	    cleanup_gather(mtag0);
	}

	else{ /* GOES_BACKWARDS(dir[j]) */
	    FORSOMEPARITY(i,s,tmp_otherparity){
		mult_adj_su3_mat_vec( &(s->link[OPP_DIR(dir[j])]),
		    &(tmp_src[i]), &(tmp_work[i]) );
	    }
	    mtag0 = start_gather_from_temp( tmp_work, sizeof(su3_vector),
	        dir[j], tmp_parity, gen_pt[0] );
	    wait_gather(mtag0);
	    FORSOMEPARITY(i,s,tmp_parity){
		 tmp_dest[i] = *(su3_vector *)gen_pt[0][i];
	    }
	    cleanup_gather(mtag0);
	}
	
	/* src for next step is dest for this one. */
	tmp_pt=tmp_src; tmp_src=tmp_dest; tmp_dest=tmp_pt;
    }  /* j=link in path */
    /* done, copy result into real dest. (tmp_src now points to result) */
    FORSOMEPARITY(i,s,parity){
        *(su3_vector *)F_PT(s,dest) = tmp_src[i];
    }
    free(tmp_src); free(tmp_dest); free(tmp_work);
  } /* end if(length>0) */
  else if( src != dest ){ /* for length=0 */
    FORSOMEPARITY(i,s,parity){
        *(su3_vector *)F_PT(s,dest) = *(su3_vector *)F_PT(s,src);
    }
  }
} /* path_transport */

/* Path transport a half_wilson_vector */
void path_transport_hwv( field_offset src, field_offset dest, int parity,
    int *dir, int length ){
    register int i;
    register site *s;
    msg_tag *mtag0;
    int j;
    half_wilson_vector *tmp_src,*tmp_dest,*tmp_work; /*source, dest and workspace*/
    half_wilson_vector *tmp_pt; /* scratch */
    int tmp_parity, tmp_otherparity; /* parity for this step */

  if( length > 0 ){
    tmp_src = (half_wilson_vector *)calloc(
      sites_on_node, sizeof(half_wilson_vector));
    tmp_dest = (half_wilson_vector *)calloc(
       sites_on_node, sizeof(half_wilson_vector));
    tmp_work = (half_wilson_vector *)calloc(
       sites_on_node, sizeof(half_wilson_vector));

    for( j=length-1; j>=0; j-- ){
	/* figure out parities for this step */
	if( j%2==0 ){
	    tmp_parity = parity;
	    switch(tmp_parity){
		case EVEN: tmp_otherparity=ODD; break;
		case ODD: tmp_otherparity=EVEN; break;
		case EVENANDODD: tmp_otherparity=EVENANDODD; break;
	    }
	}
	else { /* odd # step */
	    tmp_otherparity = parity;
	    switch(tmp_otherparity){
		case EVEN: tmp_parity=ODD; break;
		case ODD: tmp_parity=EVEN; break;
		case EVENANDODD: tmp_parity=EVENANDODD; break;
	    }
	}

	if( j==length-1 ){
	    FORSOMEPARITY(i,s,tmp_otherparity){
	        tmp_src[i] = *(half_wilson_vector *)F_PT(s,src);
	    }
	}

	if( GOES_FORWARDS(dir[j]) ) {
	    mtag0 = start_gather_from_temp( tmp_src,
		sizeof(half_wilson_vector), dir[j], tmp_parity, gen_pt[0] );
	    wait_gather(mtag0);
	    FORSOMEPARITY(i,s,tmp_parity){
		mult_su3_mat_hwvec( &(s->link[dir[j]]),
		    (half_wilson_vector *)(gen_pt[0][i]),
		    &(tmp_dest[i]) );
	    }
	    cleanup_gather(mtag0);
	}

	else{ /* GOES_BACKWARDS(dir[j]) */
	    FORSOMEPARITY(i,s,tmp_otherparity){
		mult_adj_su3_mat_hwvec( &(s->link[OPP_DIR(dir[j])]),
		    &(tmp_src[i]), &(tmp_work[i]) );
	    }
	    mtag0 = start_gather_from_temp( tmp_work,
		sizeof(half_wilson_vector), dir[j], tmp_parity, gen_pt[0] );
	    wait_gather(mtag0);
	    FORSOMEPARITY(i,s,tmp_parity){
		 tmp_dest[i] = *(half_wilson_vector *)gen_pt[0][i];
	    }
	    cleanup_gather(mtag0);
	}
	
	/* src for next step is dest for this one. */
	tmp_pt=tmp_src; tmp_src=tmp_dest; tmp_dest=tmp_pt;
    }  /* j=link in path */
    /* done, copy result into real dest. (tmp_src now points to result) */
    FORSOMEPARITY(i,s,parity){
        *(half_wilson_vector *)F_PT(s,dest) = tmp_src[i];
    }
    free(tmp_src); free(tmp_dest); free(tmp_work);
  } /* end if(length>0) */
  else if( src != dest ){ /* for length=0 */
    FORSOMEPARITY(i,s,parity){
        *(half_wilson_vector *)F_PT(s,dest) =
	  *(half_wilson_vector *)F_PT(s,src);
    }
  }
} /* path_transport_hwv */


#ifdef EO
/* Stupid dslash routine that follows all the paths.
  Should optimize by precomputing sums of paths to all the displacements
  Use dslash_fn for actions that involve only +X and +X+X+X couplings.
*/
/* D_slash routine - sets dest. on each site equal to sum of
   sources parallel transported to site, with minus sign for transport
   from negative directions.  */
void dslash_eo( field_offset src, field_offset dest, int parity ) {
   register int i;
   register site *s;
   register int ipath,otherparity;
   register double x; /* coefficient of path */

    switch(parity){
	case EVEN:	otherparity=ODD; break;
	case ODD:	otherparity=EVEN; break;
	case EVENANDODD:	otherparity=EVENANDODD; break;
    }

    /* Parallel transport by all the paths in the action.  
       Multiply by coefficient in table
    */
    FORSOMEPARITY(i,s,parity){ clearvec( (su3_vector *)F_PT(s,dest) ); }

    for( ipath=0; ipath<num_q_paths; ipath++ ){  /* loop over paths */
	    path_transport( src, F_OFFSET(tempvec[0]), parity,
		q_paths[ipath].dir, q_paths[ipath].length );
	    x=q_paths[ipath].coeff;
	    FORSOMEPARITY(i,s,parity){
		scalar_mult_add_su3_vector(  (su3_vector *)F_PT(s,dest),
		    &(s->tempvec[0]), x, (su3_vector *)F_PT(s,dest) );
	    }
    }   /* ipath */
} /* dslash_eo */
#endif /* ifdef EO */

#ifdef FN
/* Sum over paths connecting to nearest neighbor point (fat link) and to third
   nearest neighbor (longlinks) */
/* Doug Toussaint 2/4/98 */
/* modified to use t_longlinks, S. Gottlieb 7/13/01 */
/* long link calculating routine */
/* path_product() follows the path starting at step 0, and
   leaves the answer at the end of the path.  We want the answer
   at the site where the path begins.  So we look for paths with
   the opposite displacement from the displacement of the point
   that we want to transport to this site, and take the adjoint
   of the matrix at the end. clear? */
/* KS phases and APBC must be in the links. See long comment at bottom*/
void load_longlinks() {
  register int i;
  register site *s;
  int ipath,dir;
  int disp[4];
  register su3_matrix *long1;

#ifdef LLTIME
double dtime;
dtime=-dclock();
#endif
  if( phases_in != 1){
    node0_printf("BOTCH: load_longlinks needs phases in\n");
    terminate(0);
  }
  for (dir=XUP; dir<=TUP; dir++){ /* loop over longlink directions */
    /* set longlink to zero */
    FORALLSITES(i,s){
#ifdef DSLASH_TMP_LINKS
      long1 = &(t_longlink[4*i+dir]);
#else
      long1 = &(s->longlink[dir]);
#endif
      clear_su3mat( long1 );
    }

    /* loop over paths, checking for ones with total displacement 3*dir */
    for( ipath=0; ipath<num_q_paths; ipath++ ){  /* loop over paths */
	/* compute total displacement of path */
	for(i=XUP;i<=TUP;i++)disp[i]=0;
	for( i=0; i<q_paths[ipath].length; i++){
	  if( GOES_FORWARDS(q_paths[ipath].dir[i]) )
	    disp[        q_paths[ipath].dir[i]  ]++;
	  else
	    disp[OPP_DIR(q_paths[ipath].dir[i]) ]--;
	}
	for( disp[dir]+=3,i=XUP; i<=TUP; i++)if(disp[i]!=0)break;
	if( i<=TUP )continue;  /* skip if path doesn't go to right place */
/**printf("ipath = %d, found a path:  ",ipath);
for(j=0;j<q_paths[ipath].length;j++)printf("\t%d", q_paths[ipath].dir[j]);
printf("\n");**/

	path_product( q_paths[ipath].dir, q_paths[ipath].length );
	FORALLSITES(i,s){
	  su3_adjoint( &(s->tempmat1), &(s->staple) );
#ifdef DSLASH_TMP_LINKS
	  long1 = &(t_longlink[4*i+dir]);
#else
	  long1 = &(s->longlink[dir]);
#endif
          scalar_mult_add_su3_matrix( long1,
	    &(s->staple), -q_paths[ipath].coeff, long1 );
		/* minus sign in coeff. because we used backward path*/
	}
    } /* ipath */

  } /* loop over directions */

  valid_longlinks = 1;
#ifdef LLTIME
dtime += dclock();
node0_printf("LLTIME(Long):  %e\n",dtime);
#endif
}  /* load_longlinks() */

/* KS phases and APBC must be in the links. See long comment at bottom*/
void load_fatlinks() {
  register int i;
  register site *s;
  int dir;
  register su3_matrix *fat1;
#ifdef ASQ_OPTIMIZED_FATTENING
  int  nu,rho,sig ;
  double one_link ; /*needed to fix the problem with the Lepage term */
#else
  int ipath;
  int disp[4];
#endif

#ifdef LLTIME
double dtime;
dtime=-dclock();
#endif
  if( phases_in != 1){
    node0_printf("BOTCH: load_fatlinks needs phases in\n");
    terminate(0);
  }

#ifndef  ASQ_OPTIMIZED_FATTENING   /* general case code */
  for (dir=XUP; dir<=TUP; dir++){ /* loop over fatlink directions */
    /* set fatlink to zero */
    FORALLSITES(i,s){
#ifdef DSLASH_TMP_LINKS
      fat1 = &(t_fatlink[4*i+dir]);
#else
      fat1 = &(s->fatlink[dir]);
#endif
      clear_su3mat( fat1 );
    }
    
    /* loop over paths, checking for ones with total displacement 1*dir */
    for( ipath=0; ipath<num_q_paths; ipath++ ){  /* loop over paths */
	/* compute total displacement of path */
	for(i=XUP;i<=TUP;i++)disp[i]=0;
	for( i=0; i<q_paths[ipath].length; i++){
	  if( GOES_FORWARDS(q_paths[ipath].dir[i]) )
	    disp[        q_paths[ipath].dir[i]  ]++;
	  else
	    disp[OPP_DIR(q_paths[ipath].dir[i]) ]--;
	}
	for( disp[dir]+=1,i=XUP; i<=TUP; i++)if(disp[i]!=0)break;
	if( i<=TUP )continue;  /* skip if path doesn't go to right place */
/**printf("dir = %d, found a path:  ",dir);
for(j=0;j<q_paths.[ipath].length;j++)printf("\t%d", q_paths[ipath].dir[j]);
printf("\n");**/

	path_product( q_paths[ipath].dir, q_paths[ipath].length );
	FORALLSITES(i,s){
	  su3_adjoint( &(s->tempmat1), &(s->staple) );
#ifdef DSLASH_TMP_LINKS
	  fat1 = &(t_fatlink[4*i+dir]);
#else
	  fat1 = &(s->fatlink[dir]);
#endif
          scalar_mult_add_su3_matrix( fat1,
	    &(s->staple), -q_paths[ipath].coeff, fat1 );
		/* minus sign in coeff. because we used backward path*/
	}
    } /* ipath */
  } /* loop over directions */
#else	/* ASQ_OPTIMIZED_FATTENING, for Asq and Asqtad actions */
/*  Optimized fattening code for the Asq and Asqtad actions.           *
 * I assume that path 0 is the one link path 2 the 3-staple            *
 * path 3 the 5-staple path 4 the 7-staple and path 5 the Lapage term. *
 * Path 1 is the Naik term.                                            */
 
 /* to fix up the Lepage term, included by a trick below */
 one_link = (act_path_coeff[0] - 6.0*act_path_coeff[5]) ; 
 
 for (dir=XUP; dir<=TUP; dir++){
   FORALLSITES(i,s) /* Intialize fat links with c_1*U_\mu(x) */
     {
#ifdef DSLASH_TMP_LINKS
       fat1 = &(t_fatlink[4*i+dir]);
#else
       fat1 = &(s->fatlink[dir]);
#endif
       scalar_mult_su3_matrix(&(s->link[dir]), one_link ,
			      fat1 );
     }
   for(nu=XUP; nu<=TUP; nu++) if(nu!=dir)
     {
       compute_gen_staple(F_OFFSET(staple),dir,nu,F_OFFSET(link[dir]),
			  act_path_coeff[2]);
       /* The Lepage term */
       /* Note this also involves modifying c_1 (above) */
       compute_gen_staple(NULL_FP,dir,nu,F_OFFSET(staple),act_path_coeff[5]);
       for(rho=XUP; rho<=TUP; rho++) if((rho!=dir)&&(rho!=nu))
	 {
	   compute_gen_staple(F_OFFSET(tempmat1),dir,rho,F_OFFSET(staple),
			      act_path_coeff[3]);
	   for(sig=XUP; sig<=TUP; sig++)
	     if((sig!=dir)&&(sig!=nu)&&(sig!=rho))
	       {
		 compute_gen_staple(NULL_FP,dir,sig,
				    F_OFFSET(tempmat1),
				    act_path_coeff[4]);
	       } /* sig */
	 } /* rho */
     } /* nu */
 }/* dir */  
#endif

  valid_fatlinks = 1;
#ifdef LLTIME
dtime += dclock();
node0_printf("LLTIME(Fat):  %e\n",dtime);
#endif
}  /* load_fatlinks() */
#endif /* ifdef FN */


/* compare two paths, return 1 if equal, else zero */
int is_path_equal( int *path1, int* path2, int length ){
   register int i;
   for(i=0;i<length;i++)if(path1[i]!=path2[i])return(0);
   return(1);
}


/* update the  momenta with the fermion force */
/* Assumes that the conjugate gradient has been run, with the answer in
   x_off, and dslash(x_off,x_off,ODD) has been run. (fills in x_off_odd) */
/* SEE LONG COMMENTS AT END */

#ifndef ASQ_OPTIMIZED_FORCE
void eo_fermion_force( double eps, int nflavors, field_offset x_off ){
  /* note CG_solution and Dslash * solution are combined in "x_off" */
  /* New version 1/21/99.  Use forward part of Dslash to get force */
  /* see long comment at end */
  /* For each link we need x_off transported from both ends of path. */
  register int i,dir,lastdir,ipath,ilink;
  register site *s;
  int length;
  su3_matrix tmat,tmat2;
  double ferm_epsilon, coeff;
#ifdef FFTIME
  double dtime;
#endif
  msg_tag *mtag0;
  half_wilson_vector *hw_tmp0,*hw_tmp1,*tmp_pt;

#ifdef FFTIME
dtime=-dclock();
#endif
  ferm_epsilon = 2.0*(nflavors/4.0)*eps;
  hw_tmp0 = (half_wilson_vector *)
    calloc(sites_on_node, sizeof(half_wilson_vector));
  hw_tmp1 = (half_wilson_vector *)
    calloc(sites_on_node, sizeof(half_wilson_vector));
  /* Use half_wilson_vectors to store x_off transported from ends of
     path.  0 component from forward end, 1 component from back end */

  /* loop over paths, and loop over links in path */
  for( ipath=0; ipath<num_q_paths; ipath++ ){
    if(q_paths[ipath].forwback== -1)continue;	/* skip backwards dslash */
    length = q_paths[ipath].length;

    /* path transport x_off and Dslash*x_off from far end.  Sometimes
	we need them at the start point of the path, and sometimes
	one link into the path - an optimization for later */
    path_transport( x_off, F_OFFSET(tempvec[0]),
      EVENANDODD, q_paths[ipath].dir, length );
    /* use tempvec[1] for transport from starting end */
    FORALLSITES(i,s){
      hw_tmp0[i].h[0]=s->tempvec[0];
      hw_tmp0[i].h[1]=*(su3_vector *)F_PT(s,x_off);
    }

    /* A path has (length+1) points, counting the ends.  At first
	 point, no "down" direction links have their momenta "at this
	 point". At last, no "up" ... */
    for( ilink=0; ilink<=length; ilink++ ){
      if(ilink<length)dir = q_paths[ipath].dir[ilink];
      else dir=NODIR;
      coeff = ferm_epsilon*q_paths[ipath].coeff;
      if( (ilink%2)==1 )coeff = -coeff;

      /* path transport x_off and Dslash*x_off from previous point */
      /* Use "half_wilson_vector" to handle pair of vectors -
	0 component is x_off from forward end, 1 component from back end */
      /* sometimes we don't need them */
      if( (ilink>0&&ilink<length) || 
	(ilink==length && GOES_BACKWARDS(lastdir)) ){

	if( GOES_FORWARDS(lastdir) ){
	  FORALLSITES(i,s){
            mult_adj_su3_mat_hwvec( &(s->link[lastdir]),
	      &(hw_tmp0[i]), &(hw_tmp1[i]) );
	  }
	  mtag0 = start_gather_from_temp( hw_tmp1, 2*sizeof(su3_vector),
             OPP_DIR(lastdir), EVENANDODD, gen_pt[0] );
          wait_gather(mtag0);
          FORALLSITES(i,s){
	     hw_tmp0[i] = *(half_wilson_vector *)gen_pt[0][i];
	  }
          cleanup_gather(mtag0);
	}
	else{   /* GOES_BACKWARDS(lastdir) */
          mtag0 = start_gather_from_temp( hw_tmp0, 2*sizeof(su3_vector),
               OPP_DIR(lastdir), EVENANDODD, gen_pt[0] );
          wait_gather(mtag0);
          FORALLSITES(i,s){
            mult_su3_mat_hwvec( &(s->link[OPP_DIR(lastdir)]),
                (half_wilson_vector *)(gen_pt[0][i]),
		&(hw_tmp1[i]) );
          }
	  tmp_pt = hw_tmp0; hw_tmp0 = hw_tmp1; hw_tmp1 = tmp_pt;
          cleanup_gather(mtag0);
	}
      }

      /* add in contribution to the force */
      /* Put antihermitian traceless part into momentum */
      FORALLSITES(i,s){
        if( ilink<length && GOES_FORWARDS(dir) ){
          uncompress_anti_hermitian( &(s->mom[dir]), &tmat2 );
	  su3_projector( &(hw_tmp0[i].h[0]), &(hw_tmp0[i].h[1]), &tmat );
	  if( s->parity==EVEN )scalar_mult_add_su3_matrix(
	    &tmat2, &tmat,  coeff, &tmat2 );
	  else		 scalar_mult_add_su3_matrix(
	    &tmat2, &tmat, -coeff, &tmat2 );
          make_anti_hermitian( &tmat2, &(s->mom[dir]) );
        }
	if( ilink>0 && GOES_BACKWARDS(lastdir) ){
          uncompress_anti_hermitian( &(s->mom[OPP_DIR(lastdir)]), &tmat2 );
	  su3_projector( &(hw_tmp0[i].h[0]), &(hw_tmp0[i].h[1]), &tmat );
	  if( s->parity==EVEN )scalar_mult_add_su3_matrix(
	    &tmat2, &tmat, -coeff, &tmat2 );
	  else		 scalar_mult_add_su3_matrix(
	    &tmat2, &tmat,  coeff, &tmat2 );
          make_anti_hermitian( &tmat2, &(s->mom[OPP_DIR(lastdir)]) );
        }
      }
      lastdir = dir;
    } /* end loop over links in path */
  } /* end loop over paths */

  free( hw_tmp0 ); free( hw_tmp1 );

#ifdef FFTIME
dtime += dclock();
node0_printf("FFTIME:  %e\n",dtime);
/**printf("TLENGTH: %d\n",tlength);**/
#endif
} /* eo_fermion_force(version 6) */


void eo_fermion_force_3f( double eps, int nflav1, field_offset x1_off,
  int nflav2, field_offset x2_off ){
  /* note CG_solution and Dslash * solution are combined in "x_off" */
  /* New version 1/21/99.  Use forward part of Dslash to get force */
  /* 4/14/99 combine force from two different mass quarks, (eg 2+1flavors) */
  /* see long comment at end */
  /* For each link we need x_off transported from both ends of path. */
  register int i,dir,lastdir,ipath,ilink;
  register site *s;
  int length;
  su3_matrix tmat,tmat2;
  double ferm_epsilon1, ferm_epsilon2, coeff1, coeff2;
#ifdef FFTIME
  double dtime;
#endif
  msg_tag *mtag0;
  wilson_vector *w_tmp0,*w_tmp1,*tmp_pt;

#ifdef FFTIME
dtime=-dclock();
#endif
  ferm_epsilon1 = 2.0*(nflav1/4.0)*eps;
  ferm_epsilon2 = 2.0*(nflav2/4.0)*eps;
  w_tmp0 = (wilson_vector *)
    calloc(sites_on_node, sizeof(wilson_vector));
  w_tmp1 = (wilson_vector *)
    calloc(sites_on_node, sizeof(wilson_vector));
  /* Use wilson_vectors to store x_off transported from ends of
     path.  0 and 1 components from forward end, 2 and 3 components
     from back end */

  /* loop over paths, and loop over links in path */
  for( ipath=0; ipath<num_q_paths; ipath++ ){
    if(q_paths[ipath].forwback== -1)continue;	/* skip backwards dslash */
    length = q_paths[ipath].length;

    /* path transport x_off and Dslash*x_off from far end.  Sometimes
	we need them at the start point of the path, and sometimes
	one link into the path - an optimization for later */
    /**
    path_transport( x1_off, F_OFFSET(tempvec[0]),
      EVENANDODD, q_paths[ipath].dir, length );
    path_transport( x2_off, F_OFFSET(tempvec[1]),
      EVENANDODD, q_paths[ipath].dir, length );
    **/
/** WARNING!! Assumes xxx1 and xxx2 contiguous **/
if( x2_off-x1_off != sizeof(su3_vector) ){node0_printf("BOTCH\n"); exit(0);}
    path_transport_hwv( x1_off, F_OFFSET(tempvec[0]),
      EVENANDODD, q_paths[ipath].dir, length );
    /* use tempvec[2] for transport from starting end */
    FORALLSITES(i,s){
      w_tmp0[i].d[0]=s->tempvec[0];
      w_tmp0[i].d[1]=s->tempvec[1];
      w_tmp0[i].d[2]=*(su3_vector *)F_PT(s,x1_off);
      w_tmp0[i].d[3]=*(su3_vector *)F_PT(s,x2_off);
    }

    /* A path has (length+1) points, counting the ends.  At first
	 point, no "down" direction links have their momenta "at this
	 point". At last, no "up" ... */
    for( ilink=0; ilink<=length; ilink++ ){
      if(ilink<length)dir = q_paths[ipath].dir[ilink];
      else dir=NODIR;
      coeff1 = ferm_epsilon1*q_paths[ipath].coeff;
      coeff2 = ferm_epsilon2*q_paths[ipath].coeff;
      if( (ilink%2)==1 ){ coeff1 = -coeff1; coeff2 = -coeff2;}

      /* path transport x_off and Dslash*x_off from previous point */
      /* Use "wilson_vector" to handle pair of vectors -
	0 component is x_off1 from forward end, 1 component is x_off2
	from forward end,  2 and 3  components are x_off1 and x_off2 
	from back end */
      /* sometimes we don't need them */
      if( (ilink>0&&ilink<length) || 
	(ilink==length && GOES_BACKWARDS(lastdir)) ){

	if( GOES_FORWARDS(lastdir) ){
	  FORALLSITES(i,s){
            mult_adj_mat_wilson_vec( &(s->link[lastdir]),
	      &(w_tmp0[i]), &(w_tmp1[i]) );
	  }
	  mtag0 = start_gather_from_temp( w_tmp1, 4*sizeof(su3_vector),
             OPP_DIR(lastdir), EVENANDODD, gen_pt[0] );
          wait_gather(mtag0);
          FORALLSITES(i,s){
	     w_tmp0[i] = *(wilson_vector *)gen_pt[0][i];
	  }
          cleanup_gather(mtag0);
	}
	else{   /* GOES_BACKWARDS(lastdir) */
          mtag0 = start_gather_from_temp( w_tmp0, 4*sizeof(su3_vector),
               OPP_DIR(lastdir), EVENANDODD, gen_pt[0] );
          wait_gather(mtag0);
          FORALLSITES(i,s){
            mult_mat_wilson_vec( &(s->link[OPP_DIR(lastdir)]),
                (wilson_vector *)(gen_pt[0][i]),
		&(w_tmp1[i]) );
          }
	  tmp_pt = w_tmp0; w_tmp0 = w_tmp1; w_tmp1 = tmp_pt;
          cleanup_gather(mtag0);
	}
      }

      /* add in contribution to the force */
      /* Put antihermitian traceless part into momentum */
      FORALLSITES(i,s){
        if( ilink<length && GOES_FORWARDS(dir) ){
	  if( s->parity==ODD ){coeff1 *= -1.0; coeff2 *= -1.0; }
          uncompress_anti_hermitian( &(s->mom[dir]), &tmat2 );
	  su3_projector( &(w_tmp0[i].d[0]), &(w_tmp0[i].d[2]), &tmat );
	  scalar_mult_add_su3_matrix( &tmat2, &tmat,  coeff1, &tmat2 );
	  su3_projector( &(w_tmp0[i].d[1]), &(w_tmp0[i].d[3]), &tmat );
	  scalar_mult_add_su3_matrix( &tmat2, &tmat,  coeff2, &tmat2 );
          make_anti_hermitian( &tmat2, &(s->mom[dir]) );
	  if( s->parity==ODD ){coeff1 *= -1.0; coeff2 *= -1.0; }
        }
	if( ilink>0 && GOES_BACKWARDS(lastdir) ){
	  if( s->parity==EVEN ){coeff1 *= -1.0; coeff2 *= -1.0; }
          uncompress_anti_hermitian( &(s->mom[OPP_DIR(lastdir)]), &tmat2 );
	  su3_projector( &(w_tmp0[i].d[0]), &(w_tmp0[i].d[2]), &tmat );
	  scalar_mult_add_su3_matrix( &tmat2, &tmat,  coeff1, &tmat2 );
	  su3_projector( &(w_tmp0[i].d[1]), &(w_tmp0[i].d[3]), &tmat );
	  scalar_mult_add_su3_matrix( &tmat2, &tmat,  coeff2, &tmat2 );
          make_anti_hermitian( &tmat2, &(s->mom[OPP_DIR(lastdir)]) );
	  if( s->parity==EVEN ){coeff1 *= -1.0; coeff2 *= -1.0; }
        }
      }
      lastdir = dir;
    } /* end loop over links in path */
  } /* end loop over paths */

  free( w_tmp0 ); free( w_tmp1 );

#ifdef FFTIME
dtime += dclock();
node0_printf("FFTIME:  %e\n",dtime);
/**printf("TLENGTH: %d\n",tlength);**/
#endif
} /* eo_fermion_force_3f(version 6) */



#else /* ASQ_OPTIMIZED_FORCE, for Asq and Asqtad actions */
/* Optimized force code for the Asq and Asqtad actions                 *
 * I assume that path 0 is the one link path 2 the 3-staple            *
 * path 3 the 5-staple path 4 the 7-staple and path 5 the Lapage term. *
 * Path 1 is the Naik term.                                            */
#define Pmu          tempvec[0] 
#define Pnumu        tempvec[1]
#define Prhonumu     tempvec[2]
#define P7           tempvec[3]
#define P7rho        tempvec[4]              
#define P7rhonu      tempvec[5]
#define P5           tempvec[6]
#define P3           tempvec[7]
#define P5nu         tempvec[3]
#define P3mu         tempvec[3]
#define Popmu        tempvec[4]
#define Pmumumu      tempvec[4]
void eo_fermion_force( double eps, int nflavors, field_offset x_off ){
  /* note CG_solution and Dslash * solution are combined in "x_off" */
  /* New version 1/21/99.  Use forward part of Dslash to get force */
  /* see long comment at end */
  /* For each link we need x_off transported from both ends of path. */
  register int i ;
  register site *s;
  int mu,nu,rho,sig ;
  int DirectLinks[8] ;
  double ferm_epsilon, coeff;
  double OneLink, Lepage, Naik, FiveSt, ThreeSt, SevenSt ;
  su3_vector *tempvec[8] ;
  su3_vector *temp_x ;

#ifdef FFTIME
  double dtime;

  dtime=-dclock();
#endif
  ferm_epsilon = 2.0*(nflavors/4.0)*eps;
  
  /* Path coefficients times fermion epsilon */
  OneLink = act_path_coeff[0]*ferm_epsilon ; 
  Naik    = act_path_coeff[1]*ferm_epsilon ;
  ThreeSt = act_path_coeff[2]*ferm_epsilon ;
  FiveSt  = act_path_coeff[3]*ferm_epsilon ;
  SevenSt = act_path_coeff[4]*ferm_epsilon ;
  Lepage  = act_path_coeff[5]*ferm_epsilon ;
  /* *************************************** */

  /* Initialize the DirectLink flags */
  for(mu=0;mu<8;mu++)
    DirectLinks[mu] = 0 ;

  /* Allocate temporary vectors */
  for(mu=0;mu<8;mu++)
    tempvec[mu] = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));

  /*copy x_off to a temporary vector */
  temp_x = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));
  FORALLSITES(i,s) temp_x[i] = *(su3_vector *)F_PT(s,x_off) ;

  for(sig=0;sig<8;sig++)
    {
      for(mu=0;mu<8;mu++)if((mu!=sig)&&(mu!=OPP_DIR(sig)))
	{
	  u_shift_fermion(temp_x, Pmu, OPP_DIR(mu));
	  u_shift_fermion(Pmu, P3, sig);
	  if(GOES_FORWARDS(sig))
	    {
	      /* Add the force F_sig[x+mu]:         x--+             *
	       *                                   |   |             *
	       *                                   o   o             *
	       * the 1 link in the path: - (numbering starts form 0) */
	      add_force_to_mom(P3, Pmu, sig, -ThreeSt) ;
	    }
	  for(nu=0;nu<8;nu++)if((nu!=mu )&&(nu!=OPP_DIR(mu ))&&
				(nu!=sig)&&(nu!=OPP_DIR(sig)))
	    {
	      u_shift_fermion(Pmu, Pnumu, OPP_DIR(nu));
	      u_shift_fermion(Pnumu, P5, sig);
	      if(GOES_FORWARDS(sig))
		{
		  /* Add the force F_sig[x+mu+nu]:      x--+             *
		   *                                   |   |             *
		   *                                   o   o             *
		   * the 2 link in the path: + (numbering starts form 0) */
		  add_force_to_mom(P5, Pnumu, sig, FiveSt);
		}
	      for(rho=0;rho<8;rho++)if((rho!=mu )&&(rho!=OPP_DIR(mu ))&&
				       (rho!=nu )&&(rho!=OPP_DIR(nu ))&&
				       (rho!=sig)&&(rho!=OPP_DIR(sig)))
		{
		  u_shift_fermion(Pnumu, Prhonumu, OPP_DIR(rho));
		  /* Length 7 paths */
		  u_shift_fermion(Prhonumu, P7,sig);
		  if(GOES_FORWARDS(sig))
		    {
		      /* Add the force F_sig[x+mu+nu+rho]:  x--+             *
		       *                                   |   |             *
		       *                                   o   o             *
		       * the 3 link in the path: - (numbering starts form 0) */
		      add_force_to_mom(P7, Prhonumu, sig, -SevenSt ) ;
		    }
		  /*Add the force F_rho the 2(4) link in the path: +     */
		  u_shift_fermion(P7, P7rho, rho);
		  side_link_force(rho,sig,SevenSt, Pnumu, P7, Prhonumu, P7rho);
		  /* Add the P7rho vector to P5 */
		  coeff = SevenSt/FiveSt ; 
		  FORALLSITES(i,s)
		    scalar_mult_add_su3_vector(&(P5[i]),&(P7rho[i]),coeff,
					       &(P5[i]));
		}/* rho */
	      /* Length 5 paths */
	      /*Add the force F_nu the 1(3) link in the path: -     */
	      u_shift_fermion(P5,P5nu, nu);
	      side_link_force(nu,sig,-FiveSt,Pmu,P5, 
			      Pnumu,P5nu) ;
	      /* Add the P5nu vector to P3 */
	      coeff = FiveSt/ThreeSt ; 
	      FORALLSITES(i,s)
		scalar_mult_add_su3_vector(&(P3[i]),&(P5nu[i]),coeff,&(P3[i]));
	    }/* nu */

	  /* Now the Lepage term... It is the same with 5-link paths with
             nu=mu and FiveSt=Lepage. So Pnumu is really Pmumu */
	  u_shift_fermion(Pmu, Pnumu, OPP_DIR(mu));
	  u_shift_fermion(Pnumu, P5, sig);
	  if(GOES_FORWARDS(sig))
	    {
	      /* Add the force F_sig[x+mu+nu]:      x--+             *
	       *                                   |   |             *
	       *                                   o   o             *
	       * the 2 link in the path: + (numbering starts form 0) */
	      add_force_to_mom(P5, Pnumu, sig, Lepage) ;
	    }
	  /*Add the force F_nu the 1(3) link in the path: -     */
	  u_shift_fermion(P5,P5nu, mu);
	  side_link_force(mu, sig, -Lepage, Pmu, P5, Pnumu, P5nu) ;
	  /* Add the P5nu vector to P3 */
	  coeff = Lepage/ThreeSt ; 
	  FORALLSITES(i,s)
	    scalar_mult_add_su3_vector(&(P3[i]),&(P5nu[i]),coeff,&(P3[i]));

	  /* Length 3 paths (Not the Naik term) */
	  /*Add the force F_mu the 0(2) link in the path: +     */
	  if(GOES_FORWARDS(mu)) 
	    u_shift_fermion(P3,P3mu, mu );
	  /* The above shift is not needed if mu is backwards */
	  side_link_force(mu, sig, ThreeSt, temp_x, P3, Pmu, P3mu);

	  /* Finally the OneLink and the Naik term */
	  /* Check if this direction is not already done */
	  if( (!DirectLinks[mu]) ){
	    DirectLinks[mu]=1 ;
	    if(GOES_BACKWARDS(mu))/* Do only the forward terms in the Dslash */
	      {
		/* Because I have shifted with OPP_DIR(mu) Pmu is a forward *
		 * shift.                                                   */
		/* The one link */
		add_force_to_mom(Pmu, temp_x, OPP_DIR(mu), OneLink) ;
		/* For the same reason Pnumu is the forward double link */

		/* Popmu is a backward shift */
		u_shift_fermion(temp_x, Popmu, mu);
		/* The Naik */
		/* link no 1: - */
		add_force_to_mom(Pnumu, Popmu, OPP_DIR(mu), -Naik) ;
		/*Pmumumu can overwrite Popmu which is no longer needed */
		u_shift_fermion(Pnumu, Pmumumu, OPP_DIR(mu));
		/* link no 0: + */
		add_force_to_mom(Pmumumu, temp_x, OPP_DIR(mu), Naik);
	      }
	    else /* The rest of the Naik terms */
	      {
		u_shift_fermion(temp_x, Popmu, mu);
		/* link no 2: + */
		/* Pnumu is double backward shift */
		add_force_to_mom(Popmu, Pnumu, mu, Naik) ;
	      }
	  }
	}/* mu */
      /* Here we have to do together the Naik term and the one link term */
    }/*sig */

  /* Free temporary vectors */
  free(temp_x) ;
  for(mu=0;mu<8;mu++)
    free(tempvec[mu]) ;

#ifdef FFTIME
  dtime += dclock();
node0_printf("FFTIME:  %e\n",dtime);
/**printf("TLENGTH: %d\n",tlength);**/
#endif
} /* eo_fermion_force(version 6) */
#undef Pmu          
#undef Pnumu        
#undef Prhonumu     
#undef P7           
#undef P7rho        
#undef P7rhonu      
#undef P5           
#undef P3           
#undef P5nu         
#undef P3mu         
#undef Popmu        
#undef Pmumumu      

#define Pmu          hwvec[0] 
#define Pnumu        hwvec[1]
#define Prhonumu     hwvec[2]
#define P7           hwvec[3]
#define P7rho        hwvec[4]              
#define P7rhonu      hwvec[5]
#define P5           hwvec[6]
#define P3           hwvec[7]
#define P5nu         hwvec[3]
#define P3mu         hwvec[3]
#define Popmu        hwvec[4]
#define Pmumumu      hwvec[4]
void eo_fermion_force_3f( double eps, int nflav1, field_offset x1_off, 
			  int nflav2, field_offset x2_off){
  /* note CG_solution and Dslash * solution are combined in "x_off" */
  /* New version 1/21/99.  Use forward part of Dslash to get force */
  /* 4/15/99 combine force from two different mass quarks, (eg 2+1flavors) */
  /* see long comment at end */
  /* For each link we need x_off transported from both ends of path. */
  register int i ;
  register site *s;
  int mu,nu,rho,sig ;
  int DirectLinks[8] ;
  double coeff[2],ferm_epsilon ;
  double OneLink[2], Lepage[2], Naik[2], FiveSt[2], ThreeSt[2], SevenSt[2] ;
  double mNaik[2], mLepage[2], mFiveSt[2], mThreeSt[2], mSevenSt[2] ;
  half_wilson_vector *hwvec[8] ;
  half_wilson_vector *temp_x ;

#ifdef FFTIME
  double dtime;

  dtime=-dclock();
#endif
  
  /* Path coefficients times fermion epsilon */
  ferm_epsilon = 2.0*(nflav1/4.0)*eps;
  OneLink[0] = act_path_coeff[0]*ferm_epsilon ;
  Naik[0]    = act_path_coeff[1]*ferm_epsilon ; mNaik[0]    = -Naik[0] ;
  ThreeSt[0] = act_path_coeff[2]*ferm_epsilon ; mThreeSt[0] = -ThreeSt[0] ;
  FiveSt[0]  = act_path_coeff[3]*ferm_epsilon ; mFiveSt[0]  = -FiveSt[0] ;
  SevenSt[0] = act_path_coeff[4]*ferm_epsilon ; mSevenSt[0] = -SevenSt[0] ;
  Lepage[0]  = act_path_coeff[5]*ferm_epsilon ; mLepage[0]  = -Lepage[0] ;

  ferm_epsilon = 2.0*(nflav2/4.0)*eps;
  OneLink[1] = act_path_coeff[0]*ferm_epsilon ;
  Naik[1]    = act_path_coeff[1]*ferm_epsilon ; mNaik[1]    = -Naik[1] ;
  ThreeSt[1] = act_path_coeff[2]*ferm_epsilon ; mThreeSt[1] = -ThreeSt[1] ;
  FiveSt[1]  = act_path_coeff[3]*ferm_epsilon ; mFiveSt[1]  = -FiveSt[1] ;
  SevenSt[1] = act_path_coeff[4]*ferm_epsilon ; mSevenSt[1] = -SevenSt[1] ;
  Lepage[1]  = act_path_coeff[5]*ferm_epsilon ; mLepage[1]  = -Lepage[1] ;
  /* *************************************** */

  /* Initialize the DirectLink flags */
  for(mu=0;mu<8;mu++)
    DirectLinks[mu] = 0 ;

  /* Allocate temporary vectors */
  for(mu=0;mu<8;mu++)
    hwvec[mu]= 
      (half_wilson_vector *)calloc(sites_on_node, sizeof(half_wilson_vector));

  /*copy x_off to a temporary vector */
  temp_x= 
    (half_wilson_vector *)calloc(sites_on_node, sizeof(half_wilson_vector));
  FORALLSITES(i,s)
    {
      temp_x[i].h[0] = *(su3_vector *)F_PT(s,x1_off) ;
      temp_x[i].h[1] = *(su3_vector *)F_PT(s,x2_off) ;
    }

  for(sig=0;sig<8;sig++)
    {
      for(mu=0;mu<8;mu++)if((mu!=sig)&&(mu!=OPP_DIR(sig)))
	{
	  u_shift_hw_fermion(temp_x, Pmu, OPP_DIR(mu));
	  u_shift_hw_fermion(Pmu, P3, sig);
	  if(GOES_FORWARDS(sig))
	    {
	      /* Add the force F_sig[x+mu]:         x--+             *
	       *                                   |   |             *
	       *                                   o   o             *
	       * the 1 link in the path: - (numbering starts form 0) */
	      add_3f_force_to_mom(P3, Pmu, sig, mThreeSt) ;
	    }
	  for(nu=0;nu<8;nu++)if((nu!=mu )&&(nu!=OPP_DIR(mu ))&&
				(nu!=sig)&&(nu!=OPP_DIR(sig)))
	    {
	      u_shift_hw_fermion(Pmu, Pnumu, OPP_DIR(nu));
	      u_shift_hw_fermion(Pnumu, P5, sig);
	      if(GOES_FORWARDS(sig))
		{
		  /* Add the force F_sig[x+mu+nu]:      x--+             *
		   *                                   |   |             *
		   *                                   o   o             *
		   * the 2 link in the path: + (numbering starts form 0) */
		  add_3f_force_to_mom(P5, Pnumu, sig, FiveSt);
		}
	      for(rho=0;rho<8;rho++)if((rho!=mu )&&(rho!=OPP_DIR(mu ))&&
				       (rho!=nu )&&(rho!=OPP_DIR(nu ))&&
				       (rho!=sig)&&(rho!=OPP_DIR(sig)))
		{
		  u_shift_hw_fermion(Pnumu, Prhonumu, OPP_DIR(rho));
		  /* Length 7 paths */
		  u_shift_hw_fermion(Prhonumu, P7,sig);
		  if(GOES_FORWARDS(sig))
		    {
		      /* Add the force F_sig[x+mu+nu+rho]:  x--+             *
		       *                                   |   |             *
		       *                                   o   o             *
		       * the 3 link in the path: - (numbering starts form 0) */
		      add_3f_force_to_mom(P7, Prhonumu, sig, mSevenSt ) ;
		    }
		  /*Add the force F_rho the 2(4) link in the path: +     */
		  u_shift_hw_fermion(P7, P7rho, rho);
		  side_link_3f_force(rho,sig,SevenSt,Pnumu,P7,Prhonumu,P7rho);
		  /* Add the P7rho vector to P5 */
		  coeff[0] = SevenSt[0]/FiveSt[0] ; 
		  coeff[1] = SevenSt[1]/FiveSt[1] ; 
		  FORALLSITES(i,s)
		    {
		      scalar_mult_add_su3_vector(&(P5[i].h[0]),
						 &(P7rho[i].h[0]),coeff[0],
						 &(P5[i].h[0]));
		      scalar_mult_add_su3_vector(&(P5[i].h[1]),
						 &(P7rho[i].h[1]),coeff[1],
						 &(P5[i].h[1]));
		    }
		}/* rho */
	      /* Length 5 paths */
	      /*Add the force F_nu the 1(3) link in the path: -     */
	      u_shift_hw_fermion(P5,P5nu, nu);
	      side_link_3f_force(nu,sig,mFiveSt,Pmu,P5, 
			      Pnumu,P5nu) ;
	      /* Add the P5nu vector to P3 */
	      coeff[0] = FiveSt[0]/ThreeSt[0] ; 
	      coeff[1] = FiveSt[1]/ThreeSt[1] ; 
	      FORALLSITES(i,s)
		{
		  scalar_mult_add_su3_vector(&(P3[i].h[0]),
					     &(P5nu[i].h[0]), coeff[0],
					     &(P3[i].h[0]));
		  scalar_mult_add_su3_vector(&(P3[i].h[1]),
					     &(P5nu[i].h[1]), coeff[1],
					     &(P3[i].h[1]));
		}
	    }/* nu */

	  /* Now the Lepage term... It is the same with 5-link paths with
             nu=mu and FiveSt=Lepage. So Pnumu is really Pmumu */
	  u_shift_hw_fermion(Pmu, Pnumu, OPP_DIR(mu));
	  u_shift_hw_fermion(Pnumu, P5, sig);
	  if(GOES_FORWARDS(sig))
	    {
	      /* Add the force F_sig[x+mu+nu]:      x--+             *
	       *                                   |   |             *
	       *                                   o   o             *
	       * the 2 link in the path: + (numbering starts form 0) */
	      add_3f_force_to_mom(P5, Pnumu, sig, Lepage) ;
	    }
	  /*Add the force F_nu the 1(3) link in the path: -     */
	  u_shift_hw_fermion(P5,P5nu, mu);
	  side_link_3f_force(mu, sig, mLepage, Pmu, P5, Pnumu, P5nu) ;
	  /* Add the P5nu vector to P3 */
	  coeff[0] = Lepage[0]/ThreeSt[0] ; 
	  coeff[1] = Lepage[1]/ThreeSt[1] ; 
	  FORALLSITES(i,s)
	    {
	      scalar_mult_add_su3_vector(&(P3[i].h[0]),
					 &(P5nu[i].h[0]),coeff[0],
					 &(P3[i].h[0]));
	      scalar_mult_add_su3_vector(&(P3[i].h[1]),
					 &(P5nu[i].h[1]),coeff[1],
					 &(P3[i].h[1]));
	    }

	  /* Length 3 paths (Not the Naik term) */
	  /*Add the force F_mu the 0(2) link in the path: +     */
	  if(GOES_FORWARDS(mu)) 
	    u_shift_hw_fermion(P3,P3mu, mu );
	  /* The above shift is not needed if mu is backwards */
	  side_link_3f_force(mu, sig, ThreeSt, temp_x, P3, Pmu, P3mu);

	  /* Finally the OneLink and the Naik term */
	  /* Check if this direction is not already done */
	  if( (!DirectLinks[mu]) ){
	    DirectLinks[mu]=1 ;
	    if(GOES_BACKWARDS(mu))/* Do only the forward terms in the Dslash */
	      {
		/* Because I have shifted with OPP_DIR(mu) Pmu is a forward *
		 * shift.                                                   */
		/* The one link */
		add_3f_force_to_mom(Pmu, temp_x, OPP_DIR(mu), OneLink) ;
		/* For the same reason Pnumu is the forward double link */

		/* Popmu is a backward shift */
		u_shift_hw_fermion(temp_x, Popmu, mu);
		/* The Naik */
		/* link no 1: - */
		add_3f_force_to_mom(Pnumu, Popmu, OPP_DIR(mu), mNaik) ;
		/*Pmumumu can overwrite Popmu which is no longer needed */
		u_shift_hw_fermion(Pnumu, Pmumumu, OPP_DIR(mu));
		/* link no 0: + */
		add_3f_force_to_mom(Pmumumu, temp_x, OPP_DIR(mu), Naik);
	      }
	    else /* The rest of the Naik terms */
	      {
		u_shift_hw_fermion(temp_x, Popmu, mu);
		/* link no 2: + */
		/* Pnumu is double backward shift */
		add_3f_force_to_mom(Popmu, Pnumu, mu, Naik) ;
	      }
	  }
	}/* mu */
      /* Here we have to do together the Naik term and the one link term */
    }/*sig */

  /* Free temporary vectors */
  free(temp_x) ;
  for(mu=0;mu<8;mu++)
    free(hwvec[mu]) ;

#ifdef FFTIME
  dtime += dclock();
node0_printf("FFTIME:  %e\n",dtime);
/**printf("TLENGTH: %d\n",tlength);**/
#endif
} /* eo_fermion_force_3f(version 6) */
#undef Pmu          
#undef Pnumu        
#undef Prhonumu     
#undef P7           
#undef P7rho        
#undef P7rhonu      
#undef P5           
#undef P3           
#undef P5nu         
#undef P3mu         
#undef Popmu        
#undef Pmumumu      
#endif /* ASQ_OPTIMIZED_FORCE */


#ifdef  ASQ_OPTIMIZED_FATTENING   /* Asqtad action only, "_fn" executables */
#ifndef FN
BOMB THE COMPILE
#endif
void compute_gen_staple(field_offset staple, int mu, int nu, 
			field_offset link, double coef) {
  su3_matrix tmat1,tmat2;
  msg_tag *mtag0,*mtag1;
  su3_matrix *tempmat ;
  register site *s ;
  register int i ;
  register su3_matrix *fat1;

  /* Computes the staple :
                   mu
                +-------+
          nu	|	|
		|	|
		X	X
    Where the mu link can be any su3_matrix. The result is saved in staple.
    if staple==NULL_FP then the result is not saved.
    It also adds the computed staple to the fatlink[mu] with weight coef.
  */

  /* Upper staple */
  mtag0 = start_gather( link, sizeof(su3_matrix), nu, EVENANDODD, gen_pt[0] );
  mtag1 = start_gather( F_OFFSET(link[nu]), sizeof(su3_matrix), mu, 
			EVENANDODD, gen_pt[1] );
  wait_gather(mtag0);
  wait_gather(mtag1);
  
  if(staple!=NULL_FP){/* Save the staple */
    FORALLSITES(i,s){
      mult_su3_na( (su3_matrix *)gen_pt[0][i],
		   (su3_matrix *)gen_pt[1][i], &tmat1 );
      mult_su3_nn( &(s->link[nu]), &tmat1, (su3_matrix *)F_PT(s,staple) );
    }
  }
  else{ /* No need to save the staple. Add it to the fatlinks */
    FORALLSITES(i,s){
      mult_su3_na( (su3_matrix *)gen_pt[0][i],
		   (su3_matrix *)gen_pt[1][i], &tmat1 );
      mult_su3_nn( &(s->link[nu]), &tmat1, &tmat2 );
#ifdef DSLASH_TMP_LINKS
      fat1 = &(t_fatlink[4*i+mu]);
#else
      fat1 = &(s->fatlink[mu]);
#endif
      scalar_mult_add_su3_matrix(fat1, &tmat2, coef,
				 fat1) ;
    }
  }
  cleanup_gather(mtag0);
  cleanup_gather(mtag1);

  /* lower staple */
  tempmat = (su3_matrix *)calloc(sites_on_node, sizeof(su3_matrix));
  mtag0 = start_gather( F_OFFSET(link[nu]),
			sizeof(su3_matrix), mu, EVENANDODD, gen_pt[0] );
  wait_gather(mtag0);
  FORALLSITES(i,s){
    mult_su3_an( &(s->link[nu]),(su3_matrix *)F_PT(s,link), &tmat1 );
    mult_su3_nn( &(tmat1),(su3_matrix *)gen_pt[0][i], &(tempmat[i]) );
  }
  cleanup_gather(mtag0);
  mtag0 = start_gather_from_temp( tempmat, sizeof(su3_matrix),
				  OPP_DIR(nu), EVENANDODD, gen_pt[0] );
  wait_gather(mtag0);

  if(staple!=NULL_FP){/* Save the staple */
    FORALLSITES(i,s){
      add_su3_matrix( (su3_matrix *)F_PT(s,staple),(su3_matrix *)gen_pt[0][i], 
		      (su3_matrix *)F_PT(s,staple) );
#ifdef DSLASH_TMP_LINKS
      fat1 = &(t_fatlink[4*i+mu]);
#else
      fat1 = &(s->fatlink[mu]);
#endif
      scalar_mult_add_su3_matrix( fat1,
				 (su3_matrix *)F_PT(s,staple), coef, 
				 fat1 );
    }
  }
  else{ /* No need to save the staple. Add it to the fatlinks */
    FORALLSITES(i,s){
#ifdef DSLASH_TMP_LINKS
      fat1 = &(t_fatlink[4*i+mu]);
#else
      fat1 = &(s->fatlink[mu]);
#endif
      scalar_mult_add_su3_matrix( fat1,
				 (su3_matrix *)gen_pt[0][i], coef, 
				 fat1 );
    }
  }

  free(tempmat);
  cleanup_gather(mtag0);
}
#endif  /* ASQ_OPTIMIZED_FATTENING   */


#ifdef  ASQ_OPTIMIZED_FORCE
/*   Covariant shift of the src fermion field in the direction dir  *
 *  by one unit. The result is stored in dest.                       */ 
void u_shift_fermion(su3_vector *src, su3_vector *dest, int dir ) {
  su3_vector *tmpvec ; 
  msg_tag *mtag ;
  register site *s ;
  register int i ;
  
  if(GOES_FORWARDS(dir)) /* forward shift */
    {
      mtag = start_gather_from_temp(src, sizeof(su3_vector), 
				    dir, EVENANDODD, gen_pt[0]);
      wait_gather(mtag);
      FORALLSITES(i,s)
	mult_su3_mat_vec(&(s->link[dir]),(su3_vector *)(gen_pt[0][i]),
			 &(dest[i]));
      cleanup_gather(mtag);
    }
  else /* backward shift */
    {
      tmpvec = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));
      FORALLSITES(i,s)
	mult_adj_su3_mat_vec(&(s->link[OPP_DIR(dir)]),&(src[i]), &tmpvec[i]);
      mtag = start_gather_from_temp(tmpvec, sizeof(su3_vector), dir, 
				    EVENANDODD, gen_pt[0]);
      wait_gather(mtag);
      /* copy the gen_pt to the dest */
      FORALLSITES(i,s)
	dest[i] = *(su3_vector *)gen_pt[0][i];
      cleanup_gather(mtag);
      free(tmpvec) ;
    }
}

/*  Covariant shift of the src half wilson fermion field in the  *
 * direction dir   by one unit. The result is stored in dest.    */
void u_shift_hw_fermion(half_wilson_vector *src, 
			half_wilson_vector *dest, int dir ) {
  half_wilson_vector *tmpvec ; 
  msg_tag *mtag ;
  register site *s ;
  register int i ;
  
  if(GOES_FORWARDS(dir)) /* forward shift */
    {
      mtag = start_gather_from_temp(src, sizeof(half_wilson_vector), 
				    dir, EVENANDODD, gen_pt[0]);
      wait_gather(mtag);
      FORALLSITES(i,s)
	mult_su3_mat_hwvec(&(s->link[dir]),
			   (half_wilson_vector *)(gen_pt[0][i]), &(dest[i]));
      cleanup_gather(mtag);
    }
  else /* backward shift */
    {
      tmpvec = 
	(half_wilson_vector *)calloc(sites_on_node, sizeof(half_wilson_vector));
      FORALLSITES(i,s)
	mult_adj_su3_mat_hwvec(&(s->link[OPP_DIR(dir)]),&(src[i]), &tmpvec[i]);
      mtag = start_gather_from_temp(tmpvec, sizeof(half_wilson_vector), dir, 
				    EVENANDODD, gen_pt[0]);
      wait_gather(mtag);
      /* copy the gen_pt to the dest */
      FORALLSITES(i,s)
	dest[i] = *(half_wilson_vector *)gen_pt[0][i];
      cleanup_gather(mtag);
      free(tmpvec) ;
    }
}

/* Add in contribution to the force */
/* Put antihermitian traceless part into momentum */
void add_force_to_mom(su3_vector *back,su3_vector *forw,int dir,double coeff) {
  register site *s ;
  register int i ;  
  register double tmp_coeff ;

  su3_matrix tmat,tmat2;

  if(GOES_BACKWARDS(dir))
    {
      dir = OPP_DIR(dir) ; 
      coeff = -coeff ;
    }
  FORALLSITES(i,s){
    if(s->parity==ODD) 
      tmp_coeff = -coeff ;
    else
      tmp_coeff = coeff ;
    uncompress_anti_hermitian( &(s->mom[dir]), &tmat2 );
    su3_projector(&(back[i]), &(forw[i]), &tmat);
    scalar_mult_add_su3_matrix(&tmat2, &tmat,  tmp_coeff, &tmat2 );
    make_anti_hermitian( &tmat2, &(s->mom[dir]) ); 
  }
}

/* Add in contribution to the force ( 3flavor case ) */
/* Put antihermitian traceless part into momentum */
void add_3f_force_to_mom(half_wilson_vector *back,
			 half_wilson_vector *forw, int dir, double coeff[2]) {
  register site *s ;
  register int i ;  
  double tmp_coeff[2] ;

  su3_matrix tmat,tmat2;

  if(GOES_BACKWARDS(dir))
    {
      dir = OPP_DIR(dir) ; 
      coeff[0] = -coeff[0] ;
      coeff[1] = -coeff[1] ;
    }
  FORALLSITES(i,s){
    if(s->parity==ODD)
      {
	tmp_coeff[0] = -coeff[0] ;
	tmp_coeff[1] = -coeff[1] ;
      }
    else
      {
	tmp_coeff[0] = coeff[0] ;
	tmp_coeff[1] = coeff[1] ;
      }
    uncompress_anti_hermitian( &(s->mom[dir]), &tmat2 );
    su3_projector(&(back[i].h[0]), &(forw[i].h[0]), &tmat);
    scalar_mult_add_su3_matrix(&tmat2, &tmat,  tmp_coeff[0], &tmat2 );
    su3_projector(&(back[i].h[1]), &(forw[i].h[1]), &tmat);
    scalar_mult_add_su3_matrix(&tmat2, &tmat,  tmp_coeff[1], &tmat2 );
    make_anti_hermitian( &tmat2, &(s->mom[dir]) ); 
  }
}

/*  This routine is needed in order to add the force on the side link *
 * of the paths in the Asq and Asqtad actions. It gets as inputs the  *
 * direction mu of the side link and the direction nu of the Dslash   *
 * term we are dealing with. Then it takes also 4 fermion fields:     *
 * Path: the piece of the path with no hop in the nu or mu direction  *
 * Path_nu: the piece of the path with a hop in nu  but not in mu     *
 * Path_mu: is Path times the link mu                                 *
 * Path_numu: is Path_nu times the link mu                            */
void side_link_force(int mu, int nu, double coeff, 
		     su3_vector *Path   , su3_vector *Path_nu, 
		     su3_vector *Path_mu, su3_vector *Path_numu) {
  if(GOES_FORWARDS(mu))
    {
      /*                    nu           * 
       * Add the force :  +----+         *
       *               mu |    |         *
       *                  x    (x)       *
       *                  o    o         */
      if(GOES_FORWARDS(nu))
	add_force_to_mom(Path_numu, Path, mu, coeff ) ;
      else
	add_force_to_mom(Path, Path_numu, OPP_DIR(mu), -coeff );/* ? extra - */
    }
  else /*GOES_BACKWARDS(mu)*/
    {
      /* Add the force :  o    o         *
       *               mu |    |         *
       *                  x    (x)       *
       *                  +----+         *
       *                    nu           */ 
      if(GOES_FORWARDS(nu))
	add_force_to_mom(Path_nu, Path_mu, mu, -coeff) ; /* ? extra - */
      else
	add_force_to_mom(Path_mu, Path_nu, OPP_DIR(mu), coeff) ;
    }
}
 
/*  The 3 flavor version of side_link_force used *
 * to optimize fermion transports                */
void side_link_3f_force(int mu, int nu, double coeff[2], 
			half_wilson_vector *Path   , 
			half_wilson_vector *Path_nu, 
			half_wilson_vector *Path_mu, 
			half_wilson_vector *Path_numu) {
  double m_coeff[2] ;

  m_coeff[0] = -coeff[0] ;
  m_coeff[1] = -coeff[1] ;

  if(GOES_FORWARDS(mu))
    {
      /*                    nu           * 
       * Add the force :  +----+         *
       *               mu |    |         *
       *                  x    (x)       *
       *                  o    o         */
      if(GOES_FORWARDS(nu))
	add_3f_force_to_mom(Path_numu, Path, mu, coeff ) ;
      else
	add_3f_force_to_mom(Path,Path_numu,OPP_DIR(mu),m_coeff);/* ? extra - */
    }
  else /*GOES_BACKWARDS(mu)*/
    {
      /* Add the force :  o    o         *
       *               mu |    |         *
       *                  x    (x)       *
       *                  +----+         *
       *                    nu           */ 
      if(GOES_FORWARDS(nu))
	add_3f_force_to_mom(Path_nu, Path_mu, mu, m_coeff) ; /* ? extra - */
      else
	add_3f_force_to_mom(Path_mu, Path_nu, OPP_DIR(mu), coeff) ;
    }
}
#endif  /* ASQ_OPTIMIZED_FORCE   */

/* LONG COMMENTS
   Here we have combined "xxx", (offset "x_off")  which is
(M_adjoint M)^{-1} phi, with Dslash times this vector, which goes in the
odd sites of xxx.  Recall that phi is defined only on even sites.  In
computing the fermion force, we are looking at

< X |  d/dt ( Dslash_eo Dslash_oe ) | X >
=
< X | d/dt Dslash_eo | T > + < T | d/dt Dslash_oe | X >
where T = Dslash X.

The subsequent manipulations to get the coefficent of H, the momentum
matrix, in the simulation time derivative above look the same for
the two terms, except for a minus sign at the end, if we simply stick
T, which lives on odd sites, into the odd sites of X

 Each path in the action contributes terms when any link of the path
is the link for which we are computing the force.  We get a minus sign
for odd numbered links in the path, since they connect sites of the
opposite parity from what it would be for an even numbered link.
Minus signs from "going around" plaquette - ie KS phases, are supposed
to be already encoded in the path coefficients.
Minus signs from paths that go backwards are supposed to be already
encoded in the path coefficients.

Here, for example, are comments reproduced from the force routine for
the one-link plus Naik plus single-staple-fat-link action:

 The three link force has three contributions, where the link that
was differentiated is the first, second, or third link in the 3-link
path, respectively.  Diagramatically, where "O" represents the momentum,
the solid line the link corresponding to the momentum, and the dashed
lines the other links:
 

	O______________ x ............ x ...............
+
	x..............O______________x.................
+
	x..............x..............O________________
Think of this as
	< xxx | O | UUUxxx >		(  xxx, UUUX_p3 )
+
	< xxx U | O | UUxxx >		( X_m1U , UUX_p2 )
+
	< xxx U U | O | Uxxx >		( X_m2UU , UX_p1 )
where "U" indicates parallel transport, "X_p3" is xxx displaced
by +3, etc.
Note the second contribution has a relative minus sign
because it effectively contributes to the <odd|even>, or M_adjoint,
part of the force when we work on an even site. i.e., for M on
an even site, this three link path begins on an odd site.

The staple force has six contributions from each plane containing the
link direction:
Call these diagrams A-F:


	x...........x		O____________x
		    .			     .
		    .			     .
		    .			     .
		    .			     .
		    .			     .
	O___________x		x............x
	   (A)			    (B)



	x	    x		O____________x
	.	    .		.	     .
	.	    .		.	     .
	.	    .		.	     .
	.	    .		.	     .
	.	    .		.	     .
	O___________x		x	     x
	   (C)			    (D)



	x...........x		O____________x
	.			.
	.			.
	.			.
	.			.
	.			.
	O___________x		x............x
	   (E)			    (F)

As with the Naik term, diagrams C and D have a relative minus
sign because they connect sites of the other parity.

Also note an overall minus sign in the staple terms relative to the
one link term because, with the KS phase factors included, the fat
link is  "U - w3 * UUU", or the straight link MINUS w3 times the staples.

Finally, diagrams B and E get one more minus sign because the link
we are differentiating is in the opposite direction from the staple
as a whole.  You can think of this as this "U" being a correction to
a "U_adjoint", but the derivative of U is iHU and the derivative
of U_adjoint is -iHU_adjoint.

*/
/* LONG COMMENT on sign conventions
In most of the program, the KS phases and antiperiodic boundary
conditions are absorbed into the link matrices.  This greatly simplfies
multiplying by the fermion matrix.  However, it requires care in
specifying the path coefficients.  Remember that each time you
encircle a plaquette, you pick up a net minus sign from the KS phases.
Thus, when you have more than one path to the same point, you generally
have a relative minus sign for each plaquette in a surface bounded by
this path and the basic path for that displacement.

Examples:
  Fat Link:
    Positive:	X-------X

    Negative     --------
	 	|	|
		|	|
		X	X

  Naik connection, smeared
    Positive:	X-------x-------x-------X

    Negative:	---------
		|	|
		|	|
		X	x-------x-------X

    Positive:	--------x--------
		|		|
		|		|
		X		x-------X

    Negative:	--------x-------x-------x
		|			|
		|			|
		X			X
*/



/* Comment on acceptable actions.
   We construct the backwards part of dslash by reversing all the
   paths in the forwards part.  So, for example, in the p4 action
   the forwards part includes +X+Y+Y

		X
		|
		|
		X
		|
		|
	X---->--X

  so we put -X-Y-Y in the backwards part.  But this isn't the adjoint
  of U_x(0)U_y(+x)U_y(+x+y).  Since much of the code assumes that the
  backwards hop is the adjoint of the forwards (for example, in
  preventing going to 8 flavors), the code only works for actions
  where this is true.  Roughly, this means that the fat link must
  be symmetric about reflection around its midpoint.  Equivalently,
  the paths in the backwards part of Dslash are translations of the
  paths in the forwards part.  In the case of the "P4" or knight's move
  action, this means that we have to have both paths
   +X+Y+Y and +Y+Y+X to the same point, with the same coefficients.
  Alternatively, we could just use the symmetric path +Y+X+Y.
*/

invivoquark_stuff_u_shift_fermion_1452

Real time: 1.434033e+09 -- Predicted time: 1.555542e+09

Callcount: 10160

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
5081	### 1	18.4835594928402	1.745580e+05	1.706760e+05	2.2
747	### 2	1091.68869192986	1.778670e+05	1.457720e+05	18.0
1494	### 3	2060.23463563925	1.782820e+05	1.450760e+05	18.6
112	### 4	1940.04278282046	1.721280e+05	1.445440e+05	16.0
940	### 5	1283.82766135792	0.000000e+00	1.431600e+05	100.0
117	### 6	1528.94390026714	1.738960e+05	1.437440e+05	17.3
485	### 7	1193.99173118011	1.763020e+05	1.441560e+05	18.2
5971	### 8	264.910972320795	1.713880e+05	1.409000e+05	17.7
2244	### 9	131.17634454115	1.659790e+05	1.512040e+05	8.9
344	### 10	647.84990308172	1.749570e+05	1.423880e+05	18.6

                /******************  com_vanilla.c **************************************/
/* Communications routines for the SU3 program
   MIMD version 6.
   This file is communications-scheme dependent.
   Version for single processor machines.
   */
#define NOWHERE -1      /* Not an index in array of fields */

/* Modifications

   1/28/00 combined with Schroedinger functional and
           32 sublattice versions - UMH
   9/2/97  Revised to allow gathers from temporary fields.  neighbor[]
       is now list of indices, add start/restart_gather_from_temp D.T.
   8/05/97 ANSI prototyping for all routines C.D.
   10/5/96 Removed parallel I/O wrappers. Use io_ansi.c now.
   8/30/96 Added restore/save_checkpoint for compatibility
   8/05/96 Added broadcast_bytes and wrappers for system-dependent
    parallel file system calls C.D.
*/

/*
   This is a trivial version, where there is only one node, every site
   is on node ...

   initialize_machine() does any machine dependent setup at the
   very beginning.

   Some routines implement basic functions in a machine independent
   way.
   mynode() returns node number of this node.
   numnodes() returns number of nodes
   g_sync() provides a synchronization point for all nodes.
   g_floatsum() sums a floating point number over all nodes.
   g_doublesum() sums a double over all nodes.
   g_vecdoublesum() sums a vector of doubles over all nodes.
   g_complexsum() sums a single precision complex number over all nodes.
   g_dcomplexsum() sums a double precision complex number over all nodes.
   g_veccomplexsum() sums a vector of single precision complex numbers over all
	nodes.
   g_wvectorsum() sums a single precision wilson vector over all nodes.
   g_xor32() finds global exclusive or of 32-bit word
   g_floatmax() finds maximum of a floating point number over all nodes.
   g_doublemax() finds maximum of a double over all nodes.
   broadcast_float()  broadcasts a single precision number from
	node 0 to all nodes.
   broadcast_double()  broadcasts a double precision number
   broadcast_complex()  broadcasts a single precision complex number
   broadcast_dcomplex()  broadcasts a double precision complex number
   broadcast_bytes()  broadcasts a number of bytes
   send_integer() sends an integer to one other node
   receive_integer() receives an integer
   dclock() returns a double precision time, with arbitrary zero
   terminate() kills the job on all processors

   make_nn_gathers()  makes all necessary lists for communications with
   nodes containing neighbor sites.

   field_pointer_at_coordinates() returns the address of a field in
   the lattice given its coordinates.
   field_pointer_at_direction() returns the address of a field in the
   lattice at a direction from a given site.
   cleanup_field_pointer() frees the buffers that field_pointer...
   allocated.

   start_gather() starts asynchronous sends and receives required
   to gather neighbors.
   start_gather_from_temp() starts asynchronous sends and receives required
   to gather neighbors from temporary array of fields.
   wait_gather()  waits for receives to finish, insuring that the
   data has actually arrived.
   cleanup_gather() frees all the buffers that were allocated, WHICH
   MEANS THAT THE GATHERED DATA MAY SOON DISAPPEAR.
   restart_gather() repeats the internode communications of a previous
   gather.
   restart_gather_from_temp() repeats the internode communications of a
   previous gather of temporary field.

   send_field() sends a field to one other node.
   get_links() receives a field from some other node.
*/

#include "generic_includes.h"
#include <time.h>
#define MAX_GATHERS 32	/* Maximum number of gather tables */

    /* addresses of neighboring sites, NULL if off-node */
    /* neighbor[X][i] is the index of the neighbor of site lattice[i] in
       shuffle number X */
    int ** neighbor;
    /* Number of gathers (mappings) that have been set up */
    int n_gathers;

/**********************************************************************/
/* Machine initialization */
void initialize_machine(int argc, char **argv){}


/**********************************************************************/
/* Set up "comlink" structures needed by gather routines.
   make_lattice() must be called first. */
void make_nn_gathers(){

int i;
void neighbor_coords_special(
 int x,int y,int z,int t, /* coordinates of site */
 int *dirpt,              /* direction (eg XUP) */
 int fb,                  /* "forwards/backwards"  */
 int *x2p,int *y2p,int *z2p,int *t2p);
                          /* pointers to coordinates of neighbor */

    /* Allocate space for lists of pointers to neighbor sites.
       (NULL if neighbor not on this node) */
    neighbor = (int **)calloc(NDIRS, sizeof(int *));
    n_gathers=0;

    for(i=XUP;i<=TUP;i++)
	make_gather(neighbor_coords_special,&i,WANT_INVERSE,
	    ALLOW_EVEN_ODD,SWITCH_PARITY);

    /* Sort into the order we want for nearest neighbor gathers,
	so you can use XUP, XDOWN, etc. as argument in calling them. */
    sort_eight_special((void **) neighbor );
}

/**********************************************************************/
/* sort a list of eight pointers into the order we want for the
  nearest neighbor gathers:  XUP,YUP,ZUP,TUP,TDOWN,ZDOWN,YDOWN,XDOWN,
  starting from the index for the first pointer */

void sort_eight_neighborlists(int index){
  sort_eight_special((void **)&(neighbor[index]) );
}

/**********************************************************************/
/* sort a list of eight pointers into the order we want for the
  nearest neighbor gathers:  XUP,YUP,ZUP,TUP,TDOWN,ZDOWN,YDOWN,XDOWN */
void sort_eight_special(void **pt){
void *tt[8];
register int i;
    for(i=0;i<8;i++)tt[i]=pt[i];
    for(i=XUP;i<=TUP;i++){pt[i]=tt[2*i]; pt[OPP_DIR(i)]=tt[2*i+1];}
}

/**********************************************************************/
/* utility function for finding coordinates of neighbor */
/* This version for use by make_gather for nearest neighbor gathers */
void neighbor_coords_special(
 int x,int y,int z,int t, /* coordinates of site */
 int *dirpt,              /* direction (eg XUP) */
 int fb,                  /* "forwards/backwards"  */
 int *x2p,int *y2p,int *z2p,int *t2p)
                          /* pointers to coordinates of neighbor */
{
int dir;
    dir = (fb==FORWARDS) ? *dirpt : OPP_DIR(*dirpt);
    *x2p = x; *y2p = y; *z2p = z; *t2p = t;
    switch(dir){
	case XUP: *x2p = (x+1)%nx; break;
	case XDOWN: *x2p = (x+nx-1)%nx; break;
	case YUP: *y2p = (y+1)%ny; break;
	case YDOWN: *y2p = (y+ny-1)%ny; break;
	case ZUP: *z2p = (z+1)%nz; break;
	case ZDOWN: *z2p = (z+nz-1)%nz; break;
	case TUP: *t2p = (t+1)%nt; break;
	case TDOWN: *t2p = (t+nt-1)%nt; break;
	default: printf("BOTCH: bad direction\n"); exit(1);
    }
}

/**********************************************************************/
#define RECEIVE 0
#define SEND 1
/* add another gather to the list of tables */
int make_gather(
 void (*func)(int, int, int, int, int *, int, int *, int *, int *, int *),
        		/* function which defines sites to gather from */
 int *args,		/* list of arguments, to be passed to function */
 int inverse,		/* OWN_INVERSE, WANT_INVERSE, or NO_INVERSE */
 int want_even_odd,	/* ALLOW_EVEN_ODD or NO_EVEN_ODD */
 int parity_conserve)	/* {SAME,SWITCH,SCRAMBLE}_PARITY */
			/* SWITCH_PARITY for the 32 sublattice case means,
			   that for all sites in a sublattice, subl, the
			   target sublattice is the same, and different
			   from sublattice subl. */

{
register int i,k;	/* scratch */
register site *s;	/* scratch */
int dir;		/* direction */
int x,y,z,t;		/* coordinates */

#ifdef N_SUBL32
int subl, nsubl, *othersubl;	/* for 32 sublattice case */

    if( want_even_odd==ALLOW_EVEN_ODD && parity_conserve==SCRAMBLE_PARITY ){
	printf("Unsupported combination for make_gather:");
	printf(" ALLOW_EVEN_ODD and SCRAMBLE_PARITY\n");
	terminate(1);
    }
#endif

    /* we will have one or two more gathers */
    if( inverse==WANT_INVERSE ) n_gathers += 2;
    else			n_gathers += 1;
    if(n_gathers > MAX_GATHERS){
	if(this_node==0)printf("Too many gathers! change MAX_GATHERS\n");
	exit(1);
    }
    /* Allocate more space for lists of pointers to neighbor sites.
       (NULL if neighbor not on this node) */
    neighbor = (int **)realloc(neighbor, n_gathers*sizeof(int *));
    if( inverse==WANT_INVERSE) {
        neighbor[n_gathers-2] = (int *)calloc(sites_on_node, sizeof(int));
        if(neighbor[n_gathers-2]==NULL){
            printf("NODE %d: no room for neighbor vector\n",this_node);
            terminate(1);
        }
        neighbor[n_gathers-1] = (int *)calloc(sites_on_node, sizeof(int));
        if(neighbor[n_gathers-1]==NULL){
            printf("NODE %d: no room for neighbor vector\n",this_node);
            terminate(1);
        }
	dir = n_gathers-2;	/* index of gather we are working on */
    }
    else {
        neighbor[n_gathers-1] = (int *)calloc(sites_on_node, sizeof(int));
        if(neighbor[n_gathers-1]==NULL){
            printf("NODE %d: no room for neighbor vector\n",this_node);
            terminate(1);
        }
	dir = n_gathers-1;
    }

#ifndef N_SUBL32
    /* Check to see if mapping has advertised parity and inverse properties */
    /* Also check to see if it returns legal values for coordinates */
    FORALLSITES(i,s){
	/* find coordinates of neighbor who sends us data */
	func( s->x, s->y, s->z, s->t, args, FORWARDS, &x,&y,&z,&t);
	k=(x+y+z+t)%2;

	if( x<0 || y<0 || z<0  || t<0 || x>=nx || y>=ny || z>=nz || t>=nt){
	    printf("DUMMY! Your gather mapping does not stay in lattice\n");
	    printf("It mapped %d %d %d %d to %d %d %d %d\n",
		s->x,s->y,s->z,s->t,x,y,z,t);
	    terminate(1);
	}
	if( ( parity_conserve==SAME_PARITY &&
		((k==0&&s->parity==ODD)||(k==1&&s->parity==EVEN))
	    )
	  ||( parity_conserve==SWITCH_PARITY &&
		((k==0&&s->parity==EVEN)||(k==1&&s->parity==ODD))
#ifdef SCHROED_FUN
	    ) ) if( *args!=TUP && *args!=TDOWN && t!=0 && t!=(nt-1) ) {
#else
	    ) ){
#endif
	    printf("DUMMY! Your gather mapping does not obey claimed parity\n");
	    printf("It mapped %d %d %d %d to %d %d %d %d\n",
		s->x,s->y,s->z,s->t,x,y,z,t);
	    terminate(1);
	}
	if( inverse==OWN_INVERSE ){
	    int x2,y2,z2,t2;
	    func( x, y, z, t, args, FORWARDS, &x2,&y2,&z2,&t2);
	    if( s->x!=x2 || s->y!=y2 || s->z!=z2 || s->t!=t2 ){
		printf(
		    "DUMMY! Your gather mapping is not its own inverse\n");
		printf("It's square mapped %d %d %d %d to %d %d %d %d\n",
		    s->x,s->y,s->z,s->t,x2,y2,z2,t2);
		terminate(1);
	    }
	}
    }

#else	/* N_SUBL32 */

    othersubl = (int *)calloc(N_SUBL32, sizeof(int));
    if(othersubl == NULL){
	printf("NODE %d: no room for othersubl\n",this_node);
	terminate(1);
    }
    for(subl=0;subl<N_SUBL32;subl++) othersubl[subl] = NOWHERE;

    /* Check to see if mapping has advertised parity and inverse properties */
    /* Also check to see if it returns legal values for coordinates */
    FORALLSITES(i,s){
	/* find coordinates of neighbor who sends us data */
	func( s->x, s->y, s->z, s->t, args, FORWARDS, &x,&y,&z,&t);
	subl = s->parity;
	nsubl = (x%2) + 2*(y%2) + 4*(z%2) + 8*(t%2);
	nsubl += 16*((x/2+y/2+z/2+t/2)%2);

	if( want_even_odd==ALLOW_EVEN_ODD ){
	    if( othersubl[subl] == NOWHERE ){
		othersubl[subl] = nsubl;
	    }
	    else if( othersubl[subl] != nsubl ){
		printf("DUMMY! Your gather mixes up sublattices: %d vs %d\n",
		    othersubl[subl], nsubl);
		terminate(1);
	    }
	}

	if( x<0 || y<0 || z<0  || t<0 || x>=nx || y>=ny || z>=nz || t>=nt){
	    printf("DUMMY! Your gather mapping does not stay in lattice\n");
	    printf("It mapped %d %d %d %d to %d %d %d %d\n",
		s->x,s->y,s->z,s->t,x,y,z,t);
	    terminate(1);
	}
	if( parity_conserve==SAME_PARITY && subl!=nsubl ){
	    printf("DUMMY! Your gather mapping does not obey claimed parity");
	    printf(", namely SAME_PARITY\n");
	    printf("It mapped %d %d %d %d with %d to %d %d %d %d with %d\n",
		s->x,s->y,s->z,s->t,subl,x,y,z,t,nsubl);
	    terminate(1);
	}
	if( parity_conserve==SWITCH_PARITY && subl==nsubl ){
	    printf("DUMMY! Your gather mapping does not obey claimed parity");
	    printf(", namely SWITCH_PARITY\n");
	    printf("It mapped %d %d %d %d with %d to %d %d %d %d with %d\n",
		s->x,s->y,s->z,s->t,subl,x,y,z,t,nsubl);
	    terminate(1);
	}
	if( inverse==OWN_INVERSE ){
	    int x2,y2,z2,t2;
            func( x, y, z, t, args, FORWARDS, &x2,&y2,&z2,&t2);
	    if( s->x!=x2 || s->y!=y2 || s->z!=z2 || s->t!=t2 ){
	        printf(
		    "DUMMY! Your gather mapping is not its own inverse\n");
	        printf("It's square mapped %d %d %d %d to %d %d %d %d\n",
		    s->x,s->y,s->z,s->t,x2,y2,z2,t2);
	        terminate(1);
	    }
	}
    }

    free(othersubl);
#endif	/* N_SUBL32 */

    /* RECEIVE LISTS: */
    /* Fill in pointers to sites which are on this node, NOWHERE if
	they are off-node */
    FORALLSITES(i,s){
        /* find coordinates of neighbor who sends us data */
	func( s->x, s->y, s->z, s->t, args, FORWARDS, &x,&y,&z,&t);
	/* neighbor is always on node in vanilla version */
	neighbor[dir][i]= node_index(x,y,z,t);
    }

    if( inverse != WANT_INVERSE) return(dir);

    /* INVERSE GATHER */
    /* Now, if necessary, make inverse gather */
    /* In most cases, we can use the same lists as the gather, in one
	form or another.  Of course, by the time you get to here
	you know that inverse = WANT_INVERSE */
    dir++;	/* inverse gather has direction one more than original */

    /* Always set up pointers to sites on this node */
    /* scan sites in lattice */
    FORALLSITES(i,s){
	/* find coordinates of neighbor who sends us data */
	func( s->x, s->y, s->z, s->t, args, BACKWARDS, &x,&y,&z,&t);

	/* set up pointer */
	neighbor[dir][i]= node_index(x,y,z,t);
    }

    return(dir-1);
}

/**********************************************************************/
/* utility function for finding coordinates of neighbor */
void neighbor_coords(
 int x, int y, int z, int t,  /* coordinates of site */
 int dir,	              /* direction (eg XUP) */
 int *x2p, int *y2p, int *z2p, int *t2p)
                             /* pointers to coordinates of neighbor */
{
    *x2p = x; *y2p = y; *z2p = z; *t2p = t;
    switch(dir){
	case XUP: *x2p = (x+1)%nx; break;
	case XDOWN: *x2p = (x+nx-1)%nx; break;
	case YUP: *y2p = (y+1)%ny; break;
	case YDOWN: *y2p = (y+ny-1)%ny; break;
	case ZUP: *z2p = (z+1)%nz; break;
	case ZDOWN: *z2p = (z+nz-1)%nz; break;
	case TUP: *t2p = (t+1)%nt; break;
	case TDOWN: *t2p = (t+nt-1)%nt; break;
	default: printf("BOTCH: bad direction\n"); exit(1);
    }
}

/**********************************************************************/
/* Set up interrupt handlers to handle field_pointer routines
   make_lattice() must be called first. */
void start_handlers(){
  void fillfieldrequest(int type,int size,int node,int pid);
    /* nothing to do */
}

/**********************************************************************/
/* GATHER ROUTINES */
/* start_gather() returns a pointer to a list of msg_tag's, which will
   be used as input to subsequent wait_gather() and cleanup_gather() calls.
   In the single processor version, the routine will always return NULL.

   usage:  tag = start_gather( source, size, direction, parity, dest )
   example:
	msg_tag *tag;
	tag = start_gather( F_OFFSET(phi), sizeof(su3_vector), XUP,
	    EVEN, gen_pt[0] );
	  ** do other stuff **
	wait_gather(tag);
	  ** gen_pt[0][i] now contains the address of the phi
	   vector (or a copy therof) on the neighbor of site i in the
	   XUP direction for all even sites i.
	   Do whatever you want with it here.
	  **
	cleanup_gather(tag);
	  ** subsequent calls will overwrite the gathered fields. but if you
	   don't clean up, you will eventually run out of space **

    Under certain circumstances it is possible to efficiently gather
    a field that has been previously gathered.  This happens when the
    field being gathered has been modified, but the pointers (the
    destination of start_gather() ) have not been modified.  To use
    restart_gather, the original gather must have been waited for but
    not cleaned up.  The usage is:
	msg_tag *tag;
	tag = start_gather( F_OFFSET(phi), sizeof(su3_vector), XUP,
	    EVEN, gen_pt[0] );
	  ** do other stuff, but don't modify tag or gen_pt[0] **
	wait_gather(tag);
	  ** gen_pt[0][i] now contains the address of the phi
	   vector (or a copy thereof) on the neighbor of site i in the
	   XUP direction for all even sites i.
	   Do whatever you want with it here, but don't modify tag or
	   gen_pt[0].
	   Do modify the source field phi. **
	restart_gather( F_OFFSET(phi), sizeof(su3_vector), XUP,
	    EVEN, gen_pt[0], tag );
	  ** do other stuff **
	wait_gather(tag);
	  ** gen_pt[0][i] now contains the address of the modified phi.
	   The restart-wait may be repeated as often as desired.  **
	cleanup_gather(tag);
	  ** subsequent calls will overwrite the gathered fields. but if you
	   don't clean up, you will eventually run out of space **

   Internal convention for message types is (schematically):
	( f_offset(src) | sending_node | dir ) + GATHER_BASE
*/

/**********************************************************************/
#ifndef N_SUBL32

msg_tag * start_gather(
/* arguments */
 field_offset field,	/* which field? Some member of structure "site" */
 int size,		/* size in bytes of the field (eg sizeof(su3_vector))*/
 int index,		/* direction to gather from. eg XUP - index into
			   neighbor tables */
 int parity,		/* parity of sites whose neighbors we gather.
			   one of EVEN, ODD or EVENANDODD. */
 char ** dest)		/* one of the vectors of pointers */
{
/* local variables */
register int j;		/* scratch */
register site *s;	/* scratch pointer to site */

    /* set pointers in sites whose neighbors are on this node.  (If all
	neighbors are on this node, this is the only thing done.) */
    switch(parity){
	case EVEN:
	    FOREVENSITES(j,s){ 
                dest[j] = ((char *)(lattice+neighbor[index][j]))+field;
	    }
	    break;
	case ODD:
	    FORODDSITES(j,s){ 
                dest[j] = ((char *)(lattice+neighbor[index][j]))+field;
	    }
	    break;
	case EVENANDODD:
	    FORALLSITES(j,s){ 
                dest[j] = ((char *)(lattice+neighbor[index][j]))+field;
	    }
	    break;
    }

    /* return */
    return(NULL);
}

#else	/* N_SUBL32 */

msg_tag * start_gather(
/* arguments */
 field_offset field,	/* which field? Some member of structure "site" */
 int size,		/* size in bytes of the field (eg sizeof(su3_vector))*/
 int index,		/* direction to gather from. eg XUP - index into
			   neighbor tables */
 int subl,		/* subl of sites whose neighbors we gather.
			   It is EVENANDODD, if all sublattices are done. */
 char ** dest)		/* one of the vectors of pointers */
{
/* local variables */
register int j;		/* scratch */
register site *s;	/* scratch pointer to site */

    /* set pointers in sites whose neighbors are on this node.  (If all
	neighbors are on this node, this is the only thing done.) */
    if( subl == EVENANDODD ){
	FORALLSITES(j,s){
	    dest[j] = ((char *)(lattice+neighbor[index][j]))+field;
	}
    }
    else{
	FORSOMESUBLATTICE(j,s,subl){
	    dest[j] = ((char *)(lattice+neighbor[index][j]))+field;
	}
    }

    /* return */
    return(NULL);
}
#endif	/* N_SUBL32 */

/**********************************************************************/
/* Repeat a gather with the same source and destination as a
  previous gather.  The previous gather must have been waited for
  but not cleaned up.  Pointers to sites on the same node are not
  reset, and the same buffers are reused. */
void restart_gather(
/* arguments */
 field_offset field,	/* which field? Some member of structure "site" */
 int size,		/* size in bytes of the field (eg sizeof(su3_vector))*/
 int index,		/* direction to gather from. eg XUP - index into
			   neighbor tables */
 int parity,		/* parity of sites whose neighbors we gather.
			   one of EVEN, ODD or EVENANDODD. */
 char ** dest,		/* one of the vectors of pointers */
 msg_tag *mbuf)          /* previously returned by start_gather */
{
    /* Nothing to do */
}
/**********************************************************************/
#ifndef N_SUBL32

msg_tag * start_gather_from_temp(
/* arguments */
 void * field,		/* which field? Pointer returned by malloc() */
 int size,		/* size in bytes of the field (eg sizeof(su3_vector))*/
 int index,		/* direction to gather from. eg XUP - index into
			   neighbor tables */
 int parity,		/* parity of sites whose neighbors we gather.
			   one of EVEN, ODD or EVENANDODD. */
 char ** dest)		/* one of the vectors of pointers */
{
/* local variables */
register int j;		/* scratch */
register site *s;	/* scratch pointer to site */

    /* set pointers in sites whose neighbors are on this node.  (If all
	neighbors are on this node, this is the only thing done.) */
    switch(parity){
	case EVEN:
	    FOREVENSITES(j,s){ 
                dest[j] = (char *)field + neighbor[index][j]*size;
	    }
	    break;
	case ODD:
	    FORODDSITES(j,s){ 
                dest[j] = (char *)field + neighbor[index][j]*size;
	    }
	    break;
	case EVENANDODD:
	    FORALLSITES(j,s){ 
                dest[j] = (char *)field + neighbor[index][j]*size;
	    }
	    break;
    }

    /* return */
    return(NULL);
}

#else	/* N_SUBL32 */

msg_tag * start_gather_from_temp(
/* arguments */
 void * field,		/* which field? Pointer returned by malloc() */
 int size,		/* size in bytes of the field (eg sizeof(su3_vector))*/
 int index,		/* direction to gather from. eg XUP - index into
			   neighbor tables */
 int subl,		/* subl of sites whose neighbors we gather.
			   It is EVENANDODD, if all sublattices are done. */
 char ** dest)		/* one of the vectors of pointers */
{
/* local variables */
register int j;		/* scratch */
register site *s;	/* scratch pointer to site */

    /* set pointers in sites whose neighbors are on this node.  (If all
	neighbors are on this node, this is the only thing done.) */
    if( subl == EVENANDODD ){
	FORALLSITES(j,s){
	    dest[j] = (char *)field + neighbor[index][j]*size;
	}
    }
    else{
	FORSOMESUBLATTICE(j,s,subl){
	    dest[j] = (char *)field + neighbor[index][j]*size;
	}
    }

    /* return */
    return(NULL);
}
#endif	/* N_SUBL32 */

/**********************************************************************/
/* Repeat a gather with the same source and destination as a
  previous gather.  The previous gather must have been waited for
  but not cleaned up.  Pointers to sites on the same node are not
  reset, and the same buffers are reused. */
void restart_gather_from_temp(
/* arguments */
 void * field,		/* which field? Pointer returned by malloc() */
 int size,		/* size in bytes of the field (eg sizeof(su3_vector))*/
 int index,		/* direction to gather from. eg XUP - index into
			   neighbor tables */
 int parity,		/* parity of sites whose neighbors we gather.
			   one of EVEN, ODD or EVENANDODD. */
 char ** dest,		/* one of the vectors of pointers */
 msg_tag *mbuf)          /* previously returned by start_gather */
{
    /* Nothing to do */
}
/**********************************************************************/
void wait_gather(msg_tag *mbuf) {
    return;
}

/**********************************************************************/
void cleanup_gather(msg_tag *mbuf) {
    return;
}

/**********************************************************************/
/* GENERAL_GATHER ROUTINES */
/* start_general_gather() returns a pointer to a list of msg_tag's, which will
   be used as input to subsequent wait_general_gather() and
   cleanup_general_gather() calls.
   This list contains msg_tags for all receive buffers, followed by
   a msg_tag with msg_id = 0 and msg_buf = NULL, followed by msg_tags
   for all send buffers, followed by a msg_tag with id=0, buf=NULL.
   If no messages at all are required, the routine will return NULL.
   msg_buf=NULL should be a reliable indicator of no message.

   usage:  tag = start_general_gather( source, size, displacement, parity, dest)
   example:
	msg_tag *tag;
	int disp[4]; 
        disp[XUP]=1; disp[YUP]= -1; disp[ZUP] = disp[TUP] = 0;
	tag = start_general_gather( F_OFFSET(phi), sizeof(su3_vector), disp,
	    EVEN, gen_pt[0] );
	  ** do other stuff **
	wait_general_gather(tag);
	  ** gen_pt[0][i] now contains the address of the phi
	   vector (or a copy thereof) on the neighbor of site i in the
	   XUP direction for all even sites i.
	   Do whatever you want with it here.
	  **
	cleanup_general_gather(tag);
	  ** subsequent calls will overwrite the gathered fields. but if you
	   don't clean up, you will eventually run out of space **

*/
struct msg_tmp { int node,count; }; /* temporary structure for keeping
	track of messages to be sent or received */
int g_gather_flag=0;	/* flag to tell if general gather in progress */

#ifndef N_SUBL32

msg_tag * start_general_gather(
/* arguments */
 field_offset field,	/* which field? Some member of structure "site" */
 int size,		/* size in bytes of the field (eg sizeof(su3_vector))*/
 int *displacement,	/* displacement to gather from. four components */
 int parity,		/* parity of sites to which we gather.
			   one of EVEN, ODD or EVENANDODD. */
 char ** dest)		/* one of the vectors of pointers */
{
/* local variables */
register int i;		/* scratch */
register site *s;	/* scratch pointer to site */
int tx,ty,tz,tt;	/* temporary coordinates */


    /* check for gather already in progress */
    if(g_gather_flag!=0){
	fprintf(stderr,"ERROR: node %d, two general_gathers() at once!\n",
	    mynode() );
	exit(1);
    }

    /* set pointers in sites whose neighbors are on this node.  (If all
	neighbors are on this node, this is the only thing done.) Make
	list of nodes from whom we expect messages */
    FORSOMEPARITY(i,s,parity){
	if(displacement[XUP]!=0)tx = (s->x + displacement[XUP] + nx)%nx;
			else    tx = s->x;
	if(displacement[YUP]!=0)ty = (s->y + displacement[YUP] + ny)%ny;
			else    ty = s->y;
	if(displacement[ZUP]!=0)tz = (s->z + displacement[ZUP] + nz)%nz;
			else    tz = s->z;
	if(displacement[TUP]!=0)tt = (s->t + displacement[TUP] + nt)%nt;
			else    tt = s->t;
	dest[i] = F_PT( &lattice[node_index(tx,ty,tz,tt)], field );
    }

    /* mark gather in progress and return */
    g_gather_flag=1;

    return(NULL);
}

#else	/* N_SUBL32 */

msg_tag * start_general_gather(
/* arguments */
 field_offset field,	/* which field? Some member of structure "site" */
 int size,		/* size in bytes of the field (eg sizeof(su3_vector))*/
 int *displacement,	/* displacement to gather from. four components */
 int subl,		/* subl of sites whose neighbors we gather.
			   It is EVENANDODD, if all sublattices are done. */
 char ** dest)		/* one of the vectors of pointers */
{
/* local variables */
register int i;		/* scratch */
register site *s;	/* scratch pointer to site */
int tx,ty,tz,tt;	/* temporary coordinates */


    /* check for gather already in progress */
    if(g_gather_flag!=0){
	fprintf(stderr,"ERROR: node %d, two general_gathers() at once!\n",
	    mynode() );
	exit(1);
    }

    /* set pointers in sites whose neighbors are on this node.  (If all
	neighbors are on this node, this is the only thing done.) Make
	list of nodes from whom we expect messages */
    if( subl == EVENANDODD ){
	FORALLSITES(i,s){
	    if(displacement[XUP]!=0)tx = (s->x + displacement[XUP] + nx)%nx;
			else    tx = s->x;
	    if(displacement[YUP]!=0)ty = (s->y + displacement[YUP] + ny)%ny;
			else    ty = s->y;
	    if(displacement[ZUP]!=0)tz = (s->z + displacement[ZUP] + nz)%nz;
			else    tz = s->z;
	    if(displacement[TUP]!=0)tt = (s->t + displacement[TUP] + nt)%nt;
			else    tt = s->t;
	    dest[i] = (char *)( &lattice[node_index(tx,ty,tz,tt)]) + field;
	}
    }
    else{
	FORSOMESUBLATTICE(i,s,subl){
	    if(displacement[XUP]!=0)tx = (s->x + displacement[XUP] + nx)%nx;
			else    tx = s->x;
	    if(displacement[YUP]!=0)ty = (s->y + displacement[YUP] + ny)%ny;
			else    ty = s->y;
	    if(displacement[ZUP]!=0)tz = (s->z + displacement[ZUP] + nz)%nz;
			else    tz = s->z;
	    if(displacement[TUP]!=0)tt = (s->t + displacement[TUP] + nt)%nt;
			else    tt = s->t;
	    dest[i] = (char *)( &lattice[node_index(tx,ty,tz,tt)]) + field;
	}
    }

    /* mark gather in progress and return */
    g_gather_flag=1;

    return(NULL);
}
#endif	/* N_SUBL32 */

/**********************************************************************/
void wait_general_gather(msg_tag *mbuf) {
    g_gather_flag=0;
    return;
}

/**********************************************************************/
void cleanup_general_gather(msg_tag *mbuf) {
    return;
}


/**********************************************************************/
/* FIELD POINTER ROUTINES */

/* Return a pointer to a field in the lattice at some coordinates.
   If the site is on this node, just return its address.  If it is
   on another node, make a buffer, get the data, and return address of
   buffer.
   Example:  
	su3_matrix *pt;
	pt = (su3_matrix *)field_pointer_at_coordinates(
	    F_OFFSET(xlink), sizeof(su3_matrix), x,y,z,t );
	... do stuff here ...
	cleanup_field_pointer( (char *)pt );
*/
char * field_pointer_at_coordinates(
/* arguments */
 int field,	/* offset of one of the fields in lattice[] */
 int size,	/* size of the field in bytes */
 int x,int y,int z,int t)	/* coordinates of point to get field from */
{
int index;
    index=node_index(x,y,z,t);
    /* return address */
    return( F_PT( &(lattice[index]), field ) );
}

/**********************************************************************/
/* Return a pointer to a field in the lattice at some direction from
   a site on this node. (Usually the "current" site.)
   Works like field_pointer_at_coordinates. */
char * field_pointer_at_direction(
/* arguments */
 field_offset field,	/* offset of one of the fields in lattice[] */
 int size,	/* size of the field in bytes */
 site *s,	/* pointer to a site on this node */
 int direction)	/* direction of site's neighbor to get data from.
		   one of XUP,XDOWN,YUP... */
{
    return( F_PT( neighbor[direction][s-lattice], field ) );
}

/**********************************************************************/
/* free any buffers allocated by above routines.  Argument is the
address returned by field_pointer...() */
void cleanup_field_pointer(char * buf) {
}

/**********************************************************************/
/* SEND AND RECEIVE FIELD */
void send_field(char *buf, int size, int tonode) {
    printf("BOTCH: this never happens\n");
}
/**********************************************************************/
void get_field(char *buf, int size){
    printf("BOTCH: this never happens\n");
}

/**********************************************************************/
/* BASIC COMMUNICATIONS FUNCTIONS */

/* Tell what kind of machine we are on */
static char name[]="Scalar processor";
char * machine_type(){
    return(name);
}

/* Return my node number */
int mynode(){ return(0); }

/* Return number of nodes */
int numnodes(){ return(1); }

/* Synchronize all nodes */
void g_sync(){
}

/* Sum double over all nodes */
void g_floatsum( float *fpt) {
}

/* Sum double over all nodes */
void g_doublesum( double *dpt) {
}

/* Sum a vector of doubles over all nodes */
void g_vecdoublesum( double *dpt, int ndoubles) {
}

/* Sum complex over all nodes */
void g_complexsum( complex *cpt) {
}

/* Sum double_complex over all nodes */
void g_dcomplexsum( double_complex *cpt) {
}

/* Sum a vector of complex over all nodes */
void g_veccomplexsum( complex *cpt, int ncomplex){
}

/* Sum wilson_vector over all nodes */
void g_wvectorsum( wilson_vector *wvpt) {
}

/* 32 bit exclusive or over all nodes */
void g_xor32( u_int32type *pt ) {
}

/* Find maximum of double over all nodes */
void g_floatmax( float *fpt) {
}

/* Find maximum of double over all nodes */
void g_doublemax( double *dpt) {
}

/* Broadcast floating point number from node zero */
void broadcast_float(float *fpt) {
}

/* Broadcast double precision doubleing point number from node zero */
void broadcast_double(double *dpt) {
}

/* Broadcast single precision complex number from node zero */
void broadcast_complex(complex *cpt) {
}

/* Broadcast double precision complex number from node zero */
void broadcast_dcomplex(double_complex *cpt) {
}

/* Broadcast bytes from node zero */
void broadcast_bytes(char *buf,int size)  {
}

/* Send an integer to one other node */
/* This is to be called only by the node doing the sending */
void send_integer(int tonode, int *address) {
    printf("BOTCH: this never happens\n");
}

/* Receive an integer from another node */
/* Note we do not check if this was really meant for us */
void receive_integer(int *address) {
    printf("BOTCH: this never happens\n");
}

/* Double precision time */
/* This one wraps around after 36 minutes!! It gives the cpu time,
   not the wall clock time */
double dclock(){
long fine;
  /*CPU2006    fine = clock(); */
    fine = 0.;
    return( ((double)fine)/1000000.0 );
}

/* Print time stamp */
void time_stamp(char *msg){
  time_t time_stamp;
  
  if(mynode()==0){
    time(&time_stamp);
    printf("%s: %s\n",msg,ctime(&time_stamp));
  }
}

/* version of exit for multinode processes -- kill all nodes */
void terminate(int status) {
  time_stamp("termination");
  printf("Termination: status = %d\n",status);
  fflush(stdout);fflush(stderr);
  exit(status);
}

/* version of normal exit for scalar processes */
void normal_exit(int status) {
  /*CPU2006  time_stamp("exit");*/
  exit(status);
}

invivocom_vanilla_start_gather_from_temp_581

Real time: 1.784641e+08 -- Predicted time: 2.017281e+08

Callcount: 26472

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
76	### 1	24477.6146953405	7.600000e+03	6.696000e+03	11.8
65	### 2	2022.97321937322	7.760000e+03	7.020000e+03	9.5

                /******************  m_mat_an.c  (in su3.a) *****************************
*									*
* void mult_su3_an( su3_matrix *a,*b,*c )				*
* matrix multiply, first matrix is adjoint 				*
* C <-  A_adjoint*B							*
*/
#include "config.h"
#include "complex.h"
#include "su3.h"

#ifndef FAST
void mult_su3_an( su3_matrix *a, su3_matrix *b, su3_matrix *c ){
register int i,j,k;
register complex x,y;
    for(i=0;i<3;i++)for(j=0;j<3;j++){
	x.real=x.imag=0.0;
	for(k=0;k<3;k++){
	    CMULJ_( a->e[k][i] , b->e[k][j], y );
	    CSUM( x , y );
	}
	c->e[i][j] = x;
    }
}

/* "Hand coded" routines, clearer coding is up above */
#else

void mult_su3_an( su3_matrix *a, su3_matrix *b, su3_matrix *c ){
  int j;

#ifdef NATIVEDOUBLE
  register double a0r,a0i,a1r,a1i,a2r,a2i;
  register double b0r,b0i,b1r,b1i,b2r,b2i;
#else
  register double a0r,a0i,a1r,a1i,a2r,a2i;
  register double b0r,b0i,b1r,b1i,b2r,b2i;
#endif

  for(j=0;j<3;j++){
    
    a0r=a->e[0][0].real; a0i=a->e[0][0].imag;
    b0r=b->e[0][j].real; b0i=b->e[0][j].imag;
    a1r=a->e[1][0].real; a1i=a->e[1][0].imag;
    b1r=b->e[1][j].real; b1i=b->e[1][j].imag;
    a2r=a->e[2][0].real; a2i=a->e[2][0].imag;
    b2r=b->e[2][j].real; b2i=b->e[2][j].imag;
    
    c->e[0][j].real = a0r*b0r + a0i*b0i + a1r*b1r + a1i*b1i + a2r*b2r + a2i*b2i;
    c->e[0][j].imag = a0r*b0i - a0i*b0r + a1r*b1i - a1i*b1r + a2r*b2i - a2i*b2r;
  
    a0r=a->e[0][1].real; a0i=a->e[0][1].imag;
    b0r=b->e[0][j].real; b0i=b->e[0][j].imag;
    a1r=a->e[1][1].real; a1i=a->e[1][1].imag;
    b1r=b->e[1][j].real; b1i=b->e[1][j].imag;
    a2r=a->e[2][1].real; a2i=a->e[2][1].imag;
    b2r=b->e[2][j].real; b2i=b->e[2][j].imag;
    
    c->e[1][j].real = a0r*b0r + a0i*b0i + a1r*b1r + a1i*b1i + a2r*b2r + a2i*b2i;
    c->e[1][j].imag = a0r*b0i - a0i*b0r + a1r*b1i - a1i*b1r + a2r*b2i - a2i*b2r;
  
    a0r=a->e[0][2].real; a0i=a->e[0][2].imag;
    b0r=b->e[0][j].real; b0i=b->e[0][j].imag;
    a1r=a->e[1][2].real; a1i=a->e[1][2].imag;
    b1r=b->e[1][j].real; b1i=b->e[1][j].imag;
    a2r=a->e[2][2].real; a2i=a->e[2][2].imag;
    b2r=b->e[2][j].real; b2i=b->e[2][j].imag;
    
    c->e[2][j].real = a0r*b0r + a0i*b0i + a1r*b1r + a1i*b1i + a2r*b2r + a2i*b2i;
    c->e[2][j].imag = a0r*b0i - a0i*b0r + a1r*b1i - a1i*b1r + a2r*b2i - a2i*b2r;
  
    }
}

#endif	/* End of "#ifdef FAST" */

invivom_mat_an_mult_su3_an_39

Real time: 5.050794e+08 -- Predicted time: 2.255115e+09

Callcount: 3216000

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
632253	### 1	919253.191630562	5.300000e+02	1.520000e+02	71.3
2304026	### 2	815546.539081672	5.390000e+02	1.560000e+02	71.0
3084728	### 3	485037.946393254	8.370000e+02	1.480000e+02	82.3
2337338	### 4	451413.801902248	1.032000e+03	1.600000e+02	84.4
2637716	### 5	203548.460903411	9.600000e+02	1.640000e+02	82.9
3157735	### 6	162916.985807615	6.000000e+02	1.680000e+02	72.0
1859244	### 7	172645.810865939	8.780000e+02	1.440000e+02	83.5
979226	### 8	11449.7997346746	8.860000e+02	3.240000e+02	63.4
2708617	### 9	2806.6130459105	5.740000e+02	2.240000e+02	60.9

                /************************* control.c *******************************/
/* MIMD version 6 */
/* Main procedure for SU3 with dynamical fermions 			*/
/* Naik plus fat link fermions, general gauge action */

/* This version combines code for the PHI algorithm (approriate for 4
   flavors) and the R algorithm for "epsilon squared" updating of 
   1 to 4 flavors.  Compilation should occur with PHI_ALGORITHM defined
   for the former and not defined for the latter.  It also contains code
   for the hybrid Monte Carlo algorithm, for which HMC_ALGORITHM and
   PHI_ALGORITHM should be defined.  (Actually, the
   changes to control.c are minimal and the real differences will appear
   in update.c */

#define CONTROL
#include "ks_imp_includes.h"	/* definitions files and prototypes */

EXTERN  gauge_header start_lat_hdr;     /* Input gauge field header */

int main( int argc, char **argv ){
    int meascount,traj_done;
    int prompt;
    int s_iters,avs_iters,avspect_iters,
        avbcorr_iters;
    double dtime, dclock();

    initialize_machine(argc,argv);
    g_sync();
    /* set up */
    prompt = setup();
    /* loop over input sets */
    while( readin(prompt) == 0){

	/* perform warmup trajectories */
	dtime = -dclock();
	for( traj_done=0; traj_done < warms; traj_done++ ){
            update();
	}
	node0_printf("WARMUPS COMPLETED\n"); fflush(stdout);

	/* perform measuring trajectories, reunitarizing and measuring 	*/
	meascount=0;		/* number of measurements 		*/
	avspect_iters = avs_iters = avbcorr_iters = 0;
	for( traj_done=0; traj_done < trajecs; traj_done++ ){ 

	    /* do the trajectories */
	    s_iters=update();

	    /* measure every "propinterval" trajectories */
	    if( (traj_done%propinterval)==(propinterval-1) ){

	        /* call gauge_variable fermion_variable measuring routines */
		/* results are printed in output file */
	        rephase(OFF);
		g_measure( );
		rephase(ON);
		f_meas_imp(F_OFFSET(phi),F_OFFSET(xxx),mass);
#ifdef SPECTRUM 
                rephase( OFF );
                gaugefix(TUP,(double)1.8,500,(double)GAUGE_FIX_TOL,
          	   F_OFFSET(tempmat1),F_OFFSET(tempvec[0]),
			 0,NULL,NULL,0,NULL,NULL);
                rephase( ON );
#ifdef FN
		valid_fatlinks = valid_longlinks = 0;
#endif
		avspect_iters += spectrum2(mass,F_OFFSET(phi),F_OFFSET(xxx));

		/* Fix TUP Coulomb gauge - gauge links only*/
		rephase( OFF );
		gaugefix(TUP,(double)1.8,500,(double)GAUGE_FIX_TOL,
			 F_OFFSET(tempmat1),F_OFFSET(staple),
			 0,NULL,NULL,0,NULL,NULL);
		rephase( ON );
		valid_fatlinks = valid_longlinks = 0;

		/* Comment these spectrum calls out as you wish */
		/* There are no side-effects */
                avspect_iters += nl_spectrum(mass,F_OFFSET(phi),F_OFFSET(xxx),
					F_OFFSET(tempmat1),F_OFFSET(staple));
		avspect_iters += spectrum_mom(mass,mass,F_OFFSET(phi),5e-3);
		avspect_iters += spectrum_multimom(mass,mass,mass/4.,4,5e-3);
		avspect_iters += spectrum_nlpi2(mass,mass,F_OFFSET(phi),5e-3);
#ifdef FN
		if( !valid_fatlinks )load_fatlinks();
#endif
		rephase( OFF );
		hvy_pot( F_OFFSET(link[XUP]) );
		rephase( ON );
#endif
		avs_iters += s_iters;
	        ++meascount;
		fflush(stdout);
	    }
	}	/* end loop over trajectories */

	node0_printf("RUNNING COMPLETED\n"); fflush(stdout);
	if(meascount>0)  {
	    node0_printf("average cg iters for step= %e\n",
		(double)avs_iters/meascount);
#ifdef SPECTRUM
	    node0_printf("average cg iters for spectrum = %e\n",
		(double)avspect_iters/meascount);
#endif
	}

	dtime += dclock();
	if(this_node==0){
	  /*CPU2006  printf("Time = %e seconds\n",dtime);*/
	    printf("total_iters = %d\n",total_iters);
	}
	fflush(stdout);

	/* save lattice if requested */
        if( saveflag != FORGET ){
          rephase( OFF );
          save_lattice( saveflag, savefile );
          rephase( ON );
        }
    }
    return 0;
}

invivocontrol_main_32

Real time: 1.329345e+10 -- Predicted time: 2.289945e+07

Callcount: 1

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	2.289945e+07	1.330103e+10	99.8

                /****** quark_stuff.c  -- ******************/
/* MIMD version 6 */
/* quark action stuff for improved action
* D.T. 1/28/98, starting from gauge_stuff.c
* K.O. 3/99 Added optimized fattening for Asq actions
* D.T. 4/99 Combine force calculations for both mass quarks
* K.O. 4/99 Optimized force for Asq action
* S.G. 7/01, modified to use t_longlink and t_fatlink
* C.D. 10/02, consolidated quark_stuff.c and quark_stuff_tmp.c
*
* This code combines quark_stuff.c and quark_stuff_tmp.c
* with DSLASH_TMP_LINKS defined, puts links in field-major variables
* t_longlink and t_fatlink. Otherwise, puts them in the site structure.
*
* In this directory, assume all paths connect even to odd sites, etc.
* Tabulate "backwards" paths (e.g. "XDOWN" is backward path to "XUP")
* as separate parity transforms of the fundamental paths.  They will
* generally need a negative sign in Dslash.  See bottom for a long
* comment on sign conventions.
*/

/**#define FFTIME**/
/**#define LLTIME**/

#include "generic_ks_includes.h"	/* definitions files and prototypes */


#define NULL_FP -1 /* NULL field_offset to be used in the optimized version *
                    * of the load_fatlinks subroutine */

    /* Specify paths in orientation in which they appear in the
       forward part of the x component of dslash().  Rotations and
       reflections will be automatically included. Be careful
       about signs of coefficients.  See long comment at bottom. */
#include <quark_action.h>
    /* Include file specifies the basic paths */

#define GOES_FORWARDS(dir) (dir<=TUP)
#define GOES_BACKWARDS(dir) (dir>TUP)

void printpath( int *path, int length );
void path_transport( field_offset src, field_offset dest, int parity,
    int *dir, int length );
void path_transport_hwv( field_offset src, field_offset dest, int parity,
    int *dir, int length );
#ifdef  ASQ_OPTIMIZED_FATTENING
void compute_gen_staple(field_offset staple, int mu, int nu,
                        field_offset link, double coef ) ;
#endif
#ifdef  ASQ_OPTIMIZED_FORCE
#ifndef FN
BOMB THE COMPILE
#endif
void u_shift_fermion(su3_vector *src, su3_vector *dest, int dir ) ;
void add_force_to_mom(su3_vector *back, su3_vector *forw, int dir, double coef);
void side_link_force(int mu, int nu, double coeff, su3_vector *Path,
		     su3_vector *Path_nu, su3_vector *Path_mu, 
		     su3_vector *Path_numu) ;

void u_shift_hw_fermion(half_wilson_vector *src, 
			half_wilson_vector *dest, int dir ) ;
void add_3f_force_to_mom(half_wilson_vector *back,
			 half_wilson_vector *forw, int dir, double coeff[2]) ;
void side_link_3f_force(int mu, int nu, double coeff[2], 
			half_wilson_vector *Path   , 
			half_wilson_vector *Path_nu, 
			half_wilson_vector *Path_mu, 
			half_wilson_vector *Path_numu) ;
#endif

int path_num[MAX_BASIC_PATHS];	/* number of rotations/reflections for each 
					kind */
static double act_path_coeff[MAX_BASIC_PATHS]; /* actual path coefficient     *
                                               * it is equal to path_coeff   *
                                               * if not tadpole improvement  *
                                               * is specified                *
                                               * or path_coeff*u_0^(L-1) when*
                                               * tadpole improvement is      *
                                               * specified                   */

/* Array of structures, for each rotation and reflection of each kind of
	path.  */
struct {
    int dir[MAX_LENGTH];	/* directions in path */
    int length;		/* length of path */
    double coeff;	/* coefficient, including minus sign if backwards */
    double forwback;	/* +1 if in forward Dslash, -1 if in backward */
}   q_paths[MAX_NUM];
int num_q_paths;	/* number of paths in dslash */
int num_basic_paths;	/* number of paths before rotation/reflection */

int is_path_equal( int *path1, int* path2, int length );
int add_basic_path( int *vec, int length, double coeff );

/* Make table of paths in action */
void make_path_table() {

    int i,j;
#ifdef TADPOLE_IMPROVE
    int k;
#endif

    /* table of directions, 1 for each kind of path */
    /**int path_ind[MAX_BASIC_PATHS][MAX_LENGTH];**/
    /* table of coefficients in action, for each path */

    node0_printf("%s\n",quark_action_description);
    num_q_paths = 0;
    num_basic_paths = 0;


    /* add rots. and reflects to table, print out the path coefficients */
    node0_printf("path coefficients: npath  path_coeff  multiplicity\n");
    for(j=0;j<quark_action_npaths;j++) {
	double this_coeff;
	this_coeff = path_coeff[j];
#ifdef TADPOLE_IMPROVE
	for(k=1;k< path_length_in[j];k++)this_coeff /= u0;
#endif
	act_path_coeff[j] = this_coeff ;
	i = add_basic_path( path_ind[j], path_length_in[j],
	    this_coeff );
	node0_printf("                    %d      %e     %d\n",
	    j,this_coeff,i);
    }
}

/* add rotations and reflections of a path to the table.  Return
   multiplicity of paths added */
int add_basic_path( int *basic_vec, int length, double coeff ) {

    int perm[8],pp[8],ir[4];
    int j,path_num;
    int vec[MAX_LENGTH];
    int flag;

    path_num = 0;
    /* now fill the long table with all rotations and reflections
	of the fundamental path.  The path presented to us is for
        the positive x component of dslash, so if the x coordinate
        is reflected it will appear with a negative sign. */
      /* permutations */
      for(perm[0]=0;perm[0]<4;perm[0]++)
      for(perm[1]=0;perm[1]<4;perm[1]++)
      for(perm[2]=0;perm[2]<4;perm[2]++)
      for(perm[3]=0;perm[3]<4;perm[3]++){
	if(perm[0] != perm[1] && perm[0] != perm[2] 
	  && perm[0] != perm[3] && perm[1] != perm[2]
	  && perm[1] != perm[3] && perm[2] != perm[3] ) {
	  /* reflections*/
	  for(ir[0]=0;ir[0]<2;ir[0]++)
	  for(ir[1]=0;ir[1]<2;ir[1]++)
	  for(ir[2]=0;ir[2]<2;ir[2]++)
	  for(ir[3]=0;ir[3]<2;ir[3]++){
	    for(j=0;j<4;j++){
	      pp[j]=perm[j];

	      if(ir[j] == 1) pp[j]=OPP_DIR(pp[j]);
	      pp[OPP_DIR(j)]=OPP_DIR(pp[j]);
	    }
	    /* create new vector*/
	    for(j=0;j<length;j++) vec[j]=pp[basic_vec[j]];
	    for(j=length;j<MAX_LENGTH;j++) vec[j]=NODIR;

            flag=0;
	    /* check if it's a new set: */
	    for(j=0;j<num_q_paths;j++){
	      flag = is_path_equal( vec, q_paths[j].dir, MAX_LENGTH );
	      if(flag==1)break;
	    }
	    if(flag == 0 ){
	      if(num_q_paths>=MAX_NUM){
		node0_printf("OOPS: MAX_NUM too small\n");
		exit(0);
	      }
	      q_paths[num_q_paths].length=length;
	      for(j=0;j<MAX_LENGTH;j++) q_paths[num_q_paths].dir[j]=vec[j];
		/* remember to copy NODIR's, or comparison will get confused */
	      if(ir[0]==0){
		q_paths[num_q_paths].coeff =  coeff;
		q_paths[num_q_paths].forwback =  +1;
	      }
	      else{
		q_paths[num_q_paths].coeff = -coeff;
		q_paths[num_q_paths].forwback = -1;
	      }
	      num_q_paths++;
	      path_num++;
	      /**node0_printf("ADD PATH %d:  rx=%d ",num_q_paths-1,ir[0]);
		 printpath( vec, length );**/
	    }

	  } /* end reflection*/
        } /* end permutation if block */
      } /* end permutation */
    num_basic_paths++;
    return(path_num);
} /* add_basic_path */

/* parallel transport a vector to the current site along a path.
   For example, if the path is "XUP", bring in the vector from
   the +X direction to the current site. If KS phases are in lattice,
   this transport will automatically include them. 
   OK for src and dest to be the same.  OK for length=0.  */
/* NOT OPTIMIZED at the moment - do lots of extra copying.  Use temp.
   vectors rather than stuff in lattice.h */
void path_transport( field_offset src, field_offset dest, int parity,
    int *dir, int length ){
    register int i;
    register site *s;
    msg_tag *mtag0;
    int j;
    su3_vector *tmp_src,*tmp_dest,*tmp_work; /*source, dest and workspace*/
    su3_vector *tmp_pt; /* scratch */
    int tmp_parity, tmp_otherparity; /* parity for this step */

  if( length > 0 ){
    tmp_src = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));
    tmp_dest = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));
    tmp_work = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));

    for( j=length-1; j>=0; j-- ){
	/* figure out parities for this step */
	if( j%2==0 ){
	    tmp_parity = parity;
	    switch(tmp_parity){
		case EVEN: tmp_otherparity=ODD; break;
		case ODD: tmp_otherparity=EVEN; break;
		case EVENANDODD: tmp_otherparity=EVENANDODD; break;
	    }
	}
	else { /* odd # step */
	    tmp_otherparity = parity;
	    switch(tmp_otherparity){
		case EVEN: tmp_parity=ODD; break;
		case ODD: tmp_parity=EVEN; break;
		case EVENANDODD: tmp_parity=EVENANDODD; break;
	    }
	}

	if( j==length-1 ){
	    FORSOMEPARITY(i,s,tmp_otherparity){
	        tmp_src[i] = *(su3_vector *)F_PT(s,src);
	    }
	}

	if( GOES_FORWARDS(dir[j]) ) {
	    mtag0 = start_gather_from_temp( tmp_src, sizeof(su3_vector),
	        dir[j], tmp_parity, gen_pt[0] );
	    wait_gather(mtag0);
	    FORSOMEPARITY(i,s,tmp_parity){
		mult_su3_mat_vec( &(s->link[dir[j]]),
		    (su3_vector *)(gen_pt[0][i]),
		    &(tmp_dest[i]) );
	    }
	    cleanup_gather(mtag0);
	}

	else{ /* GOES_BACKWARDS(dir[j]) */
	    FORSOMEPARITY(i,s,tmp_otherparity){
		mult_adj_su3_mat_vec( &(s->link[OPP_DIR(dir[j])]),
		    &(tmp_src[i]), &(tmp_work[i]) );
	    }
	    mtag0 = start_gather_from_temp( tmp_work, sizeof(su3_vector),
	        dir[j], tmp_parity, gen_pt[0] );
	    wait_gather(mtag0);
	    FORSOMEPARITY(i,s,tmp_parity){
		 tmp_dest[i] = *(su3_vector *)gen_pt[0][i];
	    }
	    cleanup_gather(mtag0);
	}
	
	/* src for next step is dest for this one. */
	tmp_pt=tmp_src; tmp_src=tmp_dest; tmp_dest=tmp_pt;
    }  /* j=link in path */
    /* done, copy result into real dest. (tmp_src now points to result) */
    FORSOMEPARITY(i,s,parity){
        *(su3_vector *)F_PT(s,dest) = tmp_src[i];
    }
    free(tmp_src); free(tmp_dest); free(tmp_work);
  } /* end if(length>0) */
  else if( src != dest ){ /* for length=0 */
    FORSOMEPARITY(i,s,parity){
        *(su3_vector *)F_PT(s,dest) = *(su3_vector *)F_PT(s,src);
    }
  }
} /* path_transport */

/* Path transport a half_wilson_vector */
void path_transport_hwv( field_offset src, field_offset dest, int parity,
    int *dir, int length ){
    register int i;
    register site *s;
    msg_tag *mtag0;
    int j;
    half_wilson_vector *tmp_src,*tmp_dest,*tmp_work; /*source, dest and workspace*/
    half_wilson_vector *tmp_pt; /* scratch */
    int tmp_parity, tmp_otherparity; /* parity for this step */

  if( length > 0 ){
    tmp_src = (half_wilson_vector *)calloc(
      sites_on_node, sizeof(half_wilson_vector));
    tmp_dest = (half_wilson_vector *)calloc(
       sites_on_node, sizeof(half_wilson_vector));
    tmp_work = (half_wilson_vector *)calloc(
       sites_on_node, sizeof(half_wilson_vector));

    for( j=length-1; j>=0; j-- ){
	/* figure out parities for this step */
	if( j%2==0 ){
	    tmp_parity = parity;
	    switch(tmp_parity){
		case EVEN: tmp_otherparity=ODD; break;
		case ODD: tmp_otherparity=EVEN; break;
		case EVENANDODD: tmp_otherparity=EVENANDODD; break;
	    }
	}
	else { /* odd # step */
	    tmp_otherparity = parity;
	    switch(tmp_otherparity){
		case EVEN: tmp_parity=ODD; break;
		case ODD: tmp_parity=EVEN; break;
		case EVENANDODD: tmp_parity=EVENANDODD; break;
	    }
	}

	if( j==length-1 ){
	    FORSOMEPARITY(i,s,tmp_otherparity){
	        tmp_src[i] = *(half_wilson_vector *)F_PT(s,src);
	    }
	}

	if( GOES_FORWARDS(dir[j]) ) {
	    mtag0 = start_gather_from_temp( tmp_src,
		sizeof(half_wilson_vector), dir[j], tmp_parity, gen_pt[0] );
	    wait_gather(mtag0);
	    FORSOMEPARITY(i,s,tmp_parity){
		mult_su3_mat_hwvec( &(s->link[dir[j]]),
		    (half_wilson_vector *)(gen_pt[0][i]),
		    &(tmp_dest[i]) );
	    }
	    cleanup_gather(mtag0);
	}

	else{ /* GOES_BACKWARDS(dir[j]) */
	    FORSOMEPARITY(i,s,tmp_otherparity){
		mult_adj_su3_mat_hwvec( &(s->link[OPP_DIR(dir[j])]),
		    &(tmp_src[i]), &(tmp_work[i]) );
	    }
	    mtag0 = start_gather_from_temp( tmp_work,
		sizeof(half_wilson_vector), dir[j], tmp_parity, gen_pt[0] );
	    wait_gather(mtag0);
	    FORSOMEPARITY(i,s,tmp_parity){
		 tmp_dest[i] = *(half_wilson_vector *)gen_pt[0][i];
	    }
	    cleanup_gather(mtag0);
	}
	
	/* src for next step is dest for this one. */
	tmp_pt=tmp_src; tmp_src=tmp_dest; tmp_dest=tmp_pt;
    }  /* j=link in path */
    /* done, copy result into real dest. (tmp_src now points to result) */
    FORSOMEPARITY(i,s,parity){
        *(half_wilson_vector *)F_PT(s,dest) = tmp_src[i];
    }
    free(tmp_src); free(tmp_dest); free(tmp_work);
  } /* end if(length>0) */
  else if( src != dest ){ /* for length=0 */
    FORSOMEPARITY(i,s,parity){
        *(half_wilson_vector *)F_PT(s,dest) =
	  *(half_wilson_vector *)F_PT(s,src);
    }
  }
} /* path_transport_hwv */


#ifdef EO
/* Stupid dslash routine that follows all the paths.
  Should optimize by precomputing sums of paths to all the displacements
  Use dslash_fn for actions that involve only +X and +X+X+X couplings.
*/
/* D_slash routine - sets dest. on each site equal to sum of
   sources parallel transported to site, with minus sign for transport
   from negative directions.  */
void dslash_eo( field_offset src, field_offset dest, int parity ) {
   register int i;
   register site *s;
   register int ipath,otherparity;
   register double x; /* coefficient of path */

    switch(parity){
	case EVEN:	otherparity=ODD; break;
	case ODD:	otherparity=EVEN; break;
	case EVENANDODD:	otherparity=EVENANDODD; break;
    }

    /* Parallel transport by all the paths in the action.  
       Multiply by coefficient in table
    */
    FORSOMEPARITY(i,s,parity){ clearvec( (su3_vector *)F_PT(s,dest) ); }

    for( ipath=0; ipath<num_q_paths; ipath++ ){  /* loop over paths */
	    path_transport( src, F_OFFSET(tempvec[0]), parity,
		q_paths[ipath].dir, q_paths[ipath].length );
	    x=q_paths[ipath].coeff;
	    FORSOMEPARITY(i,s,parity){
		scalar_mult_add_su3_vector(  (su3_vector *)F_PT(s,dest),
		    &(s->tempvec[0]), x, (su3_vector *)F_PT(s,dest) );
	    }
    }   /* ipath */
} /* dslash_eo */
#endif /* ifdef EO */

#ifdef FN
/* Sum over paths connecting to nearest neighbor point (fat link) and to third
   nearest neighbor (longlinks) */
/* Doug Toussaint 2/4/98 */
/* modified to use t_longlinks, S. Gottlieb 7/13/01 */
/* long link calculating routine */
/* path_product() follows the path starting at step 0, and
   leaves the answer at the end of the path.  We want the answer
   at the site where the path begins.  So we look for paths with
   the opposite displacement from the displacement of the point
   that we want to transport to this site, and take the adjoint
   of the matrix at the end. clear? */
/* KS phases and APBC must be in the links. See long comment at bottom*/
void load_longlinks() {
  register int i;
  register site *s;
  int ipath,dir;
  int disp[4];
  register su3_matrix *long1;

#ifdef LLTIME
double dtime;
dtime=-dclock();
#endif
  if( phases_in != 1){
    node0_printf("BOTCH: load_longlinks needs phases in\n");
    terminate(0);
  }
  for (dir=XUP; dir<=TUP; dir++){ /* loop over longlink directions */
    /* set longlink to zero */
    FORALLSITES(i,s){
#ifdef DSLASH_TMP_LINKS
      long1 = &(t_longlink[4*i+dir]);
#else
      long1 = &(s->longlink[dir]);
#endif
      clear_su3mat( long1 );
    }

    /* loop over paths, checking for ones with total displacement 3*dir */
    for( ipath=0; ipath<num_q_paths; ipath++ ){  /* loop over paths */
	/* compute total displacement of path */
	for(i=XUP;i<=TUP;i++)disp[i]=0;
	for( i=0; i<q_paths[ipath].length; i++){
	  if( GOES_FORWARDS(q_paths[ipath].dir[i]) )
	    disp[        q_paths[ipath].dir[i]  ]++;
	  else
	    disp[OPP_DIR(q_paths[ipath].dir[i]) ]--;
	}
	for( disp[dir]+=3,i=XUP; i<=TUP; i++)if(disp[i]!=0)break;
	if( i<=TUP )continue;  /* skip if path doesn't go to right place */
/**printf("ipath = %d, found a path:  ",ipath);
for(j=0;j<q_paths[ipath].length;j++)printf("\t%d", q_paths[ipath].dir[j]);
printf("\n");**/

	path_product( q_paths[ipath].dir, q_paths[ipath].length );
	FORALLSITES(i,s){
	  su3_adjoint( &(s->tempmat1), &(s->staple) );
#ifdef DSLASH_TMP_LINKS
	  long1 = &(t_longlink[4*i+dir]);
#else
	  long1 = &(s->longlink[dir]);
#endif
          scalar_mult_add_su3_matrix( long1,
	    &(s->staple), -q_paths[ipath].coeff, long1 );
		/* minus sign in coeff. because we used backward path*/
	}
    } /* ipath */

  } /* loop over directions */

  valid_longlinks = 1;
#ifdef LLTIME
dtime += dclock();
node0_printf("LLTIME(Long):  %e\n",dtime);
#endif
}  /* load_longlinks() */

/* KS phases and APBC must be in the links. See long comment at bottom*/
void load_fatlinks() {
  register int i;
  register site *s;
  int dir;
  register su3_matrix *fat1;
#ifdef ASQ_OPTIMIZED_FATTENING
  int  nu,rho,sig ;
  double one_link ; /*needed to fix the problem with the Lepage term */
#else
  int ipath;
  int disp[4];
#endif

#ifdef LLTIME
double dtime;
dtime=-dclock();
#endif
  if( phases_in != 1){
    node0_printf("BOTCH: load_fatlinks needs phases in\n");
    terminate(0);
  }

#ifndef  ASQ_OPTIMIZED_FATTENING   /* general case code */
  for (dir=XUP; dir<=TUP; dir++){ /* loop over fatlink directions */
    /* set fatlink to zero */
    FORALLSITES(i,s){
#ifdef DSLASH_TMP_LINKS
      fat1 = &(t_fatlink[4*i+dir]);
#else
      fat1 = &(s->fatlink[dir]);
#endif
      clear_su3mat( fat1 );
    }
    
    /* loop over paths, checking for ones with total displacement 1*dir */
    for( ipath=0; ipath<num_q_paths; ipath++ ){  /* loop over paths */
	/* compute total displacement of path */
	for(i=XUP;i<=TUP;i++)disp[i]=0;
	for( i=0; i<q_paths[ipath].length; i++){
	  if( GOES_FORWARDS(q_paths[ipath].dir[i]) )
	    disp[        q_paths[ipath].dir[i]  ]++;
	  else
	    disp[OPP_DIR(q_paths[ipath].dir[i]) ]--;
	}
	for( disp[dir]+=1,i=XUP; i<=TUP; i++)if(disp[i]!=0)break;
	if( i<=TUP )continue;  /* skip if path doesn't go to right place */
/**printf("dir = %d, found a path:  ",dir);
for(j=0;j<q_paths.[ipath].length;j++)printf("\t%d", q_paths[ipath].dir[j]);
printf("\n");**/

	path_product( q_paths[ipath].dir, q_paths[ipath].length );
	FORALLSITES(i,s){
	  su3_adjoint( &(s->tempmat1), &(s->staple) );
#ifdef DSLASH_TMP_LINKS
	  fat1 = &(t_fatlink[4*i+dir]);
#else
	  fat1 = &(s->fatlink[dir]);
#endif
          scalar_mult_add_su3_matrix( fat1,
	    &(s->staple), -q_paths[ipath].coeff, fat1 );
		/* minus sign in coeff. because we used backward path*/
	}
    } /* ipath */
  } /* loop over directions */
#else	/* ASQ_OPTIMIZED_FATTENING, for Asq and Asqtad actions */
/*  Optimized fattening code for the Asq and Asqtad actions.           *
 * I assume that path 0 is the one link path 2 the 3-staple            *
 * path 3 the 5-staple path 4 the 7-staple and path 5 the Lapage term. *
 * Path 1 is the Naik term.                                            */
 
 /* to fix up the Lepage term, included by a trick below */
 one_link = (act_path_coeff[0] - 6.0*act_path_coeff[5]) ; 
 
 for (dir=XUP; dir<=TUP; dir++){
   FORALLSITES(i,s) /* Intialize fat links with c_1*U_\mu(x) */
     {
#ifdef DSLASH_TMP_LINKS
       fat1 = &(t_fatlink[4*i+dir]);
#else
       fat1 = &(s->fatlink[dir]);
#endif
       scalar_mult_su3_matrix(&(s->link[dir]), one_link ,
			      fat1 );
     }
   for(nu=XUP; nu<=TUP; nu++) if(nu!=dir)
     {
       compute_gen_staple(F_OFFSET(staple),dir,nu,F_OFFSET(link[dir]),
			  act_path_coeff[2]);
       /* The Lepage term */
       /* Note this also involves modifying c_1 (above) */
       compute_gen_staple(NULL_FP,dir,nu,F_OFFSET(staple),act_path_coeff[5]);
       for(rho=XUP; rho<=TUP; rho++) if((rho!=dir)&&(rho!=nu))
	 {
	   compute_gen_staple(F_OFFSET(tempmat1),dir,rho,F_OFFSET(staple),
			      act_path_coeff[3]);
	   for(sig=XUP; sig<=TUP; sig++)
	     if((sig!=dir)&&(sig!=nu)&&(sig!=rho))
	       {
		 compute_gen_staple(NULL_FP,dir,sig,
				    F_OFFSET(tempmat1),
				    act_path_coeff[4]);
	       } /* sig */
	 } /* rho */
     } /* nu */
 }/* dir */  
#endif

  valid_fatlinks = 1;
#ifdef LLTIME
dtime += dclock();
node0_printf("LLTIME(Fat):  %e\n",dtime);
#endif
}  /* load_fatlinks() */
#endif /* ifdef FN */


/* compare two paths, return 1 if equal, else zero */
int is_path_equal( int *path1, int* path2, int length ){
   register int i;
   for(i=0;i<length;i++)if(path1[i]!=path2[i])return(0);
   return(1);
}


/* update the  momenta with the fermion force */
/* Assumes that the conjugate gradient has been run, with the answer in
   x_off, and dslash(x_off,x_off,ODD) has been run. (fills in x_off_odd) */
/* SEE LONG COMMENTS AT END */

#ifndef ASQ_OPTIMIZED_FORCE
void eo_fermion_force( double eps, int nflavors, field_offset x_off ){
  /* note CG_solution and Dslash * solution are combined in "x_off" */
  /* New version 1/21/99.  Use forward part of Dslash to get force */
  /* see long comment at end */
  /* For each link we need x_off transported from both ends of path. */
  register int i,dir,lastdir,ipath,ilink;
  register site *s;
  int length;
  su3_matrix tmat,tmat2;
  double ferm_epsilon, coeff;
#ifdef FFTIME
  double dtime;
#endif
  msg_tag *mtag0;
  half_wilson_vector *hw_tmp0,*hw_tmp1,*tmp_pt;

#ifdef FFTIME
dtime=-dclock();
#endif
  ferm_epsilon = 2.0*(nflavors/4.0)*eps;
  hw_tmp0 = (half_wilson_vector *)
    calloc(sites_on_node, sizeof(half_wilson_vector));
  hw_tmp1 = (half_wilson_vector *)
    calloc(sites_on_node, sizeof(half_wilson_vector));
  /* Use half_wilson_vectors to store x_off transported from ends of
     path.  0 component from forward end, 1 component from back end */

  /* loop over paths, and loop over links in path */
  for( ipath=0; ipath<num_q_paths; ipath++ ){
    if(q_paths[ipath].forwback== -1)continue;	/* skip backwards dslash */
    length = q_paths[ipath].length;

    /* path transport x_off and Dslash*x_off from far end.  Sometimes
	we need them at the start point of the path, and sometimes
	one link into the path - an optimization for later */
    path_transport( x_off, F_OFFSET(tempvec[0]),
      EVENANDODD, q_paths[ipath].dir, length );
    /* use tempvec[1] for transport from starting end */
    FORALLSITES(i,s){
      hw_tmp0[i].h[0]=s->tempvec[0];
      hw_tmp0[i].h[1]=*(su3_vector *)F_PT(s,x_off);
    }

    /* A path has (length+1) points, counting the ends.  At first
	 point, no "down" direction links have their momenta "at this
	 point". At last, no "up" ... */
    for( ilink=0; ilink<=length; ilink++ ){
      if(ilink<length)dir = q_paths[ipath].dir[ilink];
      else dir=NODIR;
      coeff = ferm_epsilon*q_paths[ipath].coeff;
      if( (ilink%2)==1 )coeff = -coeff;

      /* path transport x_off and Dslash*x_off from previous point */
      /* Use "half_wilson_vector" to handle pair of vectors -
	0 component is x_off from forward end, 1 component from back end */
      /* sometimes we don't need them */
      if( (ilink>0&&ilink<length) || 
	(ilink==length && GOES_BACKWARDS(lastdir)) ){

	if( GOES_FORWARDS(lastdir) ){
	  FORALLSITES(i,s){
            mult_adj_su3_mat_hwvec( &(s->link[lastdir]),
	      &(hw_tmp0[i]), &(hw_tmp1[i]) );
	  }
	  mtag0 = start_gather_from_temp( hw_tmp1, 2*sizeof(su3_vector),
             OPP_DIR(lastdir), EVENANDODD, gen_pt[0] );
          wait_gather(mtag0);
          FORALLSITES(i,s){
	     hw_tmp0[i] = *(half_wilson_vector *)gen_pt[0][i];
	  }
          cleanup_gather(mtag0);
	}
	else{   /* GOES_BACKWARDS(lastdir) */
          mtag0 = start_gather_from_temp( hw_tmp0, 2*sizeof(su3_vector),
               OPP_DIR(lastdir), EVENANDODD, gen_pt[0] );
          wait_gather(mtag0);
          FORALLSITES(i,s){
            mult_su3_mat_hwvec( &(s->link[OPP_DIR(lastdir)]),
                (half_wilson_vector *)(gen_pt[0][i]),
		&(hw_tmp1[i]) );
          }
	  tmp_pt = hw_tmp0; hw_tmp0 = hw_tmp1; hw_tmp1 = tmp_pt;
          cleanup_gather(mtag0);
	}
      }

      /* add in contribution to the force */
      /* Put antihermitian traceless part into momentum */
      FORALLSITES(i,s){
        if( ilink<length && GOES_FORWARDS(dir) ){
          uncompress_anti_hermitian( &(s->mom[dir]), &tmat2 );
	  su3_projector( &(hw_tmp0[i].h[0]), &(hw_tmp0[i].h[1]), &tmat );
	  if( s->parity==EVEN )scalar_mult_add_su3_matrix(
	    &tmat2, &tmat,  coeff, &tmat2 );
	  else		 scalar_mult_add_su3_matrix(
	    &tmat2, &tmat, -coeff, &tmat2 );
          make_anti_hermitian( &tmat2, &(s->mom[dir]) );
        }
	if( ilink>0 && GOES_BACKWARDS(lastdir) ){
          uncompress_anti_hermitian( &(s->mom[OPP_DIR(lastdir)]), &tmat2 );
	  su3_projector( &(hw_tmp0[i].h[0]), &(hw_tmp0[i].h[1]), &tmat );
	  if( s->parity==EVEN )scalar_mult_add_su3_matrix(
	    &tmat2, &tmat, -coeff, &tmat2 );
	  else		 scalar_mult_add_su3_matrix(
	    &tmat2, &tmat,  coeff, &tmat2 );
          make_anti_hermitian( &tmat2, &(s->mom[OPP_DIR(lastdir)]) );
        }
      }
      lastdir = dir;
    } /* end loop over links in path */
  } /* end loop over paths */

  free( hw_tmp0 ); free( hw_tmp1 );

#ifdef FFTIME
dtime += dclock();
node0_printf("FFTIME:  %e\n",dtime);
/**printf("TLENGTH: %d\n",tlength);**/
#endif
} /* eo_fermion_force(version 6) */


void eo_fermion_force_3f( double eps, int nflav1, field_offset x1_off,
  int nflav2, field_offset x2_off ){
  /* note CG_solution and Dslash * solution are combined in "x_off" */
  /* New version 1/21/99.  Use forward part of Dslash to get force */
  /* 4/14/99 combine force from two different mass quarks, (eg 2+1flavors) */
  /* see long comment at end */
  /* For each link we need x_off transported from both ends of path. */
  register int i,dir,lastdir,ipath,ilink;
  register site *s;
  int length;
  su3_matrix tmat,tmat2;
  double ferm_epsilon1, ferm_epsilon2, coeff1, coeff2;
#ifdef FFTIME
  double dtime;
#endif
  msg_tag *mtag0;
  wilson_vector *w_tmp0,*w_tmp1,*tmp_pt;

#ifdef FFTIME
dtime=-dclock();
#endif
  ferm_epsilon1 = 2.0*(nflav1/4.0)*eps;
  ferm_epsilon2 = 2.0*(nflav2/4.0)*eps;
  w_tmp0 = (wilson_vector *)
    calloc(sites_on_node, sizeof(wilson_vector));
  w_tmp1 = (wilson_vector *)
    calloc(sites_on_node, sizeof(wilson_vector));
  /* Use wilson_vectors to store x_off transported from ends of
     path.  0 and 1 components from forward end, 2 and 3 components
     from back end */

  /* loop over paths, and loop over links in path */
  for( ipath=0; ipath<num_q_paths; ipath++ ){
    if(q_paths[ipath].forwback== -1)continue;	/* skip backwards dslash */
    length = q_paths[ipath].length;

    /* path transport x_off and Dslash*x_off from far end.  Sometimes
	we need them at the start point of the path, and sometimes
	one link into the path - an optimization for later */
    /**
    path_transport( x1_off, F_OFFSET(tempvec[0]),
      EVENANDODD, q_paths[ipath].dir, length );
    path_transport( x2_off, F_OFFSET(tempvec[1]),
      EVENANDODD, q_paths[ipath].dir, length );
    **/
/** WARNING!! Assumes xxx1 and xxx2 contiguous **/
if( x2_off-x1_off != sizeof(su3_vector) ){node0_printf("BOTCH\n"); exit(0);}
    path_transport_hwv( x1_off, F_OFFSET(tempvec[0]),
      EVENANDODD, q_paths[ipath].dir, length );
    /* use tempvec[2] for transport from starting end */
    FORALLSITES(i,s){
      w_tmp0[i].d[0]=s->tempvec[0];
      w_tmp0[i].d[1]=s->tempvec[1];
      w_tmp0[i].d[2]=*(su3_vector *)F_PT(s,x1_off);
      w_tmp0[i].d[3]=*(su3_vector *)F_PT(s,x2_off);
    }

    /* A path has (length+1) points, counting the ends.  At first
	 point, no "down" direction links have their momenta "at this
	 point". At last, no "up" ... */
    for( ilink=0; ilink<=length; ilink++ ){
      if(ilink<length)dir = q_paths[ipath].dir[ilink];
      else dir=NODIR;
      coeff1 = ferm_epsilon1*q_paths[ipath].coeff;
      coeff2 = ferm_epsilon2*q_paths[ipath].coeff;
      if( (ilink%2)==1 ){ coeff1 = -coeff1; coeff2 = -coeff2;}

      /* path transport x_off and Dslash*x_off from previous point */
      /* Use "wilson_vector" to handle pair of vectors -
	0 component is x_off1 from forward end, 1 component is x_off2
	from forward end,  2 and 3  components are x_off1 and x_off2 
	from back end */
      /* sometimes we don't need them */
      if( (ilink>0&&ilink<length) || 
	(ilink==length && GOES_BACKWARDS(lastdir)) ){

	if( GOES_FORWARDS(lastdir) ){
	  FORALLSITES(i,s){
            mult_adj_mat_wilson_vec( &(s->link[lastdir]),
	      &(w_tmp0[i]), &(w_tmp1[i]) );
	  }
	  mtag0 = start_gather_from_temp( w_tmp1, 4*sizeof(su3_vector),
             OPP_DIR(lastdir), EVENANDODD, gen_pt[0] );
          wait_gather(mtag0);
          FORALLSITES(i,s){
	     w_tmp0[i] = *(wilson_vector *)gen_pt[0][i];
	  }
          cleanup_gather(mtag0);
	}
	else{   /* GOES_BACKWARDS(lastdir) */
          mtag0 = start_gather_from_temp( w_tmp0, 4*sizeof(su3_vector),
               OPP_DIR(lastdir), EVENANDODD, gen_pt[0] );
          wait_gather(mtag0);
          FORALLSITES(i,s){
            mult_mat_wilson_vec( &(s->link[OPP_DIR(lastdir)]),
                (wilson_vector *)(gen_pt[0][i]),
		&(w_tmp1[i]) );
          }
	  tmp_pt = w_tmp0; w_tmp0 = w_tmp1; w_tmp1 = tmp_pt;
          cleanup_gather(mtag0);
	}
      }

      /* add in contribution to the force */
      /* Put antihermitian traceless part into momentum */
      FORALLSITES(i,s){
        if( ilink<length && GOES_FORWARDS(dir) ){
	  if( s->parity==ODD ){coeff1 *= -1.0; coeff2 *= -1.0; }
          uncompress_anti_hermitian( &(s->mom[dir]), &tmat2 );
	  su3_projector( &(w_tmp0[i].d[0]), &(w_tmp0[i].d[2]), &tmat );
	  scalar_mult_add_su3_matrix( &tmat2, &tmat,  coeff1, &tmat2 );
	  su3_projector( &(w_tmp0[i].d[1]), &(w_tmp0[i].d[3]), &tmat );
	  scalar_mult_add_su3_matrix( &tmat2, &tmat,  coeff2, &tmat2 );
          make_anti_hermitian( &tmat2, &(s->mom[dir]) );
	  if( s->parity==ODD ){coeff1 *= -1.0; coeff2 *= -1.0; }
        }
	if( ilink>0 && GOES_BACKWARDS(lastdir) ){
	  if( s->parity==EVEN ){coeff1 *= -1.0; coeff2 *= -1.0; }
          uncompress_anti_hermitian( &(s->mom[OPP_DIR(lastdir)]), &tmat2 );
	  su3_projector( &(w_tmp0[i].d[0]), &(w_tmp0[i].d[2]), &tmat );
	  scalar_mult_add_su3_matrix( &tmat2, &tmat,  coeff1, &tmat2 );
	  su3_projector( &(w_tmp0[i].d[1]), &(w_tmp0[i].d[3]), &tmat );
	  scalar_mult_add_su3_matrix( &tmat2, &tmat,  coeff2, &tmat2 );
          make_anti_hermitian( &tmat2, &(s->mom[OPP_DIR(lastdir)]) );
	  if( s->parity==EVEN ){coeff1 *= -1.0; coeff2 *= -1.0; }
        }
      }
      lastdir = dir;
    } /* end loop over links in path */
  } /* end loop over paths */

  free( w_tmp0 ); free( w_tmp1 );

#ifdef FFTIME
dtime += dclock();
node0_printf("FFTIME:  %e\n",dtime);
/**printf("TLENGTH: %d\n",tlength);**/
#endif
} /* eo_fermion_force_3f(version 6) */



#else /* ASQ_OPTIMIZED_FORCE, for Asq and Asqtad actions */
/* Optimized force code for the Asq and Asqtad actions                 *
 * I assume that path 0 is the one link path 2 the 3-staple            *
 * path 3 the 5-staple path 4 the 7-staple and path 5 the Lapage term. *
 * Path 1 is the Naik term.                                            */
#define Pmu          tempvec[0] 
#define Pnumu        tempvec[1]
#define Prhonumu     tempvec[2]
#define P7           tempvec[3]
#define P7rho        tempvec[4]              
#define P7rhonu      tempvec[5]
#define P5           tempvec[6]
#define P3           tempvec[7]
#define P5nu         tempvec[3]
#define P3mu         tempvec[3]
#define Popmu        tempvec[4]
#define Pmumumu      tempvec[4]
void eo_fermion_force( double eps, int nflavors, field_offset x_off ){
  /* note CG_solution and Dslash * solution are combined in "x_off" */
  /* New version 1/21/99.  Use forward part of Dslash to get force */
  /* see long comment at end */
  /* For each link we need x_off transported from both ends of path. */
  register int i ;
  register site *s;
  int mu,nu,rho,sig ;
  int DirectLinks[8] ;
  double ferm_epsilon, coeff;
  double OneLink, Lepage, Naik, FiveSt, ThreeSt, SevenSt ;
  su3_vector *tempvec[8] ;
  su3_vector *temp_x ;

#ifdef FFTIME
  double dtime;

  dtime=-dclock();
#endif
  ferm_epsilon = 2.0*(nflavors/4.0)*eps;
  
  /* Path coefficients times fermion epsilon */
  OneLink = act_path_coeff[0]*ferm_epsilon ; 
  Naik    = act_path_coeff[1]*ferm_epsilon ;
  ThreeSt = act_path_coeff[2]*ferm_epsilon ;
  FiveSt  = act_path_coeff[3]*ferm_epsilon ;
  SevenSt = act_path_coeff[4]*ferm_epsilon ;
  Lepage  = act_path_coeff[5]*ferm_epsilon ;
  /* *************************************** */

  /* Initialize the DirectLink flags */
  for(mu=0;mu<8;mu++)
    DirectLinks[mu] = 0 ;

  /* Allocate temporary vectors */
  for(mu=0;mu<8;mu++)
    tempvec[mu] = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));

  /*copy x_off to a temporary vector */
  temp_x = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));
  FORALLSITES(i,s) temp_x[i] = *(su3_vector *)F_PT(s,x_off) ;

  for(sig=0;sig<8;sig++)
    {
      for(mu=0;mu<8;mu++)if((mu!=sig)&&(mu!=OPP_DIR(sig)))
	{
	  u_shift_fermion(temp_x, Pmu, OPP_DIR(mu));
	  u_shift_fermion(Pmu, P3, sig);
	  if(GOES_FORWARDS(sig))
	    {
	      /* Add the force F_sig[x+mu]:         x--+             *
	       *                                   |   |             *
	       *                                   o   o             *
	       * the 1 link in the path: - (numbering starts form 0) */
	      add_force_to_mom(P3, Pmu, sig, -ThreeSt) ;
	    }
	  for(nu=0;nu<8;nu++)if((nu!=mu )&&(nu!=OPP_DIR(mu ))&&
				(nu!=sig)&&(nu!=OPP_DIR(sig)))
	    {
	      u_shift_fermion(Pmu, Pnumu, OPP_DIR(nu));
	      u_shift_fermion(Pnumu, P5, sig);
	      if(GOES_FORWARDS(sig))
		{
		  /* Add the force F_sig[x+mu+nu]:      x--+             *
		   *                                   |   |             *
		   *                                   o   o             *
		   * the 2 link in the path: + (numbering starts form 0) */
		  add_force_to_mom(P5, Pnumu, sig, FiveSt);
		}
	      for(rho=0;rho<8;rho++)if((rho!=mu )&&(rho!=OPP_DIR(mu ))&&
				       (rho!=nu )&&(rho!=OPP_DIR(nu ))&&
				       (rho!=sig)&&(rho!=OPP_DIR(sig)))
		{
		  u_shift_fermion(Pnumu, Prhonumu, OPP_DIR(rho));
		  /* Length 7 paths */
		  u_shift_fermion(Prhonumu, P7,sig);
		  if(GOES_FORWARDS(sig))
		    {
		      /* Add the force F_sig[x+mu+nu+rho]:  x--+             *
		       *                                   |   |             *
		       *                                   o   o             *
		       * the 3 link in the path: - (numbering starts form 0) */
		      add_force_to_mom(P7, Prhonumu, sig, -SevenSt ) ;
		    }
		  /*Add the force F_rho the 2(4) link in the path: +     */
		  u_shift_fermion(P7, P7rho, rho);
		  side_link_force(rho,sig,SevenSt, Pnumu, P7, Prhonumu, P7rho);
		  /* Add the P7rho vector to P5 */
		  coeff = SevenSt/FiveSt ; 
		  FORALLSITES(i,s)
		    scalar_mult_add_su3_vector(&(P5[i]),&(P7rho[i]),coeff,
					       &(P5[i]));
		}/* rho */
	      /* Length 5 paths */
	      /*Add the force F_nu the 1(3) link in the path: -     */
	      u_shift_fermion(P5,P5nu, nu);
	      side_link_force(nu,sig,-FiveSt,Pmu,P5, 
			      Pnumu,P5nu) ;
	      /* Add the P5nu vector to P3 */
	      coeff = FiveSt/ThreeSt ; 
	      FORALLSITES(i,s)
		scalar_mult_add_su3_vector(&(P3[i]),&(P5nu[i]),coeff,&(P3[i]));
	    }/* nu */

	  /* Now the Lepage term... It is the same with 5-link paths with
             nu=mu and FiveSt=Lepage. So Pnumu is really Pmumu */
	  u_shift_fermion(Pmu, Pnumu, OPP_DIR(mu));
	  u_shift_fermion(Pnumu, P5, sig);
	  if(GOES_FORWARDS(sig))
	    {
	      /* Add the force F_sig[x+mu+nu]:      x--+             *
	       *                                   |   |             *
	       *                                   o   o             *
	       * the 2 link in the path: + (numbering starts form 0) */
	      add_force_to_mom(P5, Pnumu, sig, Lepage) ;
	    }
	  /*Add the force F_nu the 1(3) link in the path: -     */
	  u_shift_fermion(P5,P5nu, mu);
	  side_link_force(mu, sig, -Lepage, Pmu, P5, Pnumu, P5nu) ;
	  /* Add the P5nu vector to P3 */
	  coeff = Lepage/ThreeSt ; 
	  FORALLSITES(i,s)
	    scalar_mult_add_su3_vector(&(P3[i]),&(P5nu[i]),coeff,&(P3[i]));

	  /* Length 3 paths (Not the Naik term) */
	  /*Add the force F_mu the 0(2) link in the path: +     */
	  if(GOES_FORWARDS(mu)) 
	    u_shift_fermion(P3,P3mu, mu );
	  /* The above shift is not needed if mu is backwards */
	  side_link_force(mu, sig, ThreeSt, temp_x, P3, Pmu, P3mu);

	  /* Finally the OneLink and the Naik term */
	  /* Check if this direction is not already done */
	  if( (!DirectLinks[mu]) ){
	    DirectLinks[mu]=1 ;
	    if(GOES_BACKWARDS(mu))/* Do only the forward terms in the Dslash */
	      {
		/* Because I have shifted with OPP_DIR(mu) Pmu is a forward *
		 * shift.                                                   */
		/* The one link */
		add_force_to_mom(Pmu, temp_x, OPP_DIR(mu), OneLink) ;
		/* For the same reason Pnumu is the forward double link */

		/* Popmu is a backward shift */
		u_shift_fermion(temp_x, Popmu, mu);
		/* The Naik */
		/* link no 1: - */
		add_force_to_mom(Pnumu, Popmu, OPP_DIR(mu), -Naik) ;
		/*Pmumumu can overwrite Popmu which is no longer needed */
		u_shift_fermion(Pnumu, Pmumumu, OPP_DIR(mu));
		/* link no 0: + */
		add_force_to_mom(Pmumumu, temp_x, OPP_DIR(mu), Naik);
	      }
	    else /* The rest of the Naik terms */
	      {
		u_shift_fermion(temp_x, Popmu, mu);
		/* link no 2: + */
		/* Pnumu is double backward shift */
		add_force_to_mom(Popmu, Pnumu, mu, Naik) ;
	      }
	  }
	}/* mu */
      /* Here we have to do together the Naik term and the one link term */
    }/*sig */

  /* Free temporary vectors */
  free(temp_x) ;
  for(mu=0;mu<8;mu++)
    free(tempvec[mu]) ;

#ifdef FFTIME
  dtime += dclock();
node0_printf("FFTIME:  %e\n",dtime);
/**printf("TLENGTH: %d\n",tlength);**/
#endif
} /* eo_fermion_force(version 6) */
#undef Pmu          
#undef Pnumu        
#undef Prhonumu     
#undef P7           
#undef P7rho        
#undef P7rhonu      
#undef P5           
#undef P3           
#undef P5nu         
#undef P3mu         
#undef Popmu        
#undef Pmumumu      

#define Pmu          hwvec[0] 
#define Pnumu        hwvec[1]
#define Prhonumu     hwvec[2]
#define P7           hwvec[3]
#define P7rho        hwvec[4]              
#define P7rhonu      hwvec[5]
#define P5           hwvec[6]
#define P3           hwvec[7]
#define P5nu         hwvec[3]
#define P3mu         hwvec[3]
#define Popmu        hwvec[4]
#define Pmumumu      hwvec[4]
void eo_fermion_force_3f( double eps, int nflav1, field_offset x1_off, 
			  int nflav2, field_offset x2_off){
  /* note CG_solution and Dslash * solution are combined in "x_off" */
  /* New version 1/21/99.  Use forward part of Dslash to get force */
  /* 4/15/99 combine force from two different mass quarks, (eg 2+1flavors) */
  /* see long comment at end */
  /* For each link we need x_off transported from both ends of path. */
  register int i ;
  register site *s;
  int mu,nu,rho,sig ;
  int DirectLinks[8] ;
  double coeff[2],ferm_epsilon ;
  double OneLink[2], Lepage[2], Naik[2], FiveSt[2], ThreeSt[2], SevenSt[2] ;
  double mNaik[2], mLepage[2], mFiveSt[2], mThreeSt[2], mSevenSt[2] ;
  half_wilson_vector *hwvec[8] ;
  half_wilson_vector *temp_x ;

#ifdef FFTIME
  double dtime;

  dtime=-dclock();
#endif
  
  /* Path coefficients times fermion epsilon */
  ferm_epsilon = 2.0*(nflav1/4.0)*eps;
  OneLink[0] = act_path_coeff[0]*ferm_epsilon ;
  Naik[0]    = act_path_coeff[1]*ferm_epsilon ; mNaik[0]    = -Naik[0] ;
  ThreeSt[0] = act_path_coeff[2]*ferm_epsilon ; mThreeSt[0] = -ThreeSt[0] ;
  FiveSt[0]  = act_path_coeff[3]*ferm_epsilon ; mFiveSt[0]  = -FiveSt[0] ;
  SevenSt[0] = act_path_coeff[4]*ferm_epsilon ; mSevenSt[0] = -SevenSt[0] ;
  Lepage[0]  = act_path_coeff[5]*ferm_epsilon ; mLepage[0]  = -Lepage[0] ;

  ferm_epsilon = 2.0*(nflav2/4.0)*eps;
  OneLink[1] = act_path_coeff[0]*ferm_epsilon ;
  Naik[1]    = act_path_coeff[1]*ferm_epsilon ; mNaik[1]    = -Naik[1] ;
  ThreeSt[1] = act_path_coeff[2]*ferm_epsilon ; mThreeSt[1] = -ThreeSt[1] ;
  FiveSt[1]  = act_path_coeff[3]*ferm_epsilon ; mFiveSt[1]  = -FiveSt[1] ;
  SevenSt[1] = act_path_coeff[4]*ferm_epsilon ; mSevenSt[1] = -SevenSt[1] ;
  Lepage[1]  = act_path_coeff[5]*ferm_epsilon ; mLepage[1]  = -Lepage[1] ;
  /* *************************************** */

  /* Initialize the DirectLink flags */
  for(mu=0;mu<8;mu++)
    DirectLinks[mu] = 0 ;

  /* Allocate temporary vectors */
  for(mu=0;mu<8;mu++)
    hwvec[mu]= 
      (half_wilson_vector *)calloc(sites_on_node, sizeof(half_wilson_vector));

  /*copy x_off to a temporary vector */
  temp_x= 
    (half_wilson_vector *)calloc(sites_on_node, sizeof(half_wilson_vector));
  FORALLSITES(i,s)
    {
      temp_x[i].h[0] = *(su3_vector *)F_PT(s,x1_off) ;
      temp_x[i].h[1] = *(su3_vector *)F_PT(s,x2_off) ;
    }

  for(sig=0;sig<8;sig++)
    {
      for(mu=0;mu<8;mu++)if((mu!=sig)&&(mu!=OPP_DIR(sig)))
	{
	  u_shift_hw_fermion(temp_x, Pmu, OPP_DIR(mu));
	  u_shift_hw_fermion(Pmu, P3, sig);
	  if(GOES_FORWARDS(sig))
	    {
	      /* Add the force F_sig[x+mu]:         x--+             *
	       *                                   |   |             *
	       *                                   o   o             *
	       * the 1 link in the path: - (numbering starts form 0) */
	      add_3f_force_to_mom(P3, Pmu, sig, mThreeSt) ;
	    }
	  for(nu=0;nu<8;nu++)if((nu!=mu )&&(nu!=OPP_DIR(mu ))&&
				(nu!=sig)&&(nu!=OPP_DIR(sig)))
	    {
	      u_shift_hw_fermion(Pmu, Pnumu, OPP_DIR(nu));
	      u_shift_hw_fermion(Pnumu, P5, sig);
	      if(GOES_FORWARDS(sig))
		{
		  /* Add the force F_sig[x+mu+nu]:      x--+             *
		   *                                   |   |             *
		   *                                   o   o             *
		   * the 2 link in the path: + (numbering starts form 0) */
		  add_3f_force_to_mom(P5, Pnumu, sig, FiveSt);
		}
	      for(rho=0;rho<8;rho++)if((rho!=mu )&&(rho!=OPP_DIR(mu ))&&
				       (rho!=nu )&&(rho!=OPP_DIR(nu ))&&
				       (rho!=sig)&&(rho!=OPP_DIR(sig)))
		{
		  u_shift_hw_fermion(Pnumu, Prhonumu, OPP_DIR(rho));
		  /* Length 7 paths */
		  u_shift_hw_fermion(Prhonumu, P7,sig);
		  if(GOES_FORWARDS(sig))
		    {
		      /* Add the force F_sig[x+mu+nu+rho]:  x--+             *
		       *                                   |   |             *
		       *                                   o   o             *
		       * the 3 link in the path: - (numbering starts form 0) */
		      add_3f_force_to_mom(P7, Prhonumu, sig, mSevenSt ) ;
		    }
		  /*Add the force F_rho the 2(4) link in the path: +     */
		  u_shift_hw_fermion(P7, P7rho, rho);
		  side_link_3f_force(rho,sig,SevenSt,Pnumu,P7,Prhonumu,P7rho);
		  /* Add the P7rho vector to P5 */
		  coeff[0] = SevenSt[0]/FiveSt[0] ; 
		  coeff[1] = SevenSt[1]/FiveSt[1] ; 
		  FORALLSITES(i,s)
		    {
		      scalar_mult_add_su3_vector(&(P5[i].h[0]),
						 &(P7rho[i].h[0]),coeff[0],
						 &(P5[i].h[0]));
		      scalar_mult_add_su3_vector(&(P5[i].h[1]),
						 &(P7rho[i].h[1]),coeff[1],
						 &(P5[i].h[1]));
		    }
		}/* rho */
	      /* Length 5 paths */
	      /*Add the force F_nu the 1(3) link in the path: -     */
	      u_shift_hw_fermion(P5,P5nu, nu);
	      side_link_3f_force(nu,sig,mFiveSt,Pmu,P5, 
			      Pnumu,P5nu) ;
	      /* Add the P5nu vector to P3 */
	      coeff[0] = FiveSt[0]/ThreeSt[0] ; 
	      coeff[1] = FiveSt[1]/ThreeSt[1] ; 
	      FORALLSITES(i,s)
		{
		  scalar_mult_add_su3_vector(&(P3[i].h[0]),
					     &(P5nu[i].h[0]), coeff[0],
					     &(P3[i].h[0]));
		  scalar_mult_add_su3_vector(&(P3[i].h[1]),
					     &(P5nu[i].h[1]), coeff[1],
					     &(P3[i].h[1]));
		}
	    }/* nu */

	  /* Now the Lepage term... It is the same with 5-link paths with
             nu=mu and FiveSt=Lepage. So Pnumu is really Pmumu */
	  u_shift_hw_fermion(Pmu, Pnumu, OPP_DIR(mu));
	  u_shift_hw_fermion(Pnumu, P5, sig);
	  if(GOES_FORWARDS(sig))
	    {
	      /* Add the force F_sig[x+mu+nu]:      x--+             *
	       *                                   |   |             *
	       *                                   o   o             *
	       * the 2 link in the path: + (numbering starts form 0) */
	      add_3f_force_to_mom(P5, Pnumu, sig, Lepage) ;
	    }
	  /*Add the force F_nu the 1(3) link in the path: -     */
	  u_shift_hw_fermion(P5,P5nu, mu);
	  side_link_3f_force(mu, sig, mLepage, Pmu, P5, Pnumu, P5nu) ;
	  /* Add the P5nu vector to P3 */
	  coeff[0] = Lepage[0]/ThreeSt[0] ; 
	  coeff[1] = Lepage[1]/ThreeSt[1] ; 
	  FORALLSITES(i,s)
	    {
	      scalar_mult_add_su3_vector(&(P3[i].h[0]),
					 &(P5nu[i].h[0]),coeff[0],
					 &(P3[i].h[0]));
	      scalar_mult_add_su3_vector(&(P3[i].h[1]),
					 &(P5nu[i].h[1]),coeff[1],
					 &(P3[i].h[1]));
	    }

	  /* Length 3 paths (Not the Naik term) */
	  /*Add the force F_mu the 0(2) link in the path: +     */
	  if(GOES_FORWARDS(mu)) 
	    u_shift_hw_fermion(P3,P3mu, mu );
	  /* The above shift is not needed if mu is backwards */
	  side_link_3f_force(mu, sig, ThreeSt, temp_x, P3, Pmu, P3mu);

	  /* Finally the OneLink and the Naik term */
	  /* Check if this direction is not already done */
	  if( (!DirectLinks[mu]) ){
	    DirectLinks[mu]=1 ;
	    if(GOES_BACKWARDS(mu))/* Do only the forward terms in the Dslash */
	      {
		/* Because I have shifted with OPP_DIR(mu) Pmu is a forward *
		 * shift.                                                   */
		/* The one link */
		add_3f_force_to_mom(Pmu, temp_x, OPP_DIR(mu), OneLink) ;
		/* For the same reason Pnumu is the forward double link */

		/* Popmu is a backward shift */
		u_shift_hw_fermion(temp_x, Popmu, mu);
		/* The Naik */
		/* link no 1: - */
		add_3f_force_to_mom(Pnumu, Popmu, OPP_DIR(mu), mNaik) ;
		/*Pmumumu can overwrite Popmu which is no longer needed */
		u_shift_hw_fermion(Pnumu, Pmumumu, OPP_DIR(mu));
		/* link no 0: + */
		add_3f_force_to_mom(Pmumumu, temp_x, OPP_DIR(mu), Naik);
	      }
	    else /* The rest of the Naik terms */
	      {
		u_shift_hw_fermion(temp_x, Popmu, mu);
		/* link no 2: + */
		/* Pnumu is double backward shift */
		add_3f_force_to_mom(Popmu, Pnumu, mu, Naik) ;
	      }
	  }
	}/* mu */
      /* Here we have to do together the Naik term and the one link term */
    }/*sig */

  /* Free temporary vectors */
  free(temp_x) ;
  for(mu=0;mu<8;mu++)
    free(hwvec[mu]) ;

#ifdef FFTIME
  dtime += dclock();
node0_printf("FFTIME:  %e\n",dtime);
/**printf("TLENGTH: %d\n",tlength);**/
#endif
} /* eo_fermion_force_3f(version 6) */
#undef Pmu          
#undef Pnumu        
#undef Prhonumu     
#undef P7           
#undef P7rho        
#undef P7rhonu      
#undef P5           
#undef P3           
#undef P5nu         
#undef P3mu         
#undef Popmu        
#undef Pmumumu      
#endif /* ASQ_OPTIMIZED_FORCE */


#ifdef  ASQ_OPTIMIZED_FATTENING   /* Asqtad action only, "_fn" executables */
#ifndef FN
BOMB THE COMPILE
#endif
void compute_gen_staple(field_offset staple, int mu, int nu, 
			field_offset link, double coef) {
  su3_matrix tmat1,tmat2;
  msg_tag *mtag0,*mtag1;
  su3_matrix *tempmat ;
  register site *s ;
  register int i ;
  register su3_matrix *fat1;

  /* Computes the staple :
                   mu
                +-------+
          nu	|	|
		|	|
		X	X
    Where the mu link can be any su3_matrix. The result is saved in staple.
    if staple==NULL_FP then the result is not saved.
    It also adds the computed staple to the fatlink[mu] with weight coef.
  */

  /* Upper staple */
  mtag0 = start_gather( link, sizeof(su3_matrix), nu, EVENANDODD, gen_pt[0] );
  mtag1 = start_gather( F_OFFSET(link[nu]), sizeof(su3_matrix), mu, 
			EVENANDODD, gen_pt[1] );
  wait_gather(mtag0);
  wait_gather(mtag1);
  
  if(staple!=NULL_FP){/* Save the staple */
    FORALLSITES(i,s){
      mult_su3_na( (su3_matrix *)gen_pt[0][i],
		   (su3_matrix *)gen_pt[1][i], &tmat1 );
      mult_su3_nn( &(s->link[nu]), &tmat1, (su3_matrix *)F_PT(s,staple) );
    }
  }
  else{ /* No need to save the staple. Add it to the fatlinks */
    FORALLSITES(i,s){
      mult_su3_na( (su3_matrix *)gen_pt[0][i],
		   (su3_matrix *)gen_pt[1][i], &tmat1 );
      mult_su3_nn( &(s->link[nu]), &tmat1, &tmat2 );
#ifdef DSLASH_TMP_LINKS
      fat1 = &(t_fatlink[4*i+mu]);
#else
      fat1 = &(s->fatlink[mu]);
#endif
      scalar_mult_add_su3_matrix(fat1, &tmat2, coef,
				 fat1) ;
    }
  }
  cleanup_gather(mtag0);
  cleanup_gather(mtag1);

  /* lower staple */
  tempmat = (su3_matrix *)calloc(sites_on_node, sizeof(su3_matrix));
  mtag0 = start_gather( F_OFFSET(link[nu]),
			sizeof(su3_matrix), mu, EVENANDODD, gen_pt[0] );
  wait_gather(mtag0);
  FORALLSITES(i,s){
    mult_su3_an( &(s->link[nu]),(su3_matrix *)F_PT(s,link), &tmat1 );
    mult_su3_nn( &(tmat1),(su3_matrix *)gen_pt[0][i], &(tempmat[i]) );
  }
  cleanup_gather(mtag0);
  mtag0 = start_gather_from_temp( tempmat, sizeof(su3_matrix),
				  OPP_DIR(nu), EVENANDODD, gen_pt[0] );
  wait_gather(mtag0);

  if(staple!=NULL_FP){/* Save the staple */
    FORALLSITES(i,s){
      add_su3_matrix( (su3_matrix *)F_PT(s,staple),(su3_matrix *)gen_pt[0][i], 
		      (su3_matrix *)F_PT(s,staple) );
#ifdef DSLASH_TMP_LINKS
      fat1 = &(t_fatlink[4*i+mu]);
#else
      fat1 = &(s->fatlink[mu]);
#endif
      scalar_mult_add_su3_matrix( fat1,
				 (su3_matrix *)F_PT(s,staple), coef, 
				 fat1 );
    }
  }
  else{ /* No need to save the staple. Add it to the fatlinks */
    FORALLSITES(i,s){
#ifdef DSLASH_TMP_LINKS
      fat1 = &(t_fatlink[4*i+mu]);
#else
      fat1 = &(s->fatlink[mu]);
#endif
      scalar_mult_add_su3_matrix( fat1,
				 (su3_matrix *)gen_pt[0][i], coef, 
				 fat1 );
    }
  }

  free(tempmat);
  cleanup_gather(mtag0);
}
#endif  /* ASQ_OPTIMIZED_FATTENING   */


#ifdef  ASQ_OPTIMIZED_FORCE
/*   Covariant shift of the src fermion field in the direction dir  *
 *  by one unit. The result is stored in dest.                       */ 
void u_shift_fermion(su3_vector *src, su3_vector *dest, int dir ) {
  su3_vector *tmpvec ; 
  msg_tag *mtag ;
  register site *s ;
  register int i ;
  
  if(GOES_FORWARDS(dir)) /* forward shift */
    {
      mtag = start_gather_from_temp(src, sizeof(su3_vector), 
				    dir, EVENANDODD, gen_pt[0]);
      wait_gather(mtag);
      FORALLSITES(i,s)
	mult_su3_mat_vec(&(s->link[dir]),(su3_vector *)(gen_pt[0][i]),
			 &(dest[i]));
      cleanup_gather(mtag);
    }
  else /* backward shift */
    {
      tmpvec = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));
      FORALLSITES(i,s)
	mult_adj_su3_mat_vec(&(s->link[OPP_DIR(dir)]),&(src[i]), &tmpvec[i]);
      mtag = start_gather_from_temp(tmpvec, sizeof(su3_vector), dir, 
				    EVENANDODD, gen_pt[0]);
      wait_gather(mtag);
      /* copy the gen_pt to the dest */
      FORALLSITES(i,s)
	dest[i] = *(su3_vector *)gen_pt[0][i];
      cleanup_gather(mtag);
      free(tmpvec) ;
    }
}

/*  Covariant shift of the src half wilson fermion field in the  *
 * direction dir   by one unit. The result is stored in dest.    */
void u_shift_hw_fermion(half_wilson_vector *src, 
			half_wilson_vector *dest, int dir ) {
  half_wilson_vector *tmpvec ; 
  msg_tag *mtag ;
  register site *s ;
  register int i ;
  
  if(GOES_FORWARDS(dir)) /* forward shift */
    {
      mtag = start_gather_from_temp(src, sizeof(half_wilson_vector), 
				    dir, EVENANDODD, gen_pt[0]);
      wait_gather(mtag);
      FORALLSITES(i,s)
	mult_su3_mat_hwvec(&(s->link[dir]),
			   (half_wilson_vector *)(gen_pt[0][i]), &(dest[i]));
      cleanup_gather(mtag);
    }
  else /* backward shift */
    {
      tmpvec = 
	(half_wilson_vector *)calloc(sites_on_node, sizeof(half_wilson_vector));
      FORALLSITES(i,s)
	mult_adj_su3_mat_hwvec(&(s->link[OPP_DIR(dir)]),&(src[i]), &tmpvec[i]);
      mtag = start_gather_from_temp(tmpvec, sizeof(half_wilson_vector), dir, 
				    EVENANDODD, gen_pt[0]);
      wait_gather(mtag);
      /* copy the gen_pt to the dest */
      FORALLSITES(i,s)
	dest[i] = *(half_wilson_vector *)gen_pt[0][i];
      cleanup_gather(mtag);
      free(tmpvec) ;
    }
}

/* Add in contribution to the force */
/* Put antihermitian traceless part into momentum */
void add_force_to_mom(su3_vector *back,su3_vector *forw,int dir,double coeff) {
  register site *s ;
  register int i ;  
  register double tmp_coeff ;

  su3_matrix tmat,tmat2;

  if(GOES_BACKWARDS(dir))
    {
      dir = OPP_DIR(dir) ; 
      coeff = -coeff ;
    }
  FORALLSITES(i,s){
    if(s->parity==ODD) 
      tmp_coeff = -coeff ;
    else
      tmp_coeff = coeff ;
    uncompress_anti_hermitian( &(s->mom[dir]), &tmat2 );
    su3_projector(&(back[i]), &(forw[i]), &tmat);
    scalar_mult_add_su3_matrix(&tmat2, &tmat,  tmp_coeff, &tmat2 );
    make_anti_hermitian( &tmat2, &(s->mom[dir]) ); 
  }
}

/* Add in contribution to the force ( 3flavor case ) */
/* Put antihermitian traceless part into momentum */
void add_3f_force_to_mom(half_wilson_vector *back,
			 half_wilson_vector *forw, int dir, double coeff[2]) {
  register site *s ;
  register int i ;  
  double tmp_coeff[2] ;

  su3_matrix tmat,tmat2;

  if(GOES_BACKWARDS(dir))
    {
      dir = OPP_DIR(dir) ; 
      coeff[0] = -coeff[0] ;
      coeff[1] = -coeff[1] ;
    }
  FORALLSITES(i,s){
    if(s->parity==ODD)
      {
	tmp_coeff[0] = -coeff[0] ;
	tmp_coeff[1] = -coeff[1] ;
      }
    else
      {
	tmp_coeff[0] = coeff[0] ;
	tmp_coeff[1] = coeff[1] ;
      }
    uncompress_anti_hermitian( &(s->mom[dir]), &tmat2 );
    su3_projector(&(back[i].h[0]), &(forw[i].h[0]), &tmat);
    scalar_mult_add_su3_matrix(&tmat2, &tmat,  tmp_coeff[0], &tmat2 );
    su3_projector(&(back[i].h[1]), &(forw[i].h[1]), &tmat);
    scalar_mult_add_su3_matrix(&tmat2, &tmat,  tmp_coeff[1], &tmat2 );
    make_anti_hermitian( &tmat2, &(s->mom[dir]) ); 
  }
}

/*  This routine is needed in order to add the force on the side link *
 * of the paths in the Asq and Asqtad actions. It gets as inputs the  *
 * direction mu of the side link and the direction nu of the Dslash   *
 * term we are dealing with. Then it takes also 4 fermion fields:     *
 * Path: the piece of the path with no hop in the nu or mu direction  *
 * Path_nu: the piece of the path with a hop in nu  but not in mu     *
 * Path_mu: is Path times the link mu                                 *
 * Path_numu: is Path_nu times the link mu                            */
void side_link_force(int mu, int nu, double coeff, 
		     su3_vector *Path   , su3_vector *Path_nu, 
		     su3_vector *Path_mu, su3_vector *Path_numu) {
  if(GOES_FORWARDS(mu))
    {
      /*                    nu           * 
       * Add the force :  +----+         *
       *               mu |    |         *
       *                  x    (x)       *
       *                  o    o         */
      if(GOES_FORWARDS(nu))
	add_force_to_mom(Path_numu, Path, mu, coeff ) ;
      else
	add_force_to_mom(Path, Path_numu, OPP_DIR(mu), -coeff );/* ? extra - */
    }
  else /*GOES_BACKWARDS(mu)*/
    {
      /* Add the force :  o    o         *
       *               mu |    |         *
       *                  x    (x)       *
       *                  +----+         *
       *                    nu           */ 
      if(GOES_FORWARDS(nu))
	add_force_to_mom(Path_nu, Path_mu, mu, -coeff) ; /* ? extra - */
      else
	add_force_to_mom(Path_mu, Path_nu, OPP_DIR(mu), coeff) ;
    }
}
 
/*  The 3 flavor version of side_link_force used *
 * to optimize fermion transports                */
void side_link_3f_force(int mu, int nu, double coeff[2], 
			half_wilson_vector *Path   , 
			half_wilson_vector *Path_nu, 
			half_wilson_vector *Path_mu, 
			half_wilson_vector *Path_numu) {
  double m_coeff[2] ;

  m_coeff[0] = -coeff[0] ;
  m_coeff[1] = -coeff[1] ;

  if(GOES_FORWARDS(mu))
    {
      /*                    nu           * 
       * Add the force :  +----+         *
       *               mu |    |         *
       *                  x    (x)       *
       *                  o    o         */
      if(GOES_FORWARDS(nu))
	add_3f_force_to_mom(Path_numu, Path, mu, coeff ) ;
      else
	add_3f_force_to_mom(Path,Path_numu,OPP_DIR(mu),m_coeff);/* ? extra - */
    }
  else /*GOES_BACKWARDS(mu)*/
    {
      /* Add the force :  o    o         *
       *               mu |    |         *
       *                  x    (x)       *
       *                  +----+         *
       *                    nu           */ 
      if(GOES_FORWARDS(nu))
	add_3f_force_to_mom(Path_nu, Path_mu, mu, m_coeff) ; /* ? extra - */
      else
	add_3f_force_to_mom(Path_mu, Path_nu, OPP_DIR(mu), coeff) ;
    }
}
#endif  /* ASQ_OPTIMIZED_FORCE   */

/* LONG COMMENTS
   Here we have combined "xxx", (offset "x_off")  which is
(M_adjoint M)^{-1} phi, with Dslash times this vector, which goes in the
odd sites of xxx.  Recall that phi is defined only on even sites.  In
computing the fermion force, we are looking at

< X |  d/dt ( Dslash_eo Dslash_oe ) | X >
=
< X | d/dt Dslash_eo | T > + < T | d/dt Dslash_oe | X >
where T = Dslash X.

The subsequent manipulations to get the coefficent of H, the momentum
matrix, in the simulation time derivative above look the same for
the two terms, except for a minus sign at the end, if we simply stick
T, which lives on odd sites, into the odd sites of X

 Each path in the action contributes terms when any link of the path
is the link for which we are computing the force.  We get a minus sign
for odd numbered links in the path, since they connect sites of the
opposite parity from what it would be for an even numbered link.
Minus signs from "going around" plaquette - ie KS phases, are supposed
to be already encoded in the path coefficients.
Minus signs from paths that go backwards are supposed to be already
encoded in the path coefficients.

Here, for example, are comments reproduced from the force routine for
the one-link plus Naik plus single-staple-fat-link action:

 The three link force has three contributions, where the link that
was differentiated is the first, second, or third link in the 3-link
path, respectively.  Diagramatically, where "O" represents the momentum,
the solid line the link corresponding to the momentum, and the dashed
lines the other links:
 

	O______________ x ............ x ...............
+
	x..............O______________x.................
+
	x..............x..............O________________
Think of this as
	< xxx | O | UUUxxx >		(  xxx, UUUX_p3 )
+
	< xxx U | O | UUxxx >		( X_m1U , UUX_p2 )
+
	< xxx U U | O | Uxxx >		( X_m2UU , UX_p1 )
where "U" indicates parallel transport, "X_p3" is xxx displaced
by +3, etc.
Note the second contribution has a relative minus sign
because it effectively contributes to the <odd|even>, or M_adjoint,
part of the force when we work on an even site. i.e., for M on
an even site, this three link path begins on an odd site.

The staple force has six contributions from each plane containing the
link direction:
Call these diagrams A-F:


	x...........x		O____________x
		    .			     .
		    .			     .
		    .			     .
		    .			     .
		    .			     .
	O___________x		x............x
	   (A)			    (B)



	x	    x		O____________x
	.	    .		.	     .
	.	    .		.	     .
	.	    .		.	     .
	.	    .		.	     .
	.	    .		.	     .
	O___________x		x	     x
	   (C)			    (D)



	x...........x		O____________x
	.			.
	.			.
	.			.
	.			.
	.			.
	O___________x		x............x
	   (E)			    (F)

As with the Naik term, diagrams C and D have a relative minus
sign because they connect sites of the other parity.

Also note an overall minus sign in the staple terms relative to the
one link term because, with the KS phase factors included, the fat
link is  "U - w3 * UUU", or the straight link MINUS w3 times the staples.

Finally, diagrams B and E get one more minus sign because the link
we are differentiating is in the opposite direction from the staple
as a whole.  You can think of this as this "U" being a correction to
a "U_adjoint", but the derivative of U is iHU and the derivative
of U_adjoint is -iHU_adjoint.

*/
/* LONG COMMENT on sign conventions
In most of the program, the KS phases and antiperiodic boundary
conditions are absorbed into the link matrices.  This greatly simplfies
multiplying by the fermion matrix.  However, it requires care in
specifying the path coefficients.  Remember that each time you
encircle a plaquette, you pick up a net minus sign from the KS phases.
Thus, when you have more than one path to the same point, you generally
have a relative minus sign for each plaquette in a surface bounded by
this path and the basic path for that displacement.

Examples:
  Fat Link:
    Positive:	X-------X

    Negative     --------
	 	|	|
		|	|
		X	X

  Naik connection, smeared
    Positive:	X-------x-------x-------X

    Negative:	---------
		|	|
		|	|
		X	x-------x-------X

    Positive:	--------x--------
		|		|
		|		|
		X		x-------X

    Negative:	--------x-------x-------x
		|			|
		|			|
		X			X
*/



/* Comment on acceptable actions.
   We construct the backwards part of dslash by reversing all the
   paths in the forwards part.  So, for example, in the p4 action
   the forwards part includes +X+Y+Y

		X
		|
		|
		X
		|
		|
	X---->--X

  so we put -X-Y-Y in the backwards part.  But this isn't the adjoint
  of U_x(0)U_y(+x)U_y(+x+y).  Since much of the code assumes that the
  backwards hop is the adjoint of the forwards (for example, in
  preventing going to 8 flavors), the code only works for actions
  where this is true.  Roughly, this means that the fat link must
  be symmetric about reflection around its midpoint.  Equivalently,
  the paths in the backwards part of Dslash are translations of the
  paths in the forwards part.  In the case of the "P4" or knight's move
  action, this means that we have to have both paths
   +X+Y+Y and +Y+Y+X to the same point, with the same coefficients.
  Alternatively, we could just use the symmetric path +Y+X+Y.
*/

invivoquark_stuff_load_fatlinks_566

Real time: 2.371161e+09 -- Predicted time: 2.412550e+09

Callcount: 22

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
14	### 1	10.0038326362174	1.107993e+08	1.072275e+08	3.2
12	### 2	12.0000792662136	1.086769e+08	1.064262e+08	2.0

                /********** update.c ****************************************************/
/* MIMD version 6 */

/*
 Update lattice.
 Improved method for 1-4 flavors:
	update U by (epsilon/2)*(1-Nf/4)
	compute PHI
	update U to epsilon/2
	compute X
	update H, full step
	update U to next time needed

 This routine does not refresh the antihermitian momenta.
 This routine begins at "integral" time, with H and U evaluated
 at same time.
*/
#include "ks_imp_includes.h"	/* definitions files and prototypes */

int update()  {
int step, iters=0;
double final_rsq;
void predict_next_xxx(double *oldtime,double *newtime,double *nexttime);
double cg_time,old_cg_time,next_cg_time;	/* simulation time for last two CG's */
#ifdef HMC_ALGORITHM
double startaction,endaction,d_action();
double xrandom;
#endif

    /* refresh the momenta */
    ranmom();

/*DEBUG*/
/**grsource_imp(F_OFFSET(phi), mass, EVENANDODD);
ks_congrad(F_OFFSET(phi),F_OFFSET(xxx),mass,niter,rsqmin,EVENANDODD,&final_rsq);
checkmul();**/
/*ENDDEBUG*/

    /* do "steps" microcanonical steps  */
    for(step=1; step <= steps; step++){
 
#ifdef PHI_ALGORITHM
        /* generate a pseudofermion configuration only at start*/
     	if(step==1){grsource_imp(F_OFFSET(phi), mass, EVEN); old_cg_time = cg_time = -1.0e6;}

#ifdef HMC_ALGORITHM
        /* find action */
        /* do conjugate gradient to get (Madj M)inverse * phi */
        if(step==1){
            /* do conjugate gradient to get (Madj M)inverse * phi */
	    iters += ks_congrad( F_OFFSET(phi), F_OFFSET(xxx), mass, 
				 niter, rsqmin, EVEN, &final_rsq );
	    cg_time = 0.0;

     	    startaction=d_action();
            /* copy link field to old_link */
	    gauge_field_copy( F_OFFSET(link[0]), F_OFFSET(old_link[0]));
        }
#endif

	/* update U's to middle of interval */
     	update_u(0.5*epsilon);

	/* save conjugate gradient solution, predict next one */
	next_cg_time = ((double)step-0.5)*epsilon;
	predict_next_xxx(&old_cg_time,&cg_time,&next_cg_time);

#else /* "R" algorithm */
       	/* first update the U's to special time interval */
       	update_u(epsilon*(0.5-nflavors/8.0));

        /* generate a pseudofermion configuration */
     	grsource_imp(F_OFFSET(phi), mass, EVEN); cg_time = -1.0e6;

	/* update U's to middle of interval */
     	update_u(epsilon*nflavors/8.0);
#endif

        /* do conjugate gradient to get (Madj M)inverse * phi */
     	iters += ks_congrad( F_OFFSET(phi), F_OFFSET(xxx), mass, 
			     niter, rsqmin, EVEN, &final_rsq );
	cg_time = ((double)step - 0.5)*epsilon;
	/* now update H by full time interval */
    	update_h(epsilon);

    	/* update U's by half time step to get to even time */
    	update_u(epsilon*0.5);

        /* reunitarize the gauge field */
	rephase( OFF );
        reunitarize();
	rephase( ON );

    }	/* end loop over microcanonical steps */

#ifdef HMC_ALGORITHM
    /* find action */
    /* do conjugate gradient to get (Madj M)inverse * phi */
    next_cg_time = steps*epsilon;
    predict_next_xxx(&old_cg_time,&cg_time,&next_cg_time);
    iters += ks_congrad( F_OFFSET(phi), F_OFFSET(xxx), mass,
			 niter, rsqmin, EVEN, &final_rsq );
    cg_time = steps*epsilon;
    endaction=d_action();
    /* decide whether to accept, if not, copy old link field back */
    /* careful - must generate only one random number for whole lattice */
    if(this_node==0)xrandom = myrand(&node_prn);
    broadcast_double(&xrandom);
    if( exp( (double)(startaction-endaction) ) < xrandom ){
	if(steps > 0)
	    gauge_field_copy( F_OFFSET(old_link[0]), F_OFFSET(link[0]) );
#ifdef FN
	valid_longlinks=valid_fatlinks=0;
#endif
	node0_printf("REJECT: delta S = %e\n", (double)(endaction-startaction));
    }
    else {
	node0_printf("ACCEPT: delta S = %e\n", (double)(endaction-startaction));
    }
#endif

    if(steps > 0)return (iters/steps);
    else return(-99);
}

#ifdef PHI_ALGORITHM
/* use linear extrapolation to predict next conjugate gradient solution */
/* only need even sites */
void predict_next_xxx(double *oldtime,double *newtime,double *nexttime)
{
register int i;
register site *s;
register double x;
su3_vector tvec;
    if( *newtime != *oldtime ) x = (*nexttime-*newtime)/(*newtime-*oldtime);
    else x = 0.0;
    if( *oldtime < 0.0 ){
        FOREVENSITES(i,s){
	    s->old_xxx = s->xxx;
        }
    }
    else  {
        FOREVENSITES(i,s){
            sub_su3_vector( &(s->xxx), &(s->old_xxx), &tvec);
	    s->old_xxx = s->xxx;
            scalar_mult_add_su3_vector( &(s->xxx), &tvec,x, &(s->xxx) );
        }
    }
    *oldtime = *newtime;
    *newtime = *nexttime;
}
#endif

invivoupdate_update_40

Real time: 1.279301e+10 -- Predicted time: 1.363889e+10

Callcount: 2

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1	### 1	1	6.607406e+09	6.371209e+09	3.5
2	### 2	1	7.031489e+09	6.446610e+09	8.3

                /****** quark_stuff.c  -- ******************/
/* MIMD version 6 */
/* quark action stuff for improved action
* D.T. 1/28/98, starting from gauge_stuff.c
* K.O. 3/99 Added optimized fattening for Asq actions
* D.T. 4/99 Combine force calculations for both mass quarks
* K.O. 4/99 Optimized force for Asq action
* S.G. 7/01, modified to use t_longlink and t_fatlink
* C.D. 10/02, consolidated quark_stuff.c and quark_stuff_tmp.c
*
* This code combines quark_stuff.c and quark_stuff_tmp.c
* with DSLASH_TMP_LINKS defined, puts links in field-major variables
* t_longlink and t_fatlink. Otherwise, puts them in the site structure.
*
* In this directory, assume all paths connect even to odd sites, etc.
* Tabulate "backwards" paths (e.g. "XDOWN" is backward path to "XUP")
* as separate parity transforms of the fundamental paths.  They will
* generally need a negative sign in Dslash.  See bottom for a long
* comment on sign conventions.
*/

/**#define FFTIME**/
/**#define LLTIME**/

#include "generic_ks_includes.h"	/* definitions files and prototypes */


#define NULL_FP -1 /* NULL field_offset to be used in the optimized version *
                    * of the load_fatlinks subroutine */

    /* Specify paths in orientation in which they appear in the
       forward part of the x component of dslash().  Rotations and
       reflections will be automatically included. Be careful
       about signs of coefficients.  See long comment at bottom. */
#include <quark_action.h>
    /* Include file specifies the basic paths */

#define GOES_FORWARDS(dir) (dir<=TUP)
#define GOES_BACKWARDS(dir) (dir>TUP)

void printpath( int *path, int length );
void path_transport( field_offset src, field_offset dest, int parity,
    int *dir, int length );
void path_transport_hwv( field_offset src, field_offset dest, int parity,
    int *dir, int length );
#ifdef  ASQ_OPTIMIZED_FATTENING
void compute_gen_staple(field_offset staple, int mu, int nu,
                        field_offset link, double coef ) ;
#endif
#ifdef  ASQ_OPTIMIZED_FORCE
#ifndef FN
BOMB THE COMPILE
#endif
void u_shift_fermion(su3_vector *src, su3_vector *dest, int dir ) ;
void add_force_to_mom(su3_vector *back, su3_vector *forw, int dir, double coef);
void side_link_force(int mu, int nu, double coeff, su3_vector *Path,
		     su3_vector *Path_nu, su3_vector *Path_mu, 
		     su3_vector *Path_numu) ;

void u_shift_hw_fermion(half_wilson_vector *src, 
			half_wilson_vector *dest, int dir ) ;
void add_3f_force_to_mom(half_wilson_vector *back,
			 half_wilson_vector *forw, int dir, double coeff[2]) ;
void side_link_3f_force(int mu, int nu, double coeff[2], 
			half_wilson_vector *Path   , 
			half_wilson_vector *Path_nu, 
			half_wilson_vector *Path_mu, 
			half_wilson_vector *Path_numu) ;
#endif

int path_num[MAX_BASIC_PATHS];	/* number of rotations/reflections for each 
					kind */
static double act_path_coeff[MAX_BASIC_PATHS]; /* actual path coefficient     *
                                               * it is equal to path_coeff   *
                                               * if not tadpole improvement  *
                                               * is specified                *
                                               * or path_coeff*u_0^(L-1) when*
                                               * tadpole improvement is      *
                                               * specified                   */

/* Array of structures, for each rotation and reflection of each kind of
	path.  */
struct {
    int dir[MAX_LENGTH];	/* directions in path */
    int length;		/* length of path */
    double coeff;	/* coefficient, including minus sign if backwards */
    double forwback;	/* +1 if in forward Dslash, -1 if in backward */
}   q_paths[MAX_NUM];
int num_q_paths;	/* number of paths in dslash */
int num_basic_paths;	/* number of paths before rotation/reflection */

int is_path_equal( int *path1, int* path2, int length );
int add_basic_path( int *vec, int length, double coeff );

/* Make table of paths in action */
void make_path_table() {

    int i,j;
#ifdef TADPOLE_IMPROVE
    int k;
#endif

    /* table of directions, 1 for each kind of path */
    /**int path_ind[MAX_BASIC_PATHS][MAX_LENGTH];**/
    /* table of coefficients in action, for each path */

    node0_printf("%s\n",quark_action_description);
    num_q_paths = 0;
    num_basic_paths = 0;


    /* add rots. and reflects to table, print out the path coefficients */
    node0_printf("path coefficients: npath  path_coeff  multiplicity\n");
    for(j=0;j<quark_action_npaths;j++) {
	double this_coeff;
	this_coeff = path_coeff[j];
#ifdef TADPOLE_IMPROVE
	for(k=1;k< path_length_in[j];k++)this_coeff /= u0;
#endif
	act_path_coeff[j] = this_coeff ;
	i = add_basic_path( path_ind[j], path_length_in[j],
	    this_coeff );
	node0_printf("                    %d      %e     %d\n",
	    j,this_coeff,i);
    }
}

/* add rotations and reflections of a path to the table.  Return
   multiplicity of paths added */
int add_basic_path( int *basic_vec, int length, double coeff ) {

    int perm[8],pp[8],ir[4];
    int j,path_num;
    int vec[MAX_LENGTH];
    int flag;

    path_num = 0;
    /* now fill the long table with all rotations and reflections
	of the fundamental path.  The path presented to us is for
        the positive x component of dslash, so if the x coordinate
        is reflected it will appear with a negative sign. */
      /* permutations */
      for(perm[0]=0;perm[0]<4;perm[0]++)
      for(perm[1]=0;perm[1]<4;perm[1]++)
      for(perm[2]=0;perm[2]<4;perm[2]++)
      for(perm[3]=0;perm[3]<4;perm[3]++){
	if(perm[0] != perm[1] && perm[0] != perm[2] 
	  && perm[0] != perm[3] && perm[1] != perm[2]
	  && perm[1] != perm[3] && perm[2] != perm[3] ) {
	  /* reflections*/
	  for(ir[0]=0;ir[0]<2;ir[0]++)
	  for(ir[1]=0;ir[1]<2;ir[1]++)
	  for(ir[2]=0;ir[2]<2;ir[2]++)
	  for(ir[3]=0;ir[3]<2;ir[3]++){
	    for(j=0;j<4;j++){
	      pp[j]=perm[j];

	      if(ir[j] == 1) pp[j]=OPP_DIR(pp[j]);
	      pp[OPP_DIR(j)]=OPP_DIR(pp[j]);
	    }
	    /* create new vector*/
	    for(j=0;j<length;j++) vec[j]=pp[basic_vec[j]];
	    for(j=length;j<MAX_LENGTH;j++) vec[j]=NODIR;

            flag=0;
	    /* check if it's a new set: */
	    for(j=0;j<num_q_paths;j++){
	      flag = is_path_equal( vec, q_paths[j].dir, MAX_LENGTH );
	      if(flag==1)break;
	    }
	    if(flag == 0 ){
	      if(num_q_paths>=MAX_NUM){
		node0_printf("OOPS: MAX_NUM too small\n");
		exit(0);
	      }
	      q_paths[num_q_paths].length=length;
	      for(j=0;j<MAX_LENGTH;j++) q_paths[num_q_paths].dir[j]=vec[j];
		/* remember to copy NODIR's, or comparison will get confused */
	      if(ir[0]==0){
		q_paths[num_q_paths].coeff =  coeff;
		q_paths[num_q_paths].forwback =  +1;
	      }
	      else{
		q_paths[num_q_paths].coeff = -coeff;
		q_paths[num_q_paths].forwback = -1;
	      }
	      num_q_paths++;
	      path_num++;
	      /**node0_printf("ADD PATH %d:  rx=%d ",num_q_paths-1,ir[0]);
		 printpath( vec, length );**/
	    }

	  } /* end reflection*/
        } /* end permutation if block */
      } /* end permutation */
    num_basic_paths++;
    return(path_num);
} /* add_basic_path */

/* parallel transport a vector to the current site along a path.
   For example, if the path is "XUP", bring in the vector from
   the +X direction to the current site. If KS phases are in lattice,
   this transport will automatically include them. 
   OK for src and dest to be the same.  OK for length=0.  */
/* NOT OPTIMIZED at the moment - do lots of extra copying.  Use temp.
   vectors rather than stuff in lattice.h */
void path_transport( field_offset src, field_offset dest, int parity,
    int *dir, int length ){
    register int i;
    register site *s;
    msg_tag *mtag0;
    int j;
    su3_vector *tmp_src,*tmp_dest,*tmp_work; /*source, dest and workspace*/
    su3_vector *tmp_pt; /* scratch */
    int tmp_parity, tmp_otherparity; /* parity for this step */

  if( length > 0 ){
    tmp_src = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));
    tmp_dest = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));
    tmp_work = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));

    for( j=length-1; j>=0; j-- ){
	/* figure out parities for this step */
	if( j%2==0 ){
	    tmp_parity = parity;
	    switch(tmp_parity){
		case EVEN: tmp_otherparity=ODD; break;
		case ODD: tmp_otherparity=EVEN; break;
		case EVENANDODD: tmp_otherparity=EVENANDODD; break;
	    }
	}
	else { /* odd # step */
	    tmp_otherparity = parity;
	    switch(tmp_otherparity){
		case EVEN: tmp_parity=ODD; break;
		case ODD: tmp_parity=EVEN; break;
		case EVENANDODD: tmp_parity=EVENANDODD; break;
	    }
	}

	if( j==length-1 ){
	    FORSOMEPARITY(i,s,tmp_otherparity){
	        tmp_src[i] = *(su3_vector *)F_PT(s,src);
	    }
	}

	if( GOES_FORWARDS(dir[j]) ) {
	    mtag0 = start_gather_from_temp( tmp_src, sizeof(su3_vector),
	        dir[j], tmp_parity, gen_pt[0] );
	    wait_gather(mtag0);
	    FORSOMEPARITY(i,s,tmp_parity){
		mult_su3_mat_vec( &(s->link[dir[j]]),
		    (su3_vector *)(gen_pt[0][i]),
		    &(tmp_dest[i]) );
	    }
	    cleanup_gather(mtag0);
	}

	else{ /* GOES_BACKWARDS(dir[j]) */
	    FORSOMEPARITY(i,s,tmp_otherparity){
		mult_adj_su3_mat_vec( &(s->link[OPP_DIR(dir[j])]),
		    &(tmp_src[i]), &(tmp_work[i]) );
	    }
	    mtag0 = start_gather_from_temp( tmp_work, sizeof(su3_vector),
	        dir[j], tmp_parity, gen_pt[0] );
	    wait_gather(mtag0);
	    FORSOMEPARITY(i,s,tmp_parity){
		 tmp_dest[i] = *(su3_vector *)gen_pt[0][i];
	    }
	    cleanup_gather(mtag0);
	}
	
	/* src for next step is dest for this one. */
	tmp_pt=tmp_src; tmp_src=tmp_dest; tmp_dest=tmp_pt;
    }  /* j=link in path */
    /* done, copy result into real dest. (tmp_src now points to result) */
    FORSOMEPARITY(i,s,parity){
        *(su3_vector *)F_PT(s,dest) = tmp_src[i];
    }
    free(tmp_src); free(tmp_dest); free(tmp_work);
  } /* end if(length>0) */
  else if( src != dest ){ /* for length=0 */
    FORSOMEPARITY(i,s,parity){
        *(su3_vector *)F_PT(s,dest) = *(su3_vector *)F_PT(s,src);
    }
  }
} /* path_transport */

/* Path transport a half_wilson_vector */
void path_transport_hwv( field_offset src, field_offset dest, int parity,
    int *dir, int length ){
    register int i;
    register site *s;
    msg_tag *mtag0;
    int j;
    half_wilson_vector *tmp_src,*tmp_dest,*tmp_work; /*source, dest and workspace*/
    half_wilson_vector *tmp_pt; /* scratch */
    int tmp_parity, tmp_otherparity; /* parity for this step */

  if( length > 0 ){
    tmp_src = (half_wilson_vector *)calloc(
      sites_on_node, sizeof(half_wilson_vector));
    tmp_dest = (half_wilson_vector *)calloc(
       sites_on_node, sizeof(half_wilson_vector));
    tmp_work = (half_wilson_vector *)calloc(
       sites_on_node, sizeof(half_wilson_vector));

    for( j=length-1; j>=0; j-- ){
	/* figure out parities for this step */
	if( j%2==0 ){
	    tmp_parity = parity;
	    switch(tmp_parity){
		case EVEN: tmp_otherparity=ODD; break;
		case ODD: tmp_otherparity=EVEN; break;
		case EVENANDODD: tmp_otherparity=EVENANDODD; break;
	    }
	}
	else { /* odd # step */
	    tmp_otherparity = parity;
	    switch(tmp_otherparity){
		case EVEN: tmp_parity=ODD; break;
		case ODD: tmp_parity=EVEN; break;
		case EVENANDODD: tmp_parity=EVENANDODD; break;
	    }
	}

	if( j==length-1 ){
	    FORSOMEPARITY(i,s,tmp_otherparity){
	        tmp_src[i] = *(half_wilson_vector *)F_PT(s,src);
	    }
	}

	if( GOES_FORWARDS(dir[j]) ) {
	    mtag0 = start_gather_from_temp( tmp_src,
		sizeof(half_wilson_vector), dir[j], tmp_parity, gen_pt[0] );
	    wait_gather(mtag0);
	    FORSOMEPARITY(i,s,tmp_parity){
		mult_su3_mat_hwvec( &(s->link[dir[j]]),
		    (half_wilson_vector *)(gen_pt[0][i]),
		    &(tmp_dest[i]) );
	    }
	    cleanup_gather(mtag0);
	}

	else{ /* GOES_BACKWARDS(dir[j]) */
	    FORSOMEPARITY(i,s,tmp_otherparity){
		mult_adj_su3_mat_hwvec( &(s->link[OPP_DIR(dir[j])]),
		    &(tmp_src[i]), &(tmp_work[i]) );
	    }
	    mtag0 = start_gather_from_temp( tmp_work,
		sizeof(half_wilson_vector), dir[j], tmp_parity, gen_pt[0] );
	    wait_gather(mtag0);
	    FORSOMEPARITY(i,s,tmp_parity){
		 tmp_dest[i] = *(half_wilson_vector *)gen_pt[0][i];
	    }
	    cleanup_gather(mtag0);
	}
	
	/* src for next step is dest for this one. */
	tmp_pt=tmp_src; tmp_src=tmp_dest; tmp_dest=tmp_pt;
    }  /* j=link in path */
    /* done, copy result into real dest. (tmp_src now points to result) */
    FORSOMEPARITY(i,s,parity){
        *(half_wilson_vector *)F_PT(s,dest) = tmp_src[i];
    }
    free(tmp_src); free(tmp_dest); free(tmp_work);
  } /* end if(length>0) */
  else if( src != dest ){ /* for length=0 */
    FORSOMEPARITY(i,s,parity){
        *(half_wilson_vector *)F_PT(s,dest) =
	  *(half_wilson_vector *)F_PT(s,src);
    }
  }
} /* path_transport_hwv */


#ifdef EO
/* Stupid dslash routine that follows all the paths.
  Should optimize by precomputing sums of paths to all the displacements
  Use dslash_fn for actions that involve only +X and +X+X+X couplings.
*/
/* D_slash routine - sets dest. on each site equal to sum of
   sources parallel transported to site, with minus sign for transport
   from negative directions.  */
void dslash_eo( field_offset src, field_offset dest, int parity ) {
   register int i;
   register site *s;
   register int ipath,otherparity;
   register double x; /* coefficient of path */

    switch(parity){
	case EVEN:	otherparity=ODD; break;
	case ODD:	otherparity=EVEN; break;
	case EVENANDODD:	otherparity=EVENANDODD; break;
    }

    /* Parallel transport by all the paths in the action.  
       Multiply by coefficient in table
    */
    FORSOMEPARITY(i,s,parity){ clearvec( (su3_vector *)F_PT(s,dest) ); }

    for( ipath=0; ipath<num_q_paths; ipath++ ){  /* loop over paths */
	    path_transport( src, F_OFFSET(tempvec[0]), parity,
		q_paths[ipath].dir, q_paths[ipath].length );
	    x=q_paths[ipath].coeff;
	    FORSOMEPARITY(i,s,parity){
		scalar_mult_add_su3_vector(  (su3_vector *)F_PT(s,dest),
		    &(s->tempvec[0]), x, (su3_vector *)F_PT(s,dest) );
	    }
    }   /* ipath */
} /* dslash_eo */
#endif /* ifdef EO */

#ifdef FN
/* Sum over paths connecting to nearest neighbor point (fat link) and to third
   nearest neighbor (longlinks) */
/* Doug Toussaint 2/4/98 */
/* modified to use t_longlinks, S. Gottlieb 7/13/01 */
/* long link calculating routine */
/* path_product() follows the path starting at step 0, and
   leaves the answer at the end of the path.  We want the answer
   at the site where the path begins.  So we look for paths with
   the opposite displacement from the displacement of the point
   that we want to transport to this site, and take the adjoint
   of the matrix at the end. clear? */
/* KS phases and APBC must be in the links. See long comment at bottom*/
void load_longlinks() {
  register int i;
  register site *s;
  int ipath,dir;
  int disp[4];
  register su3_matrix *long1;

#ifdef LLTIME
double dtime;
dtime=-dclock();
#endif
  if( phases_in != 1){
    node0_printf("BOTCH: load_longlinks needs phases in\n");
    terminate(0);
  }
  for (dir=XUP; dir<=TUP; dir++){ /* loop over longlink directions */
    /* set longlink to zero */
    FORALLSITES(i,s){
#ifdef DSLASH_TMP_LINKS
      long1 = &(t_longlink[4*i+dir]);
#else
      long1 = &(s->longlink[dir]);
#endif
      clear_su3mat( long1 );
    }

    /* loop over paths, checking for ones with total displacement 3*dir */
    for( ipath=0; ipath<num_q_paths; ipath++ ){  /* loop over paths */
	/* compute total displacement of path */
	for(i=XUP;i<=TUP;i++)disp[i]=0;
	for( i=0; i<q_paths[ipath].length; i++){
	  if( GOES_FORWARDS(q_paths[ipath].dir[i]) )
	    disp[        q_paths[ipath].dir[i]  ]++;
	  else
	    disp[OPP_DIR(q_paths[ipath].dir[i]) ]--;
	}
	for( disp[dir]+=3,i=XUP; i<=TUP; i++)if(disp[i]!=0)break;
	if( i<=TUP )continue;  /* skip if path doesn't go to right place */
/**printf("ipath = %d, found a path:  ",ipath);
for(j=0;j<q_paths[ipath].length;j++)printf("\t%d", q_paths[ipath].dir[j]);
printf("\n");**/

	path_product( q_paths[ipath].dir, q_paths[ipath].length );
	FORALLSITES(i,s){
	  su3_adjoint( &(s->tempmat1), &(s->staple) );
#ifdef DSLASH_TMP_LINKS
	  long1 = &(t_longlink[4*i+dir]);
#else
	  long1 = &(s->longlink[dir]);
#endif
          scalar_mult_add_su3_matrix( long1,
	    &(s->staple), -q_paths[ipath].coeff, long1 );
		/* minus sign in coeff. because we used backward path*/
	}
    } /* ipath */

  } /* loop over directions */

  valid_longlinks = 1;
#ifdef LLTIME
dtime += dclock();
node0_printf("LLTIME(Long):  %e\n",dtime);
#endif
}  /* load_longlinks() */

/* KS phases and APBC must be in the links. See long comment at bottom*/
void load_fatlinks() {
  register int i;
  register site *s;
  int dir;
  register su3_matrix *fat1;
#ifdef ASQ_OPTIMIZED_FATTENING
  int  nu,rho,sig ;
  double one_link ; /*needed to fix the problem with the Lepage term */
#else
  int ipath;
  int disp[4];
#endif

#ifdef LLTIME
double dtime;
dtime=-dclock();
#endif
  if( phases_in != 1){
    node0_printf("BOTCH: load_fatlinks needs phases in\n");
    terminate(0);
  }

#ifndef  ASQ_OPTIMIZED_FATTENING   /* general case code */
  for (dir=XUP; dir<=TUP; dir++){ /* loop over fatlink directions */
    /* set fatlink to zero */
    FORALLSITES(i,s){
#ifdef DSLASH_TMP_LINKS
      fat1 = &(t_fatlink[4*i+dir]);
#else
      fat1 = &(s->fatlink[dir]);
#endif
      clear_su3mat( fat1 );
    }
    
    /* loop over paths, checking for ones with total displacement 1*dir */
    for( ipath=0; ipath<num_q_paths; ipath++ ){  /* loop over paths */
	/* compute total displacement of path */
	for(i=XUP;i<=TUP;i++)disp[i]=0;
	for( i=0; i<q_paths[ipath].length; i++){
	  if( GOES_FORWARDS(q_paths[ipath].dir[i]) )
	    disp[        q_paths[ipath].dir[i]  ]++;
	  else
	    disp[OPP_DIR(q_paths[ipath].dir[i]) ]--;
	}
	for( disp[dir]+=1,i=XUP; i<=TUP; i++)if(disp[i]!=0)break;
	if( i<=TUP )continue;  /* skip if path doesn't go to right place */
/**printf("dir = %d, found a path:  ",dir);
for(j=0;j<q_paths.[ipath].length;j++)printf("\t%d", q_paths[ipath].dir[j]);
printf("\n");**/

	path_product( q_paths[ipath].dir, q_paths[ipath].length );
	FORALLSITES(i,s){
	  su3_adjoint( &(s->tempmat1), &(s->staple) );
#ifdef DSLASH_TMP_LINKS
	  fat1 = &(t_fatlink[4*i+dir]);
#else
	  fat1 = &(s->fatlink[dir]);
#endif
          scalar_mult_add_su3_matrix( fat1,
	    &(s->staple), -q_paths[ipath].coeff, fat1 );
		/* minus sign in coeff. because we used backward path*/
	}
    } /* ipath */
  } /* loop over directions */
#else	/* ASQ_OPTIMIZED_FATTENING, for Asq and Asqtad actions */
/*  Optimized fattening code for the Asq and Asqtad actions.           *
 * I assume that path 0 is the one link path 2 the 3-staple            *
 * path 3 the 5-staple path 4 the 7-staple and path 5 the Lapage term. *
 * Path 1 is the Naik term.                                            */
 
 /* to fix up the Lepage term, included by a trick below */
 one_link = (act_path_coeff[0] - 6.0*act_path_coeff[5]) ; 
 
 for (dir=XUP; dir<=TUP; dir++){
   FORALLSITES(i,s) /* Intialize fat links with c_1*U_\mu(x) */
     {
#ifdef DSLASH_TMP_LINKS
       fat1 = &(t_fatlink[4*i+dir]);
#else
       fat1 = &(s->fatlink[dir]);
#endif
       scalar_mult_su3_matrix(&(s->link[dir]), one_link ,
			      fat1 );
     }
   for(nu=XUP; nu<=TUP; nu++) if(nu!=dir)
     {
       compute_gen_staple(F_OFFSET(staple),dir,nu,F_OFFSET(link[dir]),
			  act_path_coeff[2]);
       /* The Lepage term */
       /* Note this also involves modifying c_1 (above) */
       compute_gen_staple(NULL_FP,dir,nu,F_OFFSET(staple),act_path_coeff[5]);
       for(rho=XUP; rho<=TUP; rho++) if((rho!=dir)&&(rho!=nu))
	 {
	   compute_gen_staple(F_OFFSET(tempmat1),dir,rho,F_OFFSET(staple),
			      act_path_coeff[3]);
	   for(sig=XUP; sig<=TUP; sig++)
	     if((sig!=dir)&&(sig!=nu)&&(sig!=rho))
	       {
		 compute_gen_staple(NULL_FP,dir,sig,
				    F_OFFSET(tempmat1),
				    act_path_coeff[4]);
	       } /* sig */
	 } /* rho */
     } /* nu */
 }/* dir */  
#endif

  valid_fatlinks = 1;
#ifdef LLTIME
dtime += dclock();
node0_printf("LLTIME(Fat):  %e\n",dtime);
#endif
}  /* load_fatlinks() */
#endif /* ifdef FN */


/* compare two paths, return 1 if equal, else zero */
int is_path_equal( int *path1, int* path2, int length ){
   register int i;
   for(i=0;i<length;i++)if(path1[i]!=path2[i])return(0);
   return(1);
}


/* update the  momenta with the fermion force */
/* Assumes that the conjugate gradient has been run, with the answer in
   x_off, and dslash(x_off,x_off,ODD) has been run. (fills in x_off_odd) */
/* SEE LONG COMMENTS AT END */

#ifndef ASQ_OPTIMIZED_FORCE
void eo_fermion_force( double eps, int nflavors, field_offset x_off ){
  /* note CG_solution and Dslash * solution are combined in "x_off" */
  /* New version 1/21/99.  Use forward part of Dslash to get force */
  /* see long comment at end */
  /* For each link we need x_off transported from both ends of path. */
  register int i,dir,lastdir,ipath,ilink;
  register site *s;
  int length;
  su3_matrix tmat,tmat2;
  double ferm_epsilon, coeff;
#ifdef FFTIME
  double dtime;
#endif
  msg_tag *mtag0;
  half_wilson_vector *hw_tmp0,*hw_tmp1,*tmp_pt;

#ifdef FFTIME
dtime=-dclock();
#endif
  ferm_epsilon = 2.0*(nflavors/4.0)*eps;
  hw_tmp0 = (half_wilson_vector *)
    calloc(sites_on_node, sizeof(half_wilson_vector));
  hw_tmp1 = (half_wilson_vector *)
    calloc(sites_on_node, sizeof(half_wilson_vector));
  /* Use half_wilson_vectors to store x_off transported from ends of
     path.  0 component from forward end, 1 component from back end */

  /* loop over paths, and loop over links in path */
  for( ipath=0; ipath<num_q_paths; ipath++ ){
    if(q_paths[ipath].forwback== -1)continue;	/* skip backwards dslash */
    length = q_paths[ipath].length;

    /* path transport x_off and Dslash*x_off from far end.  Sometimes
	we need them at the start point of the path, and sometimes
	one link into the path - an optimization for later */
    path_transport( x_off, F_OFFSET(tempvec[0]),
      EVENANDODD, q_paths[ipath].dir, length );
    /* use tempvec[1] for transport from starting end */
    FORALLSITES(i,s){
      hw_tmp0[i].h[0]=s->tempvec[0];
      hw_tmp0[i].h[1]=*(su3_vector *)F_PT(s,x_off);
    }

    /* A path has (length+1) points, counting the ends.  At first
	 point, no "down" direction links have their momenta "at this
	 point". At last, no "up" ... */
    for( ilink=0; ilink<=length; ilink++ ){
      if(ilink<length)dir = q_paths[ipath].dir[ilink];
      else dir=NODIR;
      coeff = ferm_epsilon*q_paths[ipath].coeff;
      if( (ilink%2)==1 )coeff = -coeff;

      /* path transport x_off and Dslash*x_off from previous point */
      /* Use "half_wilson_vector" to handle pair of vectors -
	0 component is x_off from forward end, 1 component from back end */
      /* sometimes we don't need them */
      if( (ilink>0&&ilink<length) || 
	(ilink==length && GOES_BACKWARDS(lastdir)) ){

	if( GOES_FORWARDS(lastdir) ){
	  FORALLSITES(i,s){
            mult_adj_su3_mat_hwvec( &(s->link[lastdir]),
	      &(hw_tmp0[i]), &(hw_tmp1[i]) );
	  }
	  mtag0 = start_gather_from_temp( hw_tmp1, 2*sizeof(su3_vector),
             OPP_DIR(lastdir), EVENANDODD, gen_pt[0] );
          wait_gather(mtag0);
          FORALLSITES(i,s){
	     hw_tmp0[i] = *(half_wilson_vector *)gen_pt[0][i];
	  }
          cleanup_gather(mtag0);
	}
	else{   /* GOES_BACKWARDS(lastdir) */
          mtag0 = start_gather_from_temp( hw_tmp0, 2*sizeof(su3_vector),
               OPP_DIR(lastdir), EVENANDODD, gen_pt[0] );
          wait_gather(mtag0);
          FORALLSITES(i,s){
            mult_su3_mat_hwvec( &(s->link[OPP_DIR(lastdir)]),
                (half_wilson_vector *)(gen_pt[0][i]),
		&(hw_tmp1[i]) );
          }
	  tmp_pt = hw_tmp0; hw_tmp0 = hw_tmp1; hw_tmp1 = tmp_pt;
          cleanup_gather(mtag0);
	}
      }

      /* add in contribution to the force */
      /* Put antihermitian traceless part into momentum */
      FORALLSITES(i,s){
        if( ilink<length && GOES_FORWARDS(dir) ){
          uncompress_anti_hermitian( &(s->mom[dir]), &tmat2 );
	  su3_projector( &(hw_tmp0[i].h[0]), &(hw_tmp0[i].h[1]), &tmat );
	  if( s->parity==EVEN )scalar_mult_add_su3_matrix(
	    &tmat2, &tmat,  coeff, &tmat2 );
	  else		 scalar_mult_add_su3_matrix(
	    &tmat2, &tmat, -coeff, &tmat2 );
          make_anti_hermitian( &tmat2, &(s->mom[dir]) );
        }
	if( ilink>0 && GOES_BACKWARDS(lastdir) ){
          uncompress_anti_hermitian( &(s->mom[OPP_DIR(lastdir)]), &tmat2 );
	  su3_projector( &(hw_tmp0[i].h[0]), &(hw_tmp0[i].h[1]), &tmat );
	  if( s->parity==EVEN )scalar_mult_add_su3_matrix(
	    &tmat2, &tmat, -coeff, &tmat2 );
	  else		 scalar_mult_add_su3_matrix(
	    &tmat2, &tmat,  coeff, &tmat2 );
          make_anti_hermitian( &tmat2, &(s->mom[OPP_DIR(lastdir)]) );
        }
      }
      lastdir = dir;
    } /* end loop over links in path */
  } /* end loop over paths */

  free( hw_tmp0 ); free( hw_tmp1 );

#ifdef FFTIME
dtime += dclock();
node0_printf("FFTIME:  %e\n",dtime);
/**printf("TLENGTH: %d\n",tlength);**/
#endif
} /* eo_fermion_force(version 6) */


void eo_fermion_force_3f( double eps, int nflav1, field_offset x1_off,
  int nflav2, field_offset x2_off ){
  /* note CG_solution and Dslash * solution are combined in "x_off" */
  /* New version 1/21/99.  Use forward part of Dslash to get force */
  /* 4/14/99 combine force from two different mass quarks, (eg 2+1flavors) */
  /* see long comment at end */
  /* For each link we need x_off transported from both ends of path. */
  register int i,dir,lastdir,ipath,ilink;
  register site *s;
  int length;
  su3_matrix tmat,tmat2;
  double ferm_epsilon1, ferm_epsilon2, coeff1, coeff2;
#ifdef FFTIME
  double dtime;
#endif
  msg_tag *mtag0;
  wilson_vector *w_tmp0,*w_tmp1,*tmp_pt;

#ifdef FFTIME
dtime=-dclock();
#endif
  ferm_epsilon1 = 2.0*(nflav1/4.0)*eps;
  ferm_epsilon2 = 2.0*(nflav2/4.0)*eps;
  w_tmp0 = (wilson_vector *)
    calloc(sites_on_node, sizeof(wilson_vector));
  w_tmp1 = (wilson_vector *)
    calloc(sites_on_node, sizeof(wilson_vector));
  /* Use wilson_vectors to store x_off transported from ends of
     path.  0 and 1 components from forward end, 2 and 3 components
     from back end */

  /* loop over paths, and loop over links in path */
  for( ipath=0; ipath<num_q_paths; ipath++ ){
    if(q_paths[ipath].forwback== -1)continue;	/* skip backwards dslash */
    length = q_paths[ipath].length;

    /* path transport x_off and Dslash*x_off from far end.  Sometimes
	we need them at the start point of the path, and sometimes
	one link into the path - an optimization for later */
    /**
    path_transport( x1_off, F_OFFSET(tempvec[0]),
      EVENANDODD, q_paths[ipath].dir, length );
    path_transport( x2_off, F_OFFSET(tempvec[1]),
      EVENANDODD, q_paths[ipath].dir, length );
    **/
/** WARNING!! Assumes xxx1 and xxx2 contiguous **/
if( x2_off-x1_off != sizeof(su3_vector) ){node0_printf("BOTCH\n"); exit(0);}
    path_transport_hwv( x1_off, F_OFFSET(tempvec[0]),
      EVENANDODD, q_paths[ipath].dir, length );
    /* use tempvec[2] for transport from starting end */
    FORALLSITES(i,s){
      w_tmp0[i].d[0]=s->tempvec[0];
      w_tmp0[i].d[1]=s->tempvec[1];
      w_tmp0[i].d[2]=*(su3_vector *)F_PT(s,x1_off);
      w_tmp0[i].d[3]=*(su3_vector *)F_PT(s,x2_off);
    }

    /* A path has (length+1) points, counting the ends.  At first
	 point, no "down" direction links have their momenta "at this
	 point". At last, no "up" ... */
    for( ilink=0; ilink<=length; ilink++ ){
      if(ilink<length)dir = q_paths[ipath].dir[ilink];
      else dir=NODIR;
      coeff1 = ferm_epsilon1*q_paths[ipath].coeff;
      coeff2 = ferm_epsilon2*q_paths[ipath].coeff;
      if( (ilink%2)==1 ){ coeff1 = -coeff1; coeff2 = -coeff2;}

      /* path transport x_off and Dslash*x_off from previous point */
      /* Use "wilson_vector" to handle pair of vectors -
	0 component is x_off1 from forward end, 1 component is x_off2
	from forward end,  2 and 3  components are x_off1 and x_off2 
	from back end */
      /* sometimes we don't need them */
      if( (ilink>0&&ilink<length) || 
	(ilink==length && GOES_BACKWARDS(lastdir)) ){

	if( GOES_FORWARDS(lastdir) ){
	  FORALLSITES(i,s){
            mult_adj_mat_wilson_vec( &(s->link[lastdir]),
	      &(w_tmp0[i]), &(w_tmp1[i]) );
	  }
	  mtag0 = start_gather_from_temp( w_tmp1, 4*sizeof(su3_vector),
             OPP_DIR(lastdir), EVENANDODD, gen_pt[0] );
          wait_gather(mtag0);
          FORALLSITES(i,s){
	     w_tmp0[i] = *(wilson_vector *)gen_pt[0][i];
	  }
          cleanup_gather(mtag0);
	}
	else{   /* GOES_BACKWARDS(lastdir) */
          mtag0 = start_gather_from_temp( w_tmp0, 4*sizeof(su3_vector),
               OPP_DIR(lastdir), EVENANDODD, gen_pt[0] );
          wait_gather(mtag0);
          FORALLSITES(i,s){
            mult_mat_wilson_vec( &(s->link[OPP_DIR(lastdir)]),
                (wilson_vector *)(gen_pt[0][i]),
		&(w_tmp1[i]) );
          }
	  tmp_pt = w_tmp0; w_tmp0 = w_tmp1; w_tmp1 = tmp_pt;
          cleanup_gather(mtag0);
	}
      }

      /* add in contribution to the force */
      /* Put antihermitian traceless part into momentum */
      FORALLSITES(i,s){
        if( ilink<length && GOES_FORWARDS(dir) ){
	  if( s->parity==ODD ){coeff1 *= -1.0; coeff2 *= -1.0; }
          uncompress_anti_hermitian( &(s->mom[dir]), &tmat2 );
	  su3_projector( &(w_tmp0[i].d[0]), &(w_tmp0[i].d[2]), &tmat );
	  scalar_mult_add_su3_matrix( &tmat2, &tmat,  coeff1, &tmat2 );
	  su3_projector( &(w_tmp0[i].d[1]), &(w_tmp0[i].d[3]), &tmat );
	  scalar_mult_add_su3_matrix( &tmat2, &tmat,  coeff2, &tmat2 );
          make_anti_hermitian( &tmat2, &(s->mom[dir]) );
	  if( s->parity==ODD ){coeff1 *= -1.0; coeff2 *= -1.0; }
        }
	if( ilink>0 && GOES_BACKWARDS(lastdir) ){
	  if( s->parity==EVEN ){coeff1 *= -1.0; coeff2 *= -1.0; }
          uncompress_anti_hermitian( &(s->mom[OPP_DIR(lastdir)]), &tmat2 );
	  su3_projector( &(w_tmp0[i].d[0]), &(w_tmp0[i].d[2]), &tmat );
	  scalar_mult_add_su3_matrix( &tmat2, &tmat,  coeff1, &tmat2 );
	  su3_projector( &(w_tmp0[i].d[1]), &(w_tmp0[i].d[3]), &tmat );
	  scalar_mult_add_su3_matrix( &tmat2, &tmat,  coeff2, &tmat2 );
          make_anti_hermitian( &tmat2, &(s->mom[OPP_DIR(lastdir)]) );
	  if( s->parity==EVEN ){coeff1 *= -1.0; coeff2 *= -1.0; }
        }
      }
      lastdir = dir;
    } /* end loop over links in path */
  } /* end loop over paths */

  free( w_tmp0 ); free( w_tmp1 );

#ifdef FFTIME
dtime += dclock();
node0_printf("FFTIME:  %e\n",dtime);
/**printf("TLENGTH: %d\n",tlength);**/
#endif
} /* eo_fermion_force_3f(version 6) */



#else /* ASQ_OPTIMIZED_FORCE, for Asq and Asqtad actions */
/* Optimized force code for the Asq and Asqtad actions                 *
 * I assume that path 0 is the one link path 2 the 3-staple            *
 * path 3 the 5-staple path 4 the 7-staple and path 5 the Lapage term. *
 * Path 1 is the Naik term.                                            */
#define Pmu          tempvec[0] 
#define Pnumu        tempvec[1]
#define Prhonumu     tempvec[2]
#define P7           tempvec[3]
#define P7rho        tempvec[4]              
#define P7rhonu      tempvec[5]
#define P5           tempvec[6]
#define P3           tempvec[7]
#define P5nu         tempvec[3]
#define P3mu         tempvec[3]
#define Popmu        tempvec[4]
#define Pmumumu      tempvec[4]
void eo_fermion_force( double eps, int nflavors, field_offset x_off ){
  /* note CG_solution and Dslash * solution are combined in "x_off" */
  /* New version 1/21/99.  Use forward part of Dslash to get force */
  /* see long comment at end */
  /* For each link we need x_off transported from both ends of path. */
  register int i ;
  register site *s;
  int mu,nu,rho,sig ;
  int DirectLinks[8] ;
  double ferm_epsilon, coeff;
  double OneLink, Lepage, Naik, FiveSt, ThreeSt, SevenSt ;
  su3_vector *tempvec[8] ;
  su3_vector *temp_x ;

#ifdef FFTIME
  double dtime;

  dtime=-dclock();
#endif
  ferm_epsilon = 2.0*(nflavors/4.0)*eps;
  
  /* Path coefficients times fermion epsilon */
  OneLink = act_path_coeff[0]*ferm_epsilon ; 
  Naik    = act_path_coeff[1]*ferm_epsilon ;
  ThreeSt = act_path_coeff[2]*ferm_epsilon ;
  FiveSt  = act_path_coeff[3]*ferm_epsilon ;
  SevenSt = act_path_coeff[4]*ferm_epsilon ;
  Lepage  = act_path_coeff[5]*ferm_epsilon ;
  /* *************************************** */

  /* Initialize the DirectLink flags */
  for(mu=0;mu<8;mu++)
    DirectLinks[mu] = 0 ;

  /* Allocate temporary vectors */
  for(mu=0;mu<8;mu++)
    tempvec[mu] = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));

  /*copy x_off to a temporary vector */
  temp_x = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));
  FORALLSITES(i,s) temp_x[i] = *(su3_vector *)F_PT(s,x_off) ;

  for(sig=0;sig<8;sig++)
    {
      for(mu=0;mu<8;mu++)if((mu!=sig)&&(mu!=OPP_DIR(sig)))
	{
	  u_shift_fermion(temp_x, Pmu, OPP_DIR(mu));
	  u_shift_fermion(Pmu, P3, sig);
	  if(GOES_FORWARDS(sig))
	    {
	      /* Add the force F_sig[x+mu]:         x--+             *
	       *                                   |   |             *
	       *                                   o   o             *
	       * the 1 link in the path: - (numbering starts form 0) */
	      add_force_to_mom(P3, Pmu, sig, -ThreeSt) ;
	    }
	  for(nu=0;nu<8;nu++)if((nu!=mu )&&(nu!=OPP_DIR(mu ))&&
				(nu!=sig)&&(nu!=OPP_DIR(sig)))
	    {
	      u_shift_fermion(Pmu, Pnumu, OPP_DIR(nu));
	      u_shift_fermion(Pnumu, P5, sig);
	      if(GOES_FORWARDS(sig))
		{
		  /* Add the force F_sig[x+mu+nu]:      x--+             *
		   *                                   |   |             *
		   *                                   o   o             *
		   * the 2 link in the path: + (numbering starts form 0) */
		  add_force_to_mom(P5, Pnumu, sig, FiveSt);
		}
	      for(rho=0;rho<8;rho++)if((rho!=mu )&&(rho!=OPP_DIR(mu ))&&
				       (rho!=nu )&&(rho!=OPP_DIR(nu ))&&
				       (rho!=sig)&&(rho!=OPP_DIR(sig)))
		{
		  u_shift_fermion(Pnumu, Prhonumu, OPP_DIR(rho));
		  /* Length 7 paths */
		  u_shift_fermion(Prhonumu, P7,sig);
		  if(GOES_FORWARDS(sig))
		    {
		      /* Add the force F_sig[x+mu+nu+rho]:  x--+             *
		       *                                   |   |             *
		       *                                   o   o             *
		       * the 3 link in the path: - (numbering starts form 0) */
		      add_force_to_mom(P7, Prhonumu, sig, -SevenSt ) ;
		    }
		  /*Add the force F_rho the 2(4) link in the path: +     */
		  u_shift_fermion(P7, P7rho, rho);
		  side_link_force(rho,sig,SevenSt, Pnumu, P7, Prhonumu, P7rho);
		  /* Add the P7rho vector to P5 */
		  coeff = SevenSt/FiveSt ; 
		  FORALLSITES(i,s)
		    scalar_mult_add_su3_vector(&(P5[i]),&(P7rho[i]),coeff,
					       &(P5[i]));
		}/* rho */
	      /* Length 5 paths */
	      /*Add the force F_nu the 1(3) link in the path: -     */
	      u_shift_fermion(P5,P5nu, nu);
	      side_link_force(nu,sig,-FiveSt,Pmu,P5, 
			      Pnumu,P5nu) ;
	      /* Add the P5nu vector to P3 */
	      coeff = FiveSt/ThreeSt ; 
	      FORALLSITES(i,s)
		scalar_mult_add_su3_vector(&(P3[i]),&(P5nu[i]),coeff,&(P3[i]));
	    }/* nu */

	  /* Now the Lepage term... It is the same with 5-link paths with
             nu=mu and FiveSt=Lepage. So Pnumu is really Pmumu */
	  u_shift_fermion(Pmu, Pnumu, OPP_DIR(mu));
	  u_shift_fermion(Pnumu, P5, sig);
	  if(GOES_FORWARDS(sig))
	    {
	      /* Add the force F_sig[x+mu+nu]:      x--+             *
	       *                                   |   |             *
	       *                                   o   o             *
	       * the 2 link in the path: + (numbering starts form 0) */
	      add_force_to_mom(P5, Pnumu, sig, Lepage) ;
	    }
	  /*Add the force F_nu the 1(3) link in the path: -     */
	  u_shift_fermion(P5,P5nu, mu);
	  side_link_force(mu, sig, -Lepage, Pmu, P5, Pnumu, P5nu) ;
	  /* Add the P5nu vector to P3 */
	  coeff = Lepage/ThreeSt ; 
	  FORALLSITES(i,s)
	    scalar_mult_add_su3_vector(&(P3[i]),&(P5nu[i]),coeff,&(P3[i]));

	  /* Length 3 paths (Not the Naik term) */
	  /*Add the force F_mu the 0(2) link in the path: +     */
	  if(GOES_FORWARDS(mu)) 
	    u_shift_fermion(P3,P3mu, mu );
	  /* The above shift is not needed if mu is backwards */
	  side_link_force(mu, sig, ThreeSt, temp_x, P3, Pmu, P3mu);

	  /* Finally the OneLink and the Naik term */
	  /* Check if this direction is not already done */
	  if( (!DirectLinks[mu]) ){
	    DirectLinks[mu]=1 ;
	    if(GOES_BACKWARDS(mu))/* Do only the forward terms in the Dslash */
	      {
		/* Because I have shifted with OPP_DIR(mu) Pmu is a forward *
		 * shift.                                                   */
		/* The one link */
		add_force_to_mom(Pmu, temp_x, OPP_DIR(mu), OneLink) ;
		/* For the same reason Pnumu is the forward double link */

		/* Popmu is a backward shift */
		u_shift_fermion(temp_x, Popmu, mu);
		/* The Naik */
		/* link no 1: - */
		add_force_to_mom(Pnumu, Popmu, OPP_DIR(mu), -Naik) ;
		/*Pmumumu can overwrite Popmu which is no longer needed */
		u_shift_fermion(Pnumu, Pmumumu, OPP_DIR(mu));
		/* link no 0: + */
		add_force_to_mom(Pmumumu, temp_x, OPP_DIR(mu), Naik);
	      }
	    else /* The rest of the Naik terms */
	      {
		u_shift_fermion(temp_x, Popmu, mu);
		/* link no 2: + */
		/* Pnumu is double backward shift */
		add_force_to_mom(Popmu, Pnumu, mu, Naik) ;
	      }
	  }
	}/* mu */
      /* Here we have to do together the Naik term and the one link term */
    }/*sig */

  /* Free temporary vectors */
  free(temp_x) ;
  for(mu=0;mu<8;mu++)
    free(tempvec[mu]) ;

#ifdef FFTIME
  dtime += dclock();
node0_printf("FFTIME:  %e\n",dtime);
/**printf("TLENGTH: %d\n",tlength);**/
#endif
} /* eo_fermion_force(version 6) */
#undef Pmu          
#undef Pnumu        
#undef Prhonumu     
#undef P7           
#undef P7rho        
#undef P7rhonu      
#undef P5           
#undef P3           
#undef P5nu         
#undef P3mu         
#undef Popmu        
#undef Pmumumu      

#define Pmu          hwvec[0] 
#define Pnumu        hwvec[1]
#define Prhonumu     hwvec[2]
#define P7           hwvec[3]
#define P7rho        hwvec[4]              
#define P7rhonu      hwvec[5]
#define P5           hwvec[6]
#define P3           hwvec[7]
#define P5nu         hwvec[3]
#define P3mu         hwvec[3]
#define Popmu        hwvec[4]
#define Pmumumu      hwvec[4]
void eo_fermion_force_3f( double eps, int nflav1, field_offset x1_off, 
			  int nflav2, field_offset x2_off){
  /* note CG_solution and Dslash * solution are combined in "x_off" */
  /* New version 1/21/99.  Use forward part of Dslash to get force */
  /* 4/15/99 combine force from two different mass quarks, (eg 2+1flavors) */
  /* see long comment at end */
  /* For each link we need x_off transported from both ends of path. */
  register int i ;
  register site *s;
  int mu,nu,rho,sig ;
  int DirectLinks[8] ;
  double coeff[2],ferm_epsilon ;
  double OneLink[2], Lepage[2], Naik[2], FiveSt[2], ThreeSt[2], SevenSt[2] ;
  double mNaik[2], mLepage[2], mFiveSt[2], mThreeSt[2], mSevenSt[2] ;
  half_wilson_vector *hwvec[8] ;
  half_wilson_vector *temp_x ;

#ifdef FFTIME
  double dtime;

  dtime=-dclock();
#endif
  
  /* Path coefficients times fermion epsilon */
  ferm_epsilon = 2.0*(nflav1/4.0)*eps;
  OneLink[0] = act_path_coeff[0]*ferm_epsilon ;
  Naik[0]    = act_path_coeff[1]*ferm_epsilon ; mNaik[0]    = -Naik[0] ;
  ThreeSt[0] = act_path_coeff[2]*ferm_epsilon ; mThreeSt[0] = -ThreeSt[0] ;
  FiveSt[0]  = act_path_coeff[3]*ferm_epsilon ; mFiveSt[0]  = -FiveSt[0] ;
  SevenSt[0] = act_path_coeff[4]*ferm_epsilon ; mSevenSt[0] = -SevenSt[0] ;
  Lepage[0]  = act_path_coeff[5]*ferm_epsilon ; mLepage[0]  = -Lepage[0] ;

  ferm_epsilon = 2.0*(nflav2/4.0)*eps;
  OneLink[1] = act_path_coeff[0]*ferm_epsilon ;
  Naik[1]    = act_path_coeff[1]*ferm_epsilon ; mNaik[1]    = -Naik[1] ;
  ThreeSt[1] = act_path_coeff[2]*ferm_epsilon ; mThreeSt[1] = -ThreeSt[1] ;
  FiveSt[1]  = act_path_coeff[3]*ferm_epsilon ; mFiveSt[1]  = -FiveSt[1] ;
  SevenSt[1] = act_path_coeff[4]*ferm_epsilon ; mSevenSt[1] = -SevenSt[1] ;
  Lepage[1]  = act_path_coeff[5]*ferm_epsilon ; mLepage[1]  = -Lepage[1] ;
  /* *************************************** */

  /* Initialize the DirectLink flags */
  for(mu=0;mu<8;mu++)
    DirectLinks[mu] = 0 ;

  /* Allocate temporary vectors */
  for(mu=0;mu<8;mu++)
    hwvec[mu]= 
      (half_wilson_vector *)calloc(sites_on_node, sizeof(half_wilson_vector));

  /*copy x_off to a temporary vector */
  temp_x= 
    (half_wilson_vector *)calloc(sites_on_node, sizeof(half_wilson_vector));
  FORALLSITES(i,s)
    {
      temp_x[i].h[0] = *(su3_vector *)F_PT(s,x1_off) ;
      temp_x[i].h[1] = *(su3_vector *)F_PT(s,x2_off) ;
    }

  for(sig=0;sig<8;sig++)
    {
      for(mu=0;mu<8;mu++)if((mu!=sig)&&(mu!=OPP_DIR(sig)))
	{
	  u_shift_hw_fermion(temp_x, Pmu, OPP_DIR(mu));
	  u_shift_hw_fermion(Pmu, P3, sig);
	  if(GOES_FORWARDS(sig))
	    {
	      /* Add the force F_sig[x+mu]:         x--+             *
	       *                                   |   |             *
	       *                                   o   o             *
	       * the 1 link in the path: - (numbering starts form 0) */
	      add_3f_force_to_mom(P3, Pmu, sig, mThreeSt) ;
	    }
	  for(nu=0;nu<8;nu++)if((nu!=mu )&&(nu!=OPP_DIR(mu ))&&
				(nu!=sig)&&(nu!=OPP_DIR(sig)))
	    {
	      u_shift_hw_fermion(Pmu, Pnumu, OPP_DIR(nu));
	      u_shift_hw_fermion(Pnumu, P5, sig);
	      if(GOES_FORWARDS(sig))
		{
		  /* Add the force F_sig[x+mu+nu]:      x--+             *
		   *                                   |   |             *
		   *                                   o   o             *
		   * the 2 link in the path: + (numbering starts form 0) */
		  add_3f_force_to_mom(P5, Pnumu, sig, FiveSt);
		}
	      for(rho=0;rho<8;rho++)if((rho!=mu )&&(rho!=OPP_DIR(mu ))&&
				       (rho!=nu )&&(rho!=OPP_DIR(nu ))&&
				       (rho!=sig)&&(rho!=OPP_DIR(sig)))
		{
		  u_shift_hw_fermion(Pnumu, Prhonumu, OPP_DIR(rho));
		  /* Length 7 paths */
		  u_shift_hw_fermion(Prhonumu, P7,sig);
		  if(GOES_FORWARDS(sig))
		    {
		      /* Add the force F_sig[x+mu+nu+rho]:  x--+             *
		       *                                   |   |             *
		       *                                   o   o             *
		       * the 3 link in the path: - (numbering starts form 0) */
		      add_3f_force_to_mom(P7, Prhonumu, sig, mSevenSt ) ;
		    }
		  /*Add the force F_rho the 2(4) link in the path: +     */
		  u_shift_hw_fermion(P7, P7rho, rho);
		  side_link_3f_force(rho,sig,SevenSt,Pnumu,P7,Prhonumu,P7rho);
		  /* Add the P7rho vector to P5 */
		  coeff[0] = SevenSt[0]/FiveSt[0] ; 
		  coeff[1] = SevenSt[1]/FiveSt[1] ; 
		  FORALLSITES(i,s)
		    {
		      scalar_mult_add_su3_vector(&(P5[i].h[0]),
						 &(P7rho[i].h[0]),coeff[0],
						 &(P5[i].h[0]));
		      scalar_mult_add_su3_vector(&(P5[i].h[1]),
						 &(P7rho[i].h[1]),coeff[1],
						 &(P5[i].h[1]));
		    }
		}/* rho */
	      /* Length 5 paths */
	      /*Add the force F_nu the 1(3) link in the path: -     */
	      u_shift_hw_fermion(P5,P5nu, nu);
	      side_link_3f_force(nu,sig,mFiveSt,Pmu,P5, 
			      Pnumu,P5nu) ;
	      /* Add the P5nu vector to P3 */
	      coeff[0] = FiveSt[0]/ThreeSt[0] ; 
	      coeff[1] = FiveSt[1]/ThreeSt[1] ; 
	      FORALLSITES(i,s)
		{
		  scalar_mult_add_su3_vector(&(P3[i].h[0]),
					     &(P5nu[i].h[0]), coeff[0],
					     &(P3[i].h[0]));
		  scalar_mult_add_su3_vector(&(P3[i].h[1]),
					     &(P5nu[i].h[1]), coeff[1],
					     &(P3[i].h[1]));
		}
	    }/* nu */

	  /* Now the Lepage term... It is the same with 5-link paths with
             nu=mu and FiveSt=Lepage. So Pnumu is really Pmumu */
	  u_shift_hw_fermion(Pmu, Pnumu, OPP_DIR(mu));
	  u_shift_hw_fermion(Pnumu, P5, sig);
	  if(GOES_FORWARDS(sig))
	    {
	      /* Add the force F_sig[x+mu+nu]:      x--+             *
	       *                                   |   |             *
	       *                                   o   o             *
	       * the 2 link in the path: + (numbering starts form 0) */
	      add_3f_force_to_mom(P5, Pnumu, sig, Lepage) ;
	    }
	  /*Add the force F_nu the 1(3) link in the path: -     */
	  u_shift_hw_fermion(P5,P5nu, mu);
	  side_link_3f_force(mu, sig, mLepage, Pmu, P5, Pnumu, P5nu) ;
	  /* Add the P5nu vector to P3 */
	  coeff[0] = Lepage[0]/ThreeSt[0] ; 
	  coeff[1] = Lepage[1]/ThreeSt[1] ; 
	  FORALLSITES(i,s)
	    {
	      scalar_mult_add_su3_vector(&(P3[i].h[0]),
					 &(P5nu[i].h[0]),coeff[0],
					 &(P3[i].h[0]));
	      scalar_mult_add_su3_vector(&(P3[i].h[1]),
					 &(P5nu[i].h[1]),coeff[1],
					 &(P3[i].h[1]));
	    }

	  /* Length 3 paths (Not the Naik term) */
	  /*Add the force F_mu the 0(2) link in the path: +     */
	  if(GOES_FORWARDS(mu)) 
	    u_shift_hw_fermion(P3,P3mu, mu );
	  /* The above shift is not needed if mu is backwards */
	  side_link_3f_force(mu, sig, ThreeSt, temp_x, P3, Pmu, P3mu);

	  /* Finally the OneLink and the Naik term */
	  /* Check if this direction is not already done */
	  if( (!DirectLinks[mu]) ){
	    DirectLinks[mu]=1 ;
	    if(GOES_BACKWARDS(mu))/* Do only the forward terms in the Dslash */
	      {
		/* Because I have shifted with OPP_DIR(mu) Pmu is a forward *
		 * shift.                                                   */
		/* The one link */
		add_3f_force_to_mom(Pmu, temp_x, OPP_DIR(mu), OneLink) ;
		/* For the same reason Pnumu is the forward double link */

		/* Popmu is a backward shift */
		u_shift_hw_fermion(temp_x, Popmu, mu);
		/* The Naik */
		/* link no 1: - */
		add_3f_force_to_mom(Pnumu, Popmu, OPP_DIR(mu), mNaik) ;
		/*Pmumumu can overwrite Popmu which is no longer needed */
		u_shift_hw_fermion(Pnumu, Pmumumu, OPP_DIR(mu));
		/* link no 0: + */
		add_3f_force_to_mom(Pmumumu, temp_x, OPP_DIR(mu), Naik);
	      }
	    else /* The rest of the Naik terms */
	      {
		u_shift_hw_fermion(temp_x, Popmu, mu);
		/* link no 2: + */
		/* Pnumu is double backward shift */
		add_3f_force_to_mom(Popmu, Pnumu, mu, Naik) ;
	      }
	  }
	}/* mu */
      /* Here we have to do together the Naik term and the one link term */
    }/*sig */

  /* Free temporary vectors */
  free(temp_x) ;
  for(mu=0;mu<8;mu++)
    free(hwvec[mu]) ;

#ifdef FFTIME
  dtime += dclock();
node0_printf("FFTIME:  %e\n",dtime);
/**printf("TLENGTH: %d\n",tlength);**/
#endif
} /* eo_fermion_force_3f(version 6) */
#undef Pmu          
#undef Pnumu        
#undef Prhonumu     
#undef P7           
#undef P7rho        
#undef P7rhonu      
#undef P5           
#undef P3           
#undef P5nu         
#undef P3mu         
#undef Popmu        
#undef Pmumumu      
#endif /* ASQ_OPTIMIZED_FORCE */


#ifdef  ASQ_OPTIMIZED_FATTENING   /* Asqtad action only, "_fn" executables */
#ifndef FN
BOMB THE COMPILE
#endif
void compute_gen_staple(field_offset staple, int mu, int nu, 
			field_offset link, double coef) {
  su3_matrix tmat1,tmat2;
  msg_tag *mtag0,*mtag1;
  su3_matrix *tempmat ;
  register site *s ;
  register int i ;
  register su3_matrix *fat1;

  /* Computes the staple :
                   mu
                +-------+
          nu	|	|
		|	|
		X	X
    Where the mu link can be any su3_matrix. The result is saved in staple.
    if staple==NULL_FP then the result is not saved.
    It also adds the computed staple to the fatlink[mu] with weight coef.
  */

  /* Upper staple */
  mtag0 = start_gather( link, sizeof(su3_matrix), nu, EVENANDODD, gen_pt[0] );
  mtag1 = start_gather( F_OFFSET(link[nu]), sizeof(su3_matrix), mu, 
			EVENANDODD, gen_pt[1] );
  wait_gather(mtag0);
  wait_gather(mtag1);
  
  if(staple!=NULL_FP){/* Save the staple */
    FORALLSITES(i,s){
      mult_su3_na( (su3_matrix *)gen_pt[0][i],
		   (su3_matrix *)gen_pt[1][i], &tmat1 );
      mult_su3_nn( &(s->link[nu]), &tmat1, (su3_matrix *)F_PT(s,staple) );
    }
  }
  else{ /* No need to save the staple. Add it to the fatlinks */
    FORALLSITES(i,s){
      mult_su3_na( (su3_matrix *)gen_pt[0][i],
		   (su3_matrix *)gen_pt[1][i], &tmat1 );
      mult_su3_nn( &(s->link[nu]), &tmat1, &tmat2 );
#ifdef DSLASH_TMP_LINKS
      fat1 = &(t_fatlink[4*i+mu]);
#else
      fat1 = &(s->fatlink[mu]);
#endif
      scalar_mult_add_su3_matrix(fat1, &tmat2, coef,
				 fat1) ;
    }
  }
  cleanup_gather(mtag0);
  cleanup_gather(mtag1);

  /* lower staple */
  tempmat = (su3_matrix *)calloc(sites_on_node, sizeof(su3_matrix));
  mtag0 = start_gather( F_OFFSET(link[nu]),
			sizeof(su3_matrix), mu, EVENANDODD, gen_pt[0] );
  wait_gather(mtag0);
  FORALLSITES(i,s){
    mult_su3_an( &(s->link[nu]),(su3_matrix *)F_PT(s,link), &tmat1 );
    mult_su3_nn( &(tmat1),(su3_matrix *)gen_pt[0][i], &(tempmat[i]) );
  }
  cleanup_gather(mtag0);
  mtag0 = start_gather_from_temp( tempmat, sizeof(su3_matrix),
				  OPP_DIR(nu), EVENANDODD, gen_pt[0] );
  wait_gather(mtag0);

  if(staple!=NULL_FP){/* Save the staple */
    FORALLSITES(i,s){
      add_su3_matrix( (su3_matrix *)F_PT(s,staple),(su3_matrix *)gen_pt[0][i], 
		      (su3_matrix *)F_PT(s,staple) );
#ifdef DSLASH_TMP_LINKS
      fat1 = &(t_fatlink[4*i+mu]);
#else
      fat1 = &(s->fatlink[mu]);
#endif
      scalar_mult_add_su3_matrix( fat1,
				 (su3_matrix *)F_PT(s,staple), coef, 
				 fat1 );
    }
  }
  else{ /* No need to save the staple. Add it to the fatlinks */
    FORALLSITES(i,s){
#ifdef DSLASH_TMP_LINKS
      fat1 = &(t_fatlink[4*i+mu]);
#else
      fat1 = &(s->fatlink[mu]);
#endif
      scalar_mult_add_su3_matrix( fat1,
				 (su3_matrix *)gen_pt[0][i], coef, 
				 fat1 );
    }
  }

  free(tempmat);
  cleanup_gather(mtag0);
}
#endif  /* ASQ_OPTIMIZED_FATTENING   */


#ifdef  ASQ_OPTIMIZED_FORCE
/*   Covariant shift of the src fermion field in the direction dir  *
 *  by one unit. The result is stored in dest.                       */ 
void u_shift_fermion(su3_vector *src, su3_vector *dest, int dir ) {
  su3_vector *tmpvec ; 
  msg_tag *mtag ;
  register site *s ;
  register int i ;
  
  if(GOES_FORWARDS(dir)) /* forward shift */
    {
      mtag = start_gather_from_temp(src, sizeof(su3_vector), 
				    dir, EVENANDODD, gen_pt[0]);
      wait_gather(mtag);
      FORALLSITES(i,s)
	mult_su3_mat_vec(&(s->link[dir]),(su3_vector *)(gen_pt[0][i]),
			 &(dest[i]));
      cleanup_gather(mtag);
    }
  else /* backward shift */
    {
      tmpvec = (su3_vector *)calloc(sites_on_node, sizeof(su3_vector));
      FORALLSITES(i,s)
	mult_adj_su3_mat_vec(&(s->link[OPP_DIR(dir)]),&(src[i]), &tmpvec[i]);
      mtag = start_gather_from_temp(tmpvec, sizeof(su3_vector), dir, 
				    EVENANDODD, gen_pt[0]);
      wait_gather(mtag);
      /* copy the gen_pt to the dest */
      FORALLSITES(i,s)
	dest[i] = *(su3_vector *)gen_pt[0][i];
      cleanup_gather(mtag);
      free(tmpvec) ;
    }
}

/*  Covariant shift of the src half wilson fermion field in the  *
 * direction dir   by one unit. The result is stored in dest.    */
void u_shift_hw_fermion(half_wilson_vector *src, 
			half_wilson_vector *dest, int dir ) {
  half_wilson_vector *tmpvec ; 
  msg_tag *mtag ;
  register site *s ;
  register int i ;
  
  if(GOES_FORWARDS(dir)) /* forward shift */
    {
      mtag = start_gather_from_temp(src, sizeof(half_wilson_vector), 
				    dir, EVENANDODD, gen_pt[0]);
      wait_gather(mtag);
      FORALLSITES(i,s)
	mult_su3_mat_hwvec(&(s->link[dir]),
			   (half_wilson_vector *)(gen_pt[0][i]), &(dest[i]));
      cleanup_gather(mtag);
    }
  else /* backward shift */
    {
      tmpvec = 
	(half_wilson_vector *)calloc(sites_on_node, sizeof(half_wilson_vector));
      FORALLSITES(i,s)
	mult_adj_su3_mat_hwvec(&(s->link[OPP_DIR(dir)]),&(src[i]), &tmpvec[i]);
      mtag = start_gather_from_temp(tmpvec, sizeof(half_wilson_vector), dir, 
				    EVENANDODD, gen_pt[0]);
      wait_gather(mtag);
      /* copy the gen_pt to the dest */
      FORALLSITES(i,s)
	dest[i] = *(half_wilson_vector *)gen_pt[0][i];
      cleanup_gather(mtag);
      free(tmpvec) ;
    }
}

/* Add in contribution to the force */
/* Put antihermitian traceless part into momentum */
void add_force_to_mom(su3_vector *back,su3_vector *forw,int dir,double coeff) {
  register site *s ;
  register int i ;  
  register double tmp_coeff ;

  su3_matrix tmat,tmat2;

  if(GOES_BACKWARDS(dir))
    {
      dir = OPP_DIR(dir) ; 
      coeff = -coeff ;
    }
  FORALLSITES(i,s){
    if(s->parity==ODD) 
      tmp_coeff = -coeff ;
    else
      tmp_coeff = coeff ;
    uncompress_anti_hermitian( &(s->mom[dir]), &tmat2 );
    su3_projector(&(back[i]), &(forw[i]), &tmat);
    scalar_mult_add_su3_matrix(&tmat2, &tmat,  tmp_coeff, &tmat2 );
    make_anti_hermitian( &tmat2, &(s->mom[dir]) ); 
  }
}

/* Add in contribution to the force ( 3flavor case ) */
/* Put antihermitian traceless part into momentum */
void add_3f_force_to_mom(half_wilson_vector *back,
			 half_wilson_vector *forw, int dir, double coeff[2]) {
  register site *s ;
  register int i ;  
  double tmp_coeff[2] ;

  su3_matrix tmat,tmat2;

  if(GOES_BACKWARDS(dir))
    {
      dir = OPP_DIR(dir) ; 
      coeff[0] = -coeff[0] ;
      coeff[1] = -coeff[1] ;
    }
  FORALLSITES(i,s){
    if(s->parity==ODD)
      {
	tmp_coeff[0] = -coeff[0] ;
	tmp_coeff[1] = -coeff[1] ;
      }
    else
      {
	tmp_coeff[0] = coeff[0] ;
	tmp_coeff[1] = coeff[1] ;
      }
    uncompress_anti_hermitian( &(s->mom[dir]), &tmat2 );
    su3_projector(&(back[i].h[0]), &(forw[i].h[0]), &tmat);
    scalar_mult_add_su3_matrix(&tmat2, &tmat,  tmp_coeff[0], &tmat2 );
    su3_projector(&(back[i].h[1]), &(forw[i].h[1]), &tmat);
    scalar_mult_add_su3_matrix(&tmat2, &tmat,  tmp_coeff[1], &tmat2 );
    make_anti_hermitian( &tmat2, &(s->mom[dir]) ); 
  }
}

/*  This routine is needed in order to add the force on the side link *
 * of the paths in the Asq and Asqtad actions. It gets as inputs the  *
 * direction mu of the side link and the direction nu of the Dslash   *
 * term we are dealing with. Then it takes also 4 fermion fields:     *
 * Path: the piece of the path with no hop in the nu or mu direction  *
 * Path_nu: the piece of the path with a hop in nu  but not in mu     *
 * Path_mu: is Path times the link mu                                 *
 * Path_numu: is Path_nu times the link mu                            */
void side_link_force(int mu, int nu, double coeff, 
		     su3_vector *Path   , su3_vector *Path_nu, 
		     su3_vector *Path_mu, su3_vector *Path_numu) {
  if(GOES_FORWARDS(mu))
    {
      /*                    nu           * 
       * Add the force :  +----+         *
       *               mu |    |         *
       *                  x    (x)       *
       *                  o    o         */
      if(GOES_FORWARDS(nu))
	add_force_to_mom(Path_numu, Path, mu, coeff ) ;
      else
	add_force_to_mom(Path, Path_numu, OPP_DIR(mu), -coeff );/* ? extra - */
    }
  else /*GOES_BACKWARDS(mu)*/
    {
      /* Add the force :  o    o         *
       *               mu |    |         *
       *                  x    (x)       *
       *                  +----+         *
       *                    nu           */ 
      if(GOES_FORWARDS(nu))
	add_force_to_mom(Path_nu, Path_mu, mu, -coeff) ; /* ? extra - */
      else
	add_force_to_mom(Path_mu, Path_nu, OPP_DIR(mu), coeff) ;
    }
}
 
/*  The 3 flavor version of side_link_force used *
 * to optimize fermion transports                */
void side_link_3f_force(int mu, int nu, double coeff[2], 
			half_wilson_vector *Path   , 
			half_wilson_vector *Path_nu, 
			half_wilson_vector *Path_mu, 
			half_wilson_vector *Path_numu) {
  double m_coeff[2] ;

  m_coeff[0] = -coeff[0] ;
  m_coeff[1] = -coeff[1] ;

  if(GOES_FORWARDS(mu))
    {
      /*                    nu           * 
       * Add the force :  +----+         *
       *               mu |    |         *
       *                  x    (x)       *
       *                  o    o         */
      if(GOES_FORWARDS(nu))
	add_3f_force_to_mom(Path_numu, Path, mu, coeff ) ;
      else
	add_3f_force_to_mom(Path,Path_numu,OPP_DIR(mu),m_coeff);/* ? extra - */
    }
  else /*GOES_BACKWARDS(mu)*/
    {
      /* Add the force :  o    o         *
       *               mu |    |         *
       *                  x    (x)       *
       *                  +----+         *
       *                    nu           */ 
      if(GOES_FORWARDS(nu))
	add_3f_force_to_mom(Path_nu, Path_mu, mu, m_coeff) ; /* ? extra - */
      else
	add_3f_force_to_mom(Path_mu, Path_nu, OPP_DIR(mu), coeff) ;
    }
}
#endif  /* ASQ_OPTIMIZED_FORCE   */

/* LONG COMMENTS
   Here we have combined "xxx", (offset "x_off")  which is
(M_adjoint M)^{-1} phi, with Dslash times this vector, which goes in the
odd sites of xxx.  Recall that phi is defined only on even sites.  In
computing the fermion force, we are looking at

< X |  d/dt ( Dslash_eo Dslash_oe ) | X >
=
< X | d/dt Dslash_eo | T > + < T | d/dt Dslash_oe | X >
where T = Dslash X.

The subsequent manipulations to get the coefficent of H, the momentum
matrix, in the simulation time derivative above look the same for
the two terms, except for a minus sign at the end, if we simply stick
T, which lives on odd sites, into the odd sites of X

 Each path in the action contributes terms when any link of the path
is the link for which we are computing the force.  We get a minus sign
for odd numbered links in the path, since they connect sites of the
opposite parity from what it would be for an even numbered link.
Minus signs from "going around" plaquette - ie KS phases, are supposed
to be already encoded in the path coefficients.
Minus signs from paths that go backwards are supposed to be already
encoded in the path coefficients.

Here, for example, are comments reproduced from the force routine for
the one-link plus Naik plus single-staple-fat-link action:

 The three link force has three contributions, where the link that
was differentiated is the first, second, or third link in the 3-link
path, respectively.  Diagramatically, where "O" represents the momentum,
the solid line the link corresponding to the momentum, and the dashed
lines the other links:
 

	O______________ x ............ x ...............
+
	x..............O______________x.................
+
	x..............x..............O________________
Think of this as
	< xxx | O | UUUxxx >		(  xxx, UUUX_p3 )
+
	< xxx U | O | UUxxx >		( X_m1U , UUX_p2 )
+
	< xxx U U | O | Uxxx >		( X_m2UU , UX_p1 )
where "U" indicates parallel transport, "X_p3" is xxx displaced
by +3, etc.
Note the second contribution has a relative minus sign
because it effectively contributes to the <odd|even>, or M_adjoint,
part of the force when we work on an even site. i.e., for M on
an even site, this three link path begins on an odd site.

The staple force has six contributions from each plane containing the
link direction:
Call these diagrams A-F:


	x...........x		O____________x
		    .			     .
		    .			     .
		    .			     .
		    .			     .
		    .			     .
	O___________x		x............x
	   (A)			    (B)



	x	    x		O____________x
	.	    .		.	     .
	.	    .		.	     .
	.	    .		.	     .
	.	    .		.	     .
	.	    .		.	     .
	O___________x		x	     x
	   (C)			    (D)



	x...........x		O____________x
	.			.
	.			.
	.			.
	.			.
	.			.
	O___________x		x............x
	   (E)			    (F)

As with the Naik term, diagrams C and D have a relative minus
sign because they connect sites of the other parity.

Also note an overall minus sign in the staple terms relative to the
one link term because, with the KS phase factors included, the fat
link is  "U - w3 * UUU", or the straight link MINUS w3 times the staples.

Finally, diagrams B and E get one more minus sign because the link
we are differentiating is in the opposite direction from the staple
as a whole.  You can think of this as this "U" being a correction to
a "U_adjoint", but the derivative of U is iHU and the derivative
of U_adjoint is -iHU_adjoint.

*/
/* LONG COMMENT on sign conventions
In most of the program, the KS phases and antiperiodic boundary
conditions are absorbed into the link matrices.  This greatly simplfies
multiplying by the fermion matrix.  However, it requires care in
specifying the path coefficients.  Remember that each time you
encircle a plaquette, you pick up a net minus sign from the KS phases.
Thus, when you have more than one path to the same point, you generally
have a relative minus sign for each plaquette in a surface bounded by
this path and the basic path for that displacement.

Examples:
  Fat Link:
    Positive:	X-------X

    Negative     --------
	 	|	|
		|	|
		X	X

  Naik connection, smeared
    Positive:	X-------x-------x-------X

    Negative:	---------
		|	|
		|	|
		X	x-------x-------X

    Positive:	--------x--------
		|		|
		|		|
		X		x-------X

    Negative:	--------x-------x-------x
		|			|
		|			|
		X			X
*/



/* Comment on acceptable actions.
   We construct the backwards part of dslash by reversing all the
   paths in the forwards part.  So, for example, in the p4 action
   the forwards part includes +X+Y+Y

		X
		|
		|
		X
		|
		|
	X---->--X

  so we put -X-Y-Y in the backwards part.  But this isn't the adjoint
  of U_x(0)U_y(+x)U_y(+x+y).  Since much of the code assumes that the
  backwards hop is the adjoint of the forwards (for example, in
  preventing going to 8 flavors), the code only works for actions
  where this is true.  Roughly, this means that the fat link must
  be symmetric about reflection around its midpoint.  Equivalently,
  the paths in the backwards part of Dslash are translations of the
  paths in the forwards part.  In the case of the "P4" or knight's move
  action, this means that we have to have both paths
   +X+Y+Y and +Y+Y+X to the same point, with the same coefficients.
  Alternatively, we could just use the symmetric path +Y+X+Y.
*/

invivoquark_stuff_compute_gen_staple_1406

Real time: 1.792900e+08 -- Predicted time: 0.000000e+00

Callcount: 792

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
145	### 1	66.4787883366067	0.000000e+00	2.647600e+05	100.0
294	### 2	729.493444598313	0.000000e+00	2.224120e+05	100.0

                /*****************  su3mat_copy.c  (in su3.a) ***************************
*									*
* void su3mat_copy( su3_matrix *a, su3_matrix *b )			*
* Copy an su3 matrix:  B <- A   						*
*/
#include "config.h"
#include "complex.h"
#include "su3.h"

/* Copy a su3 matrix:  b <- a   */
void su3mat_copy( su3_matrix *a, su3_matrix *b ){
register int i,j;
    for(i=0;i<3;i++)for(j=0;j<3;j++){
	b->e[i][j].real = a->e[i][j].real;
	b->e[i][j].imag = a->e[i][j].imag;
    }
}

invivosu3mat_copy_su3mat_copy_13

Real time: 4.508012e+08 -- Predicted time: 1.863417e+09

Callcount: 7984000

Invocation	Cluster	Part	Invitro (cycles)	Invivo (cycles)	Error (%)
1569622	### 1	5991501.89311295	2.110000e+02	2.800000e+01	86.7
6493442	### 2	1841489.02266925	3.160000e+02	2.400000e+01	92.4
73682	### 3	77231.2005381121	2.240000e+02	3.064000e+03	92.6

                /*******************  m_mat_nn.c  (in su3.a) ****************************
*									*
* void mult_su3_nn( su3_matrix *a,*b,*c )				*
* matrix multiply, no adjoints 						*
* C  <-  A*B								*
*/
#include "config.h"
#include "complex.h"
#include "su3.h"

#ifndef FAST
void mult_su3_nn( su3_matrix *a, su3_matrix *b, su3_matrix *c ){
register int i,j,k;
register complex x,y;
    for(i=0;i<3;i++)for(j=0;j<3;j++){
	x.real=x.imag=0.0;
	for(k=0;k<3;k++){
	    CMUL( a->e[i][k] , b->e[k][j] , y );
	    CSUM( x , y );
	}
	c->e[i][j] = x;
    }
}

/* "Hand coded" routines, clearer coding is up above */
#else
#ifdef NATIVEDOUBLE   /* RS6000 version */

void mult_su3_nn( su3_matrix *a, su3_matrix *b, su3_matrix *c ){
  int j;
  register double a0r,a0i,a1r,a1i,a2r,a2i;
  register double b0r,b0i,b1r,b1i,b2r,b2i;
  
  for(j=0;j<3;j++){
    
    a0r=a->e[0][0].real; a0i=a->e[0][0].imag;
    b0r=b->e[0][j].real; b0i=b->e[0][j].imag;
    a1r=a->e[0][1].real; a1i=a->e[0][1].imag;
    b1r=b->e[1][j].real; b1i=b->e[1][j].imag;
    a2r=a->e[0][2].real; a2i=a->e[0][2].imag;
    b2r=b->e[2][j].real; b2i=b->e[2][j].imag;
    
    c->e[0][j].real = a0r*b0r - a0i*b0i + a1r*b1r - a1i*b1i + a2r*b2r - a2i*b2i;
    c->e[0][j].imag = a0r*b0i + a0i*b0r + a1r*b1i + a1i*b1r + a2r*b2i + a2i*b2r;
    
    a0r=a->e[1][0].real; a0i=a->e[1][0].imag;
    b0r=b->e[0][j].real; b0i=b->e[0][j].imag;
    a1r=a->e[1][1].real; a1i=a->e[1][1].imag;
    b1r=b->e[1][j].real; b1i=b->e[1][j].imag;
    a2r=a->e[1][2].real; a2i=a->e[1][2].imag;
    b2r=b->e[2][j].real; b2i=b->e[2][j].imag;
    
    c->e[1][j].real = a0r*b0r - a0i*b0i + a1r*b1r - a1i*b1i + a2r*b2r - a2i*b2i;
    c->e[1][j].imag = a0r*b0i + a0i*b0r + a1r*b1i + a1i*b1r + a2r*b2i + a2i*b2r;
    
    a0r=a->e[2][0].real; a0i=a->e[2][0].imag;
    b0r=b->e[0][j].real; b0i=b->e[0][j].imag;
    a1r=a->e[2][1].real; a1i=a->e[2][1].imag;
    b1r=b->e[1][j].real; b1i=b->e[1][j].imag;
    a2r=a->e[2][2].real; a2i=a->e[2][2].imag;
    b2r=b->e[2][j].real; b2i=b->e[2][j].imag;
    
    c->e[2][j].real = a0r*b0r - a0i*b0i + a1r*b1r - a1i*b1i + a2r*b2r - a2i*b2i;
    c->e[2][j].imag = a0r*b0i + a0i*b0r + a1r*b1i + a1i*b1r + a2r*b2i + a2i*b2r;
    
  }
}
#else

void mult_su3_nn( su3_matrix *a, su3_matrix *b, su3_matrix *c ){
  int i,j;
  register double t,ar,ai,br,bi,cr,ci;
    for(i=0;i<3;i++)for(j=0;j<3;j++){

	ar=a->e[i][0].real; ai=a->e[i][0].imag;
	br=b->e[0][j].real; bi=b->e[0][j].imag;
	cr=ar*br; t=ai*bi; cr -= t;
	ci=ar*bi; t=ai*br; ci += t;

	ar=a->e[i][1].real; ai=a->e[i][1].imag;
	br=b->e[1][j].real; bi=b->e[1][j].imag;
	t=ar*br; cr += t; t=ai*bi; cr -= t;
	t=ar*bi; ci += t; t=ai*br; ci += t;

	ar=a->e[i][2].real; ai=a->e[i][2].imag;
	br=b->e[2][j].real; bi=b->e[2][j].imag;
	t=ar*br; cr += t; t=ai*bi; cr -= t;
	t=ar*bi; ci += t; t=ai*br; ci += t;

	c->e[i][j].real=cr;
	c->e[i][j].imag=ci;
    }
}
#endif	/* End of "#ifdef NATIVEDOUBLE" */
#endif	/* End of "#ifdef FAST" */

433-milc

Exec Time : 1.330863e+10 cycles

Matching : 97.2%

__invivo__quark_stuff_u_shift_fermion_1452

Real time: 1.434033e+09 -- Predicted time: 1.555542e+09

Callcount: 10160

__invivo__com_vanilla_start_gather_from_temp_581

Real time: 1.784641e+08 -- Predicted time: 2.017281e+08

Callcount: 26472

__invivo__m_mat_an_mult_su3_an_39

Real time: 5.050794e+08 -- Predicted time: 2.255115e+09

Callcount: 3216000

__invivo__control_main_32

Real time: 1.329345e+10 -- Predicted time: 2.289945e+07

Callcount: 1

__invivo__quark_stuff_load_fatlinks_566

Real time: 2.371161e+09 -- Predicted time: 2.412550e+09

Callcount: 22

__invivo__update_update_40

Real time: 1.279301e+10 -- Predicted time: 1.363889e+10

Callcount: 2

__invivo__quark_stuff_compute_gen_staple_1406

Real time: 1.792900e+08 -- Predicted time: 0.000000e+00

Callcount: 792

__invivo__su3mat_copy_su3mat_copy_13

Real time: 4.508012e+08 -- Predicted time: 1.863417e+09

Callcount: 7984000

invivoquark_stuff_u_shift_fermion_1452

invivocom_vanilla_start_gather_from_temp_581

invivom_mat_an_mult_su3_an_39

invivocontrol_main_32

invivoquark_stuff_load_fatlinks_566

invivoupdate_update_40

invivoquark_stuff_compute_gen_staple_1406

invivosu3mat_copy_su3mat_copy_13