home *** CD-ROM | disk | FTP | other *** search
/ Microsoft Programmer's Library 1.3 / Microsoft-Programers-Library-v1.3.iso / sampcode / fortran / dwhet.for < prev    next >
Encoding:
Text File  |  1988-08-11  |  12.2 KB  |  375 lines
  1. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  2. C
  3. C       "WHETSTONE INSTRUCTIONS PER SECONDS" MEASURE OF FORTRAN
  4. C       AND CPU PERFORMANCE.
  5. C
  6. C         References on Whetstones:    Computer Journal Feb 76
  7. C                                      pg 43-49 vol 19 no 1.
  8. C                                      Curnow and Wichman.
  9. C
  10. C                                      and Timing Studies using a
  11. C                                      synthetic Whetstone Benchmark
  12. C                                      S. Harbaugh & J. Forakris
  13. C
  14. C         References on FORTRAN Benchmarks:
  15. C
  16. C                                   -  Computer Languages, Jan 1986
  17. C                                   -  EDN, Oct 3, 1985, Array Processors
  18. C                                           for PCs
  19. C                                   -  Byte, Feb 1984.
  20. C                                                         
  21. C       03/03/87
  22. C          Seeing that Microsoft distributed this without
  23. C          shipping the commented version, the timing loop 
  24. C          was reworked to eliminate all do loops and to put
  25. C          in code to print the variation in the measurment to
  26. C          make it more cook-book.  The 3 loop method described
  27. C          below was eliminated because it caused confusion.
  28. C          The printout was grouped and placed at the end of the test
  29. C          so that the outputs could be checked.  
  30. C          although it is ugly code, it checks with the Ada version
  31. C          and original article.          
  32. C          Because the Whetstones are printed as a reciprical,
  33. C          you can not average Whetstones to reduce time errors,
  34. C          you must run it multiple times and accumulate time.
  35. C          (AKT)
  36. C          
  37. C
  38. C       01/01/87
  39. C          fixed second subroutine to return seconds, not centi
  40. C          seconds and used double precision variables. (AKT)
  41. C   
  42. C       12/15/86
  43. C          Modified by Microsoft, removed reading in loop 
  44. C          option, added timer routine, removed meta-commands
  45. C          on large model.  Changed default looping from 100 to 10.
  46. C
  47. C       9/24/84
  48. C 
  49. C          ADDED CODE TO THESE SO THAT IT HAS VARIABLE LOOPING
  50. C
  51. C          from DEC but DONE BY OUTSIDE CONTRACTOR, OLD STYLE CODING
  52. C          not representative of DEC coding 
  53. C   
  54. C          A. TETEWSKY, c/o 
  55. C          555 TECH SQ MS 92 
  56. C          CAMBRIDGE MASS 02139           617/258-1287
  57. C 
  58. C          benchmarking notes:   1)    insure that timer has
  59. C                                      sufficient resolution and
  60. C                                      uses elapsed CPU time
  61. C
  62. C                                2)    to be useful for mainframe
  63. C                                      comparisons, measure
  64. C                                      INTEGER*4 time and large
  65. C                                      memory model or quote
  66. C                                      both large and small model
  67. C                                      times.  It may be necessary
  68. C                                      to make the arrays in this
  69. C                                      program large enough to span
  70. C                                      a 64K byte boundary because
  71. C                                      many micro compilers will
  72. C                                      generate small model code
  73. C                                      for small arrays even with 
  74. C                                      large models.
  75. C
  76. C                                 3)   Make sure that it loops
  77. C                                      long enough to gain
  78. C                                      stability, i.e. third-second
  79. C                                      loop = first loop time.
  80. C
  81. C         research notes,
  82. C         structure and definition:
  83. C         I received this code as a black box and based on some
  84. C         study, discovered the following background.
  85. C
  86. C            n1-n10 are loop counters for 10 tests, and tests
  87. C            n1,n5, and n10 are skipped.
  88. C            computed goto's are used to skip over tests that
  89. C            are not wanted.
  90. C     
  91. C
  92. C            n1-n10 scale with I.   When I is set to 10,
  93. C            kilo whets per second = 1000/ (time for doing n1-n10),
  94. C            the definition found in the literature.
  95. C
  96. C            If I were 100, the scale factor would be 10,000.
  97. C            which explains the 10,000 discovered in this code because
  98. C            it was shipped with IMUCH wired to 100.
  99. C
  100. C            the original DEC version uses a do-loop, 
  101. C                  imuch=100
  102. C                  do 200 loop=1,3
  103. C                       i = loop*imuch
  104. C                       n1-n10 scales by I
  105. C                       ... whetstones here ...
  106. C              200 continue
  107. C 
  108. C            and it took me a while to figure out why it worked.
  109. C
  110. C            This code loops three times 
  111. C                 TIMES(1) is time for 1*I  whets
  112. C                 TIMES(2) is time for 2*I
  113. C                 TIMES(3) is time for 3*I
  114. C            and TIMES(3)-TIMES(2) =  time for 1*I. 
  115. C            As long as TIMES(3)-TIMES(2) =  TIMES(1) to
  116. C            4 digits, then the cycle counter is sufficiently
  117. C            large enough for a given clock resolution.
  118. C
  119. C            By scaling whets * (IMUCH/10), you can alter IMUCH.
  120. C            The default definition is IMUCH = 10, hence the factor
  121. C            is unity.  IMUCH should be a factor of 10.
  122. C
  123. C
  124. C            Problems I have found:
  125. C            -  the SECONDS function is a single precision number
  126. C               and as CPUs get faster, you need to loop longer
  127. C               so that significant digits are not dropped.
  128. C
  129. C
  130. C       WHETS.FOR       09/27/77     TDR
  131. C       ...WHICH IS AN IMPROVED VERSION OF:
  132. C       WHET2A.FTN      01/22/75     RBG
  133. C
  134.         DOUBLE PRECISION X1,X2,X3,X4,X,Y,Z,T,T1,T2,E1
  135.         INTEGER   J,K,L,I, N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,N11,ISAVE
  136.         COMMON    T,T1,T2,E1(4),J,K,L
  137. C
  138. C
  139.         REAL*8         BEGTIM, ENDTIM, DIFTIM
  140.         REAL*8         DT,  WHETS, WS, ERR, WERR, PERR
  141.         INTEGER*4      IO1, IO1L, KREP, MKREP
  142. C
  143.         REAL*4         SECNDS
  144.         EXTERNAL       SECNDS
  145.         
  146. C
  147. C****************************************************************
  148. C                         
  149. C
  150.         MKREP  =    2
  151.         KREP   =    0
  152.         WRITE(*,*) ' Suggest inner > 10, outer > 1 '
  153.         WRITE(*,*) ' ENTER the number of inner/outer loops  ' 
  154.         READ(*,*)  I  ,IO1 
  155.  7020   CONTINUE
  156.         WRITE(*,*) ' Starting ',IO1,' loops, inner loop = ',I  
  157. C       ***** BEGININNING OF TIMED INTERVAL *****
  158.         IO1L   = 0
  159.         BEGTIM = DBLE(SECNDS(0.0E+00) )
  160.  7010   CONTINUE
  161. C
  162. C       ... the Whetstone code here ...
  163.         ISAVE=I   
  164.         T=0.499975D00
  165.         T1=0.50025D00
  166.         T2=2.0D00
  167.         N1=0
  168.         N2=12*I
  169.         N3=14*I
  170.         N4=345*I
  171.         N5=0
  172.         N6=210*I
  173.         N7=32*I
  174.         N8=899*I
  175.         N9=616*I
  176.         N10=0
  177.         N11=93*I
  178.         N12=0
  179.         X1=1.0D0
  180.         X2=-1.0D0
  181.         X3=-1.0D0                                     
  182.         X4=-1.0D0
  183.         IF(N1)19,19,11
  184.  11     DO 18 I=1,N1,1            
  185.         X1=(X1+X2+X3-X4)*T
  186.         X2=(X1+X2-X3+X4)*T
  187.         X4=(-X1+X2+X3+X4)*T
  188.         X3=(X1-X2+X3+X4)*T
  189.  18     CONTINUE
  190.  19     CONTINUE
  191.         E1(1)=1.0D0
  192.         E1(2)=-1.0D0
  193.         E1(3)=-1.0D0
  194.         E1(4)=-1.0D0
  195.         IF(N2)29,29,21
  196.  21     DO 28 I=1,N2,1
  197.         E1(1)=(E1(1)+E1(2)+E1(3)-E1(4))*T
  198.         E1(2)=(E1(1)+E1(2)-E1(3)+E1(4))*T
  199.         E1(3)=(E1(1)-E1(2)+E1(3)+E1(4))*T
  200.         E1(4)=(-E1(1)+E1(2)+E1(3)+E1(4))*T
  201.  28     CONTINUE
  202.  29     CONTINUE
  203.         IF(N3)39,39,31
  204.  31     DO 38 I=1,N3,1
  205.  38     CALL PA(E1)
  206.  39     CONTINUE
  207.         J=1
  208.         IF(N4)49,49,41
  209.  41     DO 48 I=1,N4,1
  210.         IF(J-1)43,42,43
  211.  42     J=2
  212.         GOTO44
  213.  43     J=3
  214.  44     IF(J-2)46,46,45
  215.  45     J=0
  216.         GOTO47
  217.  46     J=1
  218.  47     IF(J-1)411,412,412
  219.  411    J=1
  220.         GOTO48
  221.  412    J=0
  222.  48     CONTINUE
  223.  49     CONTINUE
  224.         J=1
  225.         K=2
  226.         L=3
  227.         IF(N6)69,69,61
  228.  61     DO 68 I=1,N6,1
  229.         J=J*(K-J)*(L-K)
  230.         K=L*K-(L-J)*K
  231.         L=(L-K)*(K+J)
  232.         E1(L-1)=J+K+L
  233.         E1(K-1)=J*K*L
  234.  68     CONTINUE
  235.  69     CONTINUE
  236.         X=0.5D0
  237.         Y=0.5D0
  238.         IF(N7)79,79,71
  239.  71     DO 78 I=1,N7,1
  240.         X=T*DATAN(T2*DSIN(X)*DCOS(X)/(DCOS(X+Y)+DCOS(X-Y)-1.0D0))
  241.         Y=T*DATAN(T2*DSIN(Y)*DCOS(Y)/(DCOS(X+Y)+DCOS(X-Y)-1.0D0))
  242.  78     CONTINUE
  243.  79     CONTINUE
  244.         X=1.0D0
  245.         Y=1.0D0
  246.         Z=1.0D0
  247.         IF(N8)89,89,81
  248.  81     DO 88 I=1,N8,1
  249.  88     CALL P3(X,Y,Z)
  250.  89     CONTINUE
  251.         J=1
  252.         K=2
  253.         L=3
  254.         E1(1)=1.0D0
  255.         E1(2)=2.0D0
  256.         E1(3)=3.0D0
  257.         IF(N9)99,99,91
  258.  91     DO 98 I=1,N9,1
  259.  98     CALL P0
  260.  99     CONTINUE
  261.         J=2
  262.         K=3
  263.         IF(N10)109,109,101
  264.  101    DO 108 I=1,N10,1
  265.         J=J+K
  266.         K=J+K
  267.         J=J-K
  268.         K=K-J-J
  269.  108    CONTINUE
  270.  109    CONTINUE
  271.         X=0.75D0
  272.         IF(N11)119,119,111
  273.  111    DO 118 I=1,N11,1
  274.  118    X=DSQRT(DEXP(DLOG(X)/T1))
  275.  119    CONTINUE
  276.         I = ISAVE
  277. C       ... the whetstone ends here
  278. C
  279. C         ... loop counter instead of do loop ...
  280.           IO1L = IO1L + 1
  281.           IF( IO1L .LT. IO1) GOTO 7010
  282. C       ******* END of TIME INTERVALED ***********
  283. C
  284.         ENDTIM = DBLE(SECNDS(0.0E+00))
  285.         DIFTIM = ENDTIM - BEGTIM
  286. C       whets  = 1000/(TIME FOR 10 inner ITERATIONS OF PROGRAM LOOP)
  287. C       or 100 for every 1 inner count
  288.         WHETS = (100.0D+00* DBLE( FLOAT(IO1*I  ))/DIFTIM)
  289.         WRITE(*,*) ' START TIME = ',BEGTIM
  290.         WRITE(*,*) ' END   TIME = ',ENDTIM
  291.         WRITE(*,*) ' DIF   TIME = ',DIFTIM
  292. C
  293.         WRITE (*,201) WHETS
  294.   201   FORMAT(' SPEED IS: ',1PE10.3,' THOUSAND WHETSTONE',
  295.      2     ' DOUBLE PRECISION INSTRUCTIONS PER SECOND')
  296.         CALL POUT(N1,N1,N1,X1,X2,X3,X4)
  297.         CALL POUT(N2,N3,N2,E1(1),E1(2),E1(3),E1(4))
  298.         CALL POUT(N3,N2,N2,E1(1),E1(2),E1(3),E1(4))
  299.         CALL POUT(N4,J,J,X1,X2,X3,X4)
  300.         CALL POUT(N6,J,K,E1(1),E1(2),E1(3),E1(4))
  301.         CALL POUT(N7,J,K,X,X,Y,Y)
  302.         CALL POUT(N8,J,K,X,Y,Z,Z)
  303.         CALL POUT(N9,J,K,E1(1),E1(2),E1(3),E1(4))
  304.         CALL POUT(N10,J,K,X1,X2,X3,X4)
  305.         CALL POUT(N11,J,K,X,X,X,X)
  306. C
  307. C       ... repeat but double (MULTIPLY UP) inner count ...
  308.         KREP = KREP + 1
  309.         IF( KREP .LT. MKREP) THEN
  310.             DT     = DIFTIM
  311.             WT     = WHETS
  312.             I=I*MKREP
  313.             GOTO 7020
  314.         ENDIF
  315. C
  316. C       ... compute sensitivity  
  317. C
  318.         ERR =  DIFTIM - (DT*DBLE(FLOAT(MKREP)))    
  319.         WERR=  WT-WHETS          
  320.         PERR=  WERR*100.0D+00/WHETS
  321.         WRITE(*,*) ' Time ERR = ',ERR, ' seconds '
  322.         WRITE(*,*) ' Whet ERR = ',WERR,' kwhets/sec '
  323.         WRITE(*,*) ' %    ERR = ',PERR,' % whet error '
  324.         IF( DIFTIM .LT. 10.0D+00) THEN
  325.          WRITE(*,*)
  326.      1   ' TIME is less than 10 seconds, suggest larger inner loop '
  327.         ENDIF
  328. C
  329.         STOP
  330.         END
  331.         SUBROUTINE PA(E)
  332. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  333.         DOUBLE PRECISION T,T1,T2,E
  334.         COMMON T,T1,T2
  335.         DIMENSION E(4)
  336.         J=0
  337.  1      E(1)=(E(1)+E(2)+E(3)-E(4))*T
  338.         E(2)=(E(1)+E(2)-E(3)+E(4))*T
  339.         E(3)=(E(1)-E(2)+E(3)+E(4))*T
  340.         E(4)=(-E(1)+E(2)+E(3)+E(4))/T2
  341.         J=J+1
  342.         IF(J-6)1,2,2
  343.  2      CONTINUE
  344.         RETURN
  345.         END
  346.         SUBROUTINE P0
  347. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  348.         DOUBLE PRECISION T,T1,T2,E1
  349.         COMMON T,T1,T2,E1(4),J,K,L
  350.         E1(J)=E1(K)
  351.         E1(K)=E1(L)
  352.         E1(L)=E1(J)
  353.         RETURN
  354.         END
  355.         SUBROUTINE P3(X,Y,Z)
  356. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  357.         DOUBLE PRECISION T,T1,T2,X1,Y1,X,Y,Z
  358.         COMMON T,T1,T2
  359.         X1=X
  360.         Y1=Y
  361.         X1=T*(X1+Y1)
  362.         Y1=T*(X1+Y1)
  363.         Z=(X1+Y1)/T2
  364.         RETURN
  365.         END
  366.         SUBROUTINE POUT(N,J,K,X1,X2,X3,X4)
  367. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  368. C
  369. C
  370.         DOUBLE PRECISION X1,X2,X3,X4
  371.         WRITE(*,1)N,J,K,X1,X2,X3,X4
  372.  1      FORMAT('  ',3(I7,1X),4(1PD12.4,1X))
  373.         RETURN
  374.         END
  375.