/* last modified 7 Sep. 04 */

/* This QCDOC program is for SU(GROUP) lattice gauge fields.  
The lattice dimensions are in size[].  
The number of lattice dimensions is arbitrary.  
The boundaries are periodic; each dimension should be
even for the checkerboard updating used.

    HITS = metropolis hits per link 
    monte() updates lattice with conventinal metropolis algorithm
    overrelax() does an overrelaxation algorithm; somewhat faster
    at decorrelating once the energy is near the right value */

/* some timings

   8^4 SU(3)
    1  node  2.83  sec/sweep, speedup=1
    2  nodes 1.55  sec/sweep, speedup=1.82 
    4  nodes 0.84  sec/sweep, speedup=3.37
    8  nodes 0.471 sec/sweep, speedup=6.01
    16 nodes 0.266 sec/sweep, speedup=10.64

   16^4 SU(3) (needs at least 8 nodes)
     8 nodes 6.93 sec/sweep
     16 nodes 3.77 sec/sweep

   8^4 on 8 nodes: about 37% of time in getlinks()

 */

/* size and other parameters */
#define DIMENSION 4
static int size[DIMENSION]={8,8,8,8};
double beta=5.71;
#define GROUP 3 
#define HITS 10

/* DEBUG turns on more error checking in the communication routines */
#define DEBUG 1
/* maximum length message queue allowed */
#define QLEN 64
#define BLOCKSIZE (2*GROUP*GROUP*sizeof(double))

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include "include/qmemcpyoc.C"

class matrix {
public:
  double real[GROUP][GROUP];
  double imag[GROUP][GROUP];
  void conjugate();
  void project(); /* projects onto the gauge group */
  void determinant(double& rez, double& imz);
  void operator = (const double x);
  void printmatrix();
};

/* for the Metropolis table */
#define TABLELENGTH 64
matrix table[TABLELENGTH];

/* global variables */ 
int nsites;
int nlinks;
int vectorlength, maxvlength; /* work with groups of links together     */
int oldsite=-1, oldlink=-1;   /* to help makelinklist avoid duplication */
int minionnode,maxionnode;    /* first and last+1 onnode link indices   */
double rfactor=(1.0/RAND_MAX); /* for turning rand into double */
int flopcount; /* not all properly flops counted if GROUP>3 */
matrix zeromatrix;   /* convenient for resetting matrices */
matrix tempmatrix;
matrix *atemp, *btemp, *ctemp;
double *sold=NULL, *snew=NULL; /* vector of old and trial actions */
double stot=0.0, eds=0.0, acc=0.0; /* global accumulators */
/* eds will accumulate exp(delta action) where delta action is the
trial action minus the old action used in the metropolis algorithm.
When in equilibrium eds should fluctuate about unity. */
int *sitelist; /* pointer to a block of onnode sites */
int *linklist; /* an array of links for fetching     */

/* prototype things */
void prod(matrix* g1, matrix* g2, matrix* g3);
void init();
void cleanup();
void monte();
void renorm();
void overrelax();
void staple(matrix *st,int site,int link);
void metro(matrix *old,matrix *trial,double bias);
void vgroup(matrix *g);
void thirdrow(matrix* g); /* special for su(3) */
void vcopy(matrix *g1, matrix *g2);
void getlinks(matrix *g,int site,int link);
void getconjugate(matrix *g,int site,int link);
void savelinks(matrix *g,int site,int link);
void vtprod(matrix *g1,matrix *g2,double *s);
void ranmat(matrix *g);
void vprod(matrix *g1,matrix *g2,matrix *g3);
void makelinklist(int site,int link);
void newtable();
void site2coords(int site, int *x);
int  coords2site(int *x);
int  parity(int site);
int  shiftsite(int site,int link,int shift);

