Actual source code: ex2.c
petsc-3.3-p7 2013-05-11
1: static const char help[] = "STREAM benchmark for pthread implemenentations\n\n";
3: /*-----------------------------------------------------------------------*/
4: /* Program: Stream */
5: /* Revision: $Id: stream.c,v 5.9 2009/04/11 16:35:00 mccalpin Exp mccalpin $ */
6: /* Original code developed by John D. McCalpin */
7: /* Programmers: John D. McCalpin */
8: /* Joe R. Zagar */
9: /* */
10: /* This program measures memory transfer rates in MB/s for simple */
11: /* computational kernels coded in C. */
12: /*-----------------------------------------------------------------------*/
13: /* Copyright 1991-2005: John D. McCalpin */
14: /*-----------------------------------------------------------------------*/
15: /* License: */
16: /* 1. You are free to use this program and/or to redistribute */
17: /* this program. */
18: /* 2. You are free to modify this program for your own use, */
19: /* including commercial use, subject to the publication */
20: /* restrictions in item 3. */
21: /* 3. You are free to publish results obtained from running this */
22: /* program, or from works that you derive from this program, */
23: /* with the following limitations: */
24: /* 3a. In order to be referred to as "STREAM benchmark results", */
25: /* published results must be in conformance to the STREAM */
26: /* Run Rules, (briefly reviewed below) published at */
27: /* http://www.cs.virginia.edu/stream/ref.html */
28: /* and incorporated herein by reference. */
29: /* As the copyright holder, John McCalpin retains the */
30: /* right to determine conformity with the Run Rules. */
31: /* 3b. Results based on modified source code or on runs not in */
32: /* accordance with the STREAM Run Rules must be clearly */
33: /* labelled whenever they are published. Examples of */
34: /* proper labelling include: */
35: /* "tuned STREAM benchmark results" */
36: /* "based on a variant of the STREAM benchmark code" */
37: /* Other comparable, clear and reasonable labelling is */
38: /* acceptable. */
39: /* 3c. Submission of results to the STREAM benchmark web site */
40: /* is encouraged, but not required. */
41: /* 4. Use of this program or creation of derived works based on this */
42: /* program constitutes acceptance of these licensing restrictions. */
43: /* 5. Absolutely no warranty is expressed or implied. */
44: /*-----------------------------------------------------------------------*/
45: # include <stdio.h>
46: # include <math.h>
47: # include <limits.h>
48: # include <float.h>
49: # include <sys/time.h>
50: # include <petscconf.h>
51: # include <petscsys.h>
52: # include <petscthreadcomm.h>
53: /* INSTRUCTIONS:
54: *
55: * 1) Stream requires a good bit of memory to run. Adjust the
56: * value of 'N' (below) to give a 'timing calibration' of
57: * at least 20 clock-ticks. This will provide rate estimates
58: * that should be good to about 5% precision.
59: */
61: #ifndef N
62: # define N 2000000
63: #endif
64: #ifndef NTIMES
65: # define NTIMES 50
66: #endif
67: #ifndef OFFSET
68: # define OFFSET 0
69: #endif
71: /*
72: * 3) Compile the code with full optimization. Many compilers
73: * generate unreasonably bad code before the optimizer tightens
74: * things up. If the results are unreasonably good, on the
75: * other hand, the optimizer might be too smart for me!
76: *
77: * Try compiling with:
78: * cc -O stream_omp.c -o stream_omp
79: *
80: * This is known to work on Cray, SGI, IBM, and Sun machines.
81: *
82: *
83: * 4) Mail the results to mccalpin@cs.virginia.edu
84: * Be sure to include:
85: * a) computer hardware model number and software revision
86: * b) the compiler flags
87: * c) all of the output from the test case.
88: * Thanks!
89: *
90: */
92: # define HLINE "-------------------------------------------------------------\n"
94: # ifndef MIN
95: # define MIN(x,y) ((x)<(y)?(x):(y))
96: # endif
97: # ifndef MAX
98: # define MAX(x,y) ((x)>(y)?(x):(y))
99: # endif
101: #ifndef STATIC_ALLOC
102: # define STATIC_ALLOC 1
103: #endif
105: #if STATIC_ALLOC
106: static double a[N+OFFSET],
107: b[N+OFFSET],
108: c[N+OFFSET];
109: #endif
111: static double avgtime[4] = {0}, maxtime[4] = {0},
112: mintime[4] = {FLT_MAX,FLT_MAX,FLT_MAX,FLT_MAX};
114: static const char *label[4] = {"Copy: ", "Scale: ","Add: ", "Triad: "};
116: static double bytes[4] = {
117: 2 * sizeof(double) * N,
118: 2 * sizeof(double) * N,
119: 3 * sizeof(double) * N,
120: 3 * sizeof(double) * N
121: };
123: double mysecond();
124: void checkSTREAMresults();
126: #ifndef WITH_PTHREADS
127: # define WITH_PTHREADS 1
128: #endif
130: #if !STATIC_ALLOC
131: double *a, *b, *c;
132: #endif
134: #if WITH_PTHREADS
135: PetscInt nworkThreads;
136: PetscInt *trstarts;
138: PetscErrorCode tuned_STREAM_2A_Kernel(PetscInt myrank)
139: {
140: int i;
141:
142: for(i=trstarts[myrank];i<trstarts[myrank+1];i++)
143: a[i] = 2.0E0*a[i];
144:
145: return(0);
146: }
148: PetscErrorCode tuned_STREAM_Initialize_Kernel(PetscInt myrank) {
149: int i;
151: for(i=trstarts[myrank];i<trstarts[myrank+1];i++) {
152: a[i] = 1.0;
153: b[i] = 2.0;
154: c[i] = 0.0;
155: }
156: return(0);
157: }
159: PetscErrorCode tuned_STREAM_Copy_Kernel(PetscInt myrank)
160: {
161: int j;
163: for (j=trstarts[myrank]; j<trstarts[myrank+1]; j++)
164: c[j] = a[j];
165: return(0);
166: }
168: PetscErrorCode tuned_STREAM_Scale_Kernel(PetscInt myrank,double *scalarp) {
169: double scalar = *scalarp;
170: int j;
172: for (j=trstarts[myrank]; j<trstarts[myrank+1]; j++)
173: b[j] = scalar*c[j];
174: return(0);
175: }
177: PetscErrorCode tuned_STREAM_Add_Kernel(PetscInt myrank) {
178: int j;
180: for (j=trstarts[myrank]; j<trstarts[myrank+1]; j++)
181: c[j] = a[j]+b[j];
183: return(0);
184: }
186: PetscErrorCode tuned_STREAM_Triad_Kernel(PetscInt myrank,double *scalarp) {
187: double scalar = *scalarp;
188: int j;
190: for (j=trstarts[myrank]; j<trstarts[myrank+1]; j++)
191: a[j] = b[j]+scalar*c[j];
193: return(0);
194: }
195: #endif
197: int main(int argc,char *argv[])
198: {
200: int quantum, checktick();
201: int BytesPerWord;
202: int j, k;
203: double scalar=3.0, t, times[4][NTIMES];
205: PetscInitialize(&argc,&argv,0,help);
206: /* --- SETUP --- determine precision and check timing --- */
207:
208: /*printf(HLINE);
209: printf("STREAM version $Revision: 5.9 $\n");
210: printf(HLINE); */
211: BytesPerWord = sizeof(double);
212: printf("This system uses %d bytes per DOUBLE PRECISION word.\n",BytesPerWord);
213:
214: printf(HLINE);
215: #ifdef NO_LONG_LONG
216: printf("Array size = %d, Offset = %d\n" , N, OFFSET);
217: #else
218: printf("Array size = %llu, Offset = %d\n", (unsigned long long) N, OFFSET);
219: #endif
220:
221: printf("Total memory required = %.1f MB.\n",(3.0 * BytesPerWord) * ( (double) N / 1048576.0));
222: printf("Each test is run %d times, but only\n", NTIMES);
223: printf("the *best* time for each is used.\n");
224:
225: printf(HLINE);
227: #if !STATIC_ALLOC
228: a = malloc((N+OFFSET)*sizeof(double));
229: b = malloc((N+OFFSET)*sizeof(double));
230: c = malloc((N+OFFSET)*sizeof(double));
231: #endif
232:
233: #if WITH_PTHREADS
234: PetscThreadCommGetNThreads(PETSC_COMM_WORLD,&nworkThreads);
235: PetscMalloc((nworkThreads+1)*sizeof(PetscInt),&trstarts);
236: PetscInt Q,R,nloc;
237: PetscBool S;
238: Q = (N+OFFSET)/nworkThreads;
239: R = (N+OFFSET) - Q*nworkThreads;
240: trstarts[0] = 0;
241: for(j=0;j < nworkThreads;j++) {
242: S = (PetscBool)(j < R);
243: nloc = S?Q+1:Q;
244: trstarts[j+1] = trstarts[j]+nloc;
245: }
246: PetscThreadCommRunKernel(PETSC_COMM_WORLD,(PetscThreadKernel)tuned_STREAM_Initialize_Kernel,1,&scalar);
247: PetscThreadCommBarrier(PETSC_COMM_WORLD);
248: # else
249: for (j=0; j<N; j++) {
250: a[j] = 1.0;
251: b[j] = 2.0;
252: c[j] = 0.0;
253: }
254: #endif
256: /*printf(HLINE);*/
257:
258: /* Get initial value for system clock. */
259: if ( (quantum = checktick()) >= 1)
260: ;
261: /* printf("Your clock granularity/precision appears to be %d microseconds.\n", quantum); */
262: else {
263: /* printf("Your clock granularity appears to be less than one microsecond.\n"); */
264: quantum = 1;
265: }
266:
267: t = mysecond();
268:
269: #if WITH_PTHREADS
270: PetscThreadCommRunKernel(PETSC_COMM_WORLD,(PetscThreadKernel)tuned_STREAM_2A_Kernel,0);
271: PetscThreadCommBarrier(PETSC_COMM_WORLD);
272: #else
273: for (j = 0; j < N; j++)
274: a[j] = 2.0E0 * a[j];
275: #endif
276: t = 1.0E6 * (mysecond() - t);
277:
278: /* printf("Each test below will take on the order of %d microseconds.\n", (int) t );
279: printf(" (= %d clock ticks)\n", (int) (t/quantum) );
280: printf("Increase the size of the arrays if this shows that\n");
281: printf("you are not getting at least 20 clock ticks per test.\n");
282:
283: printf(HLINE);
284: */
285: /* --- MAIN LOOP --- repeat test cases NTIMES times --- */
286:
287: for (k=0; k<NTIMES; k++) {
288: times[0][k] = mysecond();
289: #if WITH_PTHREADS
290: PetscThreadCommRunKernel(PETSC_COMM_WORLD,(PetscThreadKernel)tuned_STREAM_Copy_Kernel,0);
291: PetscThreadCommBarrier(PETSC_COMM_WORLD);
292: #else
293: for (j=0; j<N; j++)
294: c[j] = a[j];
295: #endif
296: times[0][k] = mysecond() - times[0][k];
297:
298: times[1][k] = mysecond();
299: #if WITH_PTHREADS
300: PetscThreadCommRunKernel(PETSC_COMM_WORLD,(PetscThreadKernel)tuned_STREAM_Scale_Kernel,1,&scalar);
301: PetscThreadCommBarrier(PETSC_COMM_WORLD);
302: #else
303: for (j=0; j<N; j++)
304: b[j] = scalar*c[j];
305: #endif
306: times[1][k] = mysecond() - times[1][k];
307:
308: times[2][k] = mysecond();
309: #if WITH_PTHREADS
310: PetscThreadCommRunKernel(PETSC_COMM_WORLD,(PetscThreadKernel)tuned_STREAM_Add_Kernel,0);
311: PetscThreadCommBarrier(PETSC_COMM_WORLD);
312: #else
313: for (j=0; j<N; j++)
314: c[j] = a[j]+b[j];
315: #endif
316: times[2][k] = mysecond() - times[2][k];
317:
318: times[3][k] = mysecond();
319: #if WITH_PTHREADS
320: PetscThreadCommRunKernel(PETSC_COMM_WORLD,(PetscThreadKernel)tuned_STREAM_Triad_Kernel,1,&scalar);
321: PetscThreadCommBarrier(PETSC_COMM_WORLD);
322: #else
323: for (j=0; j<N; j++)
324: a[j] = b[j]+scalar*c[j];
325: #endif
326: times[3][k] = mysecond() - times[3][k];
327: }
328:
329: /* --- SUMMARY --- */
330:
331: for (k=1; k<NTIMES; k++) { /* note -- skip first iteration */
332: for (j=0; j<4; j++) {
333: avgtime[j] = avgtime[j] + times[j][k];
334: mintime[j] = MIN(mintime[j], times[j][k]);
335: maxtime[j] = MAX(maxtime[j], times[j][k]);
336: }
337: }
338:
339: printf("Function Rate (MB/s) \n");
340: for (j=0; j<4; j++) {
341: avgtime[j] = avgtime[j]/(double)(NTIMES-1);
342:
343: printf("%s%11.4f \n", label[j], 1.0E-06 * bytes[j]/mintime[j]);
344: }
345: /* printf(HLINE);*/
346: #if WITH_PTHREADS
347: PetscThreadCommBarrier(PETSC_COMM_WORLD);
348: #endif
349: /* --- Check Results --- */
350: checkSTREAMresults();
351: /* printf(HLINE);*/
352: PetscFinalize();
353: return 0;
354: }
356: # define M 20
358: int checktick() {
359: int i, minDelta, Delta;
360: double t1, t2, timesfound[M];
362: /* Collect a sequence of M unique time values from the system. */
364: for (i = 0; i < M; i++) {
365: t1 = mysecond();
366: while( ((t2=mysecond()) - t1) < 1.0E-6 )
367: ;
368: timesfound[i] = t1 = t2;
369: }
371: /*
372: * Determine the minimum difference between these M values.
373: * This result will be our estimate (in microseconds) for the
374: * clock granularity.
375: */
377: minDelta = 1000000;
378: for (i = 1; i < M; i++) {
379: Delta = (int)( 1.0E6 * (timesfound[i]-timesfound[i-1]));
380: minDelta = MIN(minDelta, MAX(Delta,0));
381: }
383: return(minDelta);
384: }
386: /* A gettimeofday routine to give access to the wall
387: clock timer on most UNIX-like systems. */
389: #include <sys/time.h>
391: double mysecond() {
392: struct timeval tp;
393: struct timezone tzp;
394: int i=0;
395:
396: i = gettimeofday(&tp,&tzp);
397: return ( (double) tp.tv_sec + (double) tp.tv_usec * 1.e-6 );
398: }
400: void checkSTREAMresults () {
401: double aj,bj,cj,scalar;
402: double asum,bsum,csum;
403: double epsilon;
404: int j,k;
406: /* reproduce initialization */
407: aj = 1.0;
408: bj = 2.0;
409: cj = 0.0;
410: /* a[] is modified during timing check */
411: aj = 2.0E0 * aj;
412: /* now execute timing loop */
413: scalar = 3.0;
414: for (k=0; k<NTIMES; k++) {
415: cj = aj;
416: bj = scalar*cj;
417: cj = aj+bj;
418: aj = bj+scalar*cj;
419: }
420: aj = aj * (double) (N);
421: bj = bj * (double) (N);
422: cj = cj * (double) (N);
423:
424: asum = 0.0;
425: bsum = 0.0;
426: csum = 0.0;
427: for (j=0; j<N; j++) {
428: asum += a[j];
429: bsum += b[j];
430: csum += c[j];
431: }
432: #ifdef VERBOSE
433: printf ("Results Comparison: \n");
434: printf (" Expected : %f %f %f \n",aj,bj,cj);
435: printf (" Observed : %f %f %f \n",asum,bsum,csum);
436: #endif
438: #ifndef abs
439: #define abs(a) ((a) >= 0 ? (a) : -(a))
440: #endif
441: epsilon = 1.e-8;
443: if (abs(aj-asum)/asum > epsilon) {
444: printf ("Failed Validation on array a[]\n");
445: printf (" Expected : %f \n",aj);
446: printf (" Observed : %f \n",asum);
447: }
448: if (abs(bj-bsum)/bsum > epsilon) {
449: printf ("Failed Validation on array b[]\n");
450: printf (" Expected : %f \n",bj);
451: printf (" Observed : %f \n",bsum);
452: }
453: if (abs(cj-csum)/csum > epsilon) {
454: printf ("Failed Validation on array c[]\n");
455: printf (" Expected : %f \n",cj);
456: printf (" Observed : %f \n",csum);
457: }
458: else {
459: ;/* printf ("Solution Validates\n"); */
460: }
461: }