#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include "hamilton.h"  /* declaration of els  */
#include "op.h" 

/* This file contains the matrix-vector multiplication routines */

#ifdef USE_MPI
#include "usempi.h"
#endif

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

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

/*
 Routine for vector producets of the positive-definite matrix 
 of the generalized eigenvalue problem.
 Here, the matrix is the identity.
*/


#ifdef USE_MPI
extern fsub OPM(integer *n, doublereal *x, doublereal *y, integer *our_comm_world)
#else
extern fsub OPM(integer *n, doublereal *x, doublereal *y)
#endif
{
  memcpy(y,x,(*n)*sizeof(doublereal));
  return SUBRETURN;  
}

/*  complex version   */

#ifdef USE_MPI
extern fsub ZOPM(integer *n, doublecomplex *x, doublecomplex *y, integer *our_comm_world)
#else
extern fsub ZOPM(integer *n, doublecomplex *x, doublecomplex *y)
#endif
{
  memcpy(y,x,(*n)*sizeof(doublecomplex));
  return SUBRETURN;  
}

/***********************************************************************/
/*
        Routine for the matrix-vector multplies
        P and Q are Given (same), R is Answer
*/

#ifdef USE_MPI

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

extern fsub OP(integer *n, doublereal *P, doublereal *Q, doublereal *R, 
	       integer *our_comm_world) 
{ 

  /*****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;
  int *recvsendarray;
  int *displs;
  element *tempsum,*given_vector;

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

  recvsendarray=malloc(comm_size*sizeof(int));
  displs=malloc((comm_size+1)*sizeof(int));
  tempsum=malloc(els.n*sizeof(element));
  given_vector=malloc(els.n*sizeof(element));

  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_INTEGER, recvsendarray, 1, MPI_INT, 
		*our_comm_world);  
  displs[0] = 0;
  for(i=1;i< comm_size; i++)
    displs[i] = displs[i-1] + recvsendarray[i-1];
  counter = MPI_Allgatherv(Q, *n, MPI_DOUBLE, given_vector, recvsendarray, 
			   displs, MPI_DOUBLE, *our_comm_world);

  /*                 END CONCATENATION PROCESS                */

  for (initcount = 0; initcount < *n; initcount++)
    R[initcount] = 0.0;
  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){ /* the order matters */
	  tempsum[counti] += els.r[countk] * given_vector[els.i[countk]];
	  countk++;
	}
      
#ifdef _CRI
#pragma _CRI ivdep
#endif
      for(; countk < PointI; ++countk){
	  tempsum[els.i[countk]] += els.r[countk] * given_vector[counti]; 
	  tempsum[counti] += els.r[countk] * given_vector[els.i[countk]];
	}     
    }

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

  free(recvsendarray);
  free(displs);
  free(tempsum);
  free(given_vector);

  return SUBRETURN;
}

#else

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

    Hamiltonian multiply used by Lanczos routines.
    This is the largest time user for very large calculations.

    The running times for matrix diagonalization shown below are 
    for HINDEX=2, MOM=1, VECTORS=1, and

       int ht[HINDEX]={0,0},np=6,kt=20,charge=1,o=0,multi=15;

    in the program "runlaso.c" 
    Using DLCQ matrix elements, IMPROVED=1 in "include.h".

    also enter values nval=5 and momentum=1 0.
    couplings: 0.03 1 0.7 0.03 -0.005 3 0.13 1 -2 1
    639674 nonzero matrix elements, 14439 states, Lanczos with DNLASO
    203*3 matrix multiplies.

    Eigenvalues:   
     1.26979644446083e+00
     1.87282569596773e+00
     3.14511664804231e+00
     3.21836991529111e+00
     3.30083925987415e+00

*/

/*
  Here are some benchmarks for the above parameters.
  The results also vary according to the version of the program
  and the compiler used.

  time:  37.7 on jael (generating matrix takes 5.53) April 2002
         (using newly remodeled ZLANDR3)

  time:  24.4 sec on barak (generating matrix takes 10.4) April 2001
         benchmark stream gives 250 246 280 233 MB/s.
         (Tyan 1834 motherboard, 667MHz Pentium III CPU)

  time:  30 sec on bvds.geneva.edu (generating matrix takes 15) April 2000
         benchmark STREAM gives 282, 279, 330, 248 MB/s.
         (Tyan 1830 motherboard, 400 MHz Pentium II CPU)

  time:  66 alk, 81 vodka, 71 rum; GNU C Power PC
         (generating the matrix takes 35.)

  time:  21 quad IBM SP2 (generating matrix takes 31)  

  time:  62 gsm CONVEX SPP (generating matrix takes 52)
         240 Mflops peak.

  time:  69 babbage; Hitachi (generating matrix takes 60)
         300 Mflops peak.

  time:  108 vpp; Fujitsu (generating matrix takes 82), however
         above -- with vectorization -- is much faster.

  time:  44 t3e.osc.edu; CRAY (generating matrix takes 25)  April 2000
         When ivdep pragma is not used, 58.
         using pragma cache_bypass always slows things down.
         600 Mflops peak, 1.2 Gbytes/sec. peak memory access.

  time:  34 T3E jaromir.psc.edu; CRAY (generating matrix takes 17) April 2000
         Used ivdep pragma below.  The benchmark STREAM gives 
         443, 450, 523, 521 MB/s.

*/

