#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include "include.h"
#if defined(USE_MPI)
#include "usempi.h"
#endif
#include "run.h"

#define MOM 0 /*  flag to do nonzero transverse momentum */
              /*  This flag should be zero for nonzero L  */
#define NLONG 15 /* number of step of L in londidudinal direction */
#define LONG 7.5 /* maximum value of L */
          /* paper II best fit: n=4 is L=10.14 */

/** This file runs the lanczos algorithm for finding the lowest
    Eigenvalues of the Hamiltonian.  It calculates the winding
    spectra for all allowed transverse momenta if nonzero winding
    flag MOM is true.                         
    
    For the heavy quark potential, a range of longitudinal
    momenta is chosen. In this case angle[1]=kmax and and angle[0]=L.
**/

/************************************+++++++*****************************

 MPI command line: 
 (the flags -stdin and -stdout should work with mpich, 
 but stdin etc. do not seem to get redirected.)
 However, the following does work:

mpirun -np 15 -v -machinefile machines.LINUX mpiwind <junk.in >&junk.tmp &

2+1 case:

mpirun -np 10 -v -machinefile machines mpiwind <junk2.in >&junk2.tmp &
mpirun -np 8 -v -machinefile machines mw <wind03.in >&wind03.tmp &

************************************+++++++*****************************/


