Actual source code: ex2.c
petsc-3.6.1 2015-08-06
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: #if defined(PETSC_HAVE_SPRNG)
50: PetscRandomSetType(ran,PETSCSPRNG);
51: #elif defined(PETSC_HAVE_RAND)
52: PetscRandomSetType(ran,PETSCRAND);
53: #endif
54: PetscRandomSetFromOptions(ran);
56: MPI_Comm_size(PETSC_COMM_WORLD, &size); /* number of nodes */
57: MPI_Comm_rank(PETSC_COMM_WORLD, &rank); /* my ranking */
59: PetscOptionsHasName(NULL, "-check_generators", &flg);
60: if (flg) {
61: PetscRandomGetValue(ran,(PetscScalar*)&r);
62: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] rval: %g\n",rank,r);
63: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
64: }
66: hinfo.n = 31;
67: hinfo.r = 0.04;
68: hinfo.dt = 1.0/12; /* a month as a period */
69: hinfo.totalNumSim = 1000;
71: PetscOptionsGetInt(NULL,"-num_of_stocks",&(hinfo.n),NULL);
72: 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);
73: PetscOptionsGetReal(NULL,"-interest_rate",&(hinfo.r),NULL);
74: PetscOptionsGetReal(NULL,"-time_interval",&(hinfo.dt),NULL);
75: PetscOptionsGetInt(NULL,"-num_of_simulations",&(hinfo.totalNumSim),NULL);
77: n = hinfo.n;
78: r = hinfo.r;
79: dt = hinfo.dt;
80: totalNumSim = hinfo.totalNumSim;
81: vol = hinfo.vol = (PetscReal*)malloc(sizeof(PetscReal)*(2*n+1));
82: St0 = hinfo.St0 = hinfo.vol + n;
83: readData(PETSC_COMM_WORLD,&hinfo);
85: numdim = n*(n+1)/2;
86: if (numdim%2 == 1) numdim++;
87: eps = (PetscReal*)malloc(sizeof(PetscReal)*numdim);
89: myNumSim = divWork(rank,totalNumSim,size);
91: x = 0;
92: for (i=0; i<myNumSim; i++) {
93: stdNormalArray(eps,numdim,ran);
94: x += basketPayoff(vol,St0,n,r,dt,eps);
95: }
97: MPI_Reduce(&x, &totalx, 1, MPIU_REAL, MPIU_SUM,0,PETSC_COMM_WORLD);
98: /* payoff = exp(-r*dt*n)*(totalx/totalNumSim);
99: 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",
100: payoff,(int)(stop - start),"parallel",size,"processors",n,(int)(1/dt),r); */
102: free(vol);
103: free(eps);
104: PetscRandomDestroy(&ran);
105: PetscFinalize();
106: return 0;
107: }
109: void stdNormalArray(PetscReal *eps, PetscInt numdim, PetscRandom ran)
110: {
111: PetscInt i;
112: PetscReal u1,u2,t;
115: for (i=0; i<numdim; i+=2) {
116: PetscRandomGetValue(ran,(PetscScalar*)&u1);CHKERRABORT(PETSC_COMM_WORLD,ierr);
117: PetscRandomGetValue(ran,(PetscScalar*)&u2);CHKERRABORT(PETSC_COMM_WORLD,ierr);
119: t = PetscSqrtReal(-2*PetscLogReal(u1));
120: eps[i] = t * PetscCosReal(2*PETSC_PI*u2);
121: eps[i+1] = t * PetscSinReal(2*PETSC_PI*u2);
122: }
123: }
126: PetscReal basketPayoff(PetscReal vol[], PetscReal St0[], PetscInt n, PetscReal r,PetscReal dt, PetscReal eps[])
127: {
128: PetscReal Stk[PETSC_MAXBSIZE], temp;
129: PetscReal payoff;
130: PetscInt maxk,i,j;
131: PetscInt pointcount=0;
133: for (i=0;i<n;i++) Stk[i] = St0[i];
135: for (i=0;i<n;i++) {
136: maxk = 0;
137: for (j=0;j<(n-i);j++) {
138: Stk[j] = mcVal(Stk[j],r,vol[j],dt,eps[pointcount++]);
139: if ((Stk[j]/St0[j]) > (Stk[maxk]/St0[maxk])) maxk = j;
140: }
141: exchange(Stk+j-1,Stk+maxk);
142: exchange(St0+j-1,St0+maxk);
143: exchange(vol+j-1,vol+maxk);
144: }
146: payoff = 0;
147: for (i=0; i<n; i++) {
148: temp = (Stk[i]/St0[i]) - 1;
149: if (temp > 0) payoff += temp;
150: }
151: return payoff;
152: }
156: PetscErrorCode readData(MPI_Comm comm,himaInfo *hinfo)
157: {
158: PetscInt i;
159: FILE *fd;
160: char temp[50];
162: PetscMPIInt rank;
163: PetscReal *v = hinfo->vol, *t = hinfo->St0;
164: PetscInt num=hinfo->n;
167: MPI_Comm_rank(comm,&rank);
168: if (!rank) {
169: PetscFOpen(PETSC_COMM_SELF,DATAFILENAME,"r",&fd);
170: for (i=0;i<num;i++) {
171: double vv,tt;
172: fscanf(fd,"%s%lf%lf",temp,&vv,&tt);
173: v[i] = vv;
174: t[i] = tt;
175: }
176: fclose(fd);
177: }
178: MPI_Bcast(v,2*num,MPIU_REAL,0,PETSC_COMM_WORLD);
179: /* 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]); */
180: return(0);
181: }
183: void exchange(PetscReal *a, PetscReal *b)
184: {
185: PetscReal t;
187: t = *a;
188: *a = *b;
189: *b = t;
190: }
192: PetscReal mcVal(PetscReal St, PetscReal r, PetscReal vol, PetscReal dt, PetscReal eps)
193: {
194: return (St * PetscExpReal((r-0.5*vol*vol)*dt + vol*PetscSqrtReal(dt)*eps));
195: }
197: PetscInt divWork(PetscMPIInt id, PetscInt num, PetscMPIInt size)
198: {
199: PetscInt numit;
201: numit = (PetscInt)(((PetscReal)num)/size);
202: numit++;
203: return numit;
204: }