/**********************************************************
 **** 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>
#define  GLOBAL(BASE,J) (((BASE/b)*nprocs*b)+pid*b+J)
#undef   OPTION_DEBUG
#define   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

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);

}

static 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,bigmax=0;
  char *zerobuffer;
  char * ddst;

  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*sizeof(int),0);
  SWARM_Barrier();
  for(i=0;i<nprocs ;i++) {
    on_thread(i) {
       if (top >=0)
         mybsize[nprocs*pid     +top   ]=b;
    }
  }
  myoffset[top]=offset;
  SWARM_Barrier();
  max=0;
  for(i=nprocs*pid     ,j=0;j<nprocs;i++,j++) {
    if (mybsize[i] > max) {
      max=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);
  SWARM_Barrier();
  for(i=nprocs -1;i>= 0;i--) {
    if (i !=0) {
     on_thread(i) {
       memcpy(zerobuffer,src,b*sizeof(char));
     }
    }
     SWARM_Barrier();
    if (i >0) {
     on_thread(i-1) {
       memcpy(&ddst[myoffset[pid]],(char *)zerobuffer,b*sizeof(char));
     }
    }
     SWARM_Barrier();
  }


  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;
 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
 } 
 SWARM_Barrier();
 PRINTMAT2("C",10+MYTHREAD,C,initnp,0);

 Cin  = (double *) malloc(b*sizeof(double));
 Cout = (double *) malloc(b*sizeof(double));
 Ci2  = (double *) malloc(b*sizeof(double));
 Cou2 = (double *) malloc(b*sizeof(double));
#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;
 } 


/*******************************************/
/***** 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();
    printf("MYT-READ=%d k=%d\n",MYTHREAD,k);
    PRINTMAT2("C-1>",MYTHREAD,Cout,b,0);
    SWARM_Barrier();
    bsp_sput(preproc,&Ci2[0],&Cou2[0],(b-MIN(k,b))*sizeof(double), MIN(k,b)*sizeof(double),TH);
    SWARM_Barrier();
    printf("MYT*READ=%d k=%d\n",MYTHREAD,k);
    PRINTMAT2("C-2>",MYTHREAD,Cou2,b,0);
    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;
 int i,j;
 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);

 exfdm(K,T,S,r,nj,sigma,dx,dividend,blocksize,&rst,&rc,nj,TH);
 fflush(stdout);
 on_one_thread{
  printf("%d %d %d  %f %f\n",nprocs,blocksize,nj,rc,rst);
 }
}


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

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

}
