benchmarks/whetstone/whetstone.c

0001 /*
0002  * C Converted Whetstone Double Precision Benchmark
0003  *      Version 1.2 22 March 1998
0004  *
0005  *  (c) Copyright 1998 Painter Engineering, Inc.
0006  *      All Rights Reserved.
0007  *
0008  *      Permission is granted to use, duplicate, and
0009  *      publish this text and program as long as it
0010  *      includes this entire comment block and limited
0011  *      rights reference.
0012  *
0013  * Converted by Rich Painter, Painter Engineering, Inc. based on the
0014  * www.netlib.org benchmark/whetstoned version obtained 16 March 1998.
0015  *
0016  * A novel approach was used here to keep the look and feel of the
0017  * FORTRAN version.  Altering the FORTRAN-based array indices,
0018  * starting at element 1, to start at element 0 for C, would require
0019  * numerous changes, including decrementing the variable indices by 1.
0020  * Instead, the array E1[] was declared 1 element larger in C.  This
0021  * allows the FORTRAN index range to function without any literal or
0022  * variable indices changes.  The array element E1[0] is simply never
0023  * used and does not alter the benchmark results.
0024  *
0025  * The major FORTRAN comment blocks were retained to minimize
0026  * differences between versions.  Modules N5 and N12, like in the
0027  * FORTRAN version, have been eliminated here.
0028  *
0029  * An optional command-line argument has been provided [-c] to
0030  * offer continuous repetition of the entire benchmark.
0031  * An optional argument for setting an alternate LOOP count is also
0032  * provided.  Define PRINTOUT to cause the POUT() function to print
0033  * outputs at various stages.  Final timing measurements should be
0034  * made with the PRINTOUT undefined.
0035  *
0036  * Questions and comments may be directed to the author at
0037  *          r.painter@ieee.org
0038  */
0039 /*
0040 C**********************************************************************
0041 C     Benchmark #2 -- Double  Precision Whetstone (A001)
0042 C
0043 C     o This is a REAL*8 version of
0044 C   the Whetstone benchmark program.
0045 C
0046 C     o DO-loop semantics are ANSI-66 compatible.
0047 C
0048 C     o Final measurements are to be made with all
0049 C   WRITE statements and FORMAT sttements removed.
0050 C
0051 C**********************************************************************
0052 */
0053
0054 /* standard C library headers required */
0055 #include <stdlib.h>
0056 #include <stdio.h>
0057 #include <string.h>
0058 #include <math.h>
0059 #ifdef __rtems__
0060 #include <rtems/test-printer.h>
0061 #define fprintf(f, ...) rtems_printf(&rtems_test_printer, __VA_ARGS__)
0062 #endif /* __rtems__ */
0063
0064 /* the following is optional depending on the timing function used */
0065 #include <sys/time.h>
0066
0067 /* map the FORTRAN math functions, etc. to the C versions */
0068 #define DSIN    sin
0069 #define DCOS    cos
0070 #define DATAN   atan
0071 #define DLOG    log
0072 #define DEXP    exp
0073 #define DSQRT   sqrt
0074 #define IF      if
0075
0076 /* function prototypes */
0077 void POUT(long N, long J, long K, double X1, double X2, double X3, double X4);
0078 void PA(double E[]);
0079 void P0(void);
0080 void P3(double X, double Y, double *Z);
0081 #define USAGE   "usage: whetdc [-c] [loops]\n"
0082
0083 /*
0084     COMMON T,T1,T2,E1(4),J,K,L
0085 */
0086 double T,T1,T2,E1[5];
0087 int J,K,L;
0088
0089 static double
0090 Time(void)
0091 {
0092     struct timeval tv;
0093
0094     gettimeofday(&tv, NULL);
0095     return (double)tv.tv_sec + (double)tv.tv_usec * 1e-6;
0096 }
0097
0098 int
0099 main(int argc, char *argv[])
0100 {
0101     /* used in the FORTRAN version */
0102     long I;
0103     long N1, N2, N3, N4, N6, N7, N8, N9, N10, N11;
0104     double X1,X2,X3,X4,X,Y,Z;
0105     long LOOP;
0106     int II, JJ;
0107
0108     /* added for this version */
0109     long loopstart;
0110     double startsec, finisec;
0111     double KIPS;
0112     int continuous;
0113
0114     loopstart = 1000;       /* see the note about LOOP below */
0115     continuous = 0;
0116
0117     II = 1;     /* start at the first arg (temp use of II here) */
0118     while (II < argc) {
0119         if (strncmp(argv[II], "-c", 2) == 0 || argv[II][0] == 'c') {
0120             continuous = 1;
0121         } else if (atol(argv[II]) > 0) {
0122             loopstart = atol(argv[II]);
0123         } else {
0124             fprintf(stderr, USAGE);
0125             return(1);
0126         }
0127         II++;
0128     }
0129
0130 LCONT:
0131 /*
0132 C
0133 C   Start benchmark timing at this point.
0134 C
0135 */
0136     startsec = Time();
0137
0138 /*
0139 C
0140 C   The actual benchmark starts here.
0141 C
0142 */
0143     T  = .499975;
0144     T1 = 0.50025;
0145     T2 = 2.0;
0146 /*
0147 C
0148 C   With loopcount LOOP=10, one million Whetstone instructions
0149 C   will be executed in EACH MAJOR LOOP..A MAJOR LOOP IS EXECUTED
0150 C   'II' TIMES TO INCREASE WALL-CLOCK TIMING ACCURACY.
0151 C
0152     LOOP = 1000;
0153 */
0154     LOOP = loopstart;
0155     II   = 1;
0156
0157     JJ = 1;
0158
0159 IILOOP:
0160     N1  = 0;
0161     N2  = 12 * LOOP;
0162     N3  = 14 * LOOP;
0163     N4  = 345 * LOOP;
0164     N6  = 210 * LOOP;
0165     N7  = 32 * LOOP;
0166     N8  = 899 * LOOP;
0167     N9  = 616 * LOOP;
0168     N10 = 0;
0169     N11 = 93 * LOOP;
0170 /*
0171 C
0172 C   Module 1: Simple identifiers
0173 C
0174 */
0175     X1  =  1.0;
0176     X2  = -1.0;
0177     X3  = -1.0;
0178     X4  = -1.0;
0179
0180     for (I = 1; I <= N1; I++) {
0181         X1 = (X1 + X2 + X3 - X4) * T;
0182         X2 = (X1 + X2 - X3 + X4) * T;
0183         X3 = (X1 - X2 + X3 + X4) * T;
0184         X4 = (-X1+ X2 + X3 + X4) * T;
0185     }
0186 #ifdef PRINTOUT
0187     IF (JJ==II)POUT(N1,N1,N1,X1,X2,X3,X4);
0188 #endif
0189
0190 /*
0191 C
0192 C   Module 2: Array elements
0193 C
0194 */
0195     E1[1] =  1.0;
0196     E1[2] = -1.0;
0197     E1[3] = -1.0;
0198     E1[4] = -1.0;
0199
0200     for (I = 1; I <= N2; I++) {
0201         E1[1] = ( E1[1] + E1[2] + E1[3] - E1[4]) * T;
0202         E1[2] = ( E1[1] + E1[2] - E1[3] + E1[4]) * T;
0203         E1[3] = ( E1[1] - E1[2] + E1[3] + E1[4]) * T;
0204         E1[4] = (-E1[1] + E1[2] + E1[3] + E1[4]) * T;
0205     }
0206
0207 #ifdef PRINTOUT
0208     IF (JJ==II)POUT(N2,N3,N2,E1[1],E1[2],E1[3],E1[4]);
0209 #endif
0210
0211 /*
0212 C
0213 C   Module 3: Array as parameter
0214 C
0215 */
0216     for (I = 1; I <= N3; I++)
0217         PA(E1);
0218
0219 #ifdef PRINTOUT
0220     IF (JJ==II)POUT(N3,N2,N2,E1[1],E1[2],E1[3],E1[4]);
0221 #endif
0222
0223 /*
0224 C
0225 C   Module 4: Conditional jumps
0226 C
0227 */
0228     J = 1;
0229     for (I = 1; I <= N4; I++) {
0230         if (J == 1)
0231             J = 2;
0232         else
0233             J = 3;
0234
0235         if (J > 2)
0236             J = 0;
0237         else
0238             J = 1;
0239
0240         if (J < 1)
0241             J = 1;
0242         else
0243             J = 0;
0244     }
0245
0246 #ifdef PRINTOUT
0247     IF (JJ==II)POUT(N4,J,J,X1,X2,X3,X4);
0248 #endif
0249
0250 /*
0251 C
0252 C   Module 5: Omitted
0253 C   Module 6: Integer arithmetic
0254 C
0255 */
0256
0257     J = 1;
0258     K = 2;
0259     L = 3;
0260
0261     for (I = 1; I <= N6; I++) {
0262         J = J * (K-J) * (L-K);
0263         K = L * K - (L-J) * K;
0264         L = (L-K) * (K+J);
0265         E1[L-1] = J + K + L;
0266         E1[K-1] = J * K * L;
0267     }
0268
0269 #ifdef PRINTOUT
0270     IF (JJ==II)POUT(N6,J,K,E1[1],E1[2],E1[3],E1[4]);
0271 #endif
0272
0273 /*
0274 C
0275 C   Module 7: Trigonometric functions
0276 C
0277 */
0278     X = 0.5;
0279     Y = 0.5;
0280
0281     for (I = 1; I <= N7; I++) {
0282         X = T * DATAN(T2*DSIN(X)*DCOS(X)/(DCOS(X+Y)+DCOS(X-Y)-1.0));
0283         Y = T * DATAN(T2*DSIN(Y)*DCOS(Y)/(DCOS(X+Y)+DCOS(X-Y)-1.0));
0284     }
0285
0286 #ifdef PRINTOUT
0287     IF (JJ==II)POUT(N7,J,K,X,X,Y,Y);
0288 #endif
0289
0290 /*
0291 C
0292 C   Module 8: Procedure calls
0293 C
0294 */
0295     X = 1.0;
0296     Y = 1.0;
0297     Z = 1.0;
0298
0299     for (I = 1; I <= N8; I++)
0300         P3(X,Y,&Z);
0301
0302 #ifdef PRINTOUT
0303     IF (JJ==II)POUT(N8,J,K,X,Y,Z,Z);
0304 #endif
0305
0306 /*
0307 C
0308 C   Module 9: Array references
0309 C
0310 */
0311     J = 1;
0312     K = 2;
0313     L = 3;
0314     E1[1] = 1.0;
0315     E1[2] = 2.0;
0316     E1[3] = 3.0;
0317
0318     for (I = 1; I <= N9; I++)
0319         P0();
0320
0321 #ifdef PRINTOUT
0322     IF (JJ==II)POUT(N9,J,K,E1[1],E1[2],E1[3],E1[4]);
0323 #endif
0324
0325 /*
0326 C
0327 C   Module 10: Integer arithmetic
0328 C
0329 */
0330     J = 2;
0331     K = 3;
0332
0333     for (I = 1; I <= N10; I++) {
0334         J = J + K;
0335         K = J + K;
0336         J = K - J;
0337         K = K - J - J;
0338     }
0339
0340 #ifdef PRINTOUT
0341     IF (JJ==II)POUT(N10,J,K,X1,X2,X3,X4);
0342 #endif
0343
0344 /*
0345 C
0346 C   Module 11: Standard functions
0347 C
0348 */
0349     X = 0.75;
0350
0351     for (I = 1; I <= N11; I++)
0352         X = DSQRT(DEXP(DLOG(X)/T1));
0353
0354 #ifdef PRINTOUT
0355     IF (JJ==II)POUT(N11,J,K,X,X,X,X);
0356 #endif
0357
0358 /*
0359 C
0360 C      THIS IS THE END OF THE MAJOR LOOP.
0361 C
0362 */
0363     if (++JJ <= II)
0364         goto IILOOP;
0365
0366 /*
0367 C
0368 C      Stop benchmark timing at this point.
0369 C
0370 */
0371     finisec = Time();
0372
0373 /*
0374 C----------------------------------------------------------------
0375 C      Performance in Whetstone KIP's per second is given by
0376 C
0377 C   (100*LOOP*II)/TIME
0378 C
0379 C      where TIME is in seconds.
0380 C--------------------------------------------------------------------
0381 */
0382     printf("\n");
0383     if (finisec-startsec <= 0) {
0384         printf("Insufficient duration- Increase the LOOP count\n");
0385         return(1);
0386     }
0387
0388     printf("Loops: %ld, Iterations: %d, Duration: %f sec.\n",
0389             LOOP, II, finisec-startsec);
0390
0391     KIPS = (100.0*LOOP*II)/(finisec-startsec);
0392     if (KIPS >= 1000.0)
0393         printf("C Converted Double Precision Whetstones: %.1f MIPS\n", KIPS/1000.0);
0394     else
0395         printf("C Converted Double Precision Whetstones: %.1f KIPS\n", KIPS);
0396
0397     if (continuous)
0398         goto LCONT;
0399
0400     return(0);
0401 }
0402
0403 void
0404 PA(double E[])
0405 {
0406     J = 0;
0407
0408 L10:
0409     E[1] = ( E[1] + E[2] + E[3] - E[4]) * T;
0410     E[2] = ( E[1] + E[2] - E[3] + E[4]) * T;
0411     E[3] = ( E[1] - E[2] + E[3] + E[4]) * T;
0412     E[4] = (-E[1] + E[2] + E[3] + E[4]) / T2;
0413     J += 1;
0414
0415     if (J < 6)
0416         goto L10;
0417 }
0418
0419 void
0420 P0(void)
0421 {
0422     E1[J] = E1[K];
0423     E1[K] = E1[L];
0424     E1[L] = E1[J];
0425 }
0426
0427 void
0428 P3(double X, double Y, double *Z)
0429 {
0430     double X1, Y1;
0431
0432     X1 = X;
0433     Y1 = Y;
0434     X1 = T * (X1 + Y1);
0435     Y1 = T * (X1 + Y1);
0436     *Z  = (X1 + Y1) / T2;
0437 }
0438
0439 #ifdef PRINTOUT
0440 void
0441 POUT(long N, long J, long K, double X1, double X2, double X3, double X4)
0442 {
0443     printf("%7ld %7ld %7ld %12.4e %12.4e %12.4e %12.4e\n",
0444                         N, J, K, X1, X2, X3, X4);
0445 }
0446 #endif