#include <stdio.h>
#include "include.h"
#include "run.h"
#include "usempi.h" /* this is for the prototypes in this file */

#if !ONE_PROCESSOR_SPECTRUM

#define MS_PRINT 0 /* flag to print out debugging information to files */

#define PRINTSLAVE 0  /* print statements associated with slaves */
#define PRINTMASTER 0  /* print statements associated with master */

static MPI_Datatype mpisectors;


void slave(void)
{ 
  integer nfig=EDIGITS; /* number of significant digits in eigenvalues */
  static sectors x, y;
  int nvals, globalkt, my_id, death=1;
  element *val, *vec, couplings[NPARAMS], angle[HINDEX+2];
#if 0
  size_t gmb;
  int gme;
#endif

#if MS_PRINT 
  FILE *fp;
  char out[10];
  int i, j;
#endif 
  
  MPI_Comm_rank(MPI_COMM_WORLD,&my_id);

#if MS_PRINT /* open debugging file with name based on rank */
  sprintf(out,"p%i.tmp",my_id);
  if((fp = fopen(out,"w")) == NULL)
    {   
      fprintf(stderr,"Can't open file %s in pusempi.c, exiting\n",out);
      goto error;
    }
#endif

  /* on CRAY, be more carefull about memory management.  */

#ifdef _CRAY
  mallopt(M_LOWFIT,1);
#endif
  
  /* Main loop of program */
      
  for(;;)
    {
      MPI_Barrier(MPI_COMM_WORLD);
      MPI_Bcast(&death, 1, MPI_INT, MASTER, MPI_COMM_WORLD);
      /* Stop program */ 
      if(death==0)
	goto exit;
      
      /* Receive coupling constants and global kt. */
	  
      else
	{
	  MPI_Bcast(couplings, NPARAMS, MPI_DOUBLE, MASTER, MPI_COMM_WORLD);
	  MPI_Bcast(&globalkt, 1, MPI_INT, MASTER, MPI_COMM_WORLD);
	  
#if MS_PRINT /* diagnostic print to debugging file*/
	  fprintf(fp,"couplings= ");
	  for(i=0; i<NPARAMS; i++)fprintf(fp,"%f ",couplings[i]);
	  fprintf(fp,"globalkt=%i\n",globalkt);
	  fprintf(fp,"About to receive Job\n");
	  fflush(fp);
#endif
	  
	  MPI_Bcast(&x, 1, mpisectors, MASTER, MPI_COMM_WORLD);
	  MPI_Bcast(&y, 1, mpisectors, MASTER, MPI_COMM_WORLD);
	  MPI_Bcast(angle, HINDEX+2, MPI_DOUBLE, MASTER, MPI_COMM_WORLD);
	  MPI_Bcast(&nvals, 1, MPI_INT, MASTER, MPI_COMM_WORLD);
	  
#if MS_PRINT /* diagnostic print to debugging file */
	  fprintf(fp,"Received job, nvals=%i\n",nvals);
	  fprintf(fp,"  ht=(");
	  for(i=0; i<HINDEX; i++)fprintf(fp," %i",x.ht[i]);
	  fprintf(fp,"), np=%i kt=%i cyclic=%i o=%i and multi=%i multi2=%i,\n",
		  x.npmax,x.kt,x.cyclic,x.o,x.multi,y.multi);
	  fprintf(fp,"  angles=");
	  for(j=0; j<HINDEX+2; j++)fprintf(fp," %g",angle[j]);
	  fprintf(fp,"\n");
	  for(j=0; j<2; j++)
	    if(x.cyclic<0||x.cyclic>1||x.o<-1||x.o>1 ||
	       y.cyclic<0||y.cyclic>1||y.o<-1||y.o>1)
	      {
		fprintf(stderr,"error in values sent, exiting\n");
		goto error;
	      }
#endif
	  
	  /*
	    construct one or two bases.  The order for makebasis() 
	    and rebasis() is important for memory fragmentation 
	    problems on the CRAY. 
	    mapcouplings should have already been run.
	  */
	  
	  
	  /* Now calculate the actual spectrum */
#if PRINTSLAVE
	  printf("\n%d Slaves -- Going into spectrum2\n", my_id);
#endif
	  spectrum2(&x, &y, &val, &vec, nvals, couplings, angle, globalkt, nfig);
#if PRINTSLAVE
	  printf("\n%d Slaves -- returning from spectrum2\n", my_id);
#endif 
	  if(val==NULL)
	    {
		  fprintf(stderr,"slave() %i detecting an error in spectrum()\n", my_id);
		  goto error;
	    }
	  /*  Now, for the J^++ sector, we need to order the 0++*
	      and 2++ according to particle number content */
	  shuffle(&x, &y, vec, val, nvals);
	  
#if MS_PRINT /* diagnositic print */
	  fprintf(fp,"  values: ");
	  for(i=0; i<nvals; i++)fprintf(fp,"%.10g ",val[i]);
	  fprintf(fp,"\n");
#endif
	  
	  /*       To prevent memory fragmentation problems on systems
		   without virtual memory paging (like the CRAY), we
		   deallocate the memory used by spectrum() in reverse order.
		   
		   The other two things that can be done  to minimize
		   memory fragmentation is to choose
		   sufficiently large values for MAX_ELS and MAX_BASE.
	  */
	  
#if defined(_CRAY)
	  freespectrum();
#endif
	  freebasis(&y);
	  freebasis(&x);
	}
      
      /* for normal case, just do next loop iteration */
#if MS_PRINT && 0 /* if everything is OK, write over old info */
      rewind(fp);
#endif
    }
  
 error:  /* On error, print out memory usage and tell the master. */
  /* The master will then send out a command that everyone quit. */
  fprintf(stderr,"Error in slave %i after receiving job, exiting.\n",my_id);
  
#if MS_PRINT
  printhels(fp);
  printelements(fp);
#endif
  
#ifdef _CRAY
#if 0 /* Choose one of the following */
  printf("MM Running malloc_stats(0), slave %i:\n",my_id);
  malloc_stats(0);
#elif 0
  printf("MM Running malloc_stats(1), slave %i:\n",my_id);
  malloc_stats(1);
#elif 1  /* the following is very long */
  printf("MM Running malloc_stats(2), slave %i:\n",my_id);
  malloc_stats(2);
#endif
#endif
  
 exit:  /* Normal exit */
#if MS_PRINT /* debugging file with name based on rank */
  fclose(fp);
#endif
  
  return;
}