/* 
   Matrix-vector multiplication routine used by LANSO

   The form shown below is chosen because
   it is the simplest.  There are a number of
   possible changes that do not affect speed in
   any substantial manner (for gcc on intel processor):
   
   Interchanging order of the loops.
   Using a temporary variable for p[i] and q[i]
   Using pointers for els.i[j] and els.r[j] (slightly slower)
   Using j=els.j[i]; top=els.j[i+1]; instead of j=*(els_j++);
   Introducing jtemp=els.j[i];
   splitting into two loops decreases speed slightly.

*/

extern fsub OP(integer *n, doublereal *s, doublereal *p, doublereal *q){
  integer i,j;
  size_t *els_j;

  if(els.n!= *n){
     fprintf(stderr,__FILE__":  dimensions don't agree in OP\n");
    exit(19);
  }
  memset(q,0,*n * sizeof(*p));

  for(i=0, els_j=els.j; i<*n; i++){
    j=*(els_j++);
    if(j<*els_j && els.i[j]==i)  /* the order matters */
      q[i] += els.r[j++]*p[i];
#ifdef _CRI
#pragma _CRI ivdep
#endif
    for(; j<*els_j; j++){
      q[els.i[j]] += els.r[j]*p[i];
      q[i] += els.r[j]*p[els.i[j]];
    }
  }
  return SUBRETURN;
}

#endif

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

#ifdef USE_MPI


void zsum(complex_element *in, complex_element *out, int *len, 
		       MPI_Datatype *dtype){
  int i;

  for(i=0; i< *len; i++){
    out[i].re += in[i].re;
    out[i].im += in[i].im;
  }
}

extern fsub ZOP(integer *n, doublecomplex *s, doublecomplex *p, 
		doublecomplex *q, integer *our_comm_world)
#else
extern fsub ZOP(integer *n, doublecomplex *s, doublecomplex *p, 
		doublecomplex *q)
#endif
{
  integer i,j;
  size_t *els_j;

#ifdef USE_MPI
  int my_id, comm_size;
  int *recvsendarray;
  MPI_Datatype mpi_complex_element;
  int *displs;
  doublecomplex *given, *result;
  MPI_Op zsum_function;
  int int_n=*n;

  MPI_Type_contiguous(2,MPI_DOUBLE,&mpi_complex_element);
  MPI_Type_commit(&mpi_complex_element);
  MPI_Op_create((MPI_User_function *)zsum,1==1,&zsum_function);
  MPI_Comm_rank(*our_comm_world, &my_id);
  MPI_Comm_size(*our_comm_world, &comm_size);

  recvsendarray=malloc(comm_size*sizeof(*recvsendarray));
  displs=malloc(comm_size*sizeof(*displs));
  MPI_Allgather(&int_n, 1, MPI_INT, recvsendarray, 1, MPI_INT, 
		*our_comm_world);  
  displs[0] = 0;
  for(i=1; i<comm_size; i++)
    displs[i] = displs[i-1] + recvsendarray[i-1];
  given=malloc(els.n*sizeof(doublecomplex));
  result=malloc(els.n*sizeof(doublecomplex));
  MPI_Allgatherv(p, *n, mpi_complex_element, given, recvsendarray, 
			   displs, mpi_complex_element, *our_comm_world);
#else
  doublecomplex *result=q, *given=p;
#endif


  memset(result,0,els.n * sizeof(*result));
  memset(q,0,*n * sizeof(*q));
  
  for(i=0, els_j=els.j; i<els.n; i++){
    j=*(els_j++);
    if(j<*els_j && els.i[j]==i){              /* the order matters */
      result[i].re += els.c[j].re*given[i].re; /* diagonal elements are real */
      result[i].im += els.c[j].re*given[i].im;
      j++;
    }
#ifdef _CRI
#pragma _CRI ivdep
#endif
    for(; j<*els_j; j++){
      result[els.i[j]].re += els.c[j].re*given[i].re - els.c[j].im*given[i].im;
      result[els.i[j]].im += els.c[j].re*given[i].im + els.c[j].im*given[i].re;
      /* signs for complex conjugate */
      result[i].re += 
	els.c[j].re*given[els.i[j]].re + els.c[j].im*given[els.i[j]].im;
      result[i].im += 
	els.c[j].re*given[els.i[j]].im - els.c[j].im*given[els.i[j]].re;
    }
  }

#ifdef USE_MPI
  MPI_Reduce_scatter(result, q, recvsendarray, mpi_complex_element, zsum_function, 
		     *our_comm_world);
  MPI_Op_free(&zsum_function);
  MPI_Type_free(&mpi_complex_element);
  MPI_Barrier(*our_comm_world);
  free(recvsendarray);
  free(displs);
  free(result);
  free(given);
#endif

  return SUBRETURN;
}

/*
  Multiplication routine used by DNLASO
  It allows the possiblility of doing several vectors at once.
*/

#ifndef USE_MPI
 
extern fsub elop(integer *n, integer *m, doublereal *p, doublereal *q){
  integer k;
  for(k=0; k<*m; k++)
    OP(n,NULL,p+*n*k,q+*n*k);
  return SUBRETURN;
}

#endif