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

#if defined(USE_MPI)
#include "usempi.h"
#endif

#if defined(USE_MPI) && !(ONE_PROCESSOR_SPECTRUM)

#define PDEBUGGER 0  /* debugging macro for OP */
#define COU 0    /* wait for COU calls before printing */


#define PRINTIT 0 /*Print values, method, time, et cetera in spectrum(...)*/
#define MYPRINT 0 /*Print various things associated with mpi */

#define DOP2 1    /*Using plandr2 if 0 does plandr*/
#define FINDVEC 1 /*Finds eigenvectors if 1, only used with plandr*/

/*
   MAX_ELS is a guess for the maximum number of nonzero
   matrix elements of the the Hamiltonian.
   If we define this carefully, memory fragmentation on the heap
   can be minimized. 
*/

#ifndef MAX_ELS
#define MAX_ELS 100000
#endif

#ifndef MAX_WORK
#define MAX_WORK 1000000
#endif

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


/* list of nonzero matrix elements used by elop and spectrum. */

static elslist els={NULL,NULL,NULL,0,0,0};

static doublereal *spectrumval=NULL, *spectrumvec=NULL;
static size_t lval=0, lvec=0;

/* arrays used by OP defined in PLANDR2 */

static doublereal *given_vector=NULL, *tempsum=NULL;
static int *recvsendarray=NULL, *displs=NULL;
static doublereal *work=NULL;
static size_t maxwork=0;

