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

#if HINDEX==2
#define PRINTIT 1 /*  Print timing and convergence information  */

/*
                 Mesons on the transverse lattice

    This program searches coupling constant space for a best fit
    to Poincare symmetries that are broken by the transverse lattice
    regulator.  

    We test for rotational invariance in the transverse plane.  
    In practice, we will choose the x^1 direction and the x^1=x^2 
    directions in order to take advandage of residual reflection 
    symmetries.

    We measure:                                                         multi

    p-truncation convergence for n=(1,0) winding mode                   15
    longitudinal potential for n=(0,0) and various L                     3
    longitudinal eigenvalue for n=(1,0) and one L                       15
    longitudinal eigenvalue for n=(1,1) and one L                       17
    winding mode potential for various n=(c,0)                          15
    winding mode potential for one n=(c,c)                              17
    spectrum for P_perp = (0,0) and (c,0)                               9, 8
                                                                        7, 10
    spectrum for P_perp = (c,c)                                         11, 12
                                                                        14, 13


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

#define FCSQUARED 1 /* Number of c^2 values in each sector and
                    direction to use in fit. This can be up to NPERPVALS */
#define FPARITY 1  /* Number of parity doublets to fit. */
#define NPERP 4 /* Number of J^{P_1 C} sectors for nonzero perp momentum 
		   this must be 4 to cover all multiplets */
#define NPERPVALS 2 /* number of states calculated in each sector */
#define NP1 1 /* Number of nonzero values of P_perp in the (c,0) direction*/
#define NP2 1 /* Number of nonzero values of P_perp in the (c,c) direction*/
              /* This can be 1.  */
#define NSTILL 10 /*Number of sectors for zero transverse momentum spectra
                    This must be 10 to cover all multiplets. */ 
#define NSTILLVALS 2  /*Number of states calculated in each sector*/
                      /*This must be large enough to use for nonzero 
                        momenta */
#define MAXTH 10 /* max number of terms in fit functions */
#define FIX_MASSES 2 /* number of masses to fix to experimental values */
#define FEXTRA 0 /* Fit criteria that would not be normally used. */  
                   /*For instance, demanding 
                     a specific value for the lattice spacing.  */

/* if these are zero, take values from lattice glueball results */

#define NLONG 0  /*Data points used to estimate longitudinal potential */
#define NWIND 0 /* Data points used to estimate winding potential */
#define PNS 0 /* number of sectors in p extrapolation */


#if 0 /* basis simon is using test */
/* This takes about 120 seconds to calculate everything */ 
#define NS 4
const static int kt[NS]={6,8,10,12}, np[NS]={4,4,4,4};
#if PNS<=4 /* check dimensions */
static int pnp1[4]={1,3,5,7}, pkt1[4]={11,11,11,11};
#endif

#elif 1 /* smaller test basis */
/* This takes about 120 seconds to calculate everything */ 
#define NS 3
const static int kt[NS]={12,14,18}, np[NS]={4,4,4};
#if PNS
static int pnp1[PNS]={1,3,5,7}, pkt1[PNS]={11,11,11,11};
#else
static int *pnp1=NULL, *pkt1=NULL;
#endif
#endif  /** Choose different possible bases **/

/* Total number of different bases */
#define N_ALL_BASES ((NSTILL+(NP1+NP2)*NPERP)*NS+PNS)
static specify_calculation all_bases[N_ALL_BASES];
static specify_calculation *s[NP1][NPERP],*t[NP2][NPERP],
  *ps1,*stillbases[NSTILL];

element par[NPARAMS];

static int wpar[NPARAMS],scalemethod;
static doublereal p1a[NP1][HINDEX]={{0.25,0.0}},
  p2a[NP2][HINDEX]={{0.25,0.25}};
static int pht1[HINDEX]={1,0};  /* winding number for p-extrapolation */

/* The error for: 
     c^2 for each state (either direction), 
     roundness of winding mode potential (two data points), and
     roundness of heavy quark potential,
     parity doublets  (fractional error).           */

#if FCSQUARED <= 4  /* check that array is big enough */
static const element c2err[NPERP][4]={{0.25,2.0,2.0,2.0}, {0.25,2.0,2.0,2.0},
		   {0.25,2.0,2.0,2.0}, {0.25,2.0,2.0,2.0}};
#endif

/* These are the parity doublets that are to be compared. */
#if FPARITY<=5  /* This is just the size of the arrays. */
/* relative error for chi^2 fit */
static element parityerr[5]={0.25,0.25,0.25,0.25,0.25};
static struct {
  int charge; 
  int multi; 
  int state;}
whichparity[5][2]={
  {{-1,7,0},{-1,3,0}}, /* components of rho */
  {{1,5,0},{1,6,0}},{{-1,5,0},{-1,6,0}}, /* lowest J_z=2 doublets */
  {{1,7,0},{1,5,0}},{{1,7,0},{1,6,0}}};
#endif

/* Error for extra fit criteria. */

#if FEXTRA<2
doublereal extraerr[1]={0.1};
#endif

/* experimental meson masses, here I include the lowest state
   in each J_z^{P_1,C} sector that is not otherwise determined
   by an O^3 rotational symmetry.  The idea is that spin
   multiplets should be fixed by FPARITY above. */
