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: }