#include <stdio.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include "include.h"
#include "run.h"

/** This file runs the lanczos algorithm for finding the lowest
    Eigenvalues of the Hamiltonian with extrapolation.  
    If MOM is set to true, then a list of nonzero
    momenta is plotted.

  The extrapolated eigenvalues for these couplings and nblock=3 and nfig=8
  are (all four sectors)  this is the "second best fit spectrum" from
  the long paper:

  16.26335596755101, 36.15032823702037, 38.78009988371912, 67.33051797970321, 
  82.3795736628269, 27.28544824668094, 50.50349556979569, 67.07730722287375, 
  73.48976280575183, 79.41572800306698, 101.6908998770983, 67.01392597234463, 
  84.4419390872399, 124.3441491973057, 108.8862003981859, 51.84687719898539, 
  86.3274343538608, 94.8972685049293, 92.1662473079009

The new best fit values are:
   16.395323 34.848537 40.969516 66.660976 82.490755
   26.534894 49.663470 65.760432 73.470178 78.317953 99.759305
   66.335519 83.438241 119.072741 111.661451
   51.806386 82.135292 97.611292 95.335727

*/

#define MOM 1 /* flag to do nonzero transverse momentum*/
              /* this is ignored for nonzero L */

#if MOM
#define NA 15 /* number of nonzero angles */
#else
#define NA 0
#endif

