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