      PROGRAM PROB14
C
C     PROBLEM 14
C
C  REFERENCE:  PROBLEMS TO TEST PARALLEL AND VECTOR LANGUAGES
C              CSD-TR 516, COMPUTER SCIENCE, PURDUE UNIVERSITY
C              JOHN R. RICE, MAY 1, 1985
C
C              REVISED BY JOHN R. RICE AND J. JING, OCT. 1, 1990
C
C
C      *************************************************
C      *      Adapted for FORTRAN D benchmarking       *
C      *    by  T. HAUPT  (haupt@sccs.npac.syr.edu)    *
C      *                                               *
C      *    Northeast Parallel Architectures Center    *
C      *   at Syracuse University, Syracuse, NY, USA   *
C      *************************************************
C
C
C       VERSION SIMD/DECmpp-1.00
C       ==================================================
C
      INTEGER KASES,JFUNK,NFUNK
      PARAMETER (KASES=4)
      INTEGER N(KASES)
      DATA N / 8192,16384,65536,262144/
      DATA NFUNK /3/
      INTEGER METH,IFUN,NK
      REAL RESULT,TRUE,A,B,ERROR

      INTEGER ITIME
      REAL TIME
C      DO  IFUN=1,NFUNK
C     remember: there is corresponding ENDDO (!)
      IFUN=2
      CALL FVALS(A,B,TRUE,IFUN)
      DO  K = 1, KASES
         NK=N(K)
       DO METH=1,3
      CALL mpTIMERstart()
       DO MANY=1,50
         CALL DOIT(NK,A,B,METH,IFUN,RESULT)
       ENDDO
      ITIME=mpTIMERelapsed()
         ERROR = RESULT - TRUE
         PRINT *, ' '
         PRINT *,'PROBLEM 14 WITH N = ',NK
         PRINT *,'METHOD',METH,' FUNCTION ',IFUN
         PRINT *,'GIVES INTEGRAL ESTIMATE =', RESULT
         PRINT *,'ERROR (ESTIMATE-TRUE VALUE) = ',ERROR
      TIME=FLOAT(ITIME)/1000
      PRINT *, 'DECmpp ELAPSED TIME:',TIME,' SECONDS'
        ENDDO
       ENDDO
C      ENDDO
      STOP
      END


      SUBROUTINE  DOIT(NK,A,B,METH,IFUN,RESULT)
      INTEGER NK,METH,IFUN
      REAL A,B,RESULT
      INTEGER NSIMP,NG
      REAL H77
      REAL, ARRAY(0:NK)    :: X,X1,X2,X3
      REAL, ARRAY(0:NK)    :: F,F1,F2,F3
      REAL H

      INTERFACE
        SUBROUTINE FUN(X,N,IFUN,F)
        INTEGER N,IFUN
        REAL X(:),F(:)
      END INTERFACE
C
      IF(METH.EQ.1) THEN
C
C        TRAPEZOIDAL RULE
C
           H = (B-A)/NK
           RESULT = 0

           X = A + H * [0:NK]
           CALL FUN(X,NK,IFUN,F)
           RESULT = (SUM(F(1:NK-1))*2.0+F(0)+F(NK))*H/2.0

      ENDIF


      IF(METH.EQ.2) THEN
C
C        SIMPSON's METHOD
C
           NSIMP = NK
           IF (MOD(NSIMP,2).EQ.1) NSIMP = NSIMP-1
           H = (B-A)/NSIMP


           X = A + H * [0:NSIMP]
           CALL FUN(X,NSIMP,IFUN,F)
           RESULT=H*(F(0)+F(NSIMP)+4.0*SUM(F(1:NSIMP-1:2))+
     *           2.0*SUM(F(2:NSIMP-2:2)))/3.0
      ENDIF

      IF(METH.EQ.3) THEN
C
C        GAUSS' METHOD
C
            NG=(NK-MOD(NK,3))/3
            H = (B-A)/NG
            H77 = .774596669241*H


c           X1(0:NG)=A+H*[0:NG]-H/2.0-H77
c           X2(0:NG)=A+H*[0:NG]-H/2.0
c           X3(0:NG)=A+H*[0:NG]-H/2.0+H77


           CALL FUN(A+H*[0:NG]-H/2.0-H77,NG,IFUN,F1)
           CALL FUN(A+H*[0:NG]-H/2.0,NG,IFUN,F2)
           CALL FUN(A+H*[0:NG]-H/2.0+H77,NG,IFUN,F3)

           RESULT = H*(5.0*(SUM(F1(1:NG))+SUM(F3(1:NG)))+
     *             8.0*SUM(F2(1:NG)))/18.0


      ENDIF

      RETURN
      END


      SUBROUTINE FUN(X,N,IFUN,F)
      INTEGER N,IFUN
      REAL X(:),F(:)

      IF (IFUN.EQ.1) F = EXP(X)
      IF (IFUN.EQ.2) F = SQRT(ABS(X-.2345))
      IF (IFUN.EQ.3) F = 1.+X*X+1./(1.+100.*X*X)
      RETURN
      END

      SUBROUTINE FVALS (A,B,TRUE,IFUN)
      IF (IFUN.EQ.1) THEN
         A = 0.
         B = 1.
         TRUE = 1.71828182845
      ENDIF
      IF (IFUN.EQ.2) THEN
         A = 0.
         B = 1.
         TRUE = .5222099422093
      ENDIF
      IF (IFUN.EQ.3) THEN
         A = -1.
         B = 2.
         TRUE = 6.29919656054
      ENDIF
      RETURN
      END


