/**********************************************************
**TrinomialTree Parallel OptionPriceValuations
**This is the testing function main() component;
**Copyright (c) 2001-2003   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  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 is offered that the code works as
**advertised or is absent of any, even damaging,
**side-effects; use this code on your own personal risk
**/
/* Test Data are available as comment at the end of this file */
#define WIN2000
#undef  WIN2000
#ifdef WIN2000 
#include "stdafx.h"
#include <winsock2.h>
#else
#include <sys/utsname.h>
#include <unistd.h>
#endif
#include <stdio.h>
#include <stdlib.h>
#include "allmbt.h"

#define GLOBAL(BASE,J) (((BASE/b)*nprocs*b)+pid*b+J)
#undef AMERICAN
#define AMERICAN
#define CALL
#undef CALL

#define PRINTMAT(X)   bspfold((void (*)(void *,void *,void *)) foldmin,(char *)&t1,(char *)&t3,sizeof(double)); bspfold((void (*)(void *,void *,void *)) foldmax,(char *)&t1,(char *)&t4,sizeof(double)); 


#ifdef   MPILIB
#define AIMPIWINMAX   128
MPI_Comm comm;
MPI_Win  win;
MPI_Win  aiwin[AIMPIWINMAX];
void*    bspmpi[AIMPIWINMAX];
int bspmpi_max = AIMPIWINMAX;
int bspmpi_first=1;
int bspmpicur_max = 0;
int mpi_pid = -1;
int mpi_nprocs = -1;
double mpi_bsp_l = 300.0;
double mpi_bsp_g = 0.5;
double mpi_basetime=0.0;
#endif

#ifdef MPILIB
int aisearch(void * address) {
 register int i;

 for (i=bspmpi_max; i >=0; i-- ){
   if (bspmpi[i] == address) {
        return(i);
   }
 }
 return(-1);
}

int aiinsert(void *address) {
 register int i;
 if (bspmpi_first == 1) {
    bspmpi_first =0;
    for(i=0;i<bspmpi_max;i++) {
          bspmpi[i]=NULL;
    }
 }
 for(i=bspmpicur_max; i>=0 ; i--){
  if (bspmpi[i] == NULL)  {
    bspmpi[i]=address;
    bspmpicur_max++;
/*
    if (AIPID()==10) fprintf(stdout,"Iroc %d i=%d %d %u\n",AIPID(),i,bspmpicur_max,address);fflush(stdout);
*/
    return (i);
  }
 }
 /* runout of space */
 return(-1);
}

int aidelete(void *address) {
 register int i;
 for(i=bspmpicur_max-1;i>=0;i--) {
   if (bspmpi[i]==address){
     bspmpi[i]= NULL;
     bspmpicur_max --;
/*
     if (AIPID()==10) fprintf(stdout,"Droc %d i=%d %d %u\n",AIPID(),i,bspmpicur_max,address);fflush(stdout);
*/
     return(i);
   }
 }
 return(-1);
}
#endif

double __my_estimate_l = (double) 500;
double __my_estimate_g = (double) 0.5;
int    __my_estimate_nprocs = (int) 1;
int    __my_first_call=(int) 0;
int    __my_estimate_brd=(int) BRD_CUTOFF;
int    __my_estimate_ppf=(int) PPF_CUTOFF;

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

void __my_estimate_bsp_values()
{
        __my_estimate_l = AIL();
        __my_estimate_g = AIG();
        __my_estimate_nprocs = AINPROCS();
        __my_estimate_brd = BRD_CUTOFF;
        __my_estimate_ppf = PPF_CUTOFF;
}

/* One superstep broadcasting to communicate command line args */
void bspbroad(int fromp, char *from, char *to, int size)
{
 int pid;
 int nprocs;
 int i;

 pid=AIPID(); /* processor id */
 nprocs=AINPROCS(); /* number of processors */
 AIREGISTER(to,size);
 AIINIT();
 AICOMMIT();
 for (i = 0; i < nprocs;i++) {
   if (pid==fromp)
        AIHPPUT(i,from,to,0,size);
 }
 AISYNC();
 AIBARRIER();
 AIDEREGISTER((void*)to);
}

/* operators for bspfold */

void foldmax(double *z, double *x, double *y)
{
  *z=(((*x)>=(*y))?(*x):(*y));
}

void foldmin(double *z, double *x, double *y)
{
  *z=(((*x)<(*y))?(*x):(*y));
}