int main() {
  int i,iter;
  double mytime;
  printf("clocks per sec = %ld\n",CLOCKS_PER_SEC);
  mycom.start();
  init();
  printf("lattice size:");
  for (i=0;i<DIMENSION-1;i++)
    printf(" %d by",size[i]);
  printf(" %d\n",size[DIMENSION-1]);
  printf("vectorlength = %d\n",vectorlength);
  printf("group=SU(%d)   beta = %6.4f\n",GROUP,beta);

  printf("test monte\n");
  for (iter=0;iter<5;iter++) {
    mytime=clock();
    for (i=0;i<5;i++) {
      monte(); 
    }
    renorm();
    mytime=clock()-mytime;
    mytime=mytime/(5*CLOCKS_PER_SEC);
    printf(" each sweep took %g seconds\n",mytime);
    printf(" running at %g megaflops\n",flopcount/(mytime*1.e6));
    mytime=(1.0E6*mytime)/(1.0*nlinks);
    printf("running at %g microseconds per link\n",mytime);
    printf("processor %d: maxqueue = %d\n",processor,mycom.maxqueue);
  } 

  printf("test overrelax\n");
  for (iter=0;iter<5;iter++) {
    mytime=clock();
    for (i=0;i<5;i++) {
      overrelax(); 
    }
    renorm();
    mytime=clock()-mytime;
    mytime=(1.0E6*mytime)/(5.0*nlinks*CLOCKS_PER_SEC);
    printf("running at %g microseconds per link\n",mytime);
  } 
  printf("processor %d: maxqueue = %d\n",processor,mycom.maxqueue);

  mycom.stop();
  cleanup();
  return 0;
}

/* lattice layout */

/* sites indexed as  x+xsize*(y+ysize*(z+zsize*t)) */
/* links indexed as  link+DIMENSION*site */
/* upper bits of remap(index) determine processor */
/* lower bits give location on the node */

/* use lookup table for faster remapping;
   if memory tight, replace with calls to remap() */
int * remapsite; 

int mask,localbits;
matrix * datastart;
#define node(index) (index>>localbits)
#define address(index) (datastart+(index&mask))
#define onnode(index) (processor==node(remaplink(index)))

void * allocate(int n, int size){
  /* allocate space for n things of given size
     divide evenly and round up to a power of two on each processor */
  int pernode;
  void * ptr;
  pernode=n;
  /* amount per processor rounded up to a power of two */
  pernode = (pernode + NumNodes() - 1) / NumNodes();        
  pernode--; /* change most significant bit if a power of two */
  localbits = 1;
  while (pernode >>= 1)               /* find most significant bit */
    localbits++;
  pernode=(1<<localbits);              /* final amount to allocate */
  mask = pernode-1;
  ptr = (matrix *) qalloc(0, pernode * size);
  if (NULL == ptr) {
    die("memory allocation failed");
  }
#if DEBUG
  printf("allocation of %d objects of size %d requested\n",n,size);
  printf("%d allocated per node; %d total actually allocated\n",
         pernode,NumNodes()*pernode);
  printf("allocated space starts at %p\n",ptr);
#endif
  return ptr;
}

int remap(int i){
  /* divide each dimension in two and spread on different processors */
  /* this uses the evenness of each dimension */
  /* for more than 2^DIMENSION processors, further division goes by 
     the highest coordinates */
  /* my old naive layout just returned i */
  int result, toppart, j, x[DIMENSION];
  div_t t;
  result=0;
  toppart=0;
  site2coords(i,x);
  for (j=DIMENSION-1;j>=0;j--){
    t=div(x[j],size[j]>>1);
    result=result*(size[j]>>1)+t.rem;
    toppart=2*toppart+t.quot;
  }
  result+=(toppart*(nsites>>DIMENSION));
  return result;
}

int remaplink(int i){
  div_t t;
  t=div(i,DIMENSION);
  return remapsite[t.quot] * DIMENSION + t.rem;
}

inline void fetch(matrix* u, int index) {
  index=remaplink(index);
  mycom.qmemcpy(processor,  u,
		node(index),address(index),
		sizeof(matrix));
}