void initmpisectors(void)
{
  
  int i,seclength[6]={HINDEX,1,1,1,1,1};
  sectors mx;
  MPI_Aint secdisp[6],sec0;
#if 0
  int master[1]={0};
  MPI_GROUP g1,g2;
#endif  
  /* make MPI datatype for the structure sectors */
  /* note that this only works for part of the structure,
     the basis etc. is not included. */
  
  MPI_Address(mx.ht,secdisp);
  MPI_Address(&mx.npmax,secdisp+1);
  MPI_Address(&mx.kt,secdisp+2);
  MPI_Address(&mx.cyclic,secdisp+3);
  MPI_Address(&mx.o,secdisp+4);
  MPI_Address(&mx.multi,secdisp+5);

  /* we find offsets relative to the "base" of the structure,
     wherever that may be (thus we make no assumptions of order. */
  MPI_Address(&mx,&sec0);
  for(i=0; i<6; i++)secdisp[i]-=sec0;
  MPI_Type_hindexed(6,seclength,secdisp,MPI_INT,&mpisectors);
  MPI_Type_commit(&mpisectors);
  /* use this in the global variable mpitype */
#if FANCY /* flag for sending format */
  mpitype[0] = mpitype[1] = mpisectors;
#endif
}

void MasterShutDown(void)
{
  int death = 0;
  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Bcast(&death, 1, MPI_INT, MASTER, MPI_COMM_WORLD);
  return;
}

void spectrum(sectors *x, sectors *y, doublereal **val, doublereal **vec, 
	    integer nval, element *couplings, element *angle, int globalkt, integer nfig)
{
  int death = 1; /* don't tell everyone to exit */
#if PRINTMASTER
  printf("\nMaster -- Got into spectrum(master)\n");
#endif

  MPI_Barrier(MPI_COMM_WORLD);

#if PRINTMASTER
  printf("\nMaster -- Passed Barrier of spectrum\n");
#endif

  MPI_Bcast(&death, 1, MPI_INT, MASTER, MPI_COMM_WORLD);
  MPI_Bcast(couplings, NPARAMS, MPI_DOUBLE, MASTER, MPI_COMM_WORLD);
  MPI_Bcast(&globalkt, 1, MPI_INT, MASTER, MPI_COMM_WORLD);
  MPI_Bcast(x, 1,  mpisectors, MASTER, MPI_COMM_WORLD);
  MPI_Bcast(y, 1, mpisectors, MASTER, MPI_COMM_WORLD);
  MPI_Bcast(angle, HINDEX + 2, MPI_DOUBLE, MASTER, MPI_COMM_WORLD);
  MPI_Bcast(&nval, 1, MPI_INT, MASTER, MPI_COMM_WORLD);

#if PRINTMASTER
  printf("\nMaster -- Going into spectrum2\n");
#endif

  spectrum2(x, y, val, vec, nval, couplings, angle, globalkt, nfig);

#if PRINTMASTER
  printf("\nMaster -- Returning from spectrum2\n");
#endif

  return;
}
#endif