Actual source code: ex2.c

petsc-3.3-p7 2013-05-11
  1: static char help[] = "Tests PetscRandom functions.\n\n";

  3: #include <stdlib.h>
  4: #include <stdio.h>
  5: #include <string.h>
  6: #include <time.h>
  7: #include <sys/types.h>

  9: #include <petscsys.h>

 11: #define PETSC_MAXBSIZE     40
 12: #define PI           3.1415926535897
 13: #define DATAFILENAME "ex2_stock.txt"

 15: struct himaInfoTag {
 16:   int           n;
 17:   double        r;
 18:   double        dt;
 19:   int           totalNumSim;
 20:   double        *St0;
 21:   double        *vol;
 22: };
 23: typedef struct himaInfoTag himaInfo;

 25: /* function protype */
 26: PetscErrorCode readData(MPI_Comm comm,himaInfo *hinfo);
 27: double mcVal(double St, double r, double vol, double dt, double eps);
 28: void exchange(double *a, double *b);
 29: double basketPayoff(double vol[], double St0[], int n, double r,double dt, double eps[]);
 30: void stdNormalArray(double *eps, int size,PetscRandom ran);
 31: unsigned long divWork(int id, unsigned long num, int np);

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

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

 42: int main(int argc, char *argv[])
 43: {
 44:     double         r,dt;
 45:     int            n;
 46:     unsigned long  i,myNumSim,totalNumSim,numdim;
 47:     /* double         payoff; */
 48:     double         *vol, *St0, x, totalx;
 49:     int            np,myid;
 50:     time_t         start,stop;
 51:     double         *eps;
 52:     himaInfo       hinfo;
 53:     PetscRandom    ran;
 55:     PetscBool      flg;

 57:     PetscInitialize(&argc,&argv,(char *)0,help);
 58: #if defined(PETSC_USE_COMPLEX)
 59:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This example does not work for scalar type complex\n");
 60: #endif
 61:     time(&start);
 62:     PetscRandomCreate(PETSC_COMM_WORLD,&ran);
 63: #if defined(PETSC_HAVE_SPRNG)
 64:     PetscRandomSetType(ran,PETSCSPRNG);
 65: #elif defined(PETSC_HAVE_RAND)
 66:     PetscRandomSetType(ran,PETSCRAND);
 67: #endif
 68:     PetscRandomSetFromOptions(ran);

 70:     MPI_Comm_size(PETSC_COMM_WORLD, &np);     /* number of nodes */
 71:     MPI_Comm_rank(PETSC_COMM_WORLD, &myid);   /* my ranking */

 73:     PetscOptionsHasName(PETSC_NULL, "-check_generators", &flg);
 74:     if (flg){
 75:       PetscRandomGetValue(ran,(PetscScalar *)&r);
 76:       PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] rval: %g\n",myid,r);
 77:       PetscSynchronizedFlush(PETSC_COMM_WORLD);
 78:     }
 79: 
 80:     hinfo.n           = 31;
 81:     hinfo.r           = 0.04;
 82:     hinfo.dt          = 1.0/12; /* a month as a period */
 83:     hinfo.totalNumSim = 1000;
 84:     PetscOptionsGetInt(PETSC_NULL,"-num_of_stocks",&(hinfo.n),PETSC_NULL);
 85:     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);
 86:     PetscOptionsGetReal(PETSC_NULL,"-interest_rate",&(hinfo.r),PETSC_NULL);
 87:     PetscOptionsGetReal(PETSC_NULL,"-time_interval",&(hinfo.dt),PETSC_NULL);
 88:     PetscOptionsGetInt(PETSC_NULL,"-num_of_simulations",&(hinfo.totalNumSim),PETSC_NULL);

 90:     n           = hinfo.n;
 91:     r           = hinfo.r;
 92:     dt          = hinfo.dt;
 93:     totalNumSim = hinfo.totalNumSim;
 94:     vol         = hinfo.vol = (double *)malloc(sizeof(double)*(2*n+1));
 95:     St0         = hinfo.St0 = hinfo.vol + n;
 96:     readData(PETSC_COMM_WORLD,&hinfo);

 98:     numdim = n*(n+1)/2;
 99:     if (numdim%2 == 1){
100:       numdim++;
101:     }
102:     eps = (double *)malloc(sizeof(double)*numdim);