inline void store(matrix* u, int index) {
  index=remaplink(index);
  mycom.qmemcpy(node(index),address(index),
		processor,  u,
		sizeof(matrix));
}

void init() {
  int i,iv;
  flopcount=0;
  zeromatrix=0.;
  /* seed random number generator */
  srand(1234*(UniqueID()+1));
#if DEBUG
  printf("first random number = %g\n",(rfactor*rand()));
#endif
  /* calculate nsites, nlinks */
  nsites=1;
  for (i=0;i<DIMENSION;i++)
    nsites*=size[i];
  nlinks=DIMENSION*nsites;
  /* make space for the lattice */
  datastart=(matrix *) allocate(nlinks,sizeof(matrix));
  /* set up table for site remapping */
  remapsite = (int *) qalloc(0, nsites*sizeof(int));
  if (NULL == remapsite) {
    die("memory allocation failed");
  }
  for (i=0;i<nsites;i++){
    remapsite[i]=remap(i);
  }
  /* find the first and last onnode links */
  minionnode=0;
  vectorlength=0;
  while (!onnode(minionnode)) minionnode++;
  for (iv=minionnode;iv<nlinks;iv++)
    if (onnode(iv)){
      vectorlength++;
      maxionnode=iv+1;
    }
  /* divide vectorlength by two parities and DIMENSION link directions */
  /* round up if odd */
  vectorlength=maxvlength=(1+vectorlength/DIMENSION)/2;
  /* allocate working arrays */
  if (0==(atemp=(matrix *) qalloc(QFAST,sizeof(matrix)*(3*vectorlength))))
    die("memory allocation problems, temp matrices");
  btemp=atemp+vectorlength;
  ctemp=btemp+vectorlength;

  if (0==(sold=(double *) qalloc(0,sizeof(double)*2*vectorlength)))
    die("memory allocation problems, sold");
  snew=sold+vectorlength;

  if (0==(sitelist=(int *) qalloc(0,sizeof(int)*2*vectorlength)))
    die("memory allocation problems, sitelist");
  linklist=sitelist+vectorlength;

  /* set initial links to identity matrix */
  *atemp=1.0;
  for (iv=0;iv<nlinks;iv++)
    if (onnode(iv))
      store(atemp,iv);
  /* initialize Metropolis table */
  newtable();

#if DEBUG
  /* test a few things */
  ranmat(atemp);
  printf("trial random matrix:\n");
  atemp[2].printmatrix();
  atemp[1]=atemp[2];
  printf("conjugate trial matrix:\n");
  (atemp[1]).conjugate();
  (atemp[1]).printmatrix();
  prod(atemp+2,atemp+1,atemp+1);
  printf("U U^dagger:\n");
  (atemp[1]).printmatrix();
  if (processor==0) {
    printf("another trial matrix for storing:\n");
    (atemp[3]).printmatrix();
    store(atemp+3,nlinks/2);
    atemp[1]=0.0;
  }
  mycom.finishcopy();
  if (processor==0) {
    printf("fetch the stored matrix\n");
    fetch(atemp+1,nlinks/2);
    mycom.finishcopy();
    (atemp[1]).printmatrix();
  }  
  else
    mycom.finishcopy();
#endif

  printf("lattice initialization done\n");
  return;
}

void cleanup(){
  if (datastart) qfree(datastart); datastart=NULL;
  if (remapsite) qfree(remapsite); remapsite=NULL;
  if (atemp) qfree(atemp); atemp=NULL;
  if (sold) qfree(sold); sold=NULL;
  if (sitelist) qfree(sitelist); sitelist=NULL;
  return;
}

void vgroup(matrix *g) {
  /* make group elements out of vectorlength matrices g. */
  int iv;
  for(iv=0;iv<vectorlength;iv++)
    g[iv].project();
  return;
}