void spectrum2(sectors *s,         /* basis, filled by initbasis(...) */ 
	      sectors *s2,        /* second basis for nonzero momentum */
	      doublereal **vals,  /* list of eigenvalues */
	      doublereal **vecin, /* list of eigenvectors */
	      integer nval,       /* number of eigenvalues */
	      element *couplings, /* coupling constants */
	      element *angle,     /* angles and momenta, length HINDEX+2 */
	      int globalkt,       /* maximum value of kt, for mapcouplings */
	      integer nfig        /* number of significant digits in answer */
	      )
{
  int i,j;
  element beta;
  integer n,ierr=0, nperm=0,bet;
  doublereal *work=NULL;
  doublereal *Bnd = NULL;

  /* maxj determines how much memory dnlaso uses.
     The minimum value is 6*nblock and is faster if
     larger.  The limiting factor should be the amount
     of memory available.  */

  integer maxj, maxmaxj = 300;
  doublereal temp;
  clock_t t1;

  /*PLANDR2 variables*/
  /*G&L variables*/ 
  integer localn;

#if DOP2
  integer lohi = -1;
#else
  doublereal endl=-1*pow(10.0,nfig), endr=0.0;
  integer ev;
#endif

  doublereal Kappa = pow(10.0,-nfig);  /* Gives desired precision  */

  integer J = 0; /*outputs j is number of lanzcos steps taken, 
		     neig is number of eigenvalues converged*/
  
  /*Lanczos block is always 1*/

  integer Msglvl = 0;  /* message level for Plandr2 (0 to 20)*/
  doublereal CondM = 1.0;
 
  integer mpicomm=MPI_COMM_WORLD;

  /*MPI variables*/  
  int my_id, comm_size;

  MPI_Comm_rank(MPI_COMM_WORLD, &my_id);
  MPI_Comm_size(MPI_COMM_WORLD, &comm_size);

  t1=clock();

#if MYPRINT
  if(my_id==MASTER)
    printf("processor %d Entering spectrum\n\n",my_id);
#endif

#if PRINTIT
  if (my_id == MASTER)
    {
  printf("    Eigenvalues, %li digits, globalkt=%i\n",
	 nfig,globalkt);
  printf("      K=%i/2, %i-particles, o=%i, multi=%i, multi2=%i (%g %g);\n",
	 s->kt,s->npmax,s->o,s->multi,s2->multi,angle[0],angle[1]);
    }
#endif


/*
  make basis, etc 

  It is more efficient if mapcouplings is called with
  the globally largest value of kt so that the improved
  matrix elements are not recalculated.
*/

  beta=mapcouplings(globalkt,couplings);

#if MYPRINT
  printf("\nAfter mapcouplings stuff\n");
#endif

  makebasis(s);

#if MYPRINT
  printf("\nAfter makebasis\n");
#endif

  if(s->o != s2->o || s->multi != s2->multi){
    if(s->cyclic && s->kt>0){
      rebasis(s,angle);
      makebasis(s2);
      rebasis(s2,angle);
      makephaselist(angle,s->kt);   
    } else {
      makebasis(s2);
      maketenslist(angle[HINDEX+1],s->kt,angle[HINDEX],beta);
    }      
  }
  
#if MYPRINT
  printf("\nAfter more makebasis stuff\n");
#endif

  /* set aside initial guess for memory for nonzero matrix elements */

  if(els.tchunk<MAX_ELS||els.tchunk<max_els)
    {
      els.tchunk=(MAX_ELS<max_els)?max_els:MAX_ELS;
      els.i=realloc(els.i,els.tchunk*sizeof(size_t));
      if(NULL==els.i)
	{
	  fprintf(stderr,"els.i not allocated\n");
	  goto error;
	}
      els.e=realloc(els.e,els.tchunk*sizeof(element));
      if(NULL==els.e)
	{
	  fprintf(stderr,"els.e not allocated\n");
	  goto error;
	}
  }


  /* Choose Hamiltonian.  If the two bases are equal then use
     real valued, else use complex case (nonzero P_1).
     Nonzero string tension must always be complex valued. */
  if(s2->o==s->o && s2->multi==s->multi)
    {
#if PRINTIT
      if(my_id == MASTER)
	{
	  printf("      using hamels");
	  fflush(stdout);
	}
#endif
      hamels(&els,s,couplings);
    }
  else
    {

#if PRINTIT
      if(my_id == MASTER)
	{
	  printf("      using hamelsc");
	  fflush(stdout);
	}
#endif 
      hamelsc(&els,s,s2,couplings);
    }

  n=els.n;

  /* Make sure that maxj is never larger than els.n (error 32 in plandr2) */
  if (n>maxmaxj) 
    maxj = maxmaxj;
  else 
    maxj = n;

  if(els.j==NULL)
    {
      fprintf(stderr,"spectrum() detecting an error in hamels or hamelsc\n");
      goto error;
    }

  
  /* remember largest number of nonzero matrix elements. */

  if(max_els<els.j[els.n])max_els=els.j[els.n];    
  
#if 0 /* print out matrix elements of the Hamiltonian */
  printels(els,25);
#elif 0  /* print out all matrix elements of the Hamiltonian. */
  {
    FILE *f1;
    char fn[]="laso1.out";
    f1=fopen(fn,"w");
    for(i=0; i<els.n+1; i++)
      fprintf(f1,"%i\n",els.j[i]);
    fclose(f1);
    fn[4]='2';
    f1=fopen(fn,"w");
    for(i=0; i<els.j[els.n]; i++)
      fprintf(f1,"%i\n",els.i[i]);
    fclose(f1);
    fn[4]='3';
    f1=fopen(fn,"w");
    for(i=0; i<els.j[els.n]; i++)
      fprintf(f1,"%.16e\n",els.e[i]);
    fclose(f1);
   }
#endif
  
#if 0  /* print out the length of the longest row */
  for(i=0, j=0; i<n; i++)
    if(els.j[i+1]-els.j[i]>j)
      j=els.j[i+1]-els.j[i];
  printf("longest row is %i\n",j);
#endif
 
  
  if(n>lval){
    if(NULL==(spectrumval=realloc(spectrumval,(lval=n)*sizeof(doublereal))))
      {
	fprintf(stderr,"spectrumval not allocated\n");
	goto error;
      }
    if(NULL==(given_vector = realloc(given_vector, lval * sizeof(doublereal))))
      {
	fprintf(stderr, "given_vector not allocated\n");
	goto error;
      }
    if(NULL==(tempsum = realloc(tempsum, lval * sizeof(doublereal))))
      {
	fprintf(stderr, "tempsum not allocated\n");
	goto error;
      }  
  }

      /*
	Set up memory allocation for plandr2 and OP
	The pointers work,val and spectrumvec are defined globally.

      */

#if 0  /* spectrumvec is not defined */
  if(lvec<n*nval)
    {
      spectrumvec=realloc(spectrumvec,(lvec=n*nval)*sizeof(doublereal));
      if(NULL==spectrumvec)
	{
	  fprintf(stderr,"spectrumvec not allocated\n");
	  goto error;
	}
    }
#endif
  
  if(n==1)
    {
      if (els.j[els.n] == 1)
	temp = els.e[0];
      else
	temp = 0.0;
      MPI_Reduce(&temp,spectrumval,1,MPI_DOUBLE,MPI_SUM,MASTER,MPI_COMM_WORLD);
#if 0
      spectrumvec[1]=1;
#endif
    }
  else 
    {
  
      /******************* choose PLANDR(2) **************************/
  
#if !DOP2 && FINDVEC
      ev = my_id + 10;
#elif !DOP2
      ev =0;
#endif

      /*MPI array allocation*/
      
      if(NULL==(recvsendarray = malloc(comm_size * sizeof(*recvsendarray))))
	{
	  fprintf(stderr, "recvsendarray not allocated\n");
	  goto error;
	}
      if(NULL == (displs = malloc(comm_size * sizeof(*displs))))
	{
	  fprintf(stderr, "displs not allocated\n");
	  goto error;
	}
      
      /*find the hunk size for each lanczos vector on each pe*/
      FindChunk(&localn, my_id, comm_size, MPI_COMM_WORLD);
      
      /* finds maximum value of work size for the work array*/
#if FINDVEC && !DOP2
      if(localn > (maxj + 1))
	bet = 5 * localn + 4 * maxj + localn + 1 + (maxj * maxj);
      else
	bet = 5 * localn + 4 * maxj + maxj + 2 + (maxj * maxj);
#else
     if(localn > (maxj + 1))
	bet = 5 * localn + 4 * maxj + localn + 1;
      else
	bet = 5 * localn + 4 * maxj + maxj + 2;
#endif  

#if 0
      work = malloc(bet * sizeof(doublereal));
      /*end of work stuff*/
#endif

      if(bet>maxwork)
	{
	  maxwork=(bet>MAX_WORK)?bet:MAX_WORK;
	  work=realloc(work,maxwork*sizeof(doublereal));
	}

/* if maxj is changed then the following statement could get us in all kinds of trouble */
      if (NULL==(Bnd = malloc(maxj * sizeof(doublereal))))
	{
	  fprintf(stderr, "Bnd not allocated\n");
	  goto error;
	}

      for (i=0; i < bet; i++)work[i] = 0.0;

      t1=clock();
#if MYPRINT
      printf("\nCalling Plandr2");
#endif
      
#if DOP2
      PLANDR2(&localn, &maxj, &nval, &lohi, &CondM, &Kappa, &J, &nperm, 
	      spectrumval, Bnd, work, &bet, &ierr, &Msglvl, &mpicomm);
#else
      PLANDR(&localn, &maxj, &nval, &CondM, &endl, &endr, &ev, &Kappa, &J, &nperm, 
	     spectrumval, Bnd, work, &ber, &ierr, &Msglvl, &mpicomm);
#endif

#if MYPRINT
      printf("\nReturning from Plandr2\n");
#endif

#if PRINTIT
      if(my_id == MASTER)
	printf(" in %f seconds, %li multiplies, maxj=%li\n",
	       (float)(clock()-t1)/(float) CLOCKS_PER_SEC,J,maxj);
#endif
      if(ierr!=0)
	{
	  fprintf(stderr, "PLANDR2 returning error = %li\n", ierr);
	  goto error;
	}
    }
  

#if PRINTIT /* Print eigenvalues */
  if (my_id == MASTER)
    {
      printf("      ");
      for(j=0, i=3; j<nval; j++)
	{
	  printf("%.10g ",spectrumval[j]);
	  if((i++)%6==0 && j+1<nval)printf("\n      ");
	}
      printf("\n");
    }
#endif

#if 0 /* debugging print statements */
  printf("ierr =% li       number of matrix accesses = %li\n",ierr,ind[0]);
  printf("\n      Eigenvalue       residual norm  "
	 "value error   vector error\n");
  for(j=0; j<nperm; j++)
    {
      printf(" %.12g",spectrumval[j]);
      for(i=1; i<4; i++)printf("  %e",spectrumval[i*nval+j]);
      printf("\n");
    } 
#endif

#if 0 /* debugging statement to print out eigenvectors */
  printf("\n First four components of the eigenvectors:\n");
  for(j=0; j<nperm; j++)
    {
      for(i=0; i<4; i++)printf("  %e",spectrumvec[i+j*n]);
      printf("\n");
    } 
#endif

  free(Bnd);
  free(recvsendarray);
  free(displs);
  free(work);
 

  /* A pointer to the list of eigenvectors in vecin */
  *vecin=spectrumvec; 
  *vals=spectrumval;
  
#if MYPRINT
  if(my_id == MASTER)
    printf("\n%d Leaving spectrum\n\n", my_id);
#endif
  
  return;

  /* Yes, there are rare occasions where a goto is OK....  */
 error: 

  fprintf(stderr,"Error in spectrum().\n");
  free(work);

  *vecin=NULL;
  *vals=NULL;
  return;

}

