/**********************************************************
 **** Option price valuations with the explicit finite
 ****  difference method. This is the file that contains
 ****  the paralle code compilable under LAM-MPI and BSPlib
 ****/
/****Copyright (c) 2002-2007  Alexandros V. Gerbessiotis. 
 **** All rights reserved.                               
 ****Permission to use, copy, modify, and distribute this software,
 ****and to incorporate it, in whole or in part, into 
 ****other software,is hereby granted without fee, provided that
 ****  (1) the above copyright notice and this permission 
 ****      notice appear in all copies of the source code, and 
 ****      the above copyright notice appear in clearly visible 
 ****      form on all supporting documentation and distribution 
 ****      media, and
 ****     (2) any redistribution of the software, in original
 ****         or modified form, be without fee and subject to
 ****         these same conditions.
 ****   No guarantee, whatsoever, is offered that the code works as
 ****   advertised or is absent of any, even damaging,
 ****   side-effects. Use this code at your own personal risk.
 ****/
#include <stdio.h>
#include "math.h"
#include "bt.h"
#include <swarm.h>
#include <sys/time.h>
#define  GLOBAL(BASE,J) (((BASE/b)*nprocs*b)+pid*b+J)
#define   OPTION_DEBUG 
#undef   OPTION_DEBUG


#define PRINTMAT2(X,Y,Z,W,B)  {printf(#X"[%d]::",Y); \
                           for(k=0;k<W;k++) { \
                             printf("%5.2f ",Z[B+k]);\
                           }\
                           printf("\n");}


/*Implements a trinomial tree algorithm using explicit finite differences.
 */

#undef AMERICAN
#define AMERICAN
#define CALL
#undef CALL



static struct timeval _tstart, _tend;
static struct timezone _tstruct;

void   avg_tstart  (void);
void   avg_tend    (void);
double avg_treturn (void);

void avg_tstart(void) {
  gettimeofday(&_tstart, &_tstruct);
}

void avg_tend(void) {
  gettimeofday(&_tend,&_tstruct);
}

double avg_treturn(void) {
  double t1,t2;
  t1=(double)_tstart.tv_sec + (double)_tstart.tv_usec/(1000*1000);
  t2=(double)_tend.tv_sec + (double)_tend.tv_usec/(1000*1000);
  return(t2-t1);
}

void check_if_null(void *x,char *s,char *t)
{
  if (NULL == x) {
     printf("%s:: Null pointer error (%s)\n",s,t);
     exit(1);
  }
}

void check_if_valid(int x,int y, char *s,char *t)
{
 if (x != y) {
  printf("%s:: Intparameter %d should be = to %d (%s)\n",
         s,x,y,t);
  exit(1);
 }
}


static void bsp_bcast(int source, char * from, char * to, int size, THREADED)
{
  int pid;
  int nprocs;
  int i;
  char *temp;
  pid= (int) MYTHREAD;
  nprocs = (int) THREADS;
  temp = (char *) SWARM_malloc(size,TH);
  SWARM_Barrier();
  on_thread(source) {
    memcpy(temp,from,size);
   }
  SWARM_Barrier();
  for(i=0;i< nprocs; i++) {
    on_thread(i) {
      memcpy(to,temp,size);
    }
  }
 SWARM_Barrier();
 SWARM_free(temp,TH);

}