void monte() {
  int color,link,hit;
  int iv,site;
  flopcount=0;
  stot=eds=acc=0.0;
  newtable();
  sitelist[0]=0; /* in case vectorlength vanishes */
  for (link=0;link<DIMENSION;link++) { /* loop over link directions     */
    for (color=0;color<2;color++) {    /* loop over checkerboard colors */
      vectorlength=0;
      /* make a list of sites with the corresponding links on node */
      for (iv=link,site=0;iv<maxionnode;iv+=DIMENSION,site++){
	if (onnode(iv)) {
          if (color==parity(site))
            sitelist[vectorlength++]=site;
	}
      }
      if (vectorlength>maxvlength) die("too large a vector");
      staple(ctemp,*sitelist,link); /* get neighborhood */
      /* get old links and calculate action */
      getlinks(atemp,*sitelist,link);
      vtprod(atemp,ctemp,sold);
      for (hit=0;hit<HITS;hit++) {
	/* get random matrices */
	ranmat(btemp);
	/* find trial elements and new action */
	vprod(atemp,btemp,btemp);
	vtprod(btemp,ctemp,snew);
	metro(atemp,btemp,beta/(1.*GROUP)); /* metropolis step */
      }
      savelinks(atemp,*sitelist,link); /* save new links */
    }
  }
  /* now accumulate averages */
  stot=mycom.dgsum(stot);
  acc=mycom.dgsum(acc);
  eds=mycom.dgsum(eds);
  stot=stot/(2.*(DIMENSION-1)*nlinks*GROUP*HITS);
  acc=acc/(nlinks*HITS);
  eds=eds/(1.*nlinks*HITS);
  printf("stot=%g, acc=%g, eds=%g\n",stot,acc,eds);
  flopcount=mycom.igsum(flopcount);
  return;
}

void overrelax() {
  matrix temporary;
  int color,link,iv,site;
  stot=eds=acc=0.0;
  sitelist[0]=0; /* in case vectorlength vanishes */
  /* loop over link directions */
  for (link=0;link<DIMENSION;link++) {
    /* loop over checkerboard colors */
    for (color=0;color<2;color++) {
      vectorlength=0;
      for (iv=link,site=0;iv<maxionnode;iv+=DIMENSION,site++){
	if (onnode(iv)) {
          if (color==parity(site))
            sitelist[vectorlength++]=iv/DIMENSION;
	}
      }
      if (vectorlength>maxvlength) die("too large a vector");

      staple(ctemp,*sitelist,link); /* get neighborhood */
      /* get old link and calculate action */
      getlinks(atemp,*sitelist,link);
      vtprod(atemp,ctemp,sold);
      /* find trial element and new action */
      vcopy(ctemp,btemp);
      vgroup(btemp);
      for(iv=0;iv<vectorlength;iv++) {
	prod(atemp+iv,btemp+iv,&temporary);	
	prod(btemp+iv,&temporary,btemp+iv);	
	btemp[iv].conjugate();
      }
      vtprod(btemp,ctemp,snew);
      metro(atemp,btemp,beta/(1.*GROUP)); /* metropolis step */
      savelinks(atemp,*sitelist,link); /* save new links */
    }
  }
  /* now accumulate averages */
  stot=mycom.dgsum(stot);
  acc=mycom.dgsum(acc);
  eds=mycom.dgsum(eds);
  stot=stot/(2.*(DIMENSION-1)*nlinks*GROUP);
  acc=acc/(nlinks);
  eds=eds/(1.*nlinks);
  printf("stot=%g, acc=%g, eds=%g\n",stot,acc,eds);
  return;
}

void getlinks(matrix *g,int site,int link) {
  /* gather alternate links into vector g
     site=starting site
     link=direction */
  int iv;
  makelinklist(site,link);
  for(iv=0;iv<vectorlength;iv++) {
    do{
      mycom.handlemessages();}
    while (2*mycom.queuesize>QLEN);
    fetch(g+iv,linklist[iv]);
  }
  mycom.finishcopy();
  return;
}