104:     myNumSim = divWork(myid,totalNumSim,np);

106:     x = 0;
107:     for (i=0;i<myNumSim;i++){
108:         stdNormalArray(eps,numdim,ran);
109:         x += basketPayoff(vol,St0,n,r,dt,eps);
110:     }

112:     MPI_Reduce(&x, &totalx, 1, MPI_DOUBLE, MPI_SUM,0,PETSC_COMM_WORLD);
113:     time(&stop);
114:     /* payoff = exp(-r*dt*n)*(totalx/totalNumSim);
115:     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",
116:      payoff,(int)(stop - start),"parallel",np,"processors",n,(int)(1/dt),r); */
117: 
118:     free(vol);
119:     free(eps);
120:     PetscRandomDestroy(&ran);
121:     PetscFinalize();
122:     return 0;
123: }

125: void stdNormalArray(double *eps, int size, PetscRandom ran)
126: {
127:   int            i;
128:   double         u1,u2,t;

131:   for (i=0;i<size;i+=2){
132:     PetscRandomGetValue(ran,(PetscScalar*)&u1);CHKERRABORT(PETSC_COMM_WORLD,ierr);
133:     PetscRandomGetValue(ran,(PetscScalar*)&u2);CHKERRABORT(PETSC_COMM_WORLD,ierr);
134: 
135:     t = sqrt(-2*log(u1));
136:     eps[i] = t * cos(2*PI*u2);
137:     eps[i+1] = t * sin(2*PI*u2);
138:   }
139: }


142: double basketPayoff(double vol[], double St0[], int n, double r,double dt, double eps[])
143: {
144:   double Stk[PETSC_MAXBSIZE], temp;
145:   double payoff;
146:   int    maxk,i,j;
147:   int    pointcount=0;
148: 
149:   for (i=0;i<n;i++) {
150:     Stk[i] = St0[i];
151:   }

153:   for (i=0;i<n;i++){
154:     maxk = 0;
155:     for (j=0;j<(n-i);j++){
156:       Stk[j] = mcVal(Stk[j],r,vol[j],dt,eps[pointcount++]);
157:       if ((Stk[j]/St0[j]) > (Stk[maxk]/St0[maxk])){
158:         maxk = j;
159:       }
160:     }
161:     exchange(Stk+j-1,Stk+maxk);
162:     exchange(St0+j-1,St0+maxk);
163:     exchange(vol+j-1,vol+maxk);
164:   }
165: 
166:   payoff = 0;
167:   for (i=0;i<n;i++){
168:     temp = (Stk[i]/St0[i]) - 1 ;
169:     if (temp > 0) payoff += temp;
170:   }
171:   return payoff;
172: }

176: PetscErrorCode readData(MPI_Comm comm,himaInfo *hinfo)
177: {
178:   int            i;
179:   FILE           *fd;
180:   char           temp[50];
182:   PetscMPIInt    rank;
183:   double         *v = hinfo->vol, *t = hinfo->St0;
184:   int            num=hinfo->n;
185: 
187:   MPI_Comm_rank(comm,&rank);
188:   if (!rank){
189:     PetscFOpen(PETSC_COMM_SELF,DATAFILENAME,"r",&fd);
190:     for (i=0;i<num;i++){
191:       fscanf(fd,"%s%lf%lf",temp,v+i,t+i);
192:     }
193:     fclose(fd);
194:   }
195:   MPI_Bcast(v,2*num,MPI_DOUBLE,0,PETSC_COMM_WORLD);
196:   //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]);
197:   return(0);
198: }

200: void exchange(double *a, double *b)
201: {
202:   double t;
203: 
204:   t = *a;
205:   *a = *b;
206:   *b = t;
207: }

209: double mcVal(double St, double r, double vol, double dt, double eps)
210: {
211:   return (St * exp((r-0.5*vol*vol)*dt + vol*sqrt(dt)*eps));
212: }

214: unsigned long divWork(int id, unsigned long num, int np)
215: {
216:   unsigned long numit;

218:   numit = (unsigned long)(((double)num)/np);
219:   numit++;
220:   return numit;
221: }