/*
  free() all memory used in spectrum 
   On the CRAY, the various quantities should be
   deallocated in reverse order of allocation.
   
   This is only of use for reducing memory fragmentation
   if freemahk() is also called.   This is because resthk() 
   calls mahk().
*/

void freespectrum(void)
{
  free(work); free(spectrumval); free(spectrumvec);
   work=NULL;  spectrumval=NULL;  spectrumvec=NULL;
    maxwork=0;  lval=0;            lvec=0; 
/* Free memory used by the structure els. */
  free(els.j);  free(els.e); free(els.i); 
  els.j=NULL;   els.e=NULL;  els.i=NULL;  
  els.jchunk=0; els.tchunk=0;            
}

/* subroutine to print out memory use of the structure els */

void printhels(FILE *fp)
{
  fprintf(fp,"MM els.e, els.i have length %i; els.j length is %i\n",
	  els.tchunk,els.jchunk);
  fprintf(fp,"MM (element)=%i and (size_t)=%i bytes, total is %i bytes\n",
	  sizeof(element),sizeof(size_t),(sizeof(element)+sizeof(size_t))*
	  els.tchunk+els.jchunk*sizeof(size_t));
}

element fdot(const element *x,const element *y,int n){
  int i;
  element t=0.0;
  for(i=0; i<n; i++)t += x[i]*y[i];
  return t;
}