void getconjugate(matrix *g,int site,int link) {
  /* gather conjugate links into vector g */
  int iv;
  getlinks(g,site,link);
  for(iv=0;iv<vectorlength;iv++) {
    g[iv].conjugate();
  }
  return;
}

void savelinks(matrix *g,int site,int link) {
  /* scatter alternate links from vector g */
  int iv;
  makelinklist(site,link);
  for(iv=0;iv<vectorlength;iv++) {
    mycom.handlemessages();
    store(g+iv,linklist[iv]);
#if DEBUG
    if (!onnode(linklist[iv]))
      printf("%d: doing an offnode store, %d\n",processor,linklist[iv]);
#endif
  }
  mycom.finishcopy();
  return;
}

void vsum(matrix *g1,matrix *g2,matrix *g3) {
  /* matrix sum of g1 and g2 to g3, vectorlength times */
  int i,j,iv;
  for(iv=0;iv<vectorlength;iv++)
    for(i=0;i<GROUP;i++)
      for(j=0;j<GROUP;j++) {
	g3[iv].real[i][j]=g1[iv].real[i][j]+g2[iv].real[i][j];
	g3[iv].imag[i][j]=g1[iv].imag[i][j]+g2[iv].imag[i][j];
      }
  flopcount+=GROUP*GROUP*vectorlength*2;
  return;
}

void vcopy(matrix *g1, matrix *g2) {
  /* matrix copy of g1 to g2, vectorlength times */
  int iv;
  for(iv=0;iv<vectorlength;iv++) {
    g2[iv]=g1[iv];
  }
  return;
}

void vtprod(matrix *g1,matrix *g2,double *s) {
  /* real trace of product g1 and g2 to s, vectorlength times */
  int i,j,iv;
  for(iv=0;iv<vectorlength;iv++) {
    s[iv]=0.0;
    for(i=0;i<GROUP;i++)
      for(j=0;j<GROUP;j++) 
	s[iv]+=g1[iv].real[i][j]*g2[iv].real[j][i]
	  -g1[iv].imag[i][j]*g2[iv].imag[j][i];   
  }
  flopcount+=GROUP*GROUP*vectorlength*4;
  return;
}

void renorm() { 
  /* put whole lattice in group; call from time to time to keep down
     drift from floating point errors */
  int iv;
  matrix g;
  for (iv=minionnode;iv<maxionnode;iv++) {
    if (onnode(iv)) {
      fetch(&g,iv);
      g.project();
      store(&g, iv);
    }
  }
  return;
}

void staple(matrix *st,int site,int link) {
  /* This subroutine calculates a vector of matrices interacting
     with with links using Wilson action.  The lattice is in l
     and the result is placed in st.   The first two matrix
     vectors are used; so st should not be there.
     links labled as
     x-----link2-----x
     link3       link1
     x-----link -----x
     link6       link4
     x-----link5-----x    */
  int iv,link1,site1,site4,site5,site2;
#define link2 link
#define link5 link
#define link3 link1
#define link4 link1
#define link6 link1
#define site3 site
#define site6 site5
  for(iv=0;iv<vectorlength;iv++) 
    st[iv]=zeromatrix;
  site1=shiftsite(site,link,1);
  /* loop over planes */
  for (link1=0;link1<DIMENSION;link1++)
    if (link1!=link) {
      site2=shiftsite(site,link3,1);
      site4=shiftsite(site1,link1,-1);
      site5=shiftsite(site,link3,-1); 
      getlinks(atemp,site1,link1);
      getconjugate(btemp,site2,link2);
      vprod(atemp,btemp,btemp);
      getconjugate(atemp,site3,link3);
      vprod(btemp,atemp,btemp);
      vsum(st,btemp,st);
      getconjugate(atemp,site4,link4);
      getconjugate(btemp,site5,link5);
      vprod(atemp,btemp,btemp);
      getlinks(atemp,site6,link6);
      vprod(btemp,atemp,btemp);
      vsum(st,btemp,st);
    }
  return;
}

