
/****Copyright (c) 2002  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
#define   OPTION_DEBUG

#define PRINTMAT(X,Y,L) fprintf(stdout,#X"[%d/%d]::",L,pid); \
 for(ind=0;ind<medsize;ind++) { \
   fprintf(stdout,"%4.2f ",Y[ind]); \
 } \
 fprintf(stdout,"\n");fflush(stdout);bsp_sync();

#define PRINTLINE1(AA,C,D,E,F,G,H) fprintf(stdout,"%d:%d "#AA"*[%d] %4.2f x [%d]%4.2f x[%d] %4.2f = %4.2f (%d, %d)\n",k-ks,pid,F,tmp       ,G,D[G],H,E[H],C[F],ks,j); fflush(stdout);

#define OUT(A,B,C,D,E,F) {fflush(stdout);fprintf(stdout,"%2d %2d %8.4f %8.4f     ***%5d,%5d\n",A,B,C,D,E,F);fflush(stdout);}

#define OUT2(A,B,C,D,E,F,G,H,I,J) {fflush(stdout);fprintf(stdout,"%2d %2d %8.4f %8.4f (%8.4f, %8.4f, %8.4f) ***%3d,%2d %2d\n",A,B,C,D,G,H,I,E,F,J);fflush(stdout);}

#define SOSx(A,Z) {fflush(stdout); if (pid==0) fprintf(stdout,"%d %d\n",A,aisearch(Z));fflush(stdout);}
#define SOS(A,Z) {;}
/*Implements a binomial tree algorithm.
 */

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