void bsp_sput(int top   ,const void * src, const void *dst, int offset,int b,THREADED)
{
  int  pid;
  int  nprocs;
  int *mybsize,*myoffset;
  int i,j;
  int max_buffer=0,bigmax=0;
  char *zerobuffer;
  char * ddst;
  int dpid,doff,db;
  
  ddst=(char *) dst;
  pid=(int) MYTHREAD;
  nprocs=(int) THREADS;
  if (offset < 0) {
    offset=0;
  }
  mybsize  = (int *) SWARM_malloc(nprocs *nprocs *sizeof(int),TH);
  myoffset = (int *) SWARM_malloc(nprocs *nprocs *sizeof(int),TH);
  memset(mybsize, nprocs*nprocs*sizeof(int),0);
  memset(myoffset,nprocs*nprocs*sizeof(int),0);
  SWARM_Barrier();
  for(i=0;i<nprocs;i++) {
    on_thread(i) {
       if (top >=0)  {
         mybsize[nprocs*pid     +top   ]=b;
         myoffset[nprocs*pid+    top   ]=offset;
       }
    }
  }
  SWARM_Barrier();

  for(i=nprocs*pid,j=0;j<nprocs;i++,j++) {
    if (mybsize[i] > max_buffer) {
      max_buffer=mybsize[i];
    }
  }
  on_one_thread{
    bigmax=0;
    for(j=0;j<nprocs*nprocs;j++) {
      if (mybsize[i] > bigmax) {
        bigmax=mybsize[i];
      }
    }
  }
  zerobuffer=(char *)SWARM_malloc(bigmax*sizeof(char),TH);

  for(i=nprocs-1;i>= 1;i--) {
      on_thread(i) {
         memcpy(zerobuffer,src,b*sizeof(char));
      }
      SWARM_Barrier();
      dpid=(i-1+nprocs)%nprocs;
      doff=myoffset[i*nprocs+dpid];
      db  =mybsize [i*nprocs+dpid];
      on_thread(dpid) {
        if (db>0) {
          memcpy(&ddst[doff],(char *)zerobuffer,db*sizeof(char));
        }
      }
      SWARM_Barrier();
  }


}


void bsp_init_buf(int **myb, int **myoff, char ***shared_buf,
                  char ***dst_adr,int nprocs, THREADED)
{
  int *mmyb, *mmyoff,i,pid;
  char **sshared, **ddst_adr;

  pid = (int) MYTHREAD;
  mmyb   = (int *) SWARM_malloc(nprocs *nprocs *sizeof(int),TH);
  mmyoff = (int *) SWARM_malloc(nprocs *nprocs *sizeof(int),TH);
  sshared = (char **) SWARM_malloc(nprocs*sizeof(char *),TH);
  ddst_adr = (char **) SWARM_malloc(nprocs*sizeof(char *),TH);
on_thread(0){
  memset(mmyb, nprocs*nprocs*sizeof(int),0);
  memset(mmyoff,nprocs*nprocs*sizeof(int),0);
  for(i=0;i<nprocs;i++) {
    sshared[i]=NULL;
/*
 *     printf("IN INIT sshared: %d : %d\n",i  ,(int) &sshared[i]);
 *     */
  }
  for(i=0;i<nprocs;i++) {
    ddst_adr[i]=NULL;
/*
 *     printf("IN INIT dst_adr: %d : %d\n",i  ,(int) &ddst_adr[i]);
 *     */
  }
}
  *shared_buf = sshared;
  *dst_adr = ddst_adr;
  *myb=mmyb;
  *myoff=mmyoff;
  return;
}

void mymemset(char *x, int n, int value)
{
  int i;
  for(i=0;i<n;i+=4) {
   x[i]=(char)value;
  }
}


void print_int(char *x, char *y,THREADED)
{
  int *z,i;
  i= (int) MYTHREAD;
  z= (int *) x;
  printf("%d: %8s %d %d %d ||",i,y,z[0],z[1],z[2]);
}


