Actual source code: ex2.c

petsc-3.7.3 2016-07-24
Report Typos and Errors
  1: static char help[] = "Tests PetscRandom functions.\n\n";

  3: #include <petscsys.h>

  5: #define PETSC_MAXBSIZE     40
  6: #define DATAFILENAME "ex2_stock.txt"

  8: struct himaInfoTag {
  9:   PetscInt    n;
 10:   PetscReal   r;
 11:   PetscReal   dt;
 12:   PetscInt    totalNumSim;
 13:   PetscReal   *St0;
 14:   PetscReal   *vol;
 15: };
 16: typedef struct himaInfoTag himaInfo;

 18: PetscErrorCode readData(MPI_Comm comm,himaInfo *hinfo);
 19: PetscReal mcVal(PetscReal St, PetscReal r, PetscReal vol, PetscReal dt, PetscReal eps);
 20: void exchange(PetscReal *a, PetscReal *b);
 21: PetscReal basketPayoff(PetscReal vol[], PetscReal St0[], PetscInt n, PetscReal r,PetscReal dt, PetscReal eps[]);
 22: void stdNormalArray(PetscReal *eps, PetscInt numdim,PetscRandom ran);
 23: PetscInt divWork(PetscMPIInt id, PetscInt num, PetscMPIInt size);

 25: /*
 26:    Contributed by Xiaoyan Zeng <zengxia@iit.edu> and Liu, Kwong Ip" <kiliu@math.hkbu.edu.hk>

 28:    Example of usage:
 29:      mpiexec -n 4 ./ex2 -num_of_stocks 30 -interest_rate 0.4 -time_interval 0.01 -num_of_simulations 10000
 30: */

 34: int main(int argc, char *argv[])
 35: {
 36:   PetscReal      r,dt;
 37:   PetscInt       n;
 38:   unsigned long  i,myNumSim,totalNumSim,numdim;
 39:   PetscReal      *vol, *St0, x, totalx;
 40:   PetscMPIInt    size,rank;
 41:   PetscReal      *eps;
 42:   himaInfo       hinfo;
 43:   PetscRandom    ran;
 45:   PetscBool      flg;

 47:   PetscInitialize(&argc,&argv,(char*)0,help);
 48:   PetscRandomCreate(PETSC_COMM_WORLD,&ran);
 49:   PetscRandomSetFromOptions(ran);

 51:   MPI_Comm_size(PETSC_COMM_WORLD, &size);       /* number of nodes */
 52:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);     /* my ranking */

 54:   PetscOptionsHasName(NULL,NULL, "-check_generators", &flg);
 55:   if (flg) {
 56:     PetscRandomGetValue(ran,(PetscScalar*)&r);
 57:     PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] rval: %g\n",rank,r);
 58:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
 59:   }

 61:   hinfo.n           = 31;
 62:   hinfo.r           = 0.04;
 63:   hinfo.dt          = 1.0/12;   /* a month as a period */
 64:   hinfo.totalNumSim = 1000;

 66:   PetscOptionsGetInt(NULL,NULL,"-num_of_stocks",&(hinfo.n),NULL);
 67:   if (hinfo.n <1 || hinfo.n > 31) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only 31 stocks listed in stock.txt. num_of_stocks %D must between 1 and 31",hinfo.n);
 68:   PetscOptionsGetReal(NULL,NULL,"-interest_rate",&(hinfo.r),NULL);
 69:   PetscOptionsGetReal(NULL,NULL,"-time_interval",&(hinfo.dt),NULL);
 70:   PetscOptionsGetInt(NULL,NULL,"-num_of_simulations",&(hinfo.totalNumSim),NULL);

 72:   n           = hinfo.n;
 73:   r           = hinfo.r;
 74:   dt          = hinfo.dt;
 75:   totalNumSim = hinfo.totalNumSim;
 76:   vol         = hinfo.vol = (PetscReal*)malloc(sizeof(PetscReal)*(2*n+1));
 77:   St0         = hinfo.St0 = hinfo.vol + n;
 78:   readData(PETSC_COMM_WORLD,&hinfo);

 80:   numdim = n*(n+1)/2;
 81:   if (numdim%2 == 1) numdim++;
 82:   eps = (PetscReal*)malloc(sizeof(PetscReal)*numdim);

 84:   myNumSim = divWork(rank,totalNumSim,size);

 86:   x = 0;
 87:   for (i=0; i<myNumSim; i++) {
 88:     stdNormalArray(eps,numdim,ran);
 89:     x += basketPayoff(vol,St0,n,r,dt,eps);
 90:   }

 92:   MPI_Reduce(&x, &totalx, 1, MPIU_REAL, MPIU_SUM,0,PETSC_COMM_WORLD);
 93:   /* payoff = exp(-r*dt*n)*(totalx/totalNumSim);
 94:   PetscPrintf(PETSC_COMM_WORLD,"Option price = $%.3f using %ds of %s computation with %d %s for %d stocks, %d trading period per year, %.2f%% interest rate\n",
 95:    payoff,(int)(stop - start),"parallel",size,"processors",n,(int)(1/dt),r); */

 97:   free(vol);
 98:   free(eps);
 99:   PetscRandomDestroy(&ran);
