1 static char help[] = "Tests PetscRandom functions.\n\n"; 2 3 #include <petscsys.h> 4 5 #define PETSC_MAXBSIZE 40 6 #define DATAFILENAME "ex2_stock.txt" 7 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; 17 18 PetscErrorCode readData(MPI_Comm,himaInfo *); 19 PetscReal mcVal(PetscReal, PetscReal, PetscReal, PetscReal, PetscReal); 20 void exchangeVal(PetscReal*, PetscReal*); 21 PetscReal basketPayoff(PetscReal[], PetscReal[], PetscInt, PetscReal,PetscReal, PetscReal[]); 22 PetscErrorCode stdNormalArray(PetscReal*, PetscInt,PetscRandom); 23 PetscInt divWork(PetscMPIInt, PetscInt, PetscMPIInt); 24 25 /* 26 Contributed by Xiaoyan Zeng <zengxia@iit.edu> and Liu, Kwong Ip" <kiliu@math.hkbu.edu.hk> 27 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 */ 31 32 int main(int argc, char *argv[]) 33 { 34 PetscReal r,dt; 35 PetscInt n; 36 unsigned long i,myNumSim,totalNumSim,numdim; 37 PetscReal *vol, *St0, x, totalx; 38 PetscMPIInt size,rank; 39 PetscReal *eps; 40 himaInfo hinfo; 41 PetscRandom ran; 42 PetscErrorCode ierr; 43 44 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 45 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&ran);CHKERRQ(ierr); 46 ierr = PetscRandomSetFromOptions(ran);CHKERRQ(ierr); 47 48 ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr); /* number of nodes */ 49 ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr); /* my ranking */ 50 51 hinfo.n = 31; 52 hinfo.r = 0.04; 53 hinfo.dt = 1.0/12; /* a month as a period */ 54 hinfo.totalNumSim = 1000; 55 56 ierr = PetscOptionsGetInt(NULL,NULL,"-num_of_stocks",&(hinfo.n),NULL);CHKERRQ(ierr); 57 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); 58 ierr = PetscOptionsGetReal(NULL,NULL,"-interest_rate",&(hinfo.r),NULL);CHKERRQ(ierr); 59 ierr = PetscOptionsGetReal(NULL,NULL,"-time_interval",&(hinfo.dt),NULL);CHKERRQ(ierr); 60 ierr = PetscOptionsGetInt(NULL,NULL,"-num_of_simulations",&(hinfo.totalNumSim),NULL);CHKERRQ(ierr); 61 62 n = hinfo.n; 63 r = hinfo.r; 64 dt = hinfo.dt; 65 totalNumSim = hinfo.totalNumSim; 66 ierr = PetscMalloc1(2*n+1,&hinfo.vol);CHKERRQ(ierr); 67 vol = hinfo.vol; 68 St0 = hinfo.St0 = hinfo.vol + n; 69 ierr = readData(PETSC_COMM_WORLD,&hinfo);CHKERRQ(ierr); 70 71 numdim = n*(n+1)/2; 72 if (numdim%2 == 1) numdim++; 73 ierr = PetscMalloc1(numdim,&eps);CHKERRQ(ierr); 74 75 myNumSim = divWork(rank,totalNumSim,size); 76 77 x = 0; 78 for (i=0; i<myNumSim; i++) { 79 ierr = stdNormalArray(eps,numdim,ran);CHKERRQ(ierr); 80 x += basketPayoff(vol,St0,n,r,dt,eps); 81 } 82 83 ierr = MPI_Reduce(&x, &totalx, 1, MPIU_REAL, MPIU_SUM,0,PETSC_COMM_WORLD);CHKERRQ(ierr); 84 /* payoff = exp(-r*dt*n)*(totalx/totalNumSim); 85 ierr = 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", 86 payoff,(int)(stop - start),"parallel",size,"processors",n,(int)(1/dt),r);CHKERRQ(ierr); */ 87 88 ierr = PetscFree(vol);CHKERRQ(ierr); 89 ierr = PetscFree(eps);CHKERRQ(ierr); 90 ierr = PetscRandomDestroy(&ran);CHKERRQ(ierr); 91 ierr = PetscFinalize(); 92 return ierr; 93 } 94 95 PetscErrorCode stdNormalArray(PetscReal *eps, PetscInt numdim, PetscRandom ran) 96 { 97 PetscInt i; 98 PetscScalar u1,u2; 99 PetscReal t; 100 PetscErrorCode ierr; 101 102 PetscFunctionBegin; 103 for (i=0; i<numdim; i+=2) { 104 ierr = PetscRandomGetValue(ran,&u1);CHKERRQ(ierr); 105 ierr = PetscRandomGetValue(ran,&u2);CHKERRQ(ierr); 106 107 t = PetscSqrtReal(-2*PetscLogReal(PetscRealPart(u1))); 108 eps[i] = t * PetscCosReal(2*PETSC_PI*PetscRealPart(u2)); 109 eps[i+1] = t * PetscSinReal(2*PETSC_PI*PetscRealPart(u2)); 110 } 111 PetscFunctionReturn(0); 112 } 113 114 115 PetscReal basketPayoff(PetscReal vol[], PetscReal St0[], PetscInt n, PetscReal r,PetscReal dt, PetscReal eps[]) 116 { 117 PetscReal Stk[PETSC_MAXBSIZE], temp; 118 PetscReal payoff; 119 PetscInt maxk,i,j; 120 PetscInt pointcount=0; 121 122 for (i=0;i<n;i++) Stk[i] = St0[i]; 123 124 for (i=0;i<n;i++) { 125 maxk = 0; 126 for (j=0;j<(n-i);j++) { 127 Stk[j] = mcVal(Stk[j],r,vol[j],dt,eps[pointcount++]); 128 if ((Stk[j]/St0[j]) > (Stk[maxk]/St0[maxk])) maxk = j; 129 } 130 exchangeVal(Stk+j-1,Stk+maxk); 131 exchangeVal(St0+j-1,St0+maxk); 132 exchangeVal(vol+j-1,vol+maxk); 133 } 134 135 payoff = 0; 136 for (i=0; i<n; i++) { 137 temp = (Stk[i]/St0[i]) - 1; 138 if (temp > 0) payoff += temp; 139 } 140 return payoff; 141 } 142 143 PetscErrorCode readData(MPI_Comm comm,himaInfo *hinfo) 144 { 145 PetscInt i; 146 FILE *fd; 147 char temp[50]; 148 PetscErrorCode ierr; 149 PetscMPIInt rank; 150 PetscReal *v = hinfo->vol, *t = hinfo->St0; 151 PetscInt num=hinfo->n; 152 153 PetscFunctionBeginUser; 154 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 155 if (!rank) { 156 ierr = PetscFOpen(PETSC_COMM_SELF,DATAFILENAME,"r",&fd);CHKERRQ(ierr); 157 for (i=0;i<num;i++) { 158 double vv,tt; 159 if (fscanf(fd,"%s%lf%lf",temp,&vv,&tt) != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file\n"); 160 v[i] = vv; 161 t[i] = tt; 162 } 163 fclose(fd); 164 } 165 ierr = MPI_Bcast(v,2*num,MPIU_REAL,0,PETSC_COMM_WORLD);CHKERRQ(ierr); 166 /* ierr = 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]); */ 167 PetscFunctionReturn(0); 168 } 169 170 void exchangeVal(PetscReal *a, PetscReal *b) 171 { 172 PetscReal t; 173 174 t = *a; 175 *a = *b; 176 *b = t; 177 } 178 179 PetscReal mcVal(PetscReal St, PetscReal r, PetscReal vol, PetscReal dt, PetscReal eps) 180 { 181 return (St * PetscExpReal((r-0.5*vol*vol)*dt + vol*PetscSqrtReal(dt)*eps)); 182 } 183 184 PetscInt divWork(PetscMPIInt id, PetscInt num, PetscMPIInt size) 185 { 186 PetscInt numit; 187 188 numit = (PetscInt)(((PetscReal)num)/size); 189 numit++; 190 return numit; 191 } 192 193 /*TEST 194 195 build: 196 requires: !comple 197 output_file: output/ex1_1.out 198 199 test: 200 nsize: 2 201 output_file: output/ex1_1.out 202 localrunfiles: ex2_stock.txt 203 204 TEST*/ 205