#if FIX_MASSES<=7   /* This is just the size of the array. */
struct {
  doublereal value;  /* mass^2 in units of the string tension */
  doublereal error;  /* mass^2 error in units of the string tension */
  int charge;
  int multi;
  int state;
} experiment[7]={{0.094,1.0,1,4,0},  /* pion */
                 {3.191,2.0,-1,3,0},  /* J_z=0 of the rho */
                 {4.991,0.6, 1,3,0},  /* a_0(983)  */
		 {7.878,1.0,-1,4,0},  /* J_z=0 of the b_1(1235) */
		 {8.200,1.0, 1,7,0},  /* J_z=1 of the a_1(1260) */
		 {9.000,1.0, 1,5,0},  /* J_z=2 of the a_2(1320) */
		 {14.75,5.0,-1,5,0}}; /* J_z=2 of the rho_2(1690) */
#endif

static integer wth,lth;

/* These are all the quantities returned by fcn */
static struct {doublereal pextrap;
  element wfitting[MAXTH]; element lfitting[MAXTH]; element rescale;
  element spectrumvals1[NP1+1][NPERP][NPERPVALS]; 
  element spectrumvals2[NP2+1][NPERP][NPERPVALS];
  element stillvals[NSTILL][NSTILLVALS]; element c2[2][NPERP][NPERPVALS];
} best;
   
#define CHI_DIM ((NP1+NP2)*NPERP*FCSQUARED+FPARITY+FIX_MASSES+FEXTRA) 
                    /*Number of dimensions of chi^2 fit*/
#define NOUT 100 /* Output file name length */

/*********************************************************************/

