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