void bsp_put(int to   ,const void * src, const void *dst, int offset,int b,THREADED)
{
  int  pid,nprocs;
  int i,j,max_buffer=0,k;
  int real_off=0,real_b=0,real_to=0;
  char **shared_buf, **dst_adr;
  char * ddst, *ssrc,*real_dest;
  int *myb,*myoff;

  ddst  =(char *) dst;
  ssrc  =(char *) src;
  pid   =(int) MYTHREAD;
  nprocs=(int) THREADS;

  if (offset < 0) {
    offset=0;
  }
  if (pid==to) {
    memcpy(&ddst[offset],ssrc,b*sizeof(char));
    return;
  }

  bsp_init_buf(&myb,&myoff,&shared_buf, &dst_adr,nprocs,TH);
  SWARM_Barrier();
  printf(" %d ---> %d\n",pid,(int)ddst);

  for(i=0;i<nprocs ;i++) {     /* pid --> to */
     if (to >=0) {
       myoff[nprocs*pid + to]=offset;
       myb[nprocs*pid + to]=b;
     }
     dst_adr[pid]=ddst;
  }
  SWARM_Barrier();

#ifdef SANITY_CHECK
#endif
on_thread(0){
 real_dest=NULL;
 for(i=0;i<nprocs;i++) {
  printf("Src %d: ",i);
  for (j=0;j<nprocs;j++) {
    printf("%2d/%2d ",myb[nprocs*i+j],myoff[nprocs*i+j]);
    if (myb[nprocs*i+j] > 0) {
        real_dest = dst_adr[j];
        real_to=j;
    }
  }
   printf(":%9d %1d\n",(int)real_dest,real_to);
 }
}
 SWARM_Barrier()  ;

  /* allocate as much memory as processor for destination proc */

  for(k=0;k<nprocs;k++) {
    on_one_thread{
     max_buffer=0;
     for(i=k,j=0;j<nprocs;i+=nprocs,j++) { /* find the max of received info to proc to*/
       if (myb[i] > max_buffer) {
         max_buffer=myb[i];
       }
     }
     printf(" To=%d max_buffer=%d\n",k,max_buffer);
     fflush(stdout);
   }
     shared_buf[k] = (char *) SWARM_malloc( max_buffer*sizeof(char),TH);
     memset(shared_buf[k],max_buffer*sizeof(char),0);
     mymemset(shared_buf[k],max_buffer*sizeof(char),pid);
   }

 on_one_thread{
  for(i=0;i<nprocs;i++){
   printf("ZORRO %d : %d %d\n",i  ,(int) &shared_buf[i  ],(int) shared_buf[i  ]);
  }
 }

 SWARM_Barrier();
   printf("SHARED BUF %d : %d %d\n",pid,(int) &shared_buf[pid],(int) shared_buf[pid]);
  SWARM_Barrier();
  /* pid from src will send to into dst  b  bytes */
  /* First pid copies from local src into shared_buf[to]  b bytes */
  /* Implicit assumption: routed info is a permutation on processors */

  printf("BEF %d : %d, ",pid,to);
  print_int(shared_buf[to],"shared_buf",TH);
  print_int(ssrc,"ssrc", TH);
  printf("\n");
  SWARM_Barrier();
  memcpy(shared_buf[to],ssrc,b*sizeof(char));
  SWARM_Barrier();
  printf("AFT %d : %d, ",pid,to);
  print_int(shared_buf[to],"shared_buf",TH);
  print_int(ssrc,"ssrc", TH);
  printf("\n");
  /* Now, processor to will collect this information and copy it into
 *    * the right buffer at the correct offset 
 *       */
#ifdef SANITY_CHECK
  for(k=0;k<nprocs;k++) {
   on_thread(k) {
     real_to = k;
     real_dest = dst_adr[real_to] ;
     real_b=0; real_off=0;
     for(i=k,j=0;j<nprocs;i+=nprocs,j++) { /* find the max of received info to proc to*/
       if (myb[i] !=0 ) {
          real_b =myb[i];
          real_off=myoff[i];
       }
     }
     memcpy(real_dest+real_off,shared_buf[real_to],real_b*sizeof(char));
    }
    SWARM_Barrier();
    }
#endif
  real_to = pid;
  real_dest= dst_adr[real_to];
  real_b=0; real_off=0;
  for(i=real_to,j=0;j<nprocs;i+=nprocs,j++) {
       if (myb[i] !=0 ) {
          real_b =myb[i];
          real_off=myoff[i];
       }
  }
  memcpy(real_dest+real_off,shared_buf[real_to],real_b*sizeof(char));
  SWARM_Barrier();

}


