Got the kernel's nose out of cmdline args
[akaros.git] / user / apps / parlib / pthread / blackscholes.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59 // Copyright (c) 2007 Intel Corp.
60
61 // Black-Scholes
62 // Analytical method for calculating European Options
63 //
64 // 
65 // Reference Source: Options, Futures, and Other Derivatives, 3rd Edition, Prentice 
66 // Hall, John C. Hull,
67
68 #include <stdio.h>
69 #include <stdlib.h>
70 #include <math.h>
71 #include <string.h>
72
73 #define ENABLE_THREADS
74
75 #ifdef ENABLE_PARSEC_HOOKS
76 #include <hooks.h>
77 #endif
78
79 // Multi-threaded pthreads header
80 #ifdef ENABLE_THREADS
81 #define MAX_THREADS 128
82 // Add the following line so that icc 9.0 is compatible with pthread lib.
83 #define __thread __threadp  
84
85 #ifdef _XOPEN_SOURCE
86 #undef _XOPEN_SOURCE
87 #define _XOPEN_SOURCE 700
88 #endif
89 #ifndef _GNU_SOURCE
90 #define _GNU_SOURCE
91 #endif
92 #ifndef __USE_XOPEN2K
93 #define __USE_XOPEN2K
94 #endif
95 #ifndef __USE_UNIX98
96 #define __USE_UNIX98
97 #endif
98 #include <pthread.h>
99 #include <time.h>
100
101 pthread_t _M4_threadsTable[MAX_THREADS];
102 pthread_mutexattr_t _M4_normalMutexAttr;
103 int _M4_numThreads = MAX_THREADS;
104
105 #undef __thread
106 pthread_barrier_t barrier;;
107 #endif
108
109 // Multi-threaded OpenMP header
110 #ifdef ENABLE_OPENMP
111 #include <omp.h>
112 #endif
113
114 // Multi-threaded header for Windows
115 #ifdef WIN32
116 #pragma warning(disable : 4305)
117 #pragma warning(disable : 4244)
118 #include <windows.h>
119 #define MAX_THREADS 128
120 #endif
121
122 //Precision to use for calculations
123 #define fptype float
124
125 #define NUM_RUNS 100
126
127 typedef struct OptionData_ {
128         fptype s;          // spot price
129         fptype strike;     // strike price
130         fptype r;          // risk-free interest rate
131         fptype divq;       // dividend rate
132         fptype v;          // volatility
133         fptype t;          // time to maturity or option expiration in years 
134                            //     (1yr = 1.0, 6mos = 0.5, 3mos = 0.25, ..., etc)  
135         const char *OptionType;  // Option type.  "P"=PUT, "C"=CALL
136         fptype divs;       // dividend vals (not used in this test)
137         fptype DGrefval;   // DerivaGem Reference Value
138 } OptionData;
139
140 //++  modified by xiongww
141
142 OptionData data_init[] = {
143     #include "optionData.txt"
144 };
145
146 OptionData *data;
147 fptype *prices;
148 int numOptions;
149
150 //--
151
152 int    * otype;
153 fptype * sptprice;
154 fptype * strike;
155 fptype * rate;
156 fptype * volatility;
157 fptype * otime;
158 int numError = 0;
159 int nThreads;
160
161 ////////////////////////////////////////////////////////////////////////////////
162 ////////////////////////////////////////////////////////////////////////////////
163 ///////////////////////////////////////////////////////////////////////////////
164 ////////////////////////////////////////////////////////////////////////////////
165 // Cumulative Normal Distribution Function
166 // See Hull, Section 11.8, P.243-244
167 #define inv_sqrt_2xPI 0.39894228040143270286
168
169 fptype CNDF ( fptype InputX ) 
170 {
171     int sign;
172
173     fptype OutputX;
174     fptype xInput;
175     fptype xNPrimeofX;
176     fptype expValues;
177     fptype xK2;
178     fptype xK2_2, xK2_3;
179     fptype xK2_4, xK2_5;
180     fptype xLocal, xLocal_1;
181     fptype xLocal_2, xLocal_3;
182
183     // Check for negative value of InputX
184     if (InputX < 0.0) {
185         InputX = -InputX;
186         sign = 1;
187     } else 
188         sign = 0;
189
190     xInput = InputX;
191  
192     // Compute NPrimeX term common to both four & six decimal accuracy calcs
193     expValues = exp(-0.5f * InputX * InputX);
194     xNPrimeofX = expValues;
195     xNPrimeofX = xNPrimeofX * inv_sqrt_2xPI;
196
197     xK2 = 0.2316419 * xInput;
198     xK2 = 1.0 + xK2;
199     xK2 = 1.0 / xK2;
200     xK2_2 = xK2 * xK2;
201     xK2_3 = xK2_2 * xK2;
202     xK2_4 = xK2_3 * xK2;
203     xK2_5 = xK2_4 * xK2;
204     
205     xLocal_1 = xK2 * 0.319381530;
206     xLocal_2 = xK2_2 * (-0.356563782);
207     xLocal_3 = xK2_3 * 1.781477937;
208     xLocal_2 = xLocal_2 + xLocal_3;
209     xLocal_3 = xK2_4 * (-1.821255978);
210     xLocal_2 = xLocal_2 + xLocal_3;
211     xLocal_3 = xK2_5 * 1.330274429;
212     xLocal_2 = xLocal_2 + xLocal_3;
213
214     xLocal_1 = xLocal_2 + xLocal_1;
215     xLocal   = xLocal_1 * xNPrimeofX;
216     xLocal   = 1.0 - xLocal;
217
218     OutputX  = xLocal;
219     
220     if (sign) {
221         OutputX = 1.0 - OutputX;
222     }
223     
224     return OutputX;
225
226
227 // For debugging
228 void print_xmm(fptype in, char* s) {
229     printf("%s: %f\n", s, in);
230 }
231
232 //////////////////////////////////////////////////////////////////////////////////////
233 //////////////////////////////////////////////////////////////////////////////////////
234 //////////////////////////////////////////////////////////////////////////////////////
235 //////////////////////////////////////////////////////////////////////////////////////
236 fptype BlkSchlsEqEuroNoDiv( fptype sptprice,
237                             fptype strike, fptype rate, fptype volatility,
238                             fptype time, int otype, float timet )    
239 {
240     fptype OptionPrice;
241
242     // local private working variables for the calculation
243     fptype xStockPrice;
244     fptype xStrikePrice;
245     fptype xRiskFreeRate;
246     fptype xVolatility;
247     fptype xTime;
248     fptype xSqrtTime;
249
250     fptype logValues;
251     fptype xLogTerm;
252     fptype xD1; 
253     fptype xD2;
254     fptype xPowerTerm;
255     fptype xDen;
256     fptype d1;
257     fptype d2;
258     fptype FutureValueX;
259     fptype NofXd1;
260     fptype NofXd2;
261     fptype NegNofXd1;
262     fptype NegNofXd2;    
263     
264     xStockPrice = sptprice;
265     xStrikePrice = strike;
266     xRiskFreeRate = rate;
267     xVolatility = volatility;
268
269     xTime = time;
270     xSqrtTime = sqrt(xTime);
271
272     logValues = log( sptprice / strike );
273         
274     xLogTerm = logValues;
275         
276     
277     xPowerTerm = xVolatility * xVolatility;
278     xPowerTerm = xPowerTerm * 0.5;
279         
280     xD1 = xRiskFreeRate + xPowerTerm;
281     xD1 = xD1 * xTime;
282     xD1 = xD1 + xLogTerm;
283
284     xDen = xVolatility * xSqrtTime;
285     xD1 = xD1 / xDen;
286     xD2 = xD1 -  xDen;
287
288     d1 = xD1;
289     d2 = xD2;
290     
291     NofXd1 = CNDF( d1 );
292     NofXd2 = CNDF( d2 );
293
294     FutureValueX = strike * ( exp( -(rate)*(time) ) );        
295     if (otype == 0) {            
296         OptionPrice = (sptprice * NofXd1) - (FutureValueX * NofXd2);
297     } else { 
298         NegNofXd1 = (1.0 - NofXd1);
299         NegNofXd2 = (1.0 - NofXd2);
300         OptionPrice = (FutureValueX * NegNofXd2) - (sptprice * NegNofXd1);
301     }
302     
303     return OptionPrice;
304 }
305
306 //////////////////////////////////////////////////////////////////////////////////////
307 //////////////////////////////////////////////////////////////////////////////////////
308 //////////////////////////////////////////////////////////////////////////////////////
309 //////////////////////////////////////////////////////////////////////////////////////
310 #ifdef WIN32
311 DWORD WINAPI bs_thread(LPVOID tid_ptr){
312 #else
313 int bs_thread(void *tid_ptr) {
314 #endif
315     int i, j;
316     fptype price;
317     fptype priceDelta;
318     int tid = *(int *)tid_ptr;
319     int start = tid * (numOptions / nThreads);
320     int end = start + (numOptions / nThreads);
321
322 //printf("Thread %d at barrier.\n", tid);
323
324 #ifdef ENABLE_THREADS
325     pthread_barrier_wait(&(barrier));
326 #endif
327
328
329     for (j=0; j<NUM_RUNS; j++) {
330 #ifdef ENABLE_OPENMP
331 #pragma omp parallel for
332         for (i=0; i<numOptions; i++) {
333 #else  //ENABLE_OPENMP
334         for (i=start; i<end; i++) {
335 #endif //ENABLE_OPENMP
336             /* Calling main function to calculate option value based on 
337              * Black & Sholes's equation.
338              */
339 //printf("Thread %d performing computation.( %d, %d ) \n", tid,j,i);
340
341             price = BlkSchlsEqEuroNoDiv( sptprice[i], strike[i],
342                                          rate[i], volatility[i], otime[i], 
343                                          otype[i], 0);
344             prices[i] = price;
345             
346 #ifdef ERR_CHK   
347             priceDelta = data[i].DGrefval - price;
348             if( fabs(priceDelta) >= 1e-4 ){
349                 printf("Error on %d. Computed=%.5f, Ref=%.5f, Delta=%.5f\n",
350                        i, price, data[i].DGrefval, priceDelta);
351                 numError ++;
352             }
353 #endif
354         }
355     }
356
357     return 0;
358 }
359
360 int main (int argc, char **argv)
361 {
362     int i;
363     int loopnum;
364     fptype * buffer;
365     int * buffer2;
366     int initOptionNum;
367 #ifdef PARSEC_VERSION
368 #define __PARSEC_STRING(x) #x
369 #define __PARSEC_XSTRING(x) __PARSEC_STRING(x)
370         printf("PARSEC Benchmark Suite Version "__PARSEC_XSTRING(PARSEC_VERSION)"\n");
371         fflush(NULL);
372 #else
373         printf("PARSEC Benchmark Suite\n");
374         fflush(NULL);
375 #endif //PARSEC_VERSION
376 #ifdef ENABLE_PARSEC_HOOKS
377    __parsec_bench_begin(__parsec_blackscholes);
378 #endif
379
380    if (argc != 3)
381         {
382                 printf("Usage:\n\t%s <nthreads> <numOptions>\n", argv[0]);
383                 exit(1);
384         }
385     nThreads = atoi(argv[1]);
386     numOptions = atoi(argv[2]);
387
388     if(nThreads > numOptions) {
389       printf("WARNING: Not enough work, reducing number of threads to match number of options.");
390       nThreads = numOptions;
391     }
392
393 #if !defined(ENABLE_THREADS) && !defined(ENABLE_OPENMP)
394     if(nThreads != 1) {
395         printf("Error: <nthreads> must be 1 (serial version)\n");
396         exit(1);
397     }
398 #endif
399
400     // alloc spaces for the option data
401     data = (OptionData*)malloc(numOptions*sizeof(OptionData));
402     prices = (fptype*)malloc(numOptions*sizeof(fptype));
403     // initialize the data array
404     initOptionNum =  ( (sizeof(data_init)) / sizeof(OptionData) );
405     for ( loopnum = 0; loopnum < numOptions; ++ loopnum )   
406     {
407         OptionData *temp = data_init + loopnum%initOptionNum;
408         data[loopnum].OptionType = (const char*)malloc(strlen(temp->OptionType)+1);
409         strcpy((char *)(data[loopnum].OptionType), temp->OptionType);
410         data[loopnum].s = temp->s;
411         data[loopnum].strike = temp->strike;
412         data[loopnum].r = temp->r;
413         data[loopnum].divq = temp->divq;
414         data[loopnum].v = temp->v;
415         data[loopnum].t = temp->t;
416         data[loopnum].divs = temp->divs;
417         data[loopnum].DGrefval = temp->DGrefval;
418     }
419 #ifdef ENABLE_THREADS
420     
421     pthread_mutexattr_init( &_M4_normalMutexAttr);
422 //    pthread_mutexattr_settype( &_M4_normalMutexAttr, PTHREAD_MUTEX_NORMAL);
423     _M4_numThreads = nThreads;
424     {
425         int _M4_i;
426         for ( _M4_i = 0; _M4_i < MAX_THREADS; _M4_i++) {
427             _M4_threadsTable[_M4_i] = pthread_self();
428         }
429     }
430 ;
431     pthread_barrier_init(&(barrier),NULL,_M4_numThreads);;
432 #endif
433     printf("Num of Options: %d\n", numOptions);
434     printf("Num of Runs: %d\n", NUM_RUNS);
435
436 #define PAD 256
437 #define LINESIZE 64
438    
439     buffer = (fptype *) malloc(5 * numOptions * sizeof(fptype) + PAD);
440     sptprice = (fptype *) (((unsigned long)buffer + PAD) & ~(LINESIZE - 1));
441     strike = sptprice + numOptions;
442     rate = strike + numOptions;
443     volatility = rate + numOptions;
444     otime = volatility + numOptions;
445     
446     buffer2 = (int *) malloc(numOptions * sizeof(fptype) + PAD);
447     otype = (int *) (((unsigned long)buffer2 + PAD) & ~(LINESIZE - 1));
448     
449     for (i=0; i<numOptions; i++) {
450         otype[i]      = (data[i].OptionType[0] == 'P') ? 1 : 0;
451         sptprice[i]   = data[i].s;
452         strike[i]     = data[i].strike;
453         rate[i]       = data[i].r;
454         volatility[i] = data[i].v;    
455         otime[i]      = data[i].t;
456     }
457     
458     printf("Size of data: %d\n", numOptions * (sizeof(OptionData) + sizeof(int)));
459
460 #ifdef ENABLE_PARSEC_HOOKS
461     __parsec_roi_begin();
462 #endif
463 #ifdef ENABLE_THREADS
464     int tids[nThreads];
465     for(i=0; i<nThreads-1; i++) {
466 //printf("Creating thread %d of %d\n", i, nThreads);
467         tids[i]=i;
468         
469     {
470         int _M4_i;
471         printf("create with arg\n");
472         for ( _M4_i = 0; _M4_i < MAX_THREADS; _M4_i++) {
473             if ( _M4_threadsTable[_M4_i] == pthread_self())    break;
474         }
475                  printf("tid=%d %d\n", _M4_i, _M4_threadsTable[_M4_i]);
476         pthread_create(&_M4_threadsTable[_M4_i],NULL,(void *(*)(void *))bs_thread,(void *)&tids[i]);
477     }
478 ;
479     }
480     tids[nThreads-1] = nThreads-1;
481     bs_thread(&tids[nThreads-1]);
482     
483     {
484         int _M4_i;
485         void *_M4_ret;
486         for ( _M4_i = 0; _M4_i < MAX_THREADS;_M4_i++) {
487             printf("tid=%d %d\n", _M4_i, _M4_threadsTable[_M4_i]);
488             if ( _M4_threadsTable[_M4_i] == pthread_self())    break;
489             pthread_join( _M4_threadsTable[_M4_i], &_M4_ret);
490         }
491     }
492 ;
493 #else//ENABLE_THREADS
494 #ifdef ENABLE_OPENMP
495     {
496         int tid=0;
497         omp_set_num_threads(nThreads);
498         bs_thread(&tid);
499     }
500 #else //ENABLE_OPENMP
501 #ifdef WIN32 
502     if (nThreads > 1)
503     {
504         HANDLE threads[MAX_THREADS];
505                 int nums[MAX_THREADS];
506                 for(i=0; i<nThreads; i++) {
507                         nums[i] = i;
508                         threads[i] = CreateThread(0, 0, bs_thread, &nums[i], 0, 0);
509                 }
510                 WaitForMultipleObjects(nThreads, threads, TRUE, INFINITE);
511     } else
512 #endif
513     {
514         int tid=0;
515         bs_thread(&tid);
516     }
517 #endif //ENABLE_OPENMP
518 #endif //ENABLE_THREADS
519 #ifdef ENABLE_PARSEC_HOOKS
520     __parsec_roi_end();
521 #endif
522
523 #ifdef ERR_CHK
524     printf("Num Errors: %d\n", numError);
525 #endif
526     for (loopnum = 0; loopnum < numOptions; ++ loopnum)
527     {
528         free((void *)(data[loopnum].OptionType));
529     }
530     free(data);
531     free(prices);
532
533 #ifdef ENABLE_PARSEC_HOOKS
534     __parsec_bench_end();
535 #endif
536
537     return 0;
538 }
539