void  tt(double K, int    T, double S, double r,
         int    N, double sigma, double dx,double dividend,
         int blocksize, double *rstock, double *roption)
{
 int i,j,k,b,ind,n;
 int odd,up;
 int pid,nprocs;
 double dt, p, discount, tmp;
 double dpd,dpu,dpm,pu,pm,pd,v,invexpdx;
 double doveru;
 int prenp,mednp,medsize,rest,initnp,nelements,base,basis,offset,
     size,fblocks,pfblocks,ks,preproc;
 double *Sinit, *C, *St,*Cin, *Cout, *Ci2,*Cou2;
 double dp,odp, oneoveru;

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

 n = 2*N;
/* 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 * ((sigma*sigma*dt+v*v*dt*dt)/(dx*dx) +v*dt/dx);
 pd  = 0.5 * ((sigma*sigma*dt+v*v*dt*dt)/(dx*dx) -v*dt/dx);
 pm  = 1-(sigma*sigma*dt+v*v*dt*dt)/(dx*dx);
 discount = exp(-r*dt);
 b=blocksize;

#ifdef SANITY_CHECK
 if (pid==0) {
 fprintf(stdout,"pu=%f pd=%f pm=%f invexpdx=%f expdx=%f\n",pu,pd,pm,invexpdx,exp(dx));fflush(stdout);
 }
#endif
/***** 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*dx)*EXP((-dx)*prenp*pid);


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


/***** 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*/
 
/***** STAGE2.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

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

 AIREG(St,medsize*sizeof(double));
 AISNC(St);
 for(i=0;i<initnp;i+=size) {
   size= MIN((nelements),(initnp-i));
   AHPPUT((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;
 }
 AISNC(St);

/***** STAGE2.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));
 Ci2  = (double *) malloc(medsize*sizeof(double));
 Cou2 = (double *) malloc(medsize*sizeof(double));
#ifdef  SANITY_CHECK
 check_if_null((void *),Cout,"bt","Cout");
 check_if_null((void *),Cin ,"bt","Cin ");
 check_if_null((void *),Cou2,"bt","Cout");
 check_if_null((void *),Ci2 ,"bt","Cin ");
#endif
 for(i=0;i<medsize;i++) {
   Cin[i]=Cout[i]=Ci2[i]=Cou2[i]=0.0;
 } 


/*******************************************/
/***** STAGE3   ::  Main loop *************/
/*******************************************/
  dpd=discount*pd;
  dpu=discount*pu;
  dpm=discount*pm;
  oneoveru= invexpdx;

k=N+2;ks=0;
/*
PRINTMAT("CF",C,k);
PRINTMAT("SF",St,k);
*/
#ifdef ALEX
#endif
  AIREG(&Cout[0],medsize*sizeof(double));
  AIREG(&Cou2[0],medsize*sizeof(double));
#ifdef BSPLIB
  AISYNC();
#endif
  for(k=N; k>0; k=k-b+1){
  /***** STAGE3.0 ::  Initsize*************/
    mednp=0;
    fblocks=0;  /*fully filled blocks*/
    pfblocks=0; /*partially or fully filled blocks*/
    fblocks= (k*2+1)/(nprocs*b);
    mednp=fblocks*b;
    pfblocks=fblocks;
    rest=(k*2+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++;
      }
    }

  /***** STAGE3.1 ::  BlockA Comp***********/
  for(base=0;base<mednp;base+=b){
    if (GLOBAL(base,0) <= (k*2)) {
      Cin  [base+b-1] = C[base];
    }
    if (GLOBAL((base),1) <= (k*2)) {
      Ci2  [base+b-1] = C[base+1];
    }
  }

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

  /***** STAGE3.2 ::  BlockB Comm***********/
  preproc = ((pid == 0)? (nprocs-1): (pid-1));
  basis   = ((pid == 0)? b         :  0     );

#ifdef MPILIB
  AISNC(&Cout[0]);
#endif
#ifdef MPILIB
  AISNC(&Cou2[0]);
#endif
  for(base=basis; base<medsize; base+=b ) {
    /* redundant  &Cin[base+b/2]      (base-basis+b/2)  b-(b/2) */
    APUT(preproc,&Cin[base+b/2],&Cout[0],(base-basis+b/2)*sizeof(double),
            (b-b/2)*sizeof(double));
    APUT(preproc,&Ci2[base+b/2],&Cou2[0],(base-basis+b/2)*sizeof(double),
            (b-b/2)*sizeof(double));
  }
#ifdef MPILIB
  AISNC(&Cou2[0]); 
#endif
  AISNC(&Cout[0]); 


   /***** STAGE3.3 ::  BlockC Comp***********/
  for(ks=1; ks<(b+1)/2 ; ks++){
    for(base=0; base<mednp; base+=b) {
      for(j=b-2*ks;j<b-2;j++){
        if (GLOBAL(base,j) <= (2*(k-ks))) {
          C[base+j]=dpu*C[base+j]+ dpm*C[base+j+1]+dpd*C[base+j+2];
          St[base+j] =St[base+j]*oneoveru;
#ifdef AMERICAN
#ifdef CALL
          C[base+j] = MAX( C[base+j], (St[base+j]-K) );
#else
          C[base+j] = MAX( C[base+j], (K-St[base+j]) );
#endif
#endif
        }
      }
      if (GLOBAL(base,(b-2)) <= (2*(k-ks))) {
        C[base+b-2]=dpu*C[base+b-2]+dpm*C[base+b-1]+dpd*Cout[base+b-ks];
        St[base+b-2]=St[base+b-2]*oneoveru;
#ifdef AMERICAN
#ifdef CALL
        C[base+b-2]=MAX( C[base+b-2], (St[base+b-2]-K) );
#else
        C[base+b-2]=MAX( C[base+b-2], (K-St[base+b-2]) );
#endif
#endif
      }
      if (GLOBAL(base,(b-1)) <= (2*(k-ks))) {
        C[base+b-1]=dpu*C[base+b-1]+dpm*Cout[base+b-ks]+dpd*Cou2[base+b-ks];
        St[base+b-1]=St[base+b-1]*oneoveru;
#ifdef AMERICAN
#ifdef CALL
        C[base+b-1]=MAX( C[base+b-1], (St[base+b-1]-K) );
#else
        C[base+b-1]=MAX( C[base+b-1], (K-St[base+b-1]) );
#endif
#endif
      }
    }
  }
  /* Store top element of Block C Computation (second line)*/
if (odd==1) {
}
  for(base=0;base < mednp; base+=b) {
     if (GLOBAL(base,1) <= (2*(k-(b+1)/2+1))) {
       Ci2[base+b/2] = C[base+1];
     }
  }
  preproc = ((pid == 0)? (nprocs-1): (pid-1));
  basis   = ((pid == 0)? b         :  0     );
#ifdef MPILIB
  AISNC(&Cou2[0]);
#endif
  for(base=basis; base<medsize; base+=b ) {
    /* redundant;as before */
    APUT(preproc,&Ci2[base+b/2],&Cou2[0],
         (base-basis+b/2)*sizeof(double),sizeof(double));
  }
  AISNC(&Cou2[0]); 

  /***** STAGE3.4a ::  BlockD Comp***********/
  for(ks=(b+1)/2; ks<=(b+1)/2 ; ks++){
    for(base=0; base<mednp; base+=b) {
      /* equivalent to j< b-2*(ks-(b+1)/2) */
      for(j=0;j<b-2;j++){
        if (GLOBAL(base,j) <= (2*(k-ks))) {
          C[base+j]=dpu*C[base+j]+ dpm*C[base+j+1]+dpd*C[base+j+2];
          St[base+j] =St[base+j]*oneoveru;
#ifdef AMERICAN
#ifdef CALL
          C[base+j] = MAX( C[base+j], (St[base+j]-K) );
#else
          C[base+j] = MAX( C[base+j], (K-St[base+j]) );
#endif
#endif
        }
        if (j==0) {
           if (GLOBAL(base,0) <= (2*(k-ks))) {
             Cin[base+b-1-ks] = C[base];
           }
        }
        if (j==1) {
           if (GLOBAL(base,1) <= (2*(k-ks))) {
             Ci2[base+b-1-ks] = C[base+1];
           }
        }
      }
      if  (GLOBAL(base,(b-2)) <= (2*(k-ks))) {
        C[base+b-2]=dpu*C[base+b-2]+dpm*C[base+b-1]+dpd*Cout[base+b-ks];
        St[base+b-2]=St[base+b-2]*oneoveru;
#ifdef AMERICAN
#ifdef CALL
        C[base+b-2]=MAX( C[base+b-2], (St[base+b-2]-K) );
#else
        C[base+b-2]=MAX( C[base+b-2], (K-St[base+b-2]) );
#endif
#endif
      }
      if (GLOBAL(base,(b-1)) <= (2*(k-ks))) {
        C[base+b-1]=dpu*C[base+b-1]+dpm*Cout[base+b-ks]+dpd*Cou2[base+b-ks];
        St[base+b-1]=St[base+b-1]*oneoveru;
#ifdef AMERICAN
#ifdef CALL
        C[base+b-1]=MAX( C[base+b-1], (St[base+b-1]-K) );
#else
        C[base+b-1]=MAX( C[base+b-1], (K-St[base+b-1]) );
#endif
#endif
      }
    }
  }


  for(ks=(b+1)/2+1; ks<b ; ks++){
/*
    up = 2*b+odd-2*ks;
*/
    up = b - 2*(ks- (b+1)/2) -odd;
    for(base=0; base<mednp; base+=b) {
      /* equivalent to j< b-2*(ks-(b+1)/2) */
      for(j=0;j<up;j++){
        if (GLOBAL(base,j) <= (2*(k-ks))) {
          C[base+j]=dpu*C[base+j]+ dpm*C[base+j+1]+dpd*C[base+j+2];
          St[base+j] =St[base+j]*oneoveru;
#ifdef AMERICAN
#ifdef CALL
          C[base+j] = MAX( C[base+j], (St[base+j]-K) );
#else
          C[base+j] = MAX( C[base+j], (K-St[base+j]) );
#endif
#endif
        }
        if (j==0) {
           if (GLOBAL(base,0) <= (2*(k-ks))) {
             Cin[base+b-1-ks] = C[base];
           }
        }
        if (j==1) {
           if (GLOBAL(base,1) <= (2*(k-ks))) {
             Ci2[base+b-1-ks] = C[base+1];
           }
        }
      }
    }
  }
   /***** STAGE3.5 ::  BlockE Comm***********/
  preproc = ((pid == 0)? (nprocs-1): (pid-1));
  basis   = ((pid == 0)? b         :  0     );
#ifdef MPILIB
  AISNC(&Cout[0]);
#endif
#ifdef MPILIB
  AISNC(&Cou2[0]);
#endif
  for(base=basis; base<medsize; base+=b ) {
   /*                                   b/2 */
    APUT(preproc,&Cin[base],&Cout[0],
          (base-basis)*sizeof(double), (b)*sizeof(double));
    APUT(preproc,&Ci2[base],&Cou2[0],
          (base-basis)*sizeof(double), (b)*sizeof(double));
  }
#ifdef MPILIB
  AISNC(&Cou2[0]); 
#endif
  AISNC(&Cout[0]); 

  /***** STAGE3.6 ::  BlockF Computation***********/
  for(ks=(b+1)/2+1; ks<b ; ks++){
/*
    up = 2*b+odd-2*ks;
*/
    up = b - 2*(ks- (b+1)/2) -odd;
    for(base=0; base<mednp; base+=b) {
      /* equivalent to j= b-2*(ks-(b+1)/2) */
      for(j=up;j<b-2;j++){
        if (GLOBAL(base,j) <= (2*(k-ks))) {
          C[base+j]=dpu*C[base+j]+ dpm*C[base+j+1]+dpd*C[base+j+2];
          St[base+j] =St[base+j]*oneoveru;
#ifdef AMERICAN
#ifdef CALL
          C[base+j] = MAX( C[base+j], (St[base+j]-K) );
#else
          C[base+j] = MAX( C[base+j], (K-St[base+j]) );
#endif
#endif
        }
      }
      if (GLOBAL(base,(b-2)) <= (2*(k-ks))) {
        C[base+b-2]=dpu*C[base+b-2]+dpm*C[base+b-1]+dpd*Cout[base+b-ks];
        St[base+b-2]=St[base+b-2]*oneoveru;
#ifdef AMERICAN
#ifdef CALL
        C[base+b-2]=MAX( C[base+b-2], (St[base+b-2]-K) );
#else
        C[base+b-2]=MAX( C[base+b-2], (K-St[base+b-2]) );
#endif
#endif
      }
      if (GLOBAL(base,(b-1)) <= (2*(k-ks))) {
        C[base+b-1]=dpu*C[base+b-1]+dpm*Cout[base+b-ks]+dpd*Cou2[base+b-ks];
        St[base+b-1]=St[base+b-1]*oneoveru;
#ifdef AMERICAN
#ifdef CALL
        C[base+b-1]=MAX( C[base+b-1], (St[base+b-1]-K) );
#else
        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 *)Cou2);
 free((void *)Ci2);
 free((void *)Cout);
 free((void *)Cin);
 free((void *)St);
 free((void *)C);
 free((void *)Sinit);


}