int MAIN(){
  clock_t t1;
  int i,k;
  integer nval=4,th[1];
#define NOUT 100 /* Output file name length */
  char out[NOUT]="monster.out",*pout; /*output file name */
  FILE *fp;
  
#if HINDEX==2  /******* values for 3+1 transverse lattice *********/
  element angle[HINDEX];
  int o=1,ht[HINDEX]={0,0},multi=9,cyclic=1;
#if MOM
  int ipt[HINDEX],multi2=8;
#endif
#define NS 4
#if 0  /* test basis */
const static int kt[NS]={10,8,10,12}, np[NS]={4,6,6,4};
#elif 1  /* BIG basis */
  const static int kt[NS]={16,16,20,26}, np[NS]={8,6,6,6};
#endif

/********************** coupling constants **************/

#if 0
  element couplings[NPARAMS]={0.1,1.0,1.0,-0.2,-0.23,10.0,-1.0,5.0,1.0},
			pextrap=-1.174339527,rescale=1.0;
#elif 0
  element couplings[NPARAMS]={0.0517766953, 1.0, 0.7584261311, -0.07623059914, 
     -0.1458991891, 8.889796764, 0.1608104465, 1.657995607, -3.268403913},
			pextrap=-1.174339527,rescale=8.099820052;
#elif 0  /* couplings from best fit in Spring 1999  */
  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, pextrap=-1.764202771;
#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,pextrap=-2.212418668;
#endif
  element *kmax=NULL;

#elif HINDEX==1  /****** values for 2+1 transverse lattice ********/
  element angle[HINDEX]={4.0},pextrap=-1.08;
#if MOM
  int ipt[HINDEX]={0},multi2=2;
#endif

#if 0 /* cutoffs used in paper I */
  int o=1,ht[HINDEX]={0},multi=1,cyclic=1;
#define NS 7
  const static int kt[NS]={20,20,20,22,24,26,28}, np[NS]={4,6,8,6,6,6,6};
  element *kmax=NULL;

#elif 0 /* Test case */
  int o=1,ht[HINDEX]={0},multi=1,cyclic=1;
#define NS 1
  const static int kt[NS]={16}, np[NS]={6};
  element *kmax=NULL;

#elif 0 /* cutoffs used in paper II */
  int o=-1,ht[HINDEX]={0},multi=1,cyclic=1;
#define NS 6
  const static int kt[NS]={18,18,20,20,24,32}, np[NS]={6,8,6,8,6,6};
  element *kmax=NULL;

#elif 0    /*444444444 four particle test case 44444444444 */
  int o=1,ht[HINDEX]={0},multi=1,cyclic=1;
#define NS 7
  const static int kt[NS]={8,16,32,64,100,150,200}, np[NS]={4,4,4,4,4,4,4};
  element *kmax=NULL;

#elif 0  /*LLLLLLLLLLL  Longitudinal string tension, n=0 (BIG basis) LLLLLLL*/
  int o=1,ht[HINDEX]={0},multi=3,cyclic=0;
#define NS 30
  const static int kt[NS]={-24,-34,-50,-100,-200,-300,
			   -24,-34,-50,-100,-200,-300,
			   -24,-34,-50,-100,-200,-300,
			   -24,-34,-50,-100,-200,-300,
			   -24,-34,-50,-100,-200,-300}, 
			    np[NS]={2,2,2,2,2,2,
				    2,2,2,2,2,2,
				    2,2,2,2,2,2,
				    2,2,2,2,2,2,
				    2,2,2,2,2,2};
  element kmax[NS]={4.0,4.0,4.0,4.0,4.0,4.0,
		    5.0,5.0,5.0,5.0,5.0,5.0,
		    7.5,7.5,7.5,7.5,7.5,7.5,
                    10.0,10.0,10.0,10.0,10.0,10.0,
		    15.0,15.0,15.0,15.0,15.0,15.0};
#elif 0  /*LLLLLLLLLLL  Longitudinal string tension, n=0 LLLLLLL*/
  int o=1,ht[HINDEX]={0},multi=3,cyclic=0;
#define NS 20
  const static int kt[NS]={-24,-34,-50,-100,-200,
			   -24,-34,-50,-100,-200,
			   -24,-34,-50,-100,-200,
			   -24,-34,-50,-100,-200}, 
			    np[NS]={2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2};
  element kmax[NS]={4.0,4.0,4.0,4.0,4.0,
		    5.0,5.0,5.0,5.0,5.0,
		    7.5,7.5,7.5,7.5,7.5,
                    10.0,10.0,10.0,10.0,10.0};
#elif 0  /*LLLLLLL Longitudinal string tension, n=1 (BIGGER basis) LLLL*/
  int o=1,ht[HINDEX]={1},multi=0,cyclic=0;
#define NS 30
  const static int kt[NS]={-51,-101,-201,-301,-401,
			   -51,-101,-201,-301,-401,
			   -51,-101,-201,-301,-401,
			   -51,-101,-201,-301,-401,
			   -51,-101,-201,-301,-401,
			   -51,-101,-201,-301,-401}, 
			    np[NS]={1,1,1,1,1,
				    1,1,1,1,1,
				    1,1,1,1,1,
				    1,1,1,1,1,
				    1,1,1,1,1,
				    1,1,1,1,1};
  element kmax[NS]={4.0,4.0,4.0,4.0,4.0,
		    5.0,5.0,5.0,5.0,5.0,
		    7.5,7.5,7.5,7.5,7.5,
                    10.0,10.0,10.0,10.0,10.0,
		    15.0,15.0,15.0,15.0,15.0,
		    25.0,25.0,25.0,25.0,25.0};
#elif 0  /*LLLLLLLLLLL  Longitudinal string tension, n=1 LLLLLLL*/
  int o=1,ht[HINDEX]={1},multi=0,cyclic=0;
#define NS 25
  const static int kt[NS]={-25,-35,-51,-101,-201,
			   -25,-35,-51,-101,-201,
			   -25,-35,-51,-101,-201,
			   -25,-35,-51,-101,-201,
			   -25,-35,-51,-101,-201}, 
			    np[NS]={1,1,1,1,1,
				    1,1,1,1,1,
				    1,1,1,1,1,
				    1,1,1,1,1,
				    1,1,1,1,1};
  element kmax[NS]={4.0,4.0,4.0,4.0,4.0,
		    5.0,5.0,5.0,5.0,5.0,
		    7.5,7.5,7.5,7.5,7.5,
                    10.0,10.0,10.0,10.0,10.0,
		    15.0,15.0,15.0,15.0,15.0};
#elif 0  /*LLLLLLLL Longitudinal string tension (slightly smaller) LLLLLLL*/
  int o=1,ht[HINDEX]={0},multi=3,cyclic=0;
#define NS 15
  const static int kt[NS]={-25,-35,-50,-75,-100,
			   -25,-35,-50,-75,-100,
			   -25,-35,-50,-75,-100}, 
			    np[NS]={2,2,2,2,2,2,2,2,2,2,2,2,2,2,2};
  element kmax[NS]={4.0,4.0,4.0,4.0,4.0,
		    5.0,5.0,5.0,5.0,5.0,
		    8.0,8.0,8.0,8.0,8.0};
#elif 1  /*LLLLLLLLLL Longitudinal string tension, smaller basis LLLLLLL*/
  int o=1,ht[HINDEX]={0},multi=3,cyclic=0;
#define NS 15
  const static int kt[NS]={-20,-28,-34,-38,-40,-20,-28,-34,-38,-40,
			   -20,-28,-34,-38,-40}, 
			    np[NS]={2,2,2,2,2,2,2,2,2,2,2,2,2,2,2};
  element kmax[NS]={4.0,4.0,4.0,4.0,4.0,
		    5.0,5.0,5.0,5.0,5.0,
		    6.0,6.0,6.0,6.0,6.0};

#elif 0  /*LLLLLLLLLL Longitudinal string tension, basis in paper II*/
  int o=1,ht[HINDEX]={0},multi=3,cyclic=0;
#define NS 6
  const static int kt[NS]={-32,-32,-32,-34,-44,-60},np[NS]={2,4,2,2,2,2};
  element kmax[NS]={4,4,5,4,5,4};
#endif

/********************** coupling constants **************/

#if 0
  element couplings[NPARAMS]={0.25,1.0,0.0,0.0,-0.5,-0.5}, rescale=1.0;
#elif 0
  element couplings[NPARAMS]={0.0517767,1.0,0.0,0.0,0.0,0.0}, rescale=1.0;
#elif 0  /* This is the best fit from paper II, old tau */
  element couplings[NPARAMS]={0.032492, 1.000000, -0.077209, -0.250795,
		       221.029589, -0.960551}, rescale=5.221990;
#elif 1  /* This is the best fit from paper II, new tau */
  element couplings[NPARAMS]={0.032492, 1.000000, -0.077209, -0.250795,
		       221.029589, 0.5, 0.5}, rescale=5.221990;
#elif 0
  /* Best fit couplings from paper I
     g^2 term set to 2 for new (more standard) convention.*/
  element couplings[NPARAMS]={0.0649839, 2., -0.201606821, -0.550384068, 
		       2.0*116.657}, rescale=3.46911443596316;
#endif

#endif /******* values for 2+1 transverse lattice *********/

  size_t j;
  element *vals,rescale2;
  static sectors s[NS],s2[NS];
#if MOM
  static element pi;
  pi= atan(1.) * 4.0;
#else
  int *dummy=NULL;
#endif

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

 #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 monster, exiting\n",out);
    exit(55);
  }
