/**********************************************************
 **** 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	"bt.h"
#include "stdio.h"
#define GLOBAL(BASE,J) (((BASE/b)*nprocs*b)+pid*b+J)

#undef   OPTION_DEBUG

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

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

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)
{
 int i,j,k,b,ind,n,in2;
 int pid,nprocs, lastproc;
 double dt, p, discount, tmp;
 double pu,pm,pd,v,invexpdx;
 int prenp,medsize,initnp,base,
     size,ks,preproc;
 double  *C, *Co,  *St,*Cin, *Cout, *Ci2,*Cou2;
 double dp,odp, oneoveru;

/***** STAGE0:: Initializations   */
 pid=AIPID(); nprocs=AINPROCS();
 lastproc= nprocs-1;

 n = 2*Nj;
/* STAGE1.0:: INIT() starts here */
 /* Computte 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);
 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;
 }

 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
 } 

 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;

#ifdef ALEX
#endif
  AIREG(&Cout[0],b*sizeof(double));
  AIREG(&Cou2[0],b*sizeof(double));
#ifdef BSPLIB
  AISYNC();
#endif
  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  ***********/

   preproc = pid-1;
   AISNC(&Cout[0]);
   if (pid !=0) {
    APUT(preproc,&Cin[0],&Cout[0],(b-MIN(k,b))*sizeof(double), MIN(k,b)*sizeof(double));
   }
   AISNC(&Cout[0]); 

   AISNC(&Cou2[0]);
   if (pid !=0) {
    APUT(preproc,&Ci2[0],&Cou2[0],(b-MIN(k,b))*sizeof(double), MIN(k,b)*sizeof(double));
   }
   AISNC(&Cou2[0]); 


  /***** 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
      }
    }
  }
 }
 if (pid ==0) {
  *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);
}