100:   PetscFinalize();
101:   return 0;
102: }

104: void stdNormalArray(PetscReal *eps, PetscInt numdim, PetscRandom ran)
105: {
106:   PetscInt       i;
107:   PetscReal      u1,u2,t;

110:   for (i=0; i<numdim; i+=2) {
111:     PetscRandomGetValue(ran,(PetscScalar*)&u1);CHKERRABORT(PETSC_COMM_WORLD,ierr);
112:     PetscRandomGetValue(ran,(PetscScalar*)&u2);CHKERRABORT(PETSC_COMM_WORLD,ierr);

114:     t        = PetscSqrtReal(-2*PetscLogReal(u1));
115:     eps[i]   = t * PetscCosReal(2*PETSC_PI*u2);
116:     eps[i+1] = t * PetscSinReal(2*PETSC_PI*u2);
117:   }
118: }


121: PetscReal basketPayoff(PetscReal vol[], PetscReal St0[], PetscInt n, PetscReal r,PetscReal dt, PetscReal eps[])
122: {
123:   PetscReal Stk[PETSC_MAXBSIZE], temp;
124:   PetscReal payoff;
125:   PetscInt  maxk,i,j;
126:   PetscInt  pointcount=0;

128:   for (i=0;i<n;i++) Stk[i] = St0[i];

130:   for (i=0;i<n;i++) {
131:     maxk = 0;
132:     for (j=0;j<(n-i);j++) {
133:       Stk[j] = mcVal(Stk[j],r,vol[j],dt,eps[pointcount++]);
134:       if ((Stk[j]/St0[j]) > (Stk[maxk]/St0[maxk])) maxk = j;
135:     }
136:     exchange(Stk+j-1,Stk+maxk);
137:     exchange(St0+j-1,St0+maxk);
138:     exchange(vol+j-1,vol+maxk);
139:   }

141:   payoff = 0;
142:   for (i=0; i<n; i++) {
143:     temp = (Stk[i]/St0[i]) - 1;
144:     if (temp > 0) payoff += temp;
145:   }
146:   return payoff;
147: }

151: PetscErrorCode readData(MPI_Comm comm,himaInfo *hinfo)
152: {
153:   PetscInt       i;
154:   FILE           *fd;
155:   char           temp[50];
157:   PetscMPIInt    rank;
158:   PetscReal      *v = hinfo->vol, *t = hinfo->St0;
159:   PetscInt       num=hinfo->n;

162:   MPI_Comm_rank(comm,&rank);
163:   if (!rank) {
164:     PetscFOpen(PETSC_COMM_SELF,DATAFILENAME,"r",&fd);
165:     for (i=0;i<num;i++) {
166:       double vv,tt;
167:       fscanf(fd,"%s%lf%lf",temp,&vv,&tt);
168:       v[i] = vv;
169:       t[i] = tt;
170:     }
171:     fclose(fd);
172:   }
173:   MPI_Bcast(v,2*num,MPIU_REAL,0,PETSC_COMM_WORLD);
174:   /* PetscPrintf(PETSC_COMM_SELF,"[%d] vol %g, ... %g; St0 %g, ... %g\n",rank,hinfo->vol[0],hinfo->vol[num-1],hinfo->St0 [0],hinfo->St0[num-1]); */
175:   return(0);
176: }

178: void exchange(PetscReal *a, PetscReal *b)
179: {
180:   PetscReal t;

182:   t  = *a;
183:   *a = *b;
184:   *b = t;
185: }

187: PetscReal mcVal(PetscReal St, PetscReal r, PetscReal vol, PetscReal dt, PetscReal eps)
188: {
189:   return (St * PetscExpReal((r-0.5*vol*vol)*dt + vol*PetscSqrtReal(dt)*eps));
190: }

192: PetscInt divWork(PetscMPIInt id, PetscInt num, PetscMPIInt size)
193: {
194:   PetscInt numit;

196:   numit = (PetscInt)(((PetscReal)num)/size);
197:   numit++;
198:   return numit;
199: }