int MAIN(int argc, char **argv){
  clock_t t1;
  element *vals,rerescale;
  /* multi must always be zero! */
#define NOUT 100 /* Output file name length */
  char out[NOUT]="wind.out",*pout;
  int i,j,k,l,nval;
  integer th[1];
  FILE *fp;

/**********************************************************/
#if HINDEX==2  /* 3+1 Dimensions */
/**********************************************************/
  
#if 0 /* longitudinal string tension, test basis */
      /* "smaller test basis" from jessica3.c */
#define NWIND 3
#define NS 4
  int kt[NWIND][NS]={{-16,-12,-12,-12},{-7,-7,-7,-11},{-8,-8,-10,-12}}, 
      np[NWIND][NS]={{2,4,2,2},{3,5,3,3},{2,4,2,2}},
      ht[NWIND][HINDEX]={{0,0},{1,0},{1,1}}, cyclic=0;
   element kmax[NS]={3,3,2,3};
   int multi[NWIND]={3,15,17}, multi2[NS]={3,15,17}, o=1;

#elif 0 /*winding spectrum, test case. */
#define NWIND 4
#define NS 4
  int cyclic=1, o=1, kt[NWIND][NS]={
    {10,10,12,16},{9,9,13,15},
    {8,8,12,16},{8,8,12,14}},
    np[NWIND][NS]={{4,6,4,4},{5,7,5,5},{6,8,6,6},{6,8,6,6}},
    ht[NWIND][HINDEX]={{2,0},{3,0},{4,0},{2,2}};
  element *kmax=NULL;
  int multi[NWIND]={15,15,15,17},*multi2=NULL;

#elif 1 /*Longitudinal potential, BIG basis, */
        /*stuff calculated in long runs.*/
#define NWIND 3
#define NS 7
  int kt[NWIND][NS]={{-34,-34,-34,-40,-50,-60,-60},
		     {-17,-17,-17,-33,-33,-55,-55},
		     {-30,-30,-30,-40,-50,-60,-60}}, 
      np[NWIND][NS]={{2,4,2,2,2,2,2},{3,5,3,3,3,3,3},{2,4,2,2,2,2,2}},
      ht[NWIND][HINDEX]={{0,0},{1,0},{1,1}}, cyclic=0;
  element kmax[NS]={4,4,5,4,5,4,6};
  int multi[NWIND]={3,15,17}, multi2[NWIND]={3,15,17}, o=1;

#elif 0 /*winding spectrum, sizes used in BIG runs. */
#define NWIND 4
#define NS 4
  int cyclic=1, o=1, kt[NWIND][NS]={
    {20,20,24,32},{19,19,25,33},
    {18,18,24,32},{16,16,24,28}},
    np[NWIND][NS]={{4,6,4,4},{5,7,5,5},{6,8,6,6},{6,8,6,6}},
    ht[NWIND][HINDEX]={{2,0},{3,0},{4,0},{2,2}};
  element *kmax=NULL;
  int multi[NWIND]={15,15,15,17},*multi2=NULL;

#elif 0 /*Longitudinal potential, but with larger bases than above */
        /*Using a limit of about 60,000 states.  
	  Global maximum of max_els is 4396898
	  Global maximum of max_base is 2058992 */
#define NWIND 9
#define NS 7
  int cyclic=0, o=1, kt[NWIND][NS]={
    {-18,-18,-18,-34,-34,-44,-44},
    {-18,-18,-18,-34,-34,-44,-44},
    {-18,-18,-18,-34,-34,-44,-44},
    {-18,-18,-18,-34,-34,-44,-44},
    {-18,-18,-18,-34,-34,-44,-44},
    {-19,-19,-19,-43,-43,-65,-65},
    {-19,-19,-19,-43,-43,-65,-65},
    {-16,-16,-16,-30,-30,-36,-36}, 
    {-16,-16,-16,-30,-30,-36,-36}}, 
   np[NWIND][NS]={{4,6,4,4,4,4,4},{4,6,4,4,4,4,4},{4,6,4,4,4,4,4},
		  {4,6,4,4,4,4,4},{4,6,4,4,4,4,4},{3,5,3,3,3,3,3},
		  {3,5,3,3,3,3,3},{4,6,4,4,4,4,4},{4,6,4,4,4,4,4}},
   ht[NWIND][HINDEX]={{0,0},{0,0},{0,0},{0,0},{0,0},
		      {1,0},{1,0},{1,1},{1,1}};
  element kmax[NS]={3.0,3.0,4.5,3.0,4.5,3.0,4.5};
  int multi[NWIND] ={3,4,5,6,7,15,16,17,18}, 
      multi2[NWIND]={3,4,5,6,7,15,16,17,18};

#elif 0 /*winding spectrum, but with larger bases than above */
        /*Using a limit of about 60,000 states. */
#define NWIND 16
#define NS 4
  int cyclic=1, o=1, kt[NWIND][NS]={
    {28,28,40,60},{23,23,25,59},{20,20,28,40},{19,19,25,33},
    {28,28,40,60},{23,23,25,59},{20,20,28,40},{19,19,25,33},
    {28,28,40,60},{18,18,24,34},
    {28,28,40,60},{18,18,24,34},
    {15,15,19,23},{19,19,31,41},{16,16,22,30},{15,15,19,25}},
  np[NWIND][NS]={{4,6,4,4},{5,7,5,5},{6,8,6,6},{7,9,7,7},
		 {4,6,4,4},{5,7,5,5},{6,8,6,6},{7,9,7,7},
		 {4,6,4,4},{6,8,6,6},
		 {4,6,4,4},{6,8,6,6},
		 {7,9,7,7},{5,7,5,5},{6,8,6,6},{7,9,7,7}},
  ht[NWIND][HINDEX]={{2,0},{3,0},{4,0},{5,0},
		     {2,0},{3,0},{4,0},{5,0},
		     {1,1},{2,2},
		     {1,1},{2,2},
		     {3,2},{2,1},{3,1},{4,1}};
  element *kmax=NULL;
  int multi[NWIND]={15,15,15,15,16,16,16,16,17,17,18,18,0,0,0,0},*multi2=NULL;

#elif 0 /*Longitudinal potential, BIG basis, stuff calculated  */
        /*In long runs. Add more values of K and K_max etc. 
          to check convergence in K and K_max.
	  Global maximum of max_els is 3138688
	  Global maximum of max_base is 1010362 */
#define NWIND 3
#define NS 30
  int kt[NWIND][NS]={{-34,-34,-34,-34,-34,-34,-34,-34,-34,-34,
                      -44,-44,-44,-44,-44,-44,-44,-44,-44,-44,
		      -58,-58,-58,-58,-58,-80,-80,-80,-80,-80},
		     {-17,-17,-17,-17,-17,-17,-17,-17,-17,-17,
		      -19,-19,-19,-19,-19,-19,-19,-19,-19,-19,
		      -43,-43,-43,-43,-65,-65,-65,-65,-65,-65},
		     {-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,
		      -36,-36,-36,-36,-36,-36,-36,-36,-36,-36,
		      -50,-50,-50,-50,-50,-80,-80,-80,-80,-80}}, 
      np[NWIND][NS]={{4,4,4,4,4,2,2,2,2,2,
		      4,4,4,4,4,2,2,2,2,2,
		      2,2,2,2,2,2,2,2,2,2},
		     {5,5,5,5,5,3,3,3,3,3,
		      5,5,5,5,5,3,3,3,3,3,
		      3,3,3,3,3,3,3,3,3,3},
		     {4,4,4,4,4,2,2,2,2,2,
		      4,4,4,4,4,2,2,2,2,2,
		      2,2,2,2,2,2,2,2,2,2}},
      ht[NWIND][HINDEX]={{0,0},{1,0},{1,1}}, cyclic=0;
  element kmax[NS]={2,4,6,8,10,2,4,6,8,10,
		    2,4,6,8,10,2,4,6,8,10,
		    2,4,6,8,10,2,4,6,8,10};
  int multi[NWIND]={3,15,17}, multi2[NWIND]={3,15,17}, o=1;

#endif

/******Couplings***********/

#if 0 /* test case (these are test couplings from file runlaso.c)*/
  element  couplings[NPARAMS]={0.03, 1.0, 0.7,0.03,-0.005, 3.0,0.13,1.0,-2.0},
               rescale=1.0,c=-1.0;
#elif 0  /* Best fit from tbig15.out */
  element  couplings[NPARAMS]={0.032492,1.000000,0.559433,-0.057116,
			-0.217486,18.311338,0.107185,2.736781,-1.135574},
               rescale=8.002241,c=-1.237609;
#elif 0  /* Best fit from cbig41.out */
  element  couplings[NPARAMS]={0.0517766953,1.0,0.7584261311,-0.07623059914,
      -0.1458991891,8.889796764,0.1608104465,1.657995607,-3.268403913},
      rescale=8.099820052,c=-1.174339527;
#elif 0  /* couplings from best fit April 1999, ttraj31.out  */
  element couplings[NPARAMS]={0.001967542671,1,0.3909839522,-0.01250282378,
     -0.06421231123,353.2492529,0.07049127522,160.6273835,
     -0.7123520719,-0.6519355869},
     rescale=11.4617271, c=-1.764202771;
#elif 0  /* couplings from best fit April 1999, intermediate, ttraj33.out */
  element couplings[NPARAMS]={0.03249196962,1,0.7317766758,-0.03327406461,
      -0.08809601245,220.8666429,0.1669097664,94.36589102,-0.9688626138,
      -0.8735016911},
     rescale=9.669785336, c=-1.24351091;
#elif 0  /* couplings from best fit April 1999, large mass, ttraj34.out */
  element couplings[NPARAMS]={0.1072401379,1,1.009994691,-0.00729168206,
     -0.1435818867,60.68327765,0.1625519907,9.60935071,-0.7883196908,
     -0.7342659374},
     rescale=6.651322009, c=-1.117951894;
#elif 1  /* couplings from best fit in summer 1999, file ttraj60.out  */
  element couplings[NPARAMS]={0.001967542671,1,0.4286235989,0.01161817889,
    -0.06758574491,436.6795354,0.09711799423,181.7956756,-0.7951513785,
    -0.7015191923},
    rescale=12.02238135,c=-2.212418668;
#endif

/**********************************************************/
#elif HINDEX==1  /* 2+1 Dimensions */
/**********************************************************/
#if 0 /* longitudinal string tension, test basis */
#define NWIND 4
#define NS 1
  int kt[NWIND][NS]={{-10},{-9},{-10},{-9}}, 
      np[NWIND][NS]={{4},{5},{6},{7}},
      ht[NWIND][HINDEX]={{0},{1},{2},{3}}, cyclic=0;
  element kmax[NS]={4};
  int multi[NWIND]={1,0,0,0}, multi2[NWIND]={1,0,0,0},o=1;

#elif 0 /* longitudinal string tension, without extrapolation  */
#define NWIND 5
#define NS 1
  int kt[NWIND][NS]={{-28},{-33},{-28},{-23},{-22}}, 
      np[NWIND][NS]={{4},{3},{4},{5},{6}},
      ht[NWIND][HINDEX]={{0},{1},{2},{3},{4}}, cyclic=0;
  element kmax[NS]={4,4,5,4,5,4};

#elif 0 /* test case, with without extrapolations, 
	      longitudinal string tension  */
#define NWIND 1
#define NS 1
  int kt[NWIND][NS]={{-18}}, 
      np[NWIND][NS]={{2}},
      ht[NWIND][HINDEX]={{0}}, cyclic=0;
  element kmax[NS]={2};

#elif 0 /* test case, with extrapolations, 
	      longitudinal string tension  */
#define NWIND 1
#define NS 5
  int kt[NWIND][NS]={{-6,-6,-8,-10,-10}}, 
      np[NWIND][NS]={{4,4,4,2,4}},
      ht[NWIND][HINDEX]={{2}}, cyclic=0;
  element kmax[NS]={4,4,5,4,5,4};

#elif 0 /*Longitudinal potential, larger basis  */
#define NWIND 5
#define NS 4
  int kt[NWIND][NS]={{-16,-16,-20,-28},{-19,-19,-23,-33},
		     {-16,-16,-20,-28},{-17,-17,-19,-23},
		     {-16,-16,-18,-22}}, 
      np[NWIND][NS]={{4,6,4,4},{3,5,3,3},
		     {4,6,4,4},{5,7,5,5},
		     {6,8,6,6}},
      ht[NWIND][HINDEX]={{0},{1},{2},{3},{4}}, cyclic=0;
   element kmax[NS]={4,4,5,4,5,4};

#elif 0 /*Longitudinal potential, basis with more particles */
#define NWIND 3
#define NS 6
  int kt[NWIND][NS]={{-22,-22,-28,-32,-40,-40},
		     {-22,-22,-28,-32,-40,-40},
		     {-19,-19,-19,-33,-33,-49}}, 
      np[NWIND][NS]={{4,6,4,4,4,4},{4,6,4,4,4,4},{3,5,3,3,3,3}},
      ht[NWIND][HINDEX]={{0},{0},{1}}, cyclic=0;
  element kmax[NS]={4,4,5,4,5,4};
  int multi[NWIND]={1,2,0}, multi2[NWIND]={1,2,0},o=1;

#elif 0 /*Longitudinal potential, basis von paper II (when applicable)*/
#define NWIND 2
#define NS 6
  int kt[NWIND][NS]={{-32,-32,-32,-34,-44,-60},
		     {-19,-19,-19,-33,-33,-49},
		     {-22,-22,-22,-34,-44,-40},
		     {-19,-19,-19,-23,-25,-31},
		     {-18,-18,-18,-22,-24,-28}}, 
      np[NWIND][NS]={{2,4,2,2,2,2},{3,5,3,3,3,3},
		     {4,6,4,4,4,4},{5,7,5,5,5,5},
		     {6,8,6,6,6,6}},
      ht[NWIND][HINDEX]={{0},{1},{2},{3},{4}}, cyclic=0;
  element kmax[NS]={4.,4.,5.,4.,5.,4.};
  int multi[NWIND]={1,0,0,0,0}, multi2[NWIND]={1,0,0,0,0},o=1;

#elif 0 /* Longitudinal potential, big basis to check convergence */
#define NWIND 2
#define NS 24
  int kt[NWIND][NS]={{-32,-32,-32,-40,-40,-40,-32,-32,-32,-40,-40,-40,
		      -60,-60,-60,-60, -80,-80,-80,-80, -120,-120,-120,-120},
		     {-19,-19,-19,-23,-23,-23,-19,-19,-19,-23,-23,-23,
		      -33,-33,-33,-33, -49,-49,-49,-49, -75,-75,-75,-75}}, 
      np[NWIND][NS]={{4,4,4,4,4,4,2,2,2,2,2,2, 
		      2,2,2,2, 2,2,2,2, 2,2,2,2},
		     {5,5,5,5,5,5,3,3,3,3,3,3, 
		      3,3,3,3, 3,3,3,3, 3,3,3,3}},
      ht[NWIND][HINDEX]={{0},{1}}, cyclic=0;
  element kmax[NS]=  {3,4,5,3,4,5,3,4,5,3,4,5, 
		      3,4,5,6, 3,4,5,6, 3,4,5,6};
  int multi[NWIND]={1,0}, multi2[NWIND]={1,0},o=1;

#elif 0 /* winding spectrum, truncations used in paper I */
#define NWIND 4
#define NS 1
  int kt[NWIND][NS]={{22},{21},{22},{21}}, np[NWIND][NS]={{6},{7},{8},{9}},
      ht[NWIND][HINDEX]={{2},{3},{4},{5}}, cyclic=1;
  element *kmax=NULL;

#elif 0 /*winding spectrum, with extrapolation */
#define NWIND 4
#define NS 4
  int kt[NWIND][NS]={{20,22,22,26},{19,21,21,25},{20,22,22,26},{19,21,21,25}}, 
      np[NWIND][NS]={{6,4,6,4},{7,6,7,6},{8,7,8,7},{9,8,9,8}},
      ht[NWIND][HINDEX]={{2},{3},{4},{5}}, cyclic=1;
  element *kmax=NULL;

#elif 0 /* winding spectrum, smaller test basis */
#define NWIND 3
#define NS 3
  static int kt[NWIND][NS]={{21,23,29},{22,24,30},{21,23,29}}, 
    np[NWIND][NS]={{5,5,3},{6,6,4},{7,7,5}},
    ht[NWIND][HINDEX]={{3},{4},{5}}, cyclic=1;
  element *kmax=NULL;

#elif 0 /* winding spectrum, Larger basis */
#define NWIND 3
#define NS 4
  static int kt[NWIND][NS]={{21,23,23,27},{20,22,22,26},{21,23,23,27}}, 
    np[NWIND][NS]={{7,5,7,5},{8,6,8,6},{9,7,9,7}},
    ht[NWIND][HINDEX]={{3},{4},{5}}, cyclic=1;
  element *kmax=NULL;
#endif

/******Couplings***********/
  element c=-1.0;
#if 0 /* test case */
  element couplings[NPARAMS]={0.1,1.0,0.0,0.0,0.0,0.0},rescale=1.0;
#elif 0  /* Best fit from paper I */
  element in[NPARAMS]= {0.032492, 1.000000, -0.100803, -0.275192, 
			58.328500, 0.0},rescale=2.6340518;
#elif 0 /* Best fit from paper I */
  element couplings[NPARAMS]={0.0649839, 2., -0.201606821, -0.5503840688, 116.657},
      rescale=3.46911443596316;
#elif 0 /* Best fit from paper II */
  element couplings[NPARAMS]={0.032492, 1.000000, -0.077209, -0.250795, 
		       221.029589, -0.960551}, rescale=5.221990;
#elif 1 /* Test case for new tau interactions from straj11.out */
  element couplings[NPARAMS]={0.032492,1.000000,0.465808,-0.515163,3.883695,
		       1.481709,-1.795479}, rescale=1.0;
#endif
/**********************************************************/
#endif /* number of HINDEX */
/**********************************************************/

  static sectors s[NWIND][NS],s2[NWIND][NS];
  element pi;
#ifdef USE_MPI
  int myrank;
#define SMAX ((NLONG+1)*NWIND*NS) /* total number of matrices to calculate */
  sectors *mpibase[SMAX][2];
  int job,njob,mpinvals[SMAX],whichvals[SMAX];
  element angle[SMAX][HINDEX+2];
  static element *mpivals=NULL;
#else 
  element angle[HINDEX];
  int ipt[HINDEX];
#endif

#ifdef BABBAGE  /* Routine to initialize Fortran on Babbage */
  if(hf_fint((char*)NULL)) exit(1);
#endif

#ifdef USE_MPI
  MPI_Init(&argc, &argv);     /* initialize MPI */
  initmpisectors();
  MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
  if(myrank!=0){
    slave();
    MPI_Finalize();   /* cleanup MPI */
    return 0;
  }
#endif
  pi= atan(1.) * 4.0;

#if 1  /* Ask for output file name*/
  printf("Output file name:  ");
  fgets(out,NOUT,stdin);
  if((pout=strchr(out,'\n'))!=NULL)*pout='\0';
  printf("%s\n",out);
#endif

  if((fp = fopen(out,"w")) == NULL){
     fprintf(stderr,__FILE__":  Can't open file %s in wind.c, exiting\n",out);
    goto cleanup;
  }

  /*  Print out the results in the output file **/
#if HINDEX==1
  fprintf(fp,"***** Spectra for the 2+1 transverse lattice *****\n");
  fprintf(fp,"Lowest %i Eigenvalue(s) for various momenta and n\n",nval);
  fprintf(fp,"Couplings:  m^2, G^2 N, la1, la_2, la_3, ");
#elif HINDEX==2
  fprintf(fp,"******* Eigenvalues for the 3+1 transverse lattice  *******\n");
  fprintf(fp,"Couplings:  m^2, G^2 N, beta, la_1, la_2, la3, la4, la5,\n");
  fprintf(fp,"            ");
#endif
#if NEWTAU
  fprintf(fp,"tau_1, tau_2\n");
#else
  fprintf(fp,"tau\n");
#endif
  fprintf(fp,"For heavy sources, K<0, divide by sqrt(G^2 N).\n");
  for(i=0; i<NWIND; i++){
    fprintf(fp,"Winding (");
    for(l=0; l<HINDEX; l++)fprintf(fp," %i",ht[i][l]);
    fprintf(fp,") spectra using (K,p) = ");
    for(j=0, k=4; j<NS; k++, j++){
      fprintf(fp,"(%i/2,%i) ",kt[i][j],np[i][j]);
      if(k%5==0 && j+1<NS)fprintf(fp,"\n ");
    }
    if(!cyclic || MOM)
      fprintf(fp,"\n with multi=%i & %i.\n",multi[i],multi2[i]);
    else
      fprintf(fp,"\n with multi=%i.\n",multi[i]);
  }
  if(MOM || (!cyclic && kt[0][0]<0))
    fprintf(fp,"Orientation symmetry indefinite");
  else
    fprintf(fp,"Orientation symmetry %i",o);
  if(!cyclic && kt[0][0]<0){
    fprintf(fp,"\nHeavy sources, K_max = ");
    for(i=0, k=8; i<NS; k++, i++){
      fprintf(fp," %g",kmax[i]);
      if(k%15==0 && i+1<NS)fprintf(fp,"\n ");
    }
  }
  fprintf(fp,"\nResult format:\n"
	  " p damping, overall scale G^2 N/sigma,");
  if(!MOM && (cyclic || kt[0][0]>0))
    fprintf(fp," orientation symmetry,");
  fprintf(fp,"\n the %i couplings,\n"
	  " winding number,  momenta (or separation)\n"
	  "   spectra in units of the string tension\n",NPARAMS);
  fflush(fp);

  /****** Begin main loop of program  **********/

  for(;;){
 
    printf("Input number of eigenvalues (0=exit):  ");
    scanf("%i",&nval);
    if(nval<1)break;

/******  Input "o", this is useful for the winding potential  *******/
#if 1   /* input parameters for basis (instead of defaults.) */
    if(!MOM && (cyclic || kt[0][0]>0)){
      printf("Input o:  ");
      scanf("%i",&o);
    }
#endif

#if 0     /* input coupling constants  */
    printf("Input p-truncation convergence:  ");
    scanf("%le",&c);
    printf("Input overall scale:  ");
    scanf("%le",&rescale);
    printf("Input %i couplings:  ",NPARAMS);
    for(i=0; i<NPARAMS; i++)scanf("%le",couplings+i);
#endif 

    fprintf(fp,"\n" PDIGITS " " PDIGITS,c,rescale);
    if(!MOM && (cyclic || kt[0][0]>0))
      fprintf(fp," %i",o);
    fprintf(fp,"\n");
    for(i=0; i<NPARAMS; i++)fprintf(fp,PDIGITS " ",couplings[i]);
    fprintf(fp,"\n");

    /* get the right form of rescale for longitudinal
       string tension */
    rerescale=cyclic?rescale:sqrt(rescale);

  /* make bases */
  for(i=0; i<NWIND; i++)
    for(j=0; j<NS; j++){
      initbasis(s[i]+j,ht[i],np[i][j],kt[i][j],cyclic,o,multi[i]);
      if(MOM || (!cyclic && kt[i][j]<0))
	initbasis(s2[i]+j,ht[i],np[i][j],kt[i][j],cyclic,-o,multi2[i]);
    }

    t1=clock();
#ifdef USE_MPI 
    /* Some global default initialization */
    for(i=0; i<SMAX; i++){
      mpinvals[i]=nval;
      for(j=0; j<HINDEX+2; j++)angle[i][j]=0.0;
    }
    job=0;
    for(i=0; i<NWIND; i++)
      for(k=0; cyclic?(k<1):(k<=NLONG); k++)
	for(j=0; j<NS; job++,j++){
	  if(job==SMAX){
	     fprintf(stderr,__FILE__":  SMAX is too small, exiting\n.");
	    goto cleanup;
	  }
	  mpibase[job][0]=s[i]+j;
	  if(MOM || !cyclic)
	    mpibase[job][1]=s2[i]+j;
	  else
	    mpibase[job][1]=s[i]+j;
	  if(!cyclic){
	    angle[job][HINDEX]=kmax[j];
	    angle[job][HINDEX+1]=(element) k*LONG/(element) NLONG;
	  }
	}
    /************ now distribute the work to the various slaves **********/
    njob=job;
    for(i=0, j=0; i<njob; i++){
      whichvals[i]=j;
      j+=mpinvals[i];
    }
    mpivals=realloc(mpivals,j*sizeof(element));
#if 0 /* debugging print */
    printf("starting subroutine master with %i jobs.\n",njob);
#endif
    master(mpibase,couplings,angle,mpinvals,whichvals,mpivals,njob);
#if 0 /* debugging print */
    printf("finished subroutine master \n");
    for(i=0; i<10; i++)
      printf("element %i, val %f, nvals=%i which=%i\n",i,mpivals[i],
	     mpinvals[i],whichvals[i]);
#endif
    /******** take results and calculate extrapolations etc. **************/
    job=0;
    for(i=0; i<NWIND; i++){
      printf("Starting winding n =");
      for(j=0; j<HINDEX;j++)
	printf(" %i",ht[i][j]);
      printf("\n");
#if MOM
      k=0;
      for(j=0; j<HINDEX; j++)ipt[j]=0;
      while(k<HINDEX){
	for(j=0; j<HINDEX; j++)
	  angle[j]=2.0*pi*(element) ipt[j]/(element) ht[i][j];
	for(j=0; j<HINDEX;j++)
	  fprintf(fp,"%i ",ht[i][j]);
	printf("  Starting angle = (");
	for(j=0; j<HINDEX; j++){
	  printf(" %i/%i",ipt[j],ht[i][j]);
	  fprintf(fp,"  %i",ipt[j]);
	}
	printf(") 2 pi \n");
	fprintf(fp,"\n");
#else
      for(k=0; cyclic?(k<1):(k<=NLONG); k++){
	printf("  %e\n",angle[job][HINDEX+1]);
	for(j=0; j<HINDEX;j++)
	  fprintf(fp,"%i ",ht[i][j]);
	fprintf(fp,"  %e\n",angle[job][HINDEX+1]);
#endif
	  
	if(!cyclic && kt[i][0]<0)
	  vals=lextrapolate1(NS,s[i],mpivals+whichvals[job],
			     th,nval,&c,kmax);
	else if(MOM)
	  vals=extrapolate1(NS,s[i],mpivals+whichvals[job],th,
			    nval,&c);
	else
	  vals=extrapolate1(NS,s[i],mpivals+whichvals[job],
			    th,nval,&c);
	job+=NS;
	for(j=0; j<nval; j++){
	  printf("    %21.14g\n",rerescale*vals[j]);
	  fprintf(fp,"    %21.14g\n",rerescale*vals[j]);
	} 
#if MOM
	/*  Using abs(ht[i][k])/2 gives the first half of the Brillion zone.
	    Removing the 1/2 gives all possible momenta   */
	for(k=0; k<HINDEX && ipt[k]==abs(ht[i][k])/2; k++);
	for(l=0; l<k; l++)ipt[l]=0;
	if(k<HINDEX)ipt[k]++;
	
#endif
      }
      fflush(fp);
      printf("Total of %f CPU seconds to run\n",
	     (float)(clock()-t1)/(float) CLOCKS_PER_SEC);
      }
      
#else  /********  Do normal calculation without MPI ************/
    for(i=0; i<NWIND; i++){
      printf("Starting winding n =");
      for(j=0; j<HINDEX;j++)
	printf("%i ",ht[i][j]);
      printf("\n");
#if MOM
      k=0;
      for(j=0; j<HINDEX; j++)ipt[j]=0;
      while(k<HINDEX){
	for(j=0; j<HINDEX; j++)
	  angle[j]=2.0*pi*(element) ipt[j]/(element) ht[i][j];
	for(j=0; j<HINDEX;j++)
	  fprintf(fp,"%i ",ht[i][j]);
	printf("  Starting angle = (");
	for(j=0; j<HINDEX; j++){
	  printf(" %i/%i",ipt[j],ht[i][j]);
	  fprintf(fp,"  %i",ipt[j]);
	}
	printf(") 2 pi \n");
	fprintf(fp,"\n");
#else
      for(k=0; cyclic?(k<1):(k<=NLONG); k++){
	angle[0]=(element) k*LONG/(element) NLONG;
	printf("  %g\n",angle[0]);
	for(j=0; j<HINDEX;j++)
	  fprintf(fp,"%i ",ht[i][j]);
	fprintf(fp,"  %g\n",angle[0]);
#endif
	
	if(!cyclic && kt[i][0]<0)
	  vals=lextrapolate(NS,s[i],s2[i],th,nval,couplings,&c,kmax,angle[0]);
	else if(MOM)
	  vals=extrapolate(NS,s[i],s2[i],th,nval,couplings,&c,ipt,angle);
	else
	  vals=extrapolate(NS,s[i],s[i],th,nval,couplings,&c,ipt,angle);
	for(j=0; j<nval; j++){
	  printf("    %21.14g\n",rerescale*vals[j]);
	  fprintf(fp,"    %21.14e\n",rerescale*vals[j]);
	} 
#if MOM
	/*  Using abs(ht[i][k])/2 gives the first half of the Brillion zone.
	    Removing the 1/2 gives all possible momenta   */
	for(k=0; k<HINDEX && ipt[k]==abs(ht[i][k])/2; k++);
	for(l=0; l<k; l++)ipt[l]=0;
	if(k<HINDEX)ipt[k]++;
	
#endif
      }
      fflush(fp);
      printf("Total of %f CPU seconds to run\n",
	     (float)(clock()-t1)/(float) CLOCKS_PER_SEC);
      }
#endif /****** end of normal calculation without MPI *********/
	
      /* free bases */

    for(i=0; i<NWIND; i++)
      for(j=0; j<NS; j++){
	freebasis(s[i]+j);
	if(MOM || (!cyclic && kt[i][j]<0))freebasis(s2[i]+j);
      }
	
  } /* end of main loop */

/*********** stop everything ***************/    

cleanup:
      fclose(fp);
#if defined(USE_MPI)
  /************* Tell slaves to shut down. **************/
#if ONE_PROCESSOR_SPECTRUM
      MPI_Comm_size(MPI_COMM_WORLD,&myrank);
  for(i=1; i<myrank; i++)
    /* only the message tag is read by the slaves. */
    MPI_Send(&myrank,1,MPI_INT,i,0,MPI_COMM_WORLD);
#else
  MasterShutDown();
#endif
  MPI_Finalize();             /* cleanup MPI */
#endif
  return 0;
}