int MAIN(int argc, char **argv){
  char out[NOUT],*pout;
  /* maxfev is the maximum number of iterations in minimizing chi^2. */
  integer mm=CHI_DIM,npar=1,maxfev,mode=1,ifail=0,info,nfev,ipvt[NPARAMS];
  /*  epsfcn should be a little less than the Lanczos accuracy. (say 10 times)
      ftol should be a little more than epsfcn
      xtol is the desired precision of the coupling constants  
      factor controls the maximum initial step size.  If the
        starting parameters are near the minimum, this must
        be decreased to maybe 1. 
        (recommended range:  0.1 to 100, default 100). */
  /*  epsfcn is the variable used to determine the step size 
      for the finite difference approximation.  Since the 
      extrapolations magnify any error in the eigenvalues,
      this should be comfortably larger than the errors in
      the lanczos routine.  However, it should also be smaller
      than xtol.                                             */
  doublereal xtol=0.01,gtol=0.0,ftol=5.0e-3,epsfcn=5.0e-5,factor=1.0;
  doublereal vpar[NPARAMS],fvec[CHI_DIM],diag[NPARAMS],fjac[CHI_DIM*NPARAMS];
  doublereal qtf[NPARAMS],wa1[NPARAMS],wa2[NPARAMS],wa3[NPARAMS],wa4[CHI_DIM];
  int count_bases=0;
  element angle[HINDEX+2];
  element chi2;
#if NSTILL==6
  int multistill[NSTILL]={3, 5, 5, 6, 6, 7},
    /*    ostill[NSTILL]={1, 1,-1, 1,-1,-1};*/
    charge_still[NSTILL]={1, 1,-1, 1,-1, 1};
#elif NSTILL==10
  int multistill[NSTILL]={3, 3, 4, 4, 5, 5, 6, 6, 7, 7},
    /* ostill[NSTILL]=   {1,-1, 1,-1, 1,-1, 1,-1, 1,-1};*/
    charge_still[NSTILL]={1,-1, 1,-1, 1,-1, 1,-1,-1, 1};
#elif NSTILL==12
  int multistill[NSTILL]={3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8},
    /*    ostill[NSTILL]={1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1};*/
    charge_still[NSTILL]={1,-1, 1,-1, 1,-1, 1,-1,-1, 1,-1, 1};
#endif
  int i,j,k,l,
    multis[NPERP]={15,15,16,16},
    chargest[NPERP]={1,-1,-1,1},
      multit[NPERP]={17,17,18,18};
  FILE *fp;

#ifdef USE_MPI
  int myrank;
#endif

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

  /* initialize parallelization via MPI when applicable.
     in this case we use a simple master/slave setup to
     distribute the construction of the Hamiltonian and the
     diagonalization of the various bases.      */

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

  /* Choose scale method and open output file */

  printf("Output file name:  ");
  fgets(out,NOUT,stdin);
  if((pout=strchr(out,'\n'))!=NULL)*pout='\0';
  printf("%s\n",out);
  if((fp = fopen(out,"w")) == NULL){
    fprintf(stderr,__FILE__":  Can't open file %s, exiting\n",out);
    goto cleanup;
  }

  printf("Input maximum number of iterations:  ");
  scanf("%li",&maxfev);

  printf("Choose method for determining the overall scale:\n"
         "    0 use scale from glueball calculations\n"
	 "    1 find scale by demanding overall minimum chi^2.\n");
  scanf("%i",&scalemethod);  
  
  /*  Print out the results in the output file **/
  fprintf(fp,"********* Eigenvalues for the 3+1 transverse lattice *********\n");
  fprintf(fp,"Couplings:\n    ");
  for(i=0; i<NPARAMS; i++){
    fprintf(fp,"  %i:%s",i,coupling_string(i));
    if(i%5==4 && i!=NPARAMS-1)fprintf(fp,"\n    ");
  }    
  fprintf(fp,"\nUse chi^2 fit with %i criteria and tolerance %g.\n",
	  CHI_DIM,xtol);
  switch(scalemethod){
  case 0:
    fprintf(fp,"Overall scale from glueball data.\n");
    break;
  case 1:
    fprintf(fp,"Overall scale from minimizing chi^2.\n");
    break;
  default:
    fprintf(stderr,__FILE__":  Incorrect scale method = %i, exiting.\n",
	    scalemethod);
    goto cleanup;
  }
#ifdef KSQUARED
  fprintf(fp,"1/K^2 term used in all extrapolations.\n");
#endif
#if FEXTRA>0
  fprintf(fp,"Includes %i nonstandard fit criteria (such as a \n"
	  " specific lattice spacing).\n",FEXTRA);
#endif
#if FTEPER>0
  fprintf(fp,"Includes a fit to %i eigenvalues from Teper.\n",FTEPER);
#endif
#if FPARITY>0
  fprintf(fp,"Parity doublets (charge,multi,state) with error for difference\n"
	  "  over average.\n");
  for(i=0; i<FPARITY; i++)
    fprintf(fp,"  (%2i,%2i,%i) and (%2i,%2i,%i), error %g\n",
	    whichparity[i][0].charge,whichparity[i][0].multi,
	    whichparity[i][0].state,
	    whichparity[i][1].charge,whichparity[i][1].multi,
	    whichparity[i][1].state,parityerr[i]);
#else
  fprintf(fp,"No parity doublets.\n");
#endif
  fprintf(fp,"Spectrum for P_perp a = (0,0),");
  for(i=0; i<NP1; i++){
    fprintf(fp," (");
    for(j=0; j<HINDEX; j++)fprintf(fp," %g",p1a[i][j]);
    fprintf(fp,")");
  }
  fprintf(fp,"\n  using %i states (multiplet, charge, c^2 error for each) =",
	  NPERPVALS);
  for(i=0; i<NPERP; i++){
    fprintf(fp,"\n  (%2i, %2i,",multis[i],chargest[i]);
    for(j=0; j<FCSQUARED; j++)
      fprintf(fp," %g",c2err[i][j]);
    fprintf(fp,")");
  }
  fprintf(fp,"\nSpectrum for P_perp a = (0,0),");
  for(i=0; i<NP2; i++){
    fprintf(fp," (");
    for(j=0; j<HINDEX; j++)fprintf(fp," %g",p2a[i][j]);
    fprintf(fp,")");
  }
  fprintf(fp,"\n  using %i states (multiplet, charge, c^2 error for each) =",
	  NPERPVALS);
  for(i=0; i<NPERP; i++){
    fprintf(fp,"\n  (%2i, %2i,",multit[i],chargest[i]);
    for(j=0; j<FCSQUARED; j++)
      fprintf(fp," %g",c2err[i][j]);
    fprintf(fp,")");
  }
#if NSTILL>0
  fprintf(fp,".\nOrdinary spectra multiplets, %i states, (charge,multiplet) =",
	  NSTILLVALS);
  for(i=0, k=5; i<NSTILL; k++, i++){
    if(k%6==0)fprintf(fp,"\n ");
    fprintf(fp," (%2i, %2i) ",charge_still[i],multistill[i]);
  }
#endif
  fprintf(fp,".\nAll spectra extrapolated using (K,p) =");
  for(i=0, k=3; i<NS; k++, i++){
    if(k%5==0)fprintf(fp,"\n ");
    fprintf(fp," (%i/2,%i) ",kt[i],np[i]);
  }
#if PNS
  fprintf(fp,"p-extrapolation using n=(");
  for(l=0; l<HINDEX; l++)fprintf(fp," %i",pht1[l]);
  fprintf(fp,") and (K,p) =");
  for(j=0, k=4; j<PNS; k++, j++){
    if(k%5==0)fprintf(fp,"\n ");
    fprintf(fp," (%i/2,%i) ",pkt1[j],pnp1[j]);
  }
#endif
  fprintf(fp,".\nResult format:\n"
	  " Fit info, # steps, chi^2, and p damping, scale G^2 N/sigma.\n"
	  " The %i couplings  (G^2 N units);\n"
	  "   which couplings -- if any -- were fit.\n"
	  " Winding potential and heavy souce potential fits.\n"
	  " The rescaled spectrum for each P_perp*a and c^2 values.\n"
#if NSTILL>0
	  " Finally come the states for the ordinary spectra.\n"
#endif
	  "\n",NPARAMS);
  fflush(fp);

  /******************* make basis structures  ****************************/
  for(k=0; k<HINDEX+2; k++)angle[k]=0;

  printf("Make %i spectrum bases,",(NP1+NP2)*NPERP*NS);
  /*  Calculate the extrapolated spectrum in the P_perp = (c,0) direction */
  for(l=0; l<NP1; l++){
    for(k=0; k<HINDEX; k++)angle[k]=p1a[l][k];
    for(j=0; j<NPERP; j++){
      s[l][j]=all_bases+count_bases;
      for(i=0; i<NS; i++)
	initspectrum(all_bases+count_bases++,TYPE_MESON,NULL,np[i],kt[i],
		     chargest[j],0,multis[j],1,angle,NPERPVALS+1,0);
    }
  }
  /*  Calculate the extrapolated spectrum in the P_perp = (c,c) direction */
  for(l=0; l<NP2; l++){
    for(k=0; k<HINDEX; k++)angle[k]=p2a[l][k];
    for(j=0; j<NPERP; j++){
      t[l][j]=all_bases+count_bases;
      for(i=0; i<NS; i++)
	initspectrum(all_bases+count_bases++,TYPE_MESON,NULL,np[i],kt[i],
		     chargest[j],0,multit[j],1,angle,NPERPVALS+1,0);
    }
  }
 
  /*  Calculate the ordinary extrapolated spectra. */
  printf(" %i ordinary spectrum bases,",NSTILL);
  for(j=0; j<NSTILL; j++){
    stillbases[j]=all_bases+count_bases;
    for(i=0; i<NS; i++)
      initspectrum(all_bases+count_bases++,TYPE_MESON,NULL,np[i],kt[i],
		   charge_still[j],0,multistill[j],0,NULL,NSTILLVALS,0);
  }

  /* Calculate the convergence in p-truncation */
  printf(" and %i n=%i bases",PNS,*pht1);
  ps1=all_bases+count_bases;
  for(j=0; j<PNS; j++)
    /* assume o=1, multi=15 is lowest state */
    initspectrum(all_bases+count_bases++,TYPE_GLUE,pht1,pnp1[j],pkt1[j],
		 0,1,15,0,NULL,1,0);  
  printf("\n");
  fflush(stdout);

  if(count_bases!=N_ALL_BASES){
    fprintf(stderr,__FILE__":%i  Number of bases don't match\n",__LINE__);
    goto cleanup;
  }


  /******************* Main loop of program ********************************/

  for(;;){
    printf("How many to fit (0 to %i, -1=stop)\n",NPARAMS-2);
    scanf("%li",&npar);
    if(npar<0)break;
#if 1 /* Input starting parameters */
    printf("Input %i starting parameters\n",NPARAMS);
    for(i=0; i<NPARAMS; i++)scanf("%lf",par+i);
#endif

    if(npar==0){
      info=0; nfev=1;
      fcn(&mm,&npar,vpar,fvec,&info);
    }
    else{
      printf("Input which %li couplings (0=m^2, ...)\n",npar);
      for(i=0;i<npar;i++)scanf("%i",wpar+i);
      for(i=0;i<npar;i++)vpar[i]=par[wpar[i]];

/*    Least squares fitting routine    
      subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,
 *                  maxfev,epsfcn,diag,mode,factor,nprint,info,nfev,
 *                 fjac,ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4)  */
      LMDIF(FCN_TO_FORTRAN,&mm,&npar,vpar,fvec,&ftol,&xtol,&gtol,
            &maxfev,&epsfcn,diag,&mode,&factor,&ifail,&info,&nfev,
            fjac,&mm,ipvt,qtf,wa1,wa2,wa3,wa4);

      if(info<1||info>4) fprintf(stderr,__FILE__":  Error in lmdif info=%li\n",info);
      if(info==0)goto cleanup;
      /*Map best values into par[] */
      for(i=0; i<npar; i++)par[wpar[i]]=vpar[i];
    }
    if(info<0){
       fprintf(stderr,__FILE__":  Error in fcn info=%li\n",info);
      if(info<=-8)goto cleanup;
    }
    chi2=fdot(fvec,fvec,mm);
    printf("error %li, %li steps chi2=%f with c_p=%f rescale=%f\n",
	   info,nfev,chi2,best.pextrap,best.rescale);
    fprintf(fp,"%li %li " PDIGITS " " PDIGITS " " PDIGITS "\n",
	    info,nfev,chi2,best.pextrap,best.rescale);
    printf("Final couplings: ");
    for(i=0; i<NPARAMS; i++){
      printf(" " PDIGITS,par[i]); fprintf(fp,PDIGITS " ",par[i]);
    }
    if(npar>0){fprintf(fp,"\n");printf("\nwhile fitting");}
    for(i=0; i<npar; i++){
      printf(" %i",wpar[i]); fprintf(fp," %i",wpar[i]);
    }
    printf("\nWinding: "); fprintf(fp,"\n");
    for(i=0; i<wth; i++){
      printf("%f ",best.wfitting[i]); fprintf(fp,"%f ",best.wfitting[i]);
    }
    printf("and long:");
    for(i=0; i<lth; i++){
      printf(" %f",best.lfitting[i]); fprintf(fp," %f",best.lfitting[i]);
    }
    printf("\n");
    fprintf(fp,"\n");
    for(i=0; i<NPERP; i++)
      for(k=0; k<NPERPVALS; k++){
	for(j=0; j<NP1+1; j++){
	  printf(" %f",best.rescale*best.spectrumvals1[j][i][k]); 
	  fprintf(fp," %f",best.rescale*best.spectrumvals1[j][i][k]);
	}
	printf(" %f\n",best.c2[0][i][k]); fprintf(fp," %f\n",best.c2[0][i][k]); 
      }
    for(i=0; i<NPERP; i++)
      for(k=0; k<NPERPVALS; k++){
	for(j=0; j<NP2+1; j++){
	  printf(" %f",best.rescale*best.spectrumvals2[j][i][k]); 
	  fprintf(fp," %f",best.rescale*best.spectrumvals2[j][i][k]);
	}
	printf(" %f\n",best.c2[1][i][k]); fprintf(fp," %f\n",best.c2[1][i][k]); 
      }
    for(i=0; i<NSTILL; i++){
      printf("charge=% i multi=%i ",charge_still[i],multistill[i]);
      for(j=0; j<NSTILLVALS; j++){
	printf("%f ",best.rescale*best.stillvals[i][j]); 
	fprintf(fp,"%f ",best.rescale*best.stillvals[i][j]);
      }
      fprintf(fp,"\n"); printf("\n");
    }
    fprintf(fp,"\n"); printf("\n");
    fflush(fp);
    fflush(stdout);
  }
  
  
cleanup:
  fclose(fp);

  /* Free up memory used.  This is mainly for debugging */
  free_spectrum();

#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=0 is read by the slaves. */
    MPI_Send(&myrank,1,MPI_INT,i,0,MPI_COMM_WORLD);
#else
  MasterShutDown();
#endif
  MPI_Finalize();             /* cleanup MPI */
#else
  printf("Maximum matrix size: n=" SIZE_MEM_PRINT ", nonzero=" SIZE_MEM_PRINT 
	 ", " SIZE_MEM_PRINT " bytes\n",
	 max_els.n,max_els.nonzero,max_els.memory);
  printf("Memory used to construct bases: " SIZE_MEM_PRINT " bytes\n",
	 reduced_basis_memory);
#endif

  return 0;
}