void metro(matrix *old,matrix *trial,double bias) {
  /* accept new for old using metropolis algorithm
     bias multiplies actions in exponential (i.e. beta/GROUP)  
     actions passed in global variables sold and snew */
  int iv;
  double temp;
  int accepted;
  for(iv=0;iv<vectorlength;iv++) {
    temp=exp(bias*(snew[iv]-sold[iv]));
    accepted=((rfactor*rand())<temp);
    if (accepted) {
      sold[iv]=snew[iv];
      old[iv]=trial[iv];
    }
    acc+=1.0*accepted; /* increment global acceptance counter */
    stot+=sold[iv]; /* and action accumulator */
    eds+=temp; /* and eds monitor */
    /* in equilibrium, eds should average to unity */
  }
  /* how should I count flops for exp and rand? */
  flopcount+=7*vectorlength;
  return;
}   

void ranmat(matrix *g) {
  /* generate a vector of random matrices */
  int iv,k;
  k=(int) (TABLELENGTH * (rfactor*rand()));
  for(iv=0;iv<vectorlength;iv++) {
    if (k>=TABLELENGTH) k=0;
    g[iv]=table[k];
    if ((rfactor*rand())<.5) g[iv].conjugate();
    k++;
  }
  flopcount+=vectorlength;
  return;
}
 
#if (GROUP==3)
void thirdrow(matrix* g) {
  /* for su(3) construct third row from first two */
  int i,j,k;
  for (i=0;i<3;i++) {
    j= (i==2) ? 0 : i+1;
    k= (i==0) ? 2 : i-1;
    g->real[2][i]=
      +g->real[0][j]*g->real[1][k]
      -g->imag[0][j]*g->imag[1][k]
      -g->real[1][j]*g->real[0][k]
      +g->imag[1][j]*g->imag[0][k];
    g->imag[2][i]=
      -g->real[0][j]*g->imag[1][k]
      -g->imag[0][j]*g->real[1][k]
      +g->real[1][j]*g->imag[0][k]
      +g->imag[1][j]*g->real[0][k];
  } 
  flopcount+=14*GROUP;
  return;
}
#endif

void prod(matrix* g1, matrix* g2, matrix* g3) {
  /* overloading * does not seem as fast as this function call */
  int i,j,k;
  matrix * t;
  if ((g1==g3)||(g2==g3)) /* use temporary */
    t=&tempmatrix;
  else
    t=g3;
  for( i=0;i<GROUP-(GROUP==3);i++) {
    for (j=0;j<GROUP;j++) {
      t->real[i][j]=
	+g1->real[i][0]*g2->real[0][j]
	-g1->imag[i][0]*g2->imag[0][j];
      t->imag[i][j]= 
	+g1->real[i][0]*g2->imag[0][j]
	+g1->imag[i][0]*g2->real[0][j];
      /* further unwrapping gains nothing substantial */
      for (k=1;k<GROUP;k++) {
	t->real[i][j]+=g1->real[i][k]*g2->real[k][j]
	  -g1->imag[i][k]*g2->imag[k][j];
	t->imag[i][j]+= g1->real[i][k]*g2->imag[k][j]
	  +g1->imag[i][k]*g2->real[k][j];
      }
    }
  }
  flopcount+=GROUP*(GROUP-(GROUP==3))*(6+(GROUP-1)*8);
# if (GROUP==3) 
  thirdrow(t);
#endif
  if (t!=g3)
    *g3=*t;
  return;
}

void vprod(matrix *g1,matrix *g2,matrix *g3) { 
  int iv; 
  for(iv=0;iv<vectorlength;iv++) 
    prod(g1+iv,g2+iv,g3+iv); 
  return; 
}

/* now for a bunch of routines to set the boundary conditions */
void site2coords(int site, int *x) {
  int k;
  div_t t;
  /* break site into site2coordss */
  while (site<0) site+=nsites;       /* make positive */
  while (site>=nsites) site-=nsites; /* and less than nsites */
  k=0;
  while (site) {
    t=div(site,size[k]);
    site=t.quot;
    x[k++]=t.rem;
  }
  while(k<DIMENSION)
    x[k++]=0;
  return;
}

