#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "include.h"
#include "run.h"

/*  Check matrix elements of the Hamiltonian and
    diagonalize using standard techniques.
    Results are written to the file "testinteraction.out."  */

#define NV 0 /* Number of eigenvectors to print out into file */
/*  Number of eigenvalues to print out */
#define NEI 10
/*  Number of matrix elements to print out. */
#define NZMAX 20
#define MOM 0 /*Flag to produce nonzero transverse or long. momentum */
/*  The following flag chooses whether a (simpler) complex representation of
    the Hamiltonian should be used.  Valid for heavy sources and nonzero MOM */
#define COMPLEX (0 && MOM)

int MAIN(){
  /*  short s=1;*/
  int i,j,nz;
  integer ierr=0,n,lwork;
  element kmax=5.0,ll=1.0,beta;
#if HINDEX==1
  int ht[HINDEX]={1}, np=3, kt=7, cyclic=1, o=1, multi=0, multi2=2; 
#if 0    /*  Test against Simon */
  element  couplings[NPARAMS]={0.1,1.0,0,0,0,-0.5,-0.5};
#elif 1 /* Best fit from paper II */
  element  couplings[NPARAMS]={0.032492,1.0,-0.077209,-0.250795,
			       221.029589,0.5,0.5};
#elif 0    /* 0.1, 2.0, -0.2, -0.23, 10.0,0.0 */
  element  couplings[NPARAMS]={ 0.1, 2.0, -0.2, -0.23, 20.0,0.0,0.0};
#else                /* 0.1 1 0  0 0 */
  element  couplings[NPARAMS]={0.0324919696232906307,1.0,0.0,0.0,0.0,-1.0};
#endif
#elif HINDEX==2
#if 1  /*  normal spectrum  */
  int ht[HINDEX]={0,0}, np=6, kt=8, cyclic=1, o=1, multi=0, multi2=0; 
#elif 0  /* nonzero winding */
  int ht[HINDEX]={0,0}, np=4, kt=4, cyclic=1, o=1, multi=9; 
#endif
#if 1 /* normal test case */
  element  couplings[NPARAMS]={0.03,1.0,0.7,0.03,-0.005,3.0,0.13,0.25,0.0};
#elif 0 /* Compare to Simon */
  element  couplings[NPARAMS]={ 0.1,1.0,0.2, 0.1,  0.15,0.0,0.15,0.0,0.0,-1.0};
#elif 0 /* beta=0.1 */
  element  couplings[NPARAMS]={0.0324919696232906307,1.0,0.7,
			       0.03,-0.005,3.0,0.13,1.0,-2.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
#endif 
  element **locate;
  sectors x={0};
  char out[]="ti.out",vout[]="vti.out";
  doublereal *vals=NULL;
  FILE *fp, *vfp;
  char jobz='v',uplo='l';
#if MOM
#if HINDEX==1
  element pt[HINDEX]={1.0};
  int ipt[HINDEX]={0};
#elif HINDEX==2
  element pt[HINDEX]={1.0,0};
  int ipt[HINDEX]={0,0};
#endif
#if COMPLEX
  doublereal *rwork=NULL;
  doublecomplex *a=NULL, *work=NULL;
#else
  sectors y={0};
#endif
#endif
#if !COMPLEX
  double *work=NULL;
  element *loci;
#endif

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

  vfp = fopen(vout,"w");
  if((fp = fopen(out,"w")) == NULL){
     fprintf(stderr,__FILE__":  Can't open file %s in testinteraction, exiting\n",out);
    exit(55);
  }
  fprintf(fp,"****** Eigenvalues for the 2+%i transverse lattice  ******\n"
             "If K<0 then abs(K) is an upper limit for momenta;\n"
             "  K_max and L have units of G^2 N.\n",HINDEX);
  fprintf(fp,"cyclic=0 means no symmetry under cyclic permutations\n");
  fprintf(fp,"o=0 means no orientation symmetry\n");
  fprintf(fp,"multi = 0    the trivial subgroup of transverse lattice\n"
             "        1(2) even(odd) full transverse parity\n"
             "        3    singlet representation of the full group\n"
             "      > 3    various other representations of the full group.\n\n");
  /* printf("double epsilon=%e\n",DBL_EPSILON);*/

  /****** Start main loop of program *****************/

  do{

  /******  Input parameters for the basis  *************/
#if 0 /* input parameters for basis (instead of defaults.) */
    printf("Input ht[%i], p truncation, kt=2*K, "
	   "cyclic, o, multiplet:\n",HINDEX);
    for(i=0;i<HINDEX;i++)scanf("%i",ht+i);
    scanf("%i %i %i %i %i",&np,&kt,&cyclic,&o,&multi);
    
    if(kt<0&&!cyclic){
      printf("Input P^+ cutoff K_max and long separation L:\n");
      scanf("%lg %lg",&kmax,&ll);
    }
#if MOM && !COMPLEX
    else {
      printf("Input second multiplet:\n");
      scanf("%i",&multi2);
    }
#endif
#endif
    fprintf(fp,"winding= ");
    for(i=0;i<HINDEX;i++)fprintf(fp,"%i ",ht[i]);
    fprintf(fp,", %i particles, K=%i/2, cyclic=%i, o=%i, multi=%i\n",
	    np,kt,cyclic,o,multi);
    if(kt<0&&!cyclic)
      fprintf(fp,"  K_max=%e L=%e\n",kmax,ll);
#if MOM && !COMPLEX
    else fprintf(fp,"  and o=%i multi2=%i\n",-o,multi2);
#endif     
    
#if 0 /******  Input coupling constants  *************/
    printf("Input %i couplings:  ",NPARAMS);
    for(i=0;i<NPARAMS;i++)scanf("%lg",couplings+i);
#endif
    
#if HINDEX==1
    fprintf(fp,"mu^2=%f G^2 N=%f la1=%f la2=%f la3=%f\n"
	    "  tau=%f\n",
	    couplings[0],couplings[1],couplings[2],
	    couplings[3],couplings[4],couplings[5]);
#elif HINDEX==2
    fprintf(fp,"mu^2=%f G^2 N^2=%f b=%f la1=%f la2=%f\n"
	    "  la3=%f la4=%f la5=%f t=%f\n",
	    couplings[0],couplings[1],couplings[2],
	    couplings[3],couplings[4],
	    couplings[5],couplings[6],couplings[7],couplings[8]);
#endif
    
    /* Run mapcouplings */
    
#if 1
    printf("beta=%f Couplings: ",
	   beta=mapcouplings(kt,couplings));
    for(i=0; i<NPARAMS; i++)printf(" %.14g",couplings[i]);
    printf("\n");
#else
    beta=mapcouplings(kt,couplings);
#endif
  
    /* Make basis or bases */
  
    initbasis(&x,ht,np,kt,cyclic,o,multi);
    makebasis(&x);
#if COMPLEX
    if(o!=0 || multi!=0){
       fprintf(stderr,__FILE__":  Must have zero o and multi for complex, exiting\n");
      exit(12);
    }
#elif MOM
    if(o==0){
       fprintf(stderr,__FILE__":  Must have nonzero o for real, exiting\n");
      exit(12);
    }
    if(kt<0 && !cyclic)
      initbasis(&y,ht,np,kt,cyclic,-o,multi);
    else 
      initbasis(&y,ht,np,kt,cyclic,-o,multi2);
    makebasis(&y);
#endif


    /* Input transverse or longitudinal momenta */
#if MOM
    if(kt<0 && !cyclic)
      maketenslist(ll,kt,kmax,beta);
    else {
      if(!iswinding(x.ht)){
	printf("Input transverse momenta (%i): ",HINDEX);
	for(i=0; i<HINDEX; i++)scanf("%lg",pt+i);
	/*printf("P_1 a = %f\n",pt);*/
	fprintf(fp,"P_1 a = (");
	for(i=0; i<HINDEX; i++)fprintf(fp,"%f,",pt[i]);
	fprintf(fp,"), only o*parity is definite\n");
	for(i=0; i<HINDEX; i++)ipt[i]=0; 
      } else {
	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("P_1 a = %i/n\n",ipt);*/
	fprintf(fp,"P_1 a = (");
	for(i=0; i<HINDEX; i++)fprintf(fp,"%i/%i,",ipt[i],ht[i]);
	fprintf(fp,").\n");
	/*      (    pi        )                                  */
	for(i=0; i<HINDEX; i++)
	  pt[i]=(2.0*(atan(1.) * 4.0)*(element) ipt[i])/((element) ht[i]);
      }
      if(x.multi==0){
	rebasis(&x,ipt,pt);
#if !COMPLEX
	if(o!=0)rebasis(&y,ipt,pt);
#endif
      }
      makephaselist(pt,kt);
    }
#endif
    
#if MOM
#if COMPLEX
    printf("Using complexham\n");
    locate=complexham(&x,couplings);
    n=x.length;
#else
    printf("Using makehamc\n");
    loci=makehamc(&x,&y,couplings);
    n=x.length+y.length;
    locate=&loci;
#endif
#else
    printf("Using makeham\n");
    loci=makeham(&x,couplings);
    n=x.length;
    locate=&loci;
#endif

    printf("First %i nonzero matrix elements (real part):\n",NZMAX);
    for(nz=0, i=0; i<n; i++)
      for(j=i; j<n && nz<NZMAX; j++){
	if(fabs(locate[0][i*n+j])>1.0e-14){
	  nz++;
	  printf("  %i %i %.14g ",i,j,locate[0][i*n+j]);
	  if(!(nz%3)||(i==n-1&&j==n-1))printf("\n");
	}
      }
    printf("\n");
#if COMPLEX
    printf("First %i nonzero matrix elements (imaginary part):\n",NZMAX);
    for(nz=0, i=0; i<n; i++)
      for(j=i; j<n && nz<NZMAX; j++){
	if(fabs(locate[1][i*n+j])>1.0e-14){
	  nz++;
	  printf(" %i %i: %.14g ",i,j,locate[1][i*n+j]);
	  if(!(nz%3))printf("\n");
	}
      }
    printf("\n");
    printf("First %i nonzero phases:\n",NZMAX);
    for(nz=0, i=0; i<n; i++)
      for(j=i; j<n && nz<NZMAX; j++){
	if(fabs(locate[1][i*n+j])>1.0e-14){
	  nz++;
	  printf(" %i %i: %.14g ",i,j,atan2(locate[1][i*n+j],
					    locate[0][i*n+j]));
	  if(!(nz%5))printf("\n");
	}
      }
    printf("\n");
#endif
    for(i=0;i<n;i++)for(j=i;j<n;j++)
      if(fabs(locate[0][j*n+i]-locate[0][i*n+j])>1.0e-12)
	printf("Asymmetric for i=%i, j=%i: %.14g %.14g\n",i,j,
	       locate[0][j*n+i],locate[0][i*n+j] );
#if COMPLEX
    for(i=0;i<n;i++)for(j=i;j<n;j++)
      if(fabs(locate[1][j*n+i]+locate[1][i*n+j])>1.0e-12)
	printf("Imaginary not Antisymmetric for i=%i, j=%i: %.14g %.14g\n",i,j,
	       locate[1][j*n+i],locate[1][i*n+j] );
#endif
    
    if(n>0){
      lwork=4*n;
      vals=realloc(vals,n*sizeof(doublereal));
#if COMPLEX
      work=realloc(work,lwork*sizeof(doublecomplex));
      rwork=realloc(rwork,3*n*sizeof(doublereal));
      a=realloc(a,n*n*sizeof(doublecomplex));
      for(i=0; i<n*n; i++){
	a[i].r=locate[0][i];
	a[i].i=locate[1][i];
      }
      printf("diagonalizing with Lapack zheev\n");
#ifdef _CRAY
      F02HAE(&jobz,&uplo,&n,a,&n,vals,rwork,work,&lwork,&ierr,1,1);
#else
      ZHEEV(&jobz,&uplo,&n,a,&n,vals,work,&lwork,rwork,&ierr,1,1);
#endif
#else
      work=realloc(work,lwork*sizeof(doublereal));
      printf("diagonalizing with lapack dsyev\n");
      DSYEV(&jobz,&uplo,&n,locate[0],&n,vals,work,&lwork,&ierr,1,1);
#endif
      if(ierr!=0) fprintf(stderr,__FILE__":  error = %li in diagonalization\n",ierr);
      printf("%i Eigenvalues:  ",NEI);
      for(i=0; i<n && i<NEI; i++){
	fprintf(fp," %.15g ",vals[i]);
	printf(" %.15g ",vals[i]);
	if(((i+1)%3)==0)fprintf(fp,"\n");
	if(((i+1)%3)==0||i==n-1||i==NEI-1)printf("\n");
      }
      for(i=0; i<NV; i++)
	for(j=0; j<n; j++)
#if COMPLEX
	  fprintf(vfp,"%g %g\n",a[n*i+j].r,a[n*i+j].i);
#else
      fprintf(vfp,"%g\n",locate[0][n*i+j]);
#endif
      fprintf(vfp,"\n");
      fflush(vfp);
    }
    printf("\n");
    fprintf(fp,"  %li states in basis\n\n",n);
    
    printf("continue? (0=no)");
    scanf("%i",&i);
  }while(i);

  return 0;
}