/****   Copyright (c) 2001, 2002  Alexandros V. Gerbessiotis
 ****
 ****   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"
#define GLOBAL(BASE,J) (((BASE/b)*nprocs*b)+pid*b+J)

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

/*Implements a binomial tree algorithm.
 */

void  bt(double K, int T, double S, double r,
         int    n, double sigma, int blocksize, double *rstock, 
         double *roption)
{
 int i,j,k,b,ind;
 int pid,nprocs;
 double dt, p, discount;
 double dxu,dxd,pu,pd,v;
 double doveru;
 int prenp,mednp,medsize,rest,initnp,nelements,base,basis,offset,
     size,fblocks,pfblocks,ks,preproc;
 double *Sinit, *C, *St,*Cin, *Cout;
 double dp,odp, oneoveru;

/***** STAGE0.0:: PID+NPROCS init. */
 pid=AIPID(); nprocs=AINPROCS();

/* STAGE1.0:: INIT() starts here */
 /* Computte Delta t, p, and discount */
 dt  = (double) T / (double) n;
 v = r - 0.5*sigma*sigma;
 dxu = sqrt(sigma*sigma*dt+v*dt*v*dt);
 dxd = -dxu;
 pu  = 0.5 + 0.5 * (v*dt /dxu);
 pd  = 1-pu;
 discount= EXP(-r*dt);

 doveru=EXP(dxd-dxu);
 b=blocksize;

/***** STAGE1.1::   */
 /* 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;
 }

/***** STAGE1.2:: */
 Sinit = (double *) malloc(initnp*sizeof(double));
 for(ind=0;ind<initnp;ind++) Sinit[ind]=0.0;
#ifdef  SANITY_CHECK
 check_if_null((void *),Sinit,"bt","Sinit");
#endif
 Sinit[0] = S*EXP(n*dxu)*EXP((dxd-dxu)*prenp*pid);


 for(i=1;i<initnp;i++) {
  Sinit[i] = Sinit[i-1] * doveru;
 }

/***** STAGE 2.1 :: InitC() */
 mednp  =((n+1)/(nprocs*b))*b;/*All procs get mednp/b blcks*/
 medsize=mednp; 
 rest   =(n+1)-(mednp*nprocs);/* 0..(rest/b) get one more*/ 
 if (rest >0) {
  if (pid <= (rest/b)) 
    medsize= mednp+b; 

  if (pid < (rest/b)) /* and for these processors */
    mednp = mednp + b ; /* mednp is updated */

  if (pid == (rest/b)) /* except the last proc */
    mednp = mednp+(rest-(rest/b)*b); 
 }
 /* mednp  :: number of nodes assigned to processor */
 /* medsize:: number of nodes rounded up to a multiple of b*/
 
/***** STAGE 2.2 :: Allocate space for copy */
 /* NOTE: instead of allocating mednp ****
  *       we ALWAYS allocate multiple of b
  *       i.e. medsize amount of memory
  */
 C = (double *) malloc(medsize*sizeof(double));
 for(ind=0;ind<medsize;ind++) C[ind]=0.0;
#ifdef  SANITY_CHECK
 check_if_null((void *),C,"bt","C");
#endif

 St = (double *) malloc(medsize*sizeof(double));
 for(ind=0;ind<medsize;ind++) St[ind]=0.0;
#ifdef  SANITY_CHECK
 check_if_null((void *),St,"bt","St");
#endif

/***** STAGE 2.3 :: Copy into allocated memory */
 
 base = pid * prenp;          /* global address*/
 offset = base - (base/b)*b;  /* offiset in a block */
 nelements=b-offset;

 AIREGISTER(St,medsize*sizeof(double));
 AISYNC();
 for(i=0;i<initnp;i+=size) {
   size= MIN((nelements),(initnp-i));
   AIHPPUT((base/b)%nprocs,&Sinit[i], &St[0],
      ((base/(b*nprocs))*b+offset)*sizeof(double), size*sizeof(double));

   base += size;
   offset=(offset+size)%b; /*or offset=(offset+size)%b*/
   nelements=b;
 }
 AISYNC();

/***** STAGE 2.4 :: Initialize option prices */
 for(i=0;i<mednp;i++) {
#ifdef CALL
   C[i] = MAX((0.0), (St[i]-K));
#else
   C[i] = MAX((0.0), (K-St[i]));
#endif
 } 
 for(i=mednp;i<medsize;i++) {
   C[i] = 0.0;
 } 


 Cin  = (double *) malloc(medsize*sizeof(double));
 Cout = (double *) malloc(medsize*sizeof(double));
#ifdef  SANITY_CHECK
 check_if_null((void *),Cout,"bt","Cout");
 check_if_null((void *),Cin ,"bt","Cin ");
#endif
 for(i=0;i<medsize;i++) {
   Cin[i]=Cout[i]=0.0;
 } 


/*******************************************/
/***** STAGE 3   ::  Main loop *************/
/*******************************************/
  dp=discount*pu;
  odp=discount*pd;
  oneoveru= 1.0/EXP(dxu);

  for(k=n; k>0; k=k-b+1){
  /***** STAGE 3.1 ::  Initsize*************/
    mednp=0;
    fblocks=0;  /*fully filled blocks*/
    pfblocks=0; /*partially or fully filled blocks*/
    fblocks= (k+1)/(nprocs*b);
    mednp=fblocks*b;
    pfblocks=fblocks;
    rest=(k+1)-mednp*nprocs;
    if (rest >0) {
      if (pid < (rest/b)) {
        mednp=mednp+b;
        fblocks++;
        pfblocks++;
      }
      if (pid == (rest/b)) {
        mednp=mednp+(rest-(rest/b)*b);
        if ((rest-(rest/b)*b)>0) pfblocks++;
      }
    }

  /***** STAGE 3.2 ::  CompBlockA***********/
  for(base=0;base<mednp;base+=b){
    if (GLOBAL(base,0) <= k) 
      Cin  [base+b-1] = C[base];
  }

  for(ks=1; ks<b; ks++) {
    for(base=0;base < mednp; base+=b) {
      for(j=0;j<b-ks;j++) {
        if (GLOBAL(base,j) <= (k-ks)) {
          C[base+j] = (dp*C[base+j] + odp*C[base+j+1]);
#ifdef AMERICAN
#ifdef CALL
          St[base+j]= St[base+j] *oneoveru;
          C[base+j] = MAX( C[base+j], (St[base+j]-K) );
#else
          St[base+j]= St[base+j] *oneoveru;
          C[base+j] = MAX( C[base+j], (K-St[base+j]) );
#endif
#endif
        }
        if (j==0)
          Cin[base+b-1-ks] = C[base];
      }
    }
  }


  /***** STAGE 3.3 ::  CommBlockH***********/
  preproc = ((pid == 0)? (nprocs-1): (pid-1));
  basis   = ((pid == 0)? b         :  0     );

  AIREGISTER(&Cout[0],medsize*sizeof(double));
  AISYNC();
  for(base=basis; base<medsize; base+=b ) {
    AIPUT(preproc,&Cin[base],&Cout[0],(base-basis)*sizeof(double),
            b*sizeof(double));
  }
  AISYNC(); 
  /***** STAGE 3.4 ::  CompBlockB***********/
  for(ks=1; ks<b ; ks++){
    for(base=0; base<mednp; base+=b) {
      for(j=b-ks;j<b-1;j++){
        if (GLOBAL(base,j) <= (k-ks)) {
          C[base+j]=dp*C[base+j]+ odp*C[base+j+1];
#ifdef AMERICAN
#ifdef CALL
          St[base+j] = St[base+j]*oneoveru;
          C[base+j] = MAX( C[base+j], (St[base+j]-K) );
#else
          St[base+j] = St[base+j]*oneoveru;
          C[base+j] = MAX( C[base+j], (K-St[base+j]) );
#endif
#endif
        }
      }
      if (GLOBAL(base,(b-1)) <= (k-ks)) {
          C[base+b-1]=dp*C[base+b-1]+ odp*Cout[base+b-ks];
#ifdef AMERICAN
#ifdef CALL
          St[base+b-1] = St[base+b-1]*oneoveru;
          C[base+b-1] = MAX( C[base+b-1], (St[base+b-1]-K) );
#else
          St[base+b-1] = St[base+b-1]*oneoveru;
          C[base+b-1] = MAX( C[base+b-1], (K-St[base+b-1]) );
#endif
#endif
       
      }
    }
  }

}

 if (pid ==0) {
  *rstock= St[0]; *roption=C[0];
 }

 free((void *)Cout);
 free((void *)Cin);
 free((void *)St);
 free((void *)C);
 free((void *)Sinit);


}