/* One superstep folding operation to compute global min/max of timing data */
void 
  bspfold(void (*oper)(void *,void*,void*),char *from,char *to,int size)
{
 int  pid; 
 int  nprocs; 
 int  i;
 char *list;
 int  lsize;

 pid=AIPID(); 
 nprocs=AINPROCS(); 
 if (1==nprocs)  /* if local, copy */
    memcpy(to,from,size*sizeof(char));
 else {    /* if global, allocate memory O(p) per processor */
    lsize = nprocs*size; 
    list=(char *) malloc(lsize*sizeof(char));
    check_if_null(list,"bspfold","list");
    AIREGISTER(list,lsize);
    AIINIT();
    AICOMMIT();
 }
 /* each proc sends its from to every other proc; results stored in list */
 for (i = 0; i < nprocs;i++) {  
   AIHPPUT(i,from,list,pid*size,size);
 }
 AISYNC();
 AIBARRIER();
 /* apply operator to local list */
 for (i = 0; i < nprocs-1;i++) {
  oper(to,(list+i*size),(list+(i+1)*size));
  memcpy((list+(i+1)*size),to,size);
 }
    /* free memory and deregister */
 
    AIDEREGISTER(list);
    free((void *)list);
}


/* MAIN function */
int	main(int argc, char *argv[])
{
 int 	i,j,k,n,nprocs,runs,debug,pid,option;
 double t3,t4,t1,a,b,c;
 char hostname[512];
 double K, S,r,u,d;
 int    T,N,blocksize;
 double result,sigma,dx,dividend;
 double rst, rc;


 AIBEGIN(AINPROCS());
 pid = AIPID(); nprocs = AINPROCS();

 /* Show hostnames of processors; */
 gethostname(hostname,sizeof(hostname));
#undef SANITY_CHECK
#ifdef SANITY_CHECK
 for(i=0;i<nprocs;i++) {
   if (i==pid) {
      printf("Processor %d of %d [%s]\n",pid,nprocs,hostname);
      fflush(stdout);
   }
 }
#endif
 AIBARRIER();

  /* Command line processing */
 n=8;debug=0;runs=1; 
 if (0 == pid ){
   if (3 != argc){
     printf("Usage: %s blocksize horizon <tdata\n\
            Integer blocksize-1 must divide integer horizon\n\
            tdata is a parameter file\n",
            argv[0]);
     fflush(stdout);
     AIABORT("exit");
   } 
   else {
    blocksize=atoi(argv[1]);
    scanf("%lf",&K);
    scanf("%d",&T);
    scanf("%lf",&S);
    scanf("%lf",&r);
    scanf("%d",&N);
    scanf("%lf",&sigma);
    scanf("%lf",&dx);
    scanf("%lf",&dividend);
    N= atoi(argv[2]);
   }
 }
   dx= sigma*sqrt(3.0*T/N);

   /* Parallel Command line processing */
   AIREGISTER(&K,sizeof(double)); AIINIT();AICOMMIT();
      bspbroad(0,(char *)&K ,(char *)&K ,sizeof(double));
   AISYNC();AIBARRIER(); AIDEREGISTER(&K    );

   AIREGISTER(&T,sizeof(int   )); AIINIT();AICOMMIT();
      bspbroad(0,(char *)&T ,(char *)&T ,sizeof(int));
   AISYNC();AIBARRIER(); AIDEREGISTER(&T    );

   AIREGISTER(&blocksize,sizeof(int   )); AIINIT();AICOMMIT();
      bspbroad(0,(char *)&blocksize,(char *)&blocksize,sizeof(int));
   AISYNC();AIBARRIER(); AIDEREGISTER(&blocksize    );

   AIREGISTER(&S,sizeof(double)); AIINIT();AICOMMIT();
      bspbroad(0,(char *)&S ,(char *)&S ,sizeof(double));
   AISYNC();AIBARRIER(); AIDEREGISTER(&S    );

   AIREGISTER(&r,sizeof(double)); AIINIT();AICOMMIT();
      bspbroad(0,(char *)&r ,(char *)&r ,sizeof(double));
   AISYNC();AIBARRIER(); AIDEREGISTER(&r    );

   AIREGISTER(&sigma,sizeof(double)); AIINIT();AICOMMIT();
      bspbroad(0,(char *)&sigma ,(char *)&sigma ,sizeof(double));
   AISYNC();AIBARRIER(); AIDEREGISTER(&sigma    );

   AIREGISTER(&dx   ,sizeof(double)); AIINIT();AICOMMIT();
      bspbroad(0,(char *)&dx    ,(char *)&dx    ,sizeof(double));
   AISYNC();AIBARRIER(); AIDEREGISTER(&dx       );

   AIREGISTER(&dividend   ,sizeof(double)); AIINIT();AICOMMIT();
      bspbroad(0,(char *)&dividend,(char *)&dividend,sizeof(double));
   AISYNC();AIBARRIER(); AIDEREGISTER(&dividend);

   AIREGISTER(&N        ,sizeof(int   )); AIINIT();AICOMMIT();
      bspbroad(0,(char *)&N,(char *)&N,sizeof(int));
   AISYNC();AIBARRIER(); AIDEREGISTER(&N            );


#ifdef SANITY_CHECK
#endif
   if ((nprocs-1)==pid) {
     fprintf(stdout, "proc=%d ::K=%f T=%d S=%f r=%f N=%d, sigma=%f b=%d dx=%f dividend=%f\n", pid,K,T,S,r,N,sigma,blocksize,dx,dividend); 
     fflush(stdout);
   }

   AIBARRIER();
   t1=0.0;
   AIBARRIER();
   t3=AITIME();
   trt(K,T,S,r,N,sigma,dx,dividend,blocksize,&rst,&rc);
   AIBARRIER();
   t4=AITIME();
   t1 += (t4-t3);
   PRINTMAT(AmerPut );
   if (0 == pid) {
    printf("RunningTime= %8.5f-%8.5f (p=%3d,L=%5.1f,g=%5.2f,b=%4d)(%14.9f %10.7f)\n",(t3),(t4),AINPROCS(),AIL(),AIG(),blocksize,rc,rst);
   }
   AIEND();
   return(0);
}