void  exfdm(double K, int    T, double S, double r,
         int    N, double sigma, double dx,double dividend,
         int blocksize, double *rstock, double *roption, int Nj,THREADED)
{
 int i,j,k,b,ind,n;
 int pid,nprocs, lastproc;
 double dt, discount ;
 double pu,pm,pd,v,invexpdx;
 int prenp,medsize,initnp,
     ks,preproc;
 double  *C, *Co,  *St,*Cin, *Cout, *Ci2,*Cou2,*Zout,*Zou2;
 double oneoveru;

/***** STAGE0:: Initializations   */
 pid=MYTHREAD; nprocs=THREADS;
 lastproc= nprocs-1;

 n = 2*Nj;

/* STAGE1.0:: INIT() starts here */
 /* Compute Delta t, p, and discount */
 dt  = (double) T / (double) N;
 v = r - dividend - 0.5*sigma*sigma;
 invexpdx=1.0/exp(dx);
 pu  = 0.5 * dt *((sigma/dx)*(sigma/dx)+v/dx);
 pm  = 1.0 - dt* (sigma/dx)*(sigma/dx) - r*dt;
 pd  = 0.5 * dt *((sigma/dx)*(sigma/dx)-v/dx);
 discount = exp(-r*dt);

/*
 printf("MY=%d K=%f T=%d S=%f r=%f N=%d s=%f dx=%f div=%f ",MYTHREAD,K,T,S,r,
 N,sigma,dx,dividend);
 printf("blocksize=%d Nj=%d\n dt=%f v=%f inv=%f pu=%f pm=%f pd=%f dis=%f",
  blocksize,Nj,dt,v,invexpdx,pu,pm,pd,discount);
*/

#ifdef OPTION_DEBUG
 printf("MYTHREAD=%d K=%f dt=%f n=%d v=%f expdx=%f iexpdx=%f (p)=(%f,%f,%f) discount=%f\n", MYTHREAD, K,dt,n,v,1.0/invexpdx,invexpdx,pu,pm,pd,discount); 
#endif

 b=blocksize;
 /* How many nodes per processor ? : initnp */
 /* How many nodes processor before: prenp  */
 if (pid != (nprocs-1))  {
   prenp = (n+1)/nprocs;
   initnp= (n+1)/nprocs;
 }
 else {
   prenp = (n+1)/nprocs;
   initnp= (n+1)- ((n+1)/nprocs)*(nprocs-1);
 }

 /* prenp for proc 0  should be set to 0 */
 if (pid ==0) {
   prenp=0;
 }

/***** STAGE1a:: Initialize Stock Prices  */
 St = (double *) malloc(initnp*sizeof(double));
#ifdef  SANITY_CHECK
 check_if_null((void *)St,"ftt","St");
#endif
 for(ind=0;ind<initnp;ind++) St[ind]=0.0;
 St[0] = S*EXP(Nj*dx)*EXP((-dx)*prenp*pid);

 for(i=1;i<initnp;i++) {
  St[i] = St[i-1] * invexpdx;
 }
#ifdef OPTION_DEBUG
PRINTMAT2("St",MYTHREAD,St,initnp,0);
#endif

 medsize=initnp;  /* Number of nodes per processor */
 
/***** STAGE1b:: Initialize Option Prices */
 C = (double *) malloc(medsize*sizeof(double));
 Co= (double *) malloc(medsize*sizeof(double));
#ifdef  SANITY_CHECK
 check_if_null((void *)C,"ftt","C");
 check_if_null((void *)Co,"ftt","Co");
#endif
 for(ind=0;ind<medsize;ind++) Co[ind]=C[ind]=0.0;

 for(i=0;i<medsize;i++) {
#ifdef CALL
   C[i] = MAX((0.0), (St[i]-K));
#else
   C[i] = MAX((0.0), (K-St[i]));
#endif
 } 
#ifdef OPTION_DEBUG
 SWARM_Barrier();
 PRINTMAT2("C",10+MYTHREAD,C,initnp,0);
#endif

 Cin  = (double *) malloc(b*sizeof(double));
 Cout = (double *) malloc(b*sizeof(double));
 Zout = (double *) SWARM_malloc(b*sizeof(double),TH);
 Ci2  = (double *) malloc(b*sizeof(double));
 Cou2 = (double *) malloc(b*sizeof(double));
 Zou2 = (double *) SWARM_malloc(b*sizeof(double),TH);
#ifdef  SANITY_CHECK
 check_if_null((void *)Cout,"ftt","Cout");
 check_if_null((void *)Cin ,"ftt","Cin ");
 check_if_null((void *)Cou2,"ftt","Cout");
 check_if_null((void *)Ci2 ,"ftt","Cin ");
#endif
 for(i=0;i<b;i++) {
   Cin[i]=Cout[i]=Ci2[i]=Cou2[i]=0.0;
 } 
 on_one_thread{
  memset(Zout,b*sizeof(double),0);
  memset(Zou2,b*sizeof(double),0);
 }

/*******************************************/
/***** STAGE2::  Main loop *****************/
/*******************************************/
  oneoveru= invexpdx;

  for(k=N; k>0; k=k-b+1){
   /***** STAGE2.1 :: Block 0 elements saved for Block 2 usage****/
   Cin[b-1] = C[0];
   Ci2[b-1] = C[1];

   /***** STAGE2.2 :: Block 1 Computation  ***********/
   for(ks=1; ks<MIN(b,(k+1)); ks++) {
    if (pid != lastproc) {
      for(j=0;j<medsize-2*ks;j++) {
        Co[j] = pu*C[j]+pm*C[j+1]+pd*C[j+2];
        St[j]= St[j]*oneoveru;
      } 
      for(j=0;j<medsize-2*ks;j++) {
#ifdef AMERICAN
#ifdef CALL
        C[j] = MAX( Co[j], (St[j]-K) );
#else
        C[j] = MAX( Co[j], (K-St[j]) );
#endif
#endif
       } 
       Cin[b-1-ks] = C[0];
       Ci2[b-1-ks] = C[1];
     }
     else {
       for(j=0;j<medsize-2;j++) {
         Co[j] = pu*C[j]+pm*C[j+1]+pd*C[j+2];
         St[j]= St[j]*oneoveru;
       }
#ifdef CALL
       Co[medsize-2]=Co[medsize-3];
       Co[medsize-1]=Co[medsize-2];
#else
       Co[medsize-2]=Co[medsize-3]+(St[medsize-3]-St[medsize-2]);
       Co[medsize-1]=Co[medsize-2]+(St[medsize-2]-St[medsize-1]);
#endif
       for(j=0;j<medsize;j++) {
#ifdef AMERICAN
#ifdef CALL
         C[j] = MAX( Co[j], (St[j]-K) );
#else
         C[j] = MAX( Co[j], (K-St[j]) );
#endif
#endif
       } 
       Cin[b-1-ks] = C[0];
       Ci2[b-1-ks] = C[1];
     }
   }

   /***** STAGE2.3 :: Block 2 Communication  ***********/

   SWARM_Barrier();
    preproc = pid-1;

    bsp_sput(preproc,&Cin[0],&Cout[0],(b-MIN(k,b))*sizeof(double), MIN(k,b)*sizeof(double), TH);
   SWARM_Barrier();

    bsp_sput(preproc,&Ci2[0],&Cou2[0],(b-MIN(k,b))*sizeof(double), MIN(k,b)*sizeof(double), TH);
    SWARM_Barrier();


  /***** STAGE2.3 ::  Block 3 Computation   ***********/
  for(ks=1; ks<MIN(b,(k+1)); ks++){
    if (pid != lastproc) {
      for(j=medsize-2*ks;j<medsize-2;j++) {
        Co[j] = pu*C[j]+pm*C[j+1]+pd*C[j+2];
        St[j]= St[j]*oneoveru;
      } 
      ind=medsize-2;
      Co[ind] = pu*C[ind]+pm*C[ind+1]+pd*Cout[b-ks];
      St[ind]= St[ind]*oneoveru;
      ind=medsize-1;
      Co[ind] = pu*C[ind]+pm*Cout[b-ks]+pd*Cou2[b-ks];
      St[ind]= St[ind]*oneoveru;
      for(j=medsize-2*ks;j<medsize;j++) {
#ifdef AMERICAN
#ifdef CALL
         C[j] = MAX( Co[j], (St[j]-K) );
#else
         C[j] = MAX( Co[j], (K-St[j]) );
#endif
#endif
      }
    }
  }
 }

 on_one_thread {
  *rstock= St[0]; *roption=C[0];
 }

 free((void *)Cou2);
 free((void *)Ci2);
 free((void *)Cout);
 free((void *)Cin);
 free((void *)St);
 free((void *)C);
 free((void *)Co);
}

