1 /*
2  * C Converted Whetstone Double Precision Benchmark
3  *		Version 1.2	22 March 1998
4  *
5  *	(c) Copyright 1998 Painter Engineering, Inc.
6  *		All Rights Reserved.
7  *
8  *		Permission is granted to use, duplicate, and
9  *		publish this text and program as long as it
10  *		includes this entire comment block and limited
11  *		rights reference.
12  *
13  * Converted by Rich Painter, Painter Engineering, Inc. based on the
14  * www.netlib.org benchmark/whetstoned version obtained 16 March 1998.
15  *
16  * A novel approach was used here to keep the look and feel of the
17  * FORTRAN version.  Altering the FORTRAN-based array indices,
18  * starting at element 1, to start at element 0 for C, would require
19  * numerous changes, including decrementing the variable indices by 1.
20  * Instead, the array E1[] was declared 1 element larger in C.  This
21  * allows the FORTRAN index range to function without any literal or
22  * variable indices changes.  The array element E1[0] is simply never
23  * used and does not alter the benchmark results.
24  *
25  * The major FORTRAN comment blocks were retained to minimize
26  * differences between versions.  Modules N5 and N12, like in the
27  * FORTRAN version, have been eliminated here.
28  *
29  * An optional command-line argument has been provided [-c] to
30  * offer continuous repetition of the entire benchmark.
31  * An optional argument for setting an alternate LOOP count is also
32  * provided.  Define PRINTOUT to cause the POUT() function to print
33  * outputs at various stages.  Final timing measurements should be
34  * made with the PRINTOUT undefined.
35  *
36  * Questions and comments may be directed to the author at
37  *			r.painter@ieee.org
38  */
39 /*
40 C**********************************************************************
41 C     Benchmark #2 -- Double  Precision Whetstone (A001)
42 C
43 C     o	This is a REAL*8 version of
44 C	the Whetstone benchmark program.
45 C
46 C     o	DO-loop semantics are ANSI-66 compatible.
47 C
48 C     o	Final measurements are to be made with all
49 C	WRITE statements and FORMAT sttements removed.
50 C
51 C**********************************************************************
52 */
53 
54 /* standard C library headers required */
55 #include <stdlib.h>
56 #include <stdio.h>
57 #include <string.h>
58 #include <math.h>
59 
60 /* the following is optional depending on the timing function used */
61 #include <time.h>
62 
63 /* map the FORTRAN math functions, etc. to the C versions */
64 #define DSIN	sin
65 #define DCOS	cos
66 #define DATAN	atan
67 #define DLOG	log
68 #define DEXP	exp
69 #define DSQRT	sqrt
70 #define IF		if
71 
72 /* function prototypes */
73 void POUT(long N, long J, long K, double X1, double X2, double X3, double X4);
74 void PA(double E[]);
75 void P0(void);
76 void P3(double X, double Y, double *Z);
77 #define USAGE	"usage: whetdc [-c] [loops]\n"
78 
79 /*
80 	COMMON T,T1,T2,E1(4),J,K,L
81 */
82 double T,T1,T2,E1[5];
83 int J,K,L;
84 
85 int
main(int argc,char * argv[])86 main(int argc, char *argv[])
87 {
88 	/* used in the FORTRAN version */
89 	long I;
90 	long N1, N2, N3, N4, N6, N7, N8, N9, N10, N11;
91 	double X1,X2,X3,X4,X,Y,Z;
92 	long LOOP;
93 	int II, JJ;
94 
95 	/* added for this version */
96 	long loopstart;
97 	long startsec, finisec;
98 	float KIPS;
99 	int continuous;
100 
101 	loopstart = 1000;		/* see the note about LOOP below */
102 	continuous = 0;
103 
104 	II = 1;		/* start at the first arg (temp use of II here) */
105 	while (II < argc) {
106 		if (strncmp(argv[II], "-c", 2) == 0 || argv[II][0] == 'c') {
107 			continuous = 1;
108 		} else if (atol(argv[II]) > 0) {
109 			loopstart = atol(argv[II]);
110 		} else {
111 			fprintf(stderr, USAGE);
112 			return(1);
113 		}
114 		II++;
115 	}
116 
117 LCONT:
118 /*
119 C
120 C	Start benchmark timing at this point.
121 C
122 */
123 	startsec = time(0);
124 
125 /*
126 C
127 C	The actual benchmark starts here.
128 C
129 */
130 	T  = .499975;
131 	T1 = 0.50025;
132 	T2 = 2.0;
133 /*
134 C
135 C	With loopcount LOOP=10, one million Whetstone instructions
136 C	will be executed in EACH MAJOR LOOP..A MAJOR LOOP IS EXECUTED
137 C	'II' TIMES TO INCREASE WALL-CLOCK TIMING ACCURACY.
138 C
139 	LOOP = 1000;
140 */
141 	LOOP = loopstart;
142 	II   = 1;
143 
144 	JJ = 1;
145 
146 IILOOP:
147 	N1  = 0;
148 	N2  = 12 * LOOP;
149 	N3  = 14 * LOOP;
150 	N4  = 345 * LOOP;
151 	N6  = 210 * LOOP;
152 	N7  = 32 * LOOP;
153 	N8  = 899 * LOOP;
154 	N9  = 616 * LOOP;
155 	N10 = 0;
156 	N11 = 93 * LOOP;
157 /*
158 C
159 C	Module 1: Simple identifiers
160 C
161 */
162 	X1  =  1.0;
163 	X2  = -1.0;
164 	X3  = -1.0;
165 	X4  = -1.0;
166 
167 	for (I = 1; I <= N1; I++) {
168 	    X1 = (X1 + X2 + X3 - X4) * T;
169 	    X2 = (X1 + X2 - X3 + X4) * T;
170 	    X3 = (X1 - X2 + X3 + X4) * T;
171 	    X4 = (-X1+ X2 + X3 + X4) * T;
172 	}
173 #ifdef PRINTOUT
174 	IF (JJ==II)POUT(N1,N1,N1,X1,X2,X3,X4);
175 #endif
176 
177 /*
178 C
179 C	Module 2: Array elements
180 C
181 */
182 	E1[1] =  1.0;
183 	E1[2] = -1.0;
184 	E1[3] = -1.0;
185 	E1[4] = -1.0;
186 
187 	for (I = 1; I <= N2; I++) {
188 	    E1[1] = ( E1[1] + E1[2] + E1[3] - E1[4]) * T;
189 	    E1[2] = ( E1[1] + E1[2] - E1[3] + E1[4]) * T;
190 	    E1[3] = ( E1[1] - E1[2] + E1[3] + E1[4]) * T;
191 	    E1[4] = (-E1[1] + E1[2] + E1[3] + E1[4]) * T;
192 	}
193 
194 #ifdef PRINTOUT
195 	IF (JJ==II)POUT(N2,N3,N2,E1[1],E1[2],E1[3],E1[4]);
196 #endif
197 
198 /*
199 C
200 C	Module 3: Array as parameter
201 C
202 */
203 	for (I = 1; I <= N3; I++)
204 		PA(E1);
205 
206 #ifdef PRINTOUT
207 	IF (JJ==II)POUT(N3,N2,N2,E1[1],E1[2],E1[3],E1[4]);
208 #endif
209 
210 /*
211 C
212 C	Module 4: Conditional jumps
213 C
214 */
215 	J = 1;
216 	for (I = 1; I <= N4; I++) {
217 		if (J == 1)
218 			J = 2;
219 		else
220 			J = 3;
221 
222 		if (J > 2)
223 			J = 0;
224 		else
225 			J = 1;
226 
227 		if (J < 1)
228 			J = 1;
229 		else
230 			J = 0;
231 	}
232 
233 #ifdef PRINTOUT
234 	IF (JJ==II)POUT(N4,J,J,X1,X2,X3,X4);
235 #endif
236 
237 /*
238 C
239 C	Module 5: Omitted
240 C 	Module 6: Integer arithmetic
241 C
242 */
243 
244 	J = 1;
245 	K = 2;
246 	L = 3;
247 
248 	for (I = 1; I <= N6; I++) {
249 	    J = J * (K-J) * (L-K);
250 	    K = L * K - (L-J) * K;
251 	    L = (L-K) * (K+J);
252 	    E1[L-1] = J + K + L;
253 	    E1[K-1] = J * K * L;
254 	}
255 
256 #ifdef PRINTOUT
257 	IF (JJ==II)POUT(N6,J,K,E1[1],E1[2],E1[3],E1[4]);
258 #endif
259 
260 /*
261 C
262 C	Module 7: Trigonometric functions
263 C
264 */
265 	X = 0.5;
266 	Y = 0.5;
267 
268 	for (I = 1; I <= N7; I++) {
269 		X = T * DATAN(T2*DSIN(X)*DCOS(X)/(DCOS(X+Y)+DCOS(X-Y)-1.0));
270 		Y = T * DATAN(T2*DSIN(Y)*DCOS(Y)/(DCOS(X+Y)+DCOS(X-Y)-1.0));
271 	}
272 
273 #ifdef PRINTOUT
274 	IF (JJ==II)POUT(N7,J,K,X,X,Y,Y);
275 #endif
276 
277 /*
278 C
279 C	Module 8: Procedure calls
280 C
281 */
282 	X = 1.0;
283 	Y = 1.0;
284 	Z = 1.0;
285 
286 	for (I = 1; I <= N8; I++)
287 		P3(X,Y,&Z);
288 
289 #ifdef PRINTOUT
290 	IF (JJ==II)POUT(N8,J,K,X,Y,Z,Z);
291 #endif
292 
293 /*
294 C
295 C	Module 9: Array references
296 C
297 */
298 	J = 1;
299 	K = 2;
300 	L = 3;
301 	E1[1] = 1.0;
302 	E1[2] = 2.0;
303 	E1[3] = 3.0;
304 
305 	for (I = 1; I <= N9; I++)
306 		P0();
307 
308 #ifdef PRINTOUT
309 	IF (JJ==II)POUT(N9,J,K,E1[1],E1[2],E1[3],E1[4]);
310 #endif
311 
312 /*
313 C
314 C	Module 10: Integer arithmetic
315 C
316 */
317 	J = 2;
318 	K = 3;
319 
320 	for (I = 1; I <= N10; I++) {
321 	    J = J + K;
322 	    K = J + K;
323 	    J = K - J;
324 	    K = K - J - J;
325 	}
326 
327 #ifdef PRINTOUT
328 	IF (JJ==II)POUT(N10,J,K,X1,X2,X3,X4);
329 #endif
330 
331 /*
332 C
333 C	Module 11: Standard functions
334 C
335 */
336 	X = 0.75;
337 
338 	for (I = 1; I <= N11; I++)
339 		X = DSQRT(DEXP(DLOG(X)/T1));
340 
341 #ifdef PRINTOUT
342 	IF (JJ==II)POUT(N11,J,K,X,X,X,X);
343 #endif
344 
345 /*
346 C
347 C      THIS IS THE END OF THE MAJOR LOOP.
348 C
349 */
350 	if (++JJ <= II)
351 		goto IILOOP;
352 
353 /*
354 C
355 C      Stop benchmark timing at this point.
356 C
357 */
358 	finisec = time(0);
359 
360 /*
361 C----------------------------------------------------------------
362 C      Performance in Whetstone KIP's per second is given by
363 C
364 C	(100*LOOP*II)/TIME
365 C
366 C      where TIME is in seconds.
367 C--------------------------------------------------------------------
368 */
369 	printf("\n");
370 	if (finisec-startsec <= 0) {
371 		printf("Insufficient duration- Increase the LOOP count\n");
372 		return(1);
373 	}
374 
375 	printf("Loops: %ld, Iterations: %d, Duration: %ld sec.\n",
376 			LOOP, II, finisec-startsec);
377 
378 	KIPS = (100.0*LOOP*II)/(float)(finisec-startsec);
379 	if (KIPS >= 1000.0)
380 		printf("C Converted Double Precision Whetstones: %.1f MIPS\n", KIPS/1000.0);
381 	else
382 		printf("C Converted Double Precision Whetstones: %.1f KIPS\n", KIPS);
383 
384 	if (continuous)
385 		goto LCONT;
386 
387 	return(0);
388 }
389 
390 void
PA(double E[])391 PA(double E[])
392 {
393 	J = 0;
394 
395 L10:
396 	E[1] = ( E[1] + E[2] + E[3] - E[4]) * T;
397 	E[2] = ( E[1] + E[2] - E[3] + E[4]) * T;
398 	E[3] = ( E[1] - E[2] + E[3] + E[4]) * T;
399 	E[4] = (-E[1] + E[2] + E[3] + E[4]) / T2;
400 	J += 1;
401 
402 	if (J < 6)
403 		goto L10;
404 }
405 
406 void
P0(void)407 P0(void)
408 {
409 	E1[J] = E1[K];
410 	E1[K] = E1[L];
411 	E1[L] = E1[J];
412 }
413 
414 void
P3(double X,double Y,double * Z)415 P3(double X, double Y, double *Z)
416 {
417 	double X1, Y1;
418 
419 	X1 = X;
420 	Y1 = Y;
421 	X1 = T * (X1 + Y1);
422 	Y1 = T * (X1 + Y1);
423 	*Z  = (X1 + Y1) / T2;
424 }
425 
426 #ifdef PRINTOUT
427 void
POUT(long N,long J,long K,double X1,double X2,double X3,double X4)428 POUT(long N, long J, long K, double X1, double X2, double X3, double X4)
429 {
430 	printf("%7ld %7ld %7ld %12.4e %12.4e %12.4e %12.4e\n",
431 						N, J, K, X1, X2, X3, X4);
432 }
433 #endif
434