/*** If tdata is
 *   100.0
 *1
 *100.0
 *0.06
 *6
 *0.2
 *0.2
 *0.03
 *
 *and you run  EXECUTABLE 256 8192 <tdata
 **The option prices that will be reported are OR SHOULD BE as follows.
 **EuropeanCallOption=    9.134966234 (S= 100.000000)
 **EuropeanPutOption=    6.266866238 (S= 100.000000)
 **AmericanCallOption=    9.134974973 (S= 100.000000)
 **AmericanPutOption=    6.620313528 (S= 100.000000)
 * The following values should also be reported by the sequential program toption.c/tomain.c
 * found in directory Seq along with the data file tdata
 */
/**********************************************************
#define GLOBAL(BASE,J) (((BASE/b)*nprocs*b)+pid*b+J)
#define   OPTION_DEBUG
#undef   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) {;}

/*trt: GeneralizedTrinomialTree Parallel option price valuations
 * K : strike price
 * T : time horizon
 * S : stock/asset price
 * r : interest rate
 * n : periods in horizon
 * sigma : sigma
 * dx:
 *dividend:
 * blocksize : BSP algorithmblocksize forthe problem(read algorithm)
 * rstock:stockpriceat origin of tree (value makes sense at proc 0)
 * roption: That's what you want! (at processor 0).
 */


void  trt(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;

 AIOREGISTER(St,medsize*sizeof(double));
 AIOINIT((char *)St);AIOCOMMIT();
 for(i=0;i<initnp;i+=size) {
   size= MIN((nelements),(initnp-i));
   AIOHPPUT((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;
 }
 AIOSYNC(St);
 AIOBARRIER();

/***** 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);
*/
  AIOREGISTER(&Cout[0],medsize*sizeof(double));
  AIOREGISTER(&Cou2[0],medsize*sizeof(double));
  AIOCOMMIT();
  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     );

  AIOSYNC(&Cout[0]);
  AIOSYNC(&Cou2[0]);
  for(base=basis; base<medsize; base+=b ) {
    /* redundant  &Cin[base+b/2]      (base-basis+b/2)  b-(b/2) */
    AIOPUT(preproc,&Cin[base+b/2],&Cout[0],(base-basis+b/2)*sizeof(double),
            (b-b/2)*sizeof(double));
    AIOPUT(preproc,&Ci2[base+b/2],&Cou2[0],(base-basis+b/2)*sizeof(double),
            (b-b/2)*sizeof(double));
  }
  AIOSYNC(&Cou2[0]); 
  AIOSYNC(&Cout[0]); 
AIOBARRIER();


   /***** 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     );
  AIOSYNC(&Cou2[0]);
  for(base=basis; base<medsize; base+=b ) {
    /* redundant;as before */
    AIOPUT(preproc,&Ci2[base+b/2],&Cou2[0],
         (base-basis+b/2)*sizeof(double),sizeof(double));
  }
  AIOSYNC(&Cou2[0]); 
AIOBARRIER();

  /***** 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     );
  AIOSYNC(&Cout[0]);
  AIOSYNC(&Cou2[0]);
  for(base=basis; base<medsize; base+=b ) {
   /*                                   b/2 */
    AIOPUT(preproc,&Cin[base],&Cout[0],
          (base-basis)*sizeof(double), (b)*sizeof(double));
    AIOPUT(preproc,&Ci2[base],&Cou2[0],
          (base-basis)*sizeof(double), (b)*sizeof(double));
  }
  AIOSYNC(&Cou2[0]); 
  AIOSYNC(&Cout[0]); 
AIOBARRIER();

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

}