/*****************************************************************

  The following routine calculates the fit criteria for
  a given set of fit parameters.  It returns the vector
  fvecc which is then squared to give the chi^2.

  The values that are calculated are reused from previous
  iterations if the relevant couplings are unchanged.
  This includes:

         pextrap, wfitting, lfitting, wvalue1, lvalue1,
	 lvalue2, spectrumvals1, spectrumvals2, stillvals
         
  Information about the spectra are returned in the structure
  best.  The members of the structure best 
  are from the lowest chi^2 iteration of fcn (which is
  not necessarily the last iteration of fcn.)

  On error, fcn sets *iflag to a negative integer.
  *flag<=-8 are "fatal" errors.

*******************************************************************/


extern fsub fcn(integer *m, integer *npar, doublereal *xp, 
	       doublereal *fvecc, integer *iflag){
  /* fitc is the index for the fit criterion */
  int i,k,j,l,fitc;
  integer th[1];
  element *vals=NULL; /* this is realloc'ed by the extrapolation routines */
  element sum,temp,rescale,c2[2][NPERP][NPERPVALS];
  static element oldpar[NPARAMS],bestpar[NPARAMS],bestchi2;
  static doublereal spectrumvals1[NP1+1][NPERP][NPERPVALS],
    spectrumvals2[NP2+1][NPERP][NPERPVALS],
    stillvals[NSTILL][NSTILLVALS];
  static element pextrap;
  static element wfitting[MAXTH]={0.181946,-0.805984,0.689809};
  static element lfitting[MAXTH]={0.092962,0.211746,-0.151139};
  
  clock_t tif;

  int job;
  element *all_vals[N_ALL_BASES];
#if defined(USE_MPI) && ONE_PROCESSOR_SPECTRUM
#else
  full_basis basis={0};
  eigensystem eigen=INITIAL_EIGENSYSTEM;
#endif
  tif=clock();

  /* bestpar is used later to see if unfitted variables have changed*/
  for(i=0; i<*npar; i++)bestpar[wpar[i]]=par[wpar[i]]=xp[i];

#if PRINTIT
  printf("  Parameters: ");
  for(i=0; i<NPARAMS; i++)printf(" %.10g",par[i]);
  printf("\n");
#endif


  /* allocate memory to store eigenvalues */

  for(i=0; i<N_ALL_BASES; i++)
    all_vals[i]=malloc(all_bases[i].nval*sizeof(element));


#if defined(USE_MPI) && ONE_PROCESSOR_SPECTRUM

  /************ now distribute the work to the various slaves **********/

#if 0 /* debugging print */
  printf("starting subroutine master with %i jobs.\n",njob);
#endif
  if(master(mpibase,par,angle,mpinvals,all_vals,mpivals,njob)){
     fprintf(stderr,__FILE__":  master() is returning an error, exiting\n");
    *iflag=-9;
    return SUBRETURN;
  }
#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],all_vals[i]);
#endif
#else  /****** calculate everything without using MPI **************/

  for(job=0; job<N_ALL_BASES; job++){
    spectrum(all_bases+job,par,&eigen,&basis);
    for(i=0; i<all_bases[job].nval; i++)
      all_vals[job][i]=eigen.values[i];
  }
#endif
  /******** take results and calculate extrapolations etc. **************/

  /* Calculate the convergence in p-truncation */
#if PNS
  pextrap=pextrapolate(PNS,ps1,all_vals+(ps1-all_bases));
#else
  pextrap=-1.764202771;
#endif

  /*  Calculate the transverse string tension.  The wht matrix
      must be reformatted into a form compatable with
      the FORTRAN routince WWEIGHTS.    */
#if NWIND
  wweights(weights,&wth,nwind,wht0);
  PSEUDOINVERSE(weights,&wth,&nwind);
  for(j=0; j<wth; j++)wfitting[j]=0.0;
  for(i=0; i<nwind; i++){
    extrapolate(&vals,WNS,ws0[i],all_vals+(ws0[i]-all_bases),th,1,&pextrap);
    /* winding n=(2,0) eigenvalue saved */ 
    if(dot(wht0[i],wht0[i],HINDEX)==4)wvalue20=vals[0];
    for(j=0; j<wth; j++)
      wfitting[j] += vals[0]*weights[i+j*nwind];
  }
#else
  wth=3;
#endif


  /*  Calculate the longitudinal string tension  */
#if NLONG
  lweights(weights,&lth,nlong,mlong0);
  PSEUDOINVERSE(weights,&lth,&nlong);
  for(j=0; j<lth; j++)lfitting[j]=0.0;
  for(i=0; i<nlong; i++){
    lextrapolate(&vals,LNS,ls0[i],all_vals+(ls0[i]-all_bases),
		  th,1,&pextrap,lkmax0);
    for(j=0; j<lth; j++)
      lfitting[j]+=vals[0]*weights[i+j*nlong];
    }
#else
  lth=3;
#endif

  /*  Calculate the ordinary extrapolated spectra. */
  for(i=0; i<NSTILL; i++){
    extrapolate(&vals,NS,stillbases[i],all_vals+(stillbases[i]-all_bases),
		 th,NSTILLVALS,&pextrap);
    for(k=0; k<NSTILLVALS; k++)stillvals[i][k]=vals[k];
  }
  
  /*  Calculate the extrapolated spectrum in the P_perp = (c,0) direction */
  /*  Sort states before adding (and calculate 1 extra value) */
  for(l=0; l<NP1; l++)
    for(i=0; i<NPERP; i++){
      extrapolate(&vals,NS,s[l][i],all_vals+(s[l][i]-all_bases),
		   th,NPERPVALS+1,&pextrap);
      /* sort and add states: we assume NPERPVALS+1 eigenvalues */
      k=0;
      addterms(spectrumvals1[l+1][i],&k,NPERPVALS,vals,NPERPVALS+1);
#if 0 /* debugging print */
      printf("  l=%i, spectrum %i P_perp=(%f,%f)\n",l,i,p1a[l][0],p1a[l][1]);
#endif
    } 
  
  /*  Calculate the extrapolated spectrum in the P_perp = (c,c) direction */
  /*  Sort states before adding (and calculate 1 extra value) */
  for(l=0; l<NP2; l++)
    for(i=0; i<NPERP; i++){
      extrapolate(&vals,NS,t[l][i],all_vals+(t[l][i]-all_bases),
		   th,NPERPVALS+1,&pextrap);
      /* sort and add states: we assume NPERPVALS+1 eigenvalues */
      k=0;
      addterms(spectrumvals2[l+1][i],&k,NPERPVALS,vals,NPERPVALS+1);
#if 0 /* debugging print */
      printf("  l=%i spectrum %i P_perp=(%f,%f)\n",l,i,p2a[l][0],p2a[l][1]);
#endif
    }

  /* Free up memory used.  This is mainly for debugging */
  free_eigensystem(&eigen);
  free_spectrum();
  free_basis(&basis);
  for(i=0; i<N_ALL_BASES; i++)
    free(all_vals[i]);


  /*  Put ordinary spectrum values into nonzero P_perp array*/
  /*  Here we sort the states before adding.  */
  /*  The ordering of the states must be compatable with */
  /*  the ordering obtained when calculating nonzero p_perp */
  for(i=0; i<NPERP; i++){
    for(j=0, k=l=0; j<NSTILL; j++){
      if(stillbases[j]->charge!=s[0][i]->charge)continue;
      /* assume that the two spin 1 multiplets are degenerate */
      if(((stillbases[j]->multi==3 || stillbases[j]->multi==5 ||
	    stillbases[j]->multi==7) && s[0][i]->multi==15) ||
	  ((stillbases[j]->multi==4 || stillbases[j]->multi==6 ||
	    stillbases[j]->multi==7) && s[0][i]->multi==16))
	addterms(spectrumvals1[0][i],&k,NPERPVALS,stillvals[j],NSTILLVALS);
     if(((stillbases[j]->multi==3 || stillbases[j]->multi==6 ||
	    stillbases[j]->multi==7) && t[0][i]->multi==17) ||
	  ((stillbases[j]->multi==4 || stillbases[j]->multi==5 ||
	    stillbases[j]->multi==7) && t[0][i]->multi==18))
	addterms(spectrumvals2[0][i],&l,NPERPVALS,stillvals[j],NSTILLVALS);
    }
    if(k!=NPERPVALS || l!=NPERPVALS){
       fprintf(stderr,__FILE__":  Incorrect multiplet mapping from states "
	       "in ordinary\n"
	      "spectrum to those in nonzero P_perp, exiting\n");
       fprintf(stderr,"          k=%i, l=%i and NPERPVAL=%i\n",
	       k,l,NPERPVALS);
      *iflag=-12;
      return SUBRETURN;
    }	
  }
    
#if 1 /*debugging print of the entire spectrum */
  for(k=0; k<NPERPVALS; k++)
  for(i=0; i<NPERP; i++){
    for(j=0; j<NP1+1; j++)
      printf(" %g",spectrumvals1[j][i][k]);
    printf(",  ");
    for(j=0; j<NP1+1; j++)
      printf(" %g",spectrumvals2[j][i][k]);
    printf("\n");
  } 
#endif
    
  /*  copy couplings into oldparams */
  for(i=0; i<NPARAMS; i++)oldpar[i]=par[i];
  
  /* There are several possible methods for setting
     the overall scale of the Hamiltonian  */
  
  switch(scalemethod){         
  case 0: /* Get scale from glueball results */
    rescale=11.4617271;
    break;
  case 1:
 /* Find scale by minimizing overall chi^2 */
    /* This does not include the consistancy of the L nonzero,
       n=1 datum because it is too complicated */
    /* initial guess is from glueball data */
    rescale=11.4617271;
    for(j=0; ; j++){
      element sum=0.0, zz=0.0;
      for(j=0; j<NPERP; j++)
	for(i=0; i<FCSQUARED && i<NPERPVALS; i++){
	  c2[0][j][i]=wfitting[0]*pow(rescale,2)*
	    (spectrumvals1[1][j][i]-spectrumvals1[0][j][i])/
	    fdot(p1a[0],p1a[0],HINDEX);
	  temp=2.0*c2[0][j][i]/rescale;
	  sum += temp*(c2[0][j][i]-1.0)/pow(c2err[j][i],2);
	  zz +=  pow(temp/c2err[j][i],2);
	  
	  c2[1][j][i]=wfitting[0]*pow(rescale,2)*
	    (spectrumvals2[1][j][i]-spectrumvals2[0][j][i])/
	    fdot(p2a[0],p2a[0],HINDEX);
	  temp=2.0*c2[1][j][i]/rescale;
	  sum += temp*(c2[1][j][i]-1.0)/pow(c2err[j][i],2);
	  zz +=  pow(temp/c2err[j][i],2);
	}
#if FTEPER==1
      sum += stillvals[0][0]*(stillvals[0][0]*rescale-teper[0])/
	pow(tepererr[0],2);
      zz += pow(stillvals[0][0]/tepererr[0],2);
#endif
  /* Add fit to the experimental value for the mesons */
  for(j=0; j<FIX_MASSES; j++)
    for(i=0; i<NSTILL; i++)
      if(stillbases[i]->multi==experiment[j].multi && 
	 stillbases[i]->charge==experiment[j].charge){
	sum += stillvals[i][experiment[j].state]*
	  (rescale*stillvals[i][experiment[j].state]-
		       experiment[j].value)/pow(experiment[j].error,2);
	zz += pow(stillvals[i][experiment[j].state]/experiment[j].error,2);
	break;
      }
#if FSCALE
      sum += lfitting[0]*(rescale*lfitting[0]-1.0)/pow(lerr0,2);
      zz += pow(lfitting[0]/lerr0,2);
#endif
      temp=rescale;
      /* Here I include a damping factor which will
	 slow convergence but hopefully stablize it. */
      rescale -= 0.25*sum/zz;
#if 0 /* debugging print */
      if(j<100)printf("rescale:  j=%i rescale=%.16g\n",j,rescale);
#endif
      /* converge as well as we can, assuming round-off errors */
      if(fabs(temp-rescale)<=(2.0+2.0*FCSQUARED)*ELEMENT_EPSILON)break;
      if(j==1000){
	 fprintf(stderr,__FILE__":  Scale determination not converged, "
		 "rescale=%g\n",rescale);
	*iflag=-2;
	break;
      }
    }
    break;
  default:
    fprintf(stderr,__FILE__":  in fcn, scale method=%i is wrong.\n",
	    scalemethod);
    *iflag=-14;
    rescale=0.123456789;  /*Find scale by using glueball results */
  }
  /* If rescale is negative for some reason, 
     set it to something large and positive.        */
  if(rescale<0.0){
    rescale=50.0;
     fprintf(stderr,__FILE__":  rescale is negative, set recale=%f.\n",
	     rescale);
  }

  /*  Calculate the desired dispersion using the overall scale 
      and lattice spacing from glueball calculations.
      Find error of speed of light squared
      with coupling g^2 N and "a" as determined above.
      One can weight the value for our value of sigma. */

  fitc=0;
  temp=wfitting[0]*pow(rescale,2)/fdot(p1a[0],p1a[0],HINDEX);
  for(i=0; i<NPERP; i++)
    for(j=0; j<NPERPVALS; j++)
    c2[0][i][j] = (spectrumvals1[1][i][j]-spectrumvals1[0][i][j])*temp;
  temp=wfitting[0]*pow(rescale,2.0)/fdot(p2a[0],p2a[0],HINDEX);
  for(i=0; i<NPERP; i++)
    for(j=0; j<NPERPVALS; j++)
    c2[1][i][j] = (spectrumvals2[1][i][j]-spectrumvals2[0][i][j])*temp;
  for(k=0; k<2; k++)
    for(i=0; i<NPERP; i++)
      for(j=0; j<FCSQUARED &&  j<NPERPVALS; j++)
      fvecc[fitc++] = (c2[k][i][j]-1.0)/c2err[i][j];
  
  /* If the 0^++ or the n=(2,0) is tachyonic, modify the 
     0^++ c^2 errors to reflect this */

  for(i=0; i<NSTILL; i++)
    if(stillvals[i][0]<0.0){
      fprintf(stderr,__FILE__":  Lowest state gone tachyonc.\n");
      fvecc[i*NSTILLVALS]=100.0*fabs(stillvals[i][0]);
    }

  /* Add parity doublets to fitting criterion using fractional
     difference.   */

  for(k=0; k<FPARITY; k++){
    for(i=0; i<NSTILL; i++)
      if(stillbases[i]->multi==whichparity[k][0].multi && 
	 stillbases[i]->charge==whichparity[k][0].charge)break;
    for(j=0; j<NSTILL; j++)
      if(stillbases[j]->multi==whichparity[k][1].multi && 
	 stillbases[j]->charge==whichparity[k][1].charge)break;
    if(i<NSTILL && j<NSTILL)
      fvecc[fitc++]=(stillvals[i][whichparity[k][0].state]-
		     stillvals[j][whichparity[k][1].state])/
	(0.5*(stillvals[i][whichparity[k][0].state]+
	      stillvals[j][whichparity[k][1].state])*parityerr[k]);
    else {
      fprintf(stderr,__FILE__":  Incorrect multiplets in ordinary spectra, "
	      " exiting\n");
      *iflag=-15;
      return SUBRETURN;
    }
  }

  /* Add fit to the experimental value for the mesons */
  for(j=0; j<FIX_MASSES; j++)
    for(i=0; i<NSTILL; i++)
      if(stillbases[i]->multi==experiment[j].multi && 
	 stillbases[i]->charge==experiment[j].charge){
	fvecc[fitc++]=(rescale*stillvals[i][experiment[j].state]-
		       experiment[j].value)/experiment[j].error;
	break;
      }

    /* Demand that the lattice spacing be a specific value. */
#if FEXTRA==1
    fvecc[fitc++]=(1.0+2.0*sqrt(par[0]/par[1])-wfitting[0]*rescale)/
      extraerr[0];
#endif
    
  /* after including all fit criteria, make sure we 
     have actually filled the fvecc vector. */
  if(fitc!=*m){
    fprintf(stderr,__FILE__":  Fit criteria %i wrong, should be %li\n",
	    fitc,*m);
    *iflag=-75;
    return SUBRETURN;
  }
    
    /* check if this is the best chi^2 so far,
       if so, copy values into the structure best.
       bestpar makes sure that the unfitted variables are
       the same; that is the previous best values are
       from the same fit. */
  sum=fdot(fvecc,fvecc,*m);
  for(i=0; i<NPARAMS && fabs(bestpar[i]-par[i])<=2*ELEMENT_EPSILON; i++);
  if(*iflag>=0 && (i<NPARAMS || sum<bestchi2)){
    bestchi2=sum;
    for(i=0; i<NPARAMS; i++)bestpar[i]=par[i];
    best.pextrap=pextrap; 
    best.rescale=rescale;
    for(i=0; i<NPERP; i++)
      for(j=0; j<NPERPVALS; j++)
	for(k=0; k<2; k++)
	  best.c2[k][i][j]=c2[k][i][j];
    for(i=0; i<NPERP; i++)
      for(k=0; k<NPERPVALS; k++){
	for(j=0; j<NP1+1; j++)
	  best.spectrumvals1[j][i][k]=spectrumvals1[j][i][k];
	for(j=0; j<NP2+1; j++)
	  best.spectrumvals2[j][i][k]=spectrumvals2[j][i][k];
      }
    for(i=0; i<NSTILL; i++)
      for(j=0; j<NSTILLVALS; j++)
	best.stillvals[i][j]=stillvals[i][j];
    for(i=0; i<wth; i++)best.wfitting[i]=wfitting[i];
    for(i=0; i<lth; i++)best.lfitting[i]=lfitting[i];
  }

#if PRINTIT
  printf("  fvecc values =");
  for(i=0, k=3; i<*m; i++){
    printf(" "PDIGITS,fvecc[i]);
    if((k++)%6==0 && i+1<*m)printf("\n  ");
  }
  printf("\n  Chi^2="PDIGITS" in %f CPU seconds.\n\n",
	 fdot(fvecc,fvecc,*m),
	 (float)(clock()-tif)/(float) CLOCKS_PER_SEC);
#endif
#if PRINTIT 
  fflush(stdout);
#endif

  free(vals); /* vals is realloc()ed by extrapolation routines */

  return SUBRETURN; /* normal return from subroutine */

}

#endif /* HINDEX==2 */