int coords2site(int *x) {
  /* locate site from site2coords */
  int site,k;
  site=0;
  for (k=DIMENSION-1;k>=0;k--) {
    site*=size[k];
    site+=x[k];
  }
#if DEBUG
  if (site>=nsites) die("bad site index");
#endif
  return site;
}

void makelinklist(int site,int link) { 
  /* construct linklist[] to be a vector of links for fetching
     if site is shifted from sitelist[0], make the corresponding
     shift for the links, taking into account boundary conditions */ 
  int iv;
  while (site<0) site+=nsites; 
  while (site>=nsites) site-=nsites; 
  if ((site==oldsite)&&(link==oldlink)) {
    return; /* reuse old list */
  }
  oldsite=site;
  oldlink=link;
  if (site==sitelist[0]){ 
    for(iv=0;iv<vectorlength;iv++) 
      linklist[iv]=link+DIMENSION*sitelist[iv];
  }
  else {
    int x[DIMENSION],xs[DIMENSION];
    site2coords(site,xs);
    site2coords(sitelist[0],x); 
    int i; 
    for (i=0;i<DIMENSION;i++) {
      xs[i]=xs[i]-x[i]; /* shift to site from sitelist[0] */
      while (xs[i]<0) xs[i]+=size[i]; /* make positive */
    } 
    for(iv=0;iv<vectorlength;iv++){
      site2coords(sitelist[iv],x); 
      for (i=0;i<DIMENSION;i++) {
	x[i]+=xs[i]; /* shift the site */
	while (x[i]>=size[i]) x[i]-=size[i]; /* periodic boundaries */
      }
      linklist[iv]=link+DIMENSION*coords2site(x); 
    }
  }
  return;
}

int shiftsite(int site,int link,int shift) { 
/* returns site shifted in direction of link by amount shift */ 
  int x[DIMENSION];
  site2coords(site,x); 
  x[link]+=shift; 
  while (x[link]<0)
    x[link]+=size[link]; 
  while (x[link]>=size[link]) 
    x[link]-=size[link];
  return coords2site(x); 
}

int parity(int site) {
  /* returns 1 if an odd site */
  int k,result,x[DIMENSION];
  site2coords(site,x);
  result=0;
  for (k=0;k<DIMENSION;k++){
    result^=x[k];
  }
  return result&1;
} 

void newtable() {
 int i,j,iv;
  for (iv=0;iv<TABLELENGTH;iv++) {
    table[iv]=beta/GROUP; /* bias towards the identity */
    for(i=0;i<GROUP;i++)for(j=0;j<GROUP;j++) {
      table[iv].real[i][j]+=(rfactor*rand())-0.5;
      table[iv].imag[i][j]+=(rfactor*rand())-0.5;
    } 
    table[iv].project();
    if ((rfactor*rand())<.5) table[iv].conjugate();
  }
  flopcount+=TABLELENGTH*(1+GROUP*GROUP*4);
  return;
}

/* implement the matrix class functions */
void matrix::operator= (const double x) {
  int i,j;
  for(i=0;i<GROUP;i++)
    for(j=0;j<GROUP;j++) {
      real[i][j]= ((i==j) ? x : 0.);
      imag[i][j]=0.;
    }
  return;
}

void matrix::printmatrix() {
  for (int i=0; i<GROUP;i++) {
    for(int j=0;j<GROUP;j++) {
      printf("\t(%6.4g, %6.4g)",real[i][j],imag[i][j]);
    }
    printf("\n");
  }
  printf("\n");
  return;
}

void matrix::conjugate() {
  /* converts a matrix to its conjugate */
  double temp;
  for(int i=0;i<GROUP;i++){
    imag[i][i]=-imag[i][i];
    for(int j=0;j<i;j++) {
      temp=real[i][j];
      real[i][j]=real[j][i];
      real[j][i]=temp;
      temp=imag[i][j];
      imag[i][j]=-imag[j][i];
      imag[j][i]=-temp;
    }
  }
  return;
}