#if HINDEX==1
  fprintf(fp,"****** Eigenvalues for the 2+1 transverse lattice  ******\n");
  fprintf(fp,"Couplings:  m^2, G^2 N, lambda_1/a, lambda_2/a,"
	  " lambda_3/a tau\n");
#elif HINDEX==2
  fprintf(fp,"****** Eigenvalues for the 3+1 transverse lattice  ******\n"
	     "Couplings: m^2 G^2 N box la1 la2 la3 la4\n la5 tau\n");
#endif
  for(i=0; i<NPARAMS; i++)fprintf(fp," %f",couplings[i]);
  fprintf(fp,"\nOverall scale G^2 N/sigma = %e\n",rescale);
  fprintf(fp,"Winding number (");
  for(i=0;i<HINDEX;i++)fprintf(fp," %i",ht[i]);
  fprintf(fp,"), o=%i, multiplet=%i\n",o,multi);
#if MOM
  if(cyclic)
    fprintf(fp," along with o=%i and multiplet=%i\n",-o,multi2);
#endif
  if(cyclic){
    fprintf(fp,"Extrapolation using (K,p) =");
    for(i=0, k=3; i<NS; i++){
      fprintf(fp," (%i/2,%i) ",kt[i],np[i]);
      if((k++)%5==0)fprintf(fp,"\n ");
    }
  } else {
    fprintf(fp,"Extrapolation using (K,p,kmax) =");
    for(i=0, k=3; i<NS; i++){
      fprintf(fp," (%i/2,%i,%g) ",kt[i],np[i],kmax[i]);
      if((k++)%5==0)fprintf(fp,"\n ");
    }
  }
  fprintf(fp,"\nEigenvalues are in units of sigma (or sqrt{sigma})\n\n");
  fflush(fp);
  
  /****** Begin main loop of program  **********/

  for(;;){
    
    printf("Input number of eigenvalues (0=exit):  ");
    scanf("%li",&nval);
    if(nval<1)break;
    
/**** Input "multi" and "o," this is useful for doing several spectra  *****/
    
#if MOM && 1 /* input multiplets  (instead of defaults.) */
    if(cyclic && kt[0]>0){
      printf("Input multi, o, multi2:  ");
      scanf("%i %i %i",&multi,&o,&multi2);
    }
#endif

    /***********Constructing bases ********************/
    /* At the end of a previous loop, this was deallocated */
    
    for(i=0; i<NS; i++){
      initbasis(s+i,ht,np[i],kt[i],cyclic,o,multi);
      if(!cyclic)
	initbasis(s2+i,ht,np[i],kt[i],cyclic,-o,multi);
      else if(MOM)
	initbasis(s2+i,ht,np[i],kt[i],cyclic,-o,multi2);
    }
    
#if 0            /* input coupling constants  */
    printf("Input %i couplings:  ",NPARAMS);
    for(i=0; i<NPARAMS; i++)scanf("%le",couplings+i);
    printf("\nInput overall scale (rescale):  ");
    scanf("%le",&rescale);
    printf("\n");
#endif 
    
    for(k=0; k<=NA; k++){
      t1=clock();
      if(!cyclic)
	vals=lextrapolate(NS,s,s2,th,nval,couplings,&pextrap,kmax,angle[0]);
      else {  
#if MOM
	
/**************0***** Set angle ******** ***********/

 	for(i=0;i<HINDEX;i++)angle[i]=0.0;  
	/* 2*pi is the full Brillion zone */
#if 1  /* quadratic spacing of the momenta */
	angle[0]=((element) k*k)/((element) NA*NA)*pi;
#else  /* linear spacing of the momenta  */
	angle[0]=(element) k/(element) NA*pi;
#endif
        if(multi==11||multi==12||multi==13||multi==14){
	   fprintf(stderr,__FILE__":  Using ad hoc method to set angles\n");
	  for(i=1;i<HINDEX;i++)angle[i]=angle[0]; 
	}   
	
	printf("Starting extrapolate angle=");
	for(i=0;i<HINDEX;i++){
	  printf("%f ",angle[i]);
	  fprintf(fp,"%f ",angle[i]);
	}
	printf("\n");
	fprintf(fp,"\n");

	vals=extrapolate(NS,s,s2,th,nval,couplings,&pextrap,ipt,angle);
#else
	vals=extrapolate(NS,s,s,th,nval,couplings,&pextrap,dummy,angle);
#endif
      }
      printf("Total of %f CPU seconds to run extrapolation\n",
	     (float)(clock()-t1)/(float) CLOCKS_PER_SEC);

/****************** Print out results *************/

      if(!cyclic)
	rescale2=sqrt(rescale);
      else 
	rescale2=rescale;
      printf(" c = %e\n",pextrap);
      for(j=0; j<nval; j++){
	printf("     %21.14e\n",vals[j]*rescale2);
	fprintf(fp,"%21.14e\n",vals[j]*rescale2);
      } 
      fprintf(fp,"\n");
      fflush(fp);
    }

  /***********Deconstructing bases ********************/

    for(i=0; i<NS; i++){
      freebasis(s+i);
      if(!cyclic || MOM)freebasis(s2+i);
    }

  }

  fclose(fp);
  return 0;
}