#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include "spectrum.h"
#include "improved.h"  /* needed for global_kt */
#ifdef USE_MPI
#include "usempi.h"
#endif

#define MOM 1     /* Flag for nonzero P_1 or L  */
#define VECTORS 0 /* Flag for finding/printing eigenvectors */

/*
  This file runs the lanczos algorithm for finding the lowest
  Eigenvalues of the Hamiltonian.                        
*/

#ifdef __uxp__   /*  Fujitsu vpp */

extern int fcn(integer *m, integer *npar, doublereal *xp, 
	       doublereal *fvecc, integer *iflag){
  return 0;
}
#endif

int MAIN(int argc, char **argv){
  int i;
  enum state_type type=TYPE_SOURCE;
  specify_calculation x;
  eigensystem eigen=INITIAL_EIGENSYSTEM;
  full_basis basis={0};
  clock_t t1;
  time_t t2,tf;
  element angle[HINDEX+2];

#if HINDEX==2

#if MOM
  /* Profile examples from the end of 1998 use o=-1,multi=9,multi2=8 */
  /*                          6    20        1           15         */
  int     ht[HINDEX]={1,0},np=3,kt=8,charge=0,o=0,multi=15;
  element kmax=3.5,ll=1.265;
#else  /*   1                    24 6   8       1      8   */
  int     ht[HINDEX]={2,0},np=5,kt=4,charge=0,o=0,multi=0;
#endif

#if 0 /* meson test case, gauge only */
  element  couplings[NPARAMS]={0.03250809*0, 1.0, 0, 0, 0,
       0, 0, 0, 0, 0,  -0.323*0, 0,  0.362*0, 0};
#elif 0  /* meson test case, all couplings */
  element  couplings[NPARAMS]={0.03250809, 1.0, 0.688, -0.038, -0.09,
       374.0, 0.158, 143.0, 999, 999, -0.323, 0.162, 0.362, 0};
#elif 0  /* glueball test couplings */
  element  couplings[NPARAMS]={0.03, 1.0, 0.7,0.03,-0.005, 3.0,0.13,1.0,
                               -2.0,1.0,1,1,1,0};
#elif 1  /* May 2002, junk.in */
  element  couplings[NPARAMS]={0.001967542671,1,0.3909839522,-0.01250282378,
			       -0.06421231123,353.2492529,0.07049127522,
			       160.6273835,-0.7123520719,-0.6519355869,
			       0,0,0,0};
    /* In older results, lambda_5 is scaled differently, use 0.25 or 1.0 */
#elif 0  /* older test couplings */
  element  couplings[NPARAMS]={0.1, 1.0, 0.2,0.1,0.15,0.0,0.15,0.0,0.0,0.0};
#elif 0 /* couplings used in first 3+1 paper */
  element  couplings[NPARAMS]={0.051777,1.0,0.759751,-0.075265,-0.139860,
			       8.559161,0.157577,1.847000,-3.702121};
#endif

#elif HINDEX==1
#if MOM
  int   ht[HINDEX]={0},np=8,kt=28,charge=1,o=0,multi=1;
  element kmax=3.0,ll=3.0;
#else /*  1             0     24    8       1      8   */
  int   ht[HINDEX]={0},kt=28,np=8,charge=0,o=1,multi=1;
#endif
#if 1
  element couplings[NPARAMS]={0.1,2.0,-0.2,-0.23,20.0,0.0};
#else
  element couplings[NPARAMS]={0.0324919696232906307,1.0,0.0,0.0,0.0,-2.0};
#endif
#endif
  size_t j;
  integer nval=5;
#if MOM
  int ipt[HINDEX];
#if 0
  element transmom[HINDEX];
#endif
#endif

#ifdef USE_MPI
  int myrank;
  MPI_Init(&argc, &argv);     /* initialize MPI */
  mpi_group=MPI_COMM_WORLD;   /* communicator used by spectrum(...) */
  init_mpi_types();
  MPI_Comm_rank(MPI_COMM_WORLD,&myrank); 
#else
  int myrank=0;  /* dummy variable for the non-parallel case */
#endif

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

  
#if MOM
  for(i=0;i<HINDEX;i++)ipt[i]=0;
#endif

  /******  Input parameters for the basis  *************/
  if(myrank==0){
#if 0   /* input parameters for basis (instead of defaults.) */
    printf("Input ht[%i], p truncation, kt=2*K, "
	   "charge, o, multiplet:\n",HINDEX);
    for(i=0;i<HINDEX;i++)scanf("%i",ht+i);
    scanf("%i %i %i %i %i",&np,&kt,&charge,&o,&multi);
#if MOM
    if(type=TYPE_SOURCE){
      printf("Input P^+ cutoff kmax and long separation ll:\n");
      scanf("%lg %lg",&kmax,&ll);
    }
#endif
#endif
    
    /* print out basis properties */
    printf("winding=");
    for(i=0;i<HINDEX;i++)printf(" %i",ht[i]);
    printf(", %i particles, K=%i/2, charge=%i, o=%i, multi=%i",
	   np,kt,charge,o,multi);
#if MOM
    if(type==TYPE_SOURCE)
      printf("\n  K_max=%e L=%e",kmax,ll);
#endif
    printf("\n");
    
    printf("Input number of eigenvalues:  ");
    scanf("%li",&nval);
    printf("\n");
    printf("Couplings:");
    for(i=0;i<NPARAMS;i++) printf(" %.14g",couplings[i]);
    printf("\n");
  }
   
  /************ Start main loop of program ****************/

  do{

    if(myrank==0){
#if 0 /******  Input coupling constants  *************/
      printf("Input %i couplings:  \n",NPARAMS);
      for(i=0;i<NPARAMS;i++) scanf("%le",couplings+i);
      printf("Couplings:");
      for(i=0;i<NPARAMS;i++) printf(" %.14g",couplings[i]);
      printf("\n");
#endif

#if 0 /****** Input basis parameters *******************/
      printf("Input ht[%i], p truncation, kt=2*K, "
	     "charge, o, multiplet:\n",HINDEX);
      for(i=0;i<HINDEX;i++)scanf("%i",ht+i);
      scanf("%i %i %i %i %i",&np,&kt,&charge,&o,&multi);
#endif
      
#if MOM
      if(type==TYPE_SOURCE){
	angle[HINDEX]=kmax;
	angle[HINDEX+1]=ll;
      } else {    
	/*  check for nonzero winding number */
	if(!iswinding(ht)){
	  printf("Input transverse momenta (%i): ",HINDEX);
	  for(i=0; i<HINDEX; i++){
	    scanf("%lg",angle+i);
#if 0
	    scanf("%lg",transmom+i);
#endif
	  }
	  printf("\n");
	} else {
	  if(multi!=0) fprintf(stderr,__FILE__":  Need multi=0 for nonzero winding\n");
	  printf("Input integer transverse momenta; to (");
	  for(i=0; i<HINDEX; i++)printf("%i,",ht[i]-1);
	  printf("):");
	  for(i=0; i<HINDEX; i++)scanf("%i",ipt+i);
	  printf("\n");
	  for(i=0; i<HINDEX; i++){
#if 0
	    transmom[i]=2.0*(atan(1.) * 4.0)*(element) ipt[i]/(element) ht[i];
#endif
	    angle[i]=2.0*(atan(1.) * 4.0)*(element) ipt[i]/(element) ht[i];
	  }
	}
      }
#endif

      initspectrum(&x,type,ht,np,kt,charge,o,multi,
		   MOM,angle,nval,VECTORS);
    }

    /* The following should probably be built into initspectrum()  */
#ifdef USE_MPI
    MPI_Bcast(&x,1, mpi_specify_calculation,0,MPI_COMM_WORLD);
    MPI_Bcast(&global_kt,1,MPI_INT,0,MPI_COMM_WORLD);
    MPI_Bcast(&couplings,NPARAMS,MPI_DOUBLE,0,MPI_COMM_WORLD);
#endif

    if(myrank==0)printf("Starting calculation of spectrum.\n");
    t1=clock();
    time(&t2);
    spectrum(&x,couplings,&eigen,&basis);
    time(&tf);
    if(myrank==0)printf("%.2f CPU seconds (%li wall) to construct "
			"Hamiltonian and run Lanczos\n",
	   (clock()-t1)/(float) CLOCKS_PER_SEC,tf-t2);
    
    if(eigen.values==NULL){
       fprintf(stderr,__FILE__":  Error in Spectrum, exiting\n");
      break;
    }
    
#if !VECTORS  /* just print out eigenvalues */
    if(myrank==0)  
      for(j=0; j<nval; j++){
	printf("    %21.14e\n",eigen.values[j]);
      } 
#else /* print out some eigenvector info too */
    for(j=0; j<nval; j++){
      doublereal x=0,np[NPMAX+1],temp;
      for(i=0; i<=basis.nptrunc; i++)np[i]=0.0;
      for(i=0; i<basis.length; i++){
	if(eigen.vectors.r!=NULL){
	  x+=temp=pow(eigen.vectors.r[j*basis.length+i],2);
	}
	else {
	  x+=temp=pow(eigen.vectors.c[j*basis.length+i].re,2)+
	    pow(eigen.vectors.c[j*basis.length+i].im,2);
	}
	np[basis.np[i]]+=temp;
      }
      printf("    %21.14e, norm=%g",eigen.values[j],x);
      for(i=0; i<=basis.nptrunc; i++)
	if(np[i]>1.0e-8)printf(", %i: %g",i,np[i]);
      printf("\n");
    } 
#endif

    if(myrank==0){
      printf("continue? (0=no)");
      fflush(stdout);
      scanf("%i",&i);
    }
#ifdef USE_MPI
    MPI_Bcast(&i,1,MPI_INT,0,MPI_COMM_WORLD);
#endif
  }while(i);

  /* print out memory usage (before freeing memory) */
  if(myrank==0)
    print_memory_usage(stdout,&basis);

  /* Free up memory used.  This is mainly for debugging */
  free_eigensystem(&eigen);
  free_spectrum();
  free_basis(&basis);

#ifdef USE_MPI
  /* MasterShutDown();*/
  MPI_Finalize();             /* cleanup MPI */
#endif
  return 0;
}