void matrix::project() {
  /* projects a matrix onto the group SU(GROUP) */
  int i,j,k,nmax;
  nmax=GROUP-(GROUP<4); /* 2 and 3 are special */
  /* loop over rows */
  for (i=0;i<nmax;i++) {
    /* normalize i'th row */
    double temp=0.0;
    for (j=0;j<GROUP;j++)
      temp+=real[i][j]*real[i][j]+imag[i][j]*imag[i][j];
    temp=1/sqrt(temp);   
    for (j=0;j<GROUP;j++) {
      real[i][j]*=temp;
      imag[i][j]*=temp;
    }
    /* orthogonalize remaining rows */
    for (k=i+1;k<nmax;k++) {
      double adotbr=0.0,adotbi=0.0;
      for (j=0;j<GROUP;j++) {
	adotbr+=real[i][j]*real[k][j]+imag[i][j]*imag[k][j];
	adotbi+=real[i][j]*imag[k][j]-imag[i][j]*real[k][j];
      }
      for (j=0;j<GROUP;j++) {
	real[k][j]-= (adotbr*real[i][j]-adotbi*imag[i][j]);
	imag[k][j]-= (adotbr*imag[i][j]+adotbi*real[i][j]);
      } 
    } /* end of k loop */
  } /* end of i loop */
  flopcount+=nmax*(GROUP*6+1+nmax*GROUP*16/2);
  /* remove determinant, treating group=2 or 3 as special cases */
#if GROUP==3
  thirdrow(this);
#elif GROUP==2
  real[1][0]=-real[0][1];
  real[1][1]= real[0][0];
  imag[1][0]= imag[0][1];
  imag[1][1]=-imag[0][0];
#else
  /* remove the determinant from the first row */
  double x,y,w;
  determinant(x,y);
  for (i=0;i<GROUP;i++) {
    w=real[0][i]*x
      +imag[0][i]*y;
    imag[0][i]=imag[0][i]*x
      -real[0][i]*y;
    real[0][i]=w;
  }
#endif
  return;
}

#if GROUP>3
/* don't need this for su2 or su3 */
void matrix::determinant(double& detr, double& deti) {
  /* subroutine to calculate determinant of matrix */
  /* could perhaps improve later by doing group=2 and 3 by hand,
     but as of now those cases don't call this routine */
  int i,j,k,im,km;
  double temp;
  matrix copy;
  copy=*this;
  im=GROUP-1;
  for (i=0;i<im;i++) {
    km=i+1;
    /* find magnitude of i'th diagonal element */
    temp=1./(copy.real[i][i]*copy.real[i][i]
	     +copy.imag[i][i]*copy.imag[i][i]);
    for (k=km;k<GROUP;k++) {
      /* subtract part of row i from row k to make [k][i] element vanish */
      /* inner product of the rows */
      detr=(copy.real[k][i]*copy.real[i][i]
	    +copy.imag[k][i]*copy.imag[i][i])*temp;
      deti=(copy.imag[k][i]*copy.real[i][i]
	    -copy.real[k][i]*copy.imag[i][i])*temp;
      for (j=km;j<GROUP;j++) {
	copy.real[k][j]+= -detr*copy.real[i][j]
	  +deti*copy.imag[i][j];
	copy.imag[k][j]+=  -detr*copy.imag[i][j]
	  -deti*copy.real[i][j];
      }
    }  
  }  
  /* multiply diagonal elements */
  detr=copy.real[0][0];
  deti=copy.imag[0][0];
  for (i=1;i<GROUP;i++) {
    temp=detr*copy.real[i][i]
      -deti*copy.imag[i][i];
    deti=deti*copy.real[i][i]
      +detr*copy.imag[i][i];
    detr=temp;
  }
  return;
}
#endif