static void main_efdm (THREADED)
{
 int nk,T,nj,blocksize;
 double K,S,r,sigma,dx,dividend;
 FILE *fp ;
 double rst,rc;
 int pid;
 int nprocs;

 pid=MYTHREAD;
 nprocs=THREADS;

 on_one_thread{
    fp=fopen("efdata1","r");
    fscanf(fp,"%d",&blocksize); /* blocksize-1 must divide n */
    fscanf(fp,"%d",&nj);   /* This is the horizon not nk: disregard nk nj=N */
    fscanf(fp,"%lf",&K);
    fscanf(fp,"%d",&T);
    fscanf(fp,"%lf",&S);
    fscanf(fp,"%lf",&r);
    fscanf(fp,"%d",&nk);
    fscanf(fp,"%lf",&sigma); 
    fscanf(fp,"%lf",&dx); 
    fscanf(fp,"%lf",&dividend); 
    dx= sigma*sqrt(3.0*T/nj);
/*
    printf("%d %d %d %d \n",blocksize, nk,nj,T);
    printf("%f %f %f %f ",K,S,r,sigma);
    printf("%f %f \n",dx,dividend);

    printf("THREAD::: %d %d \n",MYTHREAD,THREADS);
*/
 }

 bsp_bcast(0,(char *) &blocksize,(char *)&blocksize,sizeof(int),TH);
 bsp_bcast(0,(char *) &nj       ,(char *)&nj       ,sizeof(int),TH);
 bsp_bcast(0,(char *) &K        ,(char *)&K        ,sizeof(double),TH);
 bsp_bcast(0,(char *) &T        ,(char *)&T        ,sizeof(int),TH);
 bsp_bcast(0,(char *) &nk       ,(char *)&nk       ,sizeof(int),TH);
 bsp_bcast(0,(char *) &S        ,(char *)&S        ,sizeof(double),TH);
 bsp_bcast(0,(char *) &r        ,(char *)&r        ,sizeof(double),TH);
 bsp_bcast(0,(char *) &sigma    ,(char *)&sigma    ,sizeof(double),TH);
 bsp_bcast(0,(char *) &dx       ,(char *)&dx       ,sizeof(double),TH);
 bsp_bcast(0,(char *) &dividend ,(char *)&dividend ,sizeof(double),TH);

/*
 printf("pid=%d b=%d nj=%d nj=%d T=%d  K=%f S=%f r=%f sigma=%f  dx=%f div=%f\n",MYTHREAD,
  blocksize,nj,nj,T,  K,S,r,sigma,  dx, dividend);
*/


  SWARM_Barrier();
  avg_tstart();   /* t1=clock(); */
 exfdm(K,T,S,r,nj,sigma,dx,dividend,blocksize,&rst,&rc,nj,TH);
  SWARM_Barrier();
  avg_tend();  /* t2=clock(); */
  printf("%3d of %3d %18.8f secs\n",pid,nprocs,avg_treturn()/1.0);
  SWARM_Barrier();
  on_one_thread{
  printf("%3d %6d %6d  %18.8f %18.8f %18.8f secs\n",nprocs,blocksize,nj,rc,rst,avg_treturn()/1.0);
  }
}


int	main(int argc, char *argv[])
{

  SWARM_Init(&argc,&argv);
  SWARM_Run((void*) main_efdm);
  SWARM_Finalize();
  return(0);

}