int dot(const int *x,const int *y,int n){
  int i,t=0;
  for(i=0; i<n; i++)t += x[i]*y[i];
  return t;
}

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

/*  This is the subroutine that stores vectors for dnlaso.
    See the documentation file "laso.doc"      */

#define DEBUGIO 0 /* flag to print out debug statements */



/*******************G & L PLANDR2 functions***************/

extern fsub OP(integer *n, doublereal *P, doublereal *Q, doublereal *R, 
	       integer *our_comm_world) 
/*P and Q are Given (same), R is Answer*/
{ 

#if PDEBUGGER 
  static int plandrcounter = 0;
#endif

  /*****Give a local copy of the given vector to all pe's*****/
  int my_id, comm_size,int_n=*n; 
  integer counter, i;
  integer counti, countk, initcount, PointI;

#if PDEBUGGER
  plandrcounter++;  
  if (plandrcounter >= COU)
    {
      printf("\n     Plandr has called OP");
      printf("\n     Doing MPI rank and size stuff");
    }
#endif

  MPI_Comm_rank(*our_comm_world, &my_id);
  MPI_Comm_size(*our_comm_world, &comm_size);

#if PDEBUGGER
  if(plandrcounter >= COU)
    printf("\nmyid=%d call number=%d\n", my_id, plandrcounter);
  printf("els.n=%i %p %p\n",els.n,given_vector[0],given_vector,
	 tempsum);
#endif

  for(counter=0; counter< els.n; counter++)
    {
      given_vector[counter] = 0.0;
      tempsum[counter] = 0.0;
    } 
  

  /*   CONCATENATES EACH PROCESSOR'S SMALL GIVEN INTO LARGE GIVEN FOR EVERYBODY'S USE      */

  MPI_Allgather(&int_n, 1, MPI_INT, recvsendarray, 1, MPI_INT, *our_comm_world);
  
#if PDEBUGGER
  if(plandrcounter >= COU)
    printf("\n%d %d Leaving first MPI routine Allgather", my_id, plandrcounter);
#endif
  
  displs[0] = 0;
  for(i=1;i< comm_size; i++)
    displs[i] = displs[i-1] + recvsendarray[i-1];
  
#if PDEBUGGER && 0
  if(plandrcounter >= COU)
    printf("\n%d %d Calling MPI routine Allgatherv", my_id, plandrcounter);
#endif

#if PDEBUGGER
  MPI_Errhandler_set(MPI_COMM_WORLD, MPI_ERRORS_RETURN);
#endif

  counter = MPI_Allgatherv(Q, *n, MPI_DOUBLE, given_vector, recvsendarray, displs, MPI_DOUBLE, 
		*our_comm_world);

  /*                 END CONCATENATION PROCESS                */

#if PDEBUGGER
  if(plandrcounter >= COU)
    printf("\n%d %d Returning from MPI routine Allgatherv", my_id, plandrcounter);
#endif  
  
  for (initcount = 0; initcount < *n; initcount++)
    R[initcount] = 0;
  
#if PDEBUGGER
  if(plandrcounter >= COU)
    printf("\n%d %d initcount initialized", my_id, plandrcounter);
#endif  

  for(counti = 0; counti < els.n; ++counti)
    {
      countk = els.j[counti];
      PointI = els.j[counti + 1];

      /* see if this row is not empty and if it has a diagonal */
      /* Both conditions are needed for the last element of the matrix */
      if (countk<PointI && els.i[countk] == counti)
	{
	  tempsum[counti] += els.e[countk] * given_vector[els.i[countk]];
	  countk++;
	}
      
#ifdef _CRI
#pragma _CRI ivdep
#endif
      for(; countk < PointI; ++countk)
	{
	  tempsum[els.i[countk]] += els.e[countk] * given_vector[counti]; 
	  tempsum[counti]+= els.e[countk] * given_vector[els.i[countk]];
	}     
    }

#if PDEBUGGER
  if(plandrcounter >= COU)
    printf("\n%d %d Matrix vector multiply finished successfully\n\n", my_id, plandrcounter);
#endif 

  MPI_Reduce_scatter(tempsum, R, recvsendarray, MPI_DOUBLE, MPI_SUM, *our_comm_world);

  MPI_Barrier(*our_comm_world);

#if PDEBUGGER 
  if(plandrcounter >= COU)
    printf("\n%d %d Returning to plandr2\n", my_id, plandrcounter);
#endif

  return SUBRETURN;
}

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


extern fsub OPM(integer *n, doublereal *x, doublereal *y, integer *our_comm_world)
{
  integer count;
  for(count = 0; count < *n; count++)
    y[count] = x[count];
  return SUBRETURN;
  
}

void FindChunk(integer *Size, int PeID, int Comm_Size, MPI_Comm Our_Comm)
{
  int tempsize;  
  *Size = (els.n * (PeID + 1))/Comm_Size - (els.n * PeID)/Comm_Size;
  
  MPI_Reduce(Size, &tempsize, 1, MPI_INT, MPI_SUM, 0, Our_Comm);
  
  if(PeID == 0 && els.n != tempsize)
    fprintf(stderr,"The tempsizes don't add up! tempsize=%i n=%i\n",
	    tempsize,els.n);
}
#endif