      PROGRAM PROB15
C
C     PROBLEM 15
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
      PARAMETER (KASES=2)
      INTEGER N(KASES)
      DATA N /8192,16384/
      INTEGER NK,NPTS
      PARAMETER(NPTS=10)
      REAL ERRMAX(NPTS,2),DECAY(NPTS,2)
      INTEGER ITIME
      REAL TIME

      DO 50 K = 1, KASES

      CALL mpTIMERstart()
      DO MANY=1,10
      NK=N(K)
      CALL DOIT(NK,ERRMAX,DECAY)
      ENDDO
      ITIME=mpTIMERelapsed()

      DO 170 KASE=1,2
         IF (KASE.EQ.1) THEN
            PRINT 110,NK
  110       FORMAT (4X,'PROBLEM 15 WITH NK=',I8/)
            PRINT 120
  120       FORMAT (9X,'EQUALLY SPACED POINTS'/)
            PRINT 130
  130       FORMAT (4X,'N      MAX ERROR      DECAY EXPONENT')
         ELSE
            PRINT 140
  140       FORMAT (/9X,'CHEBYSHEV SPACED POINTS'/)
            PRINT 130
         ENDIF
         DO 160 I = 2, NPTS
            PRINT 150,I,ERRMAX(I,KASE),DECAY(I,KASE)
  150       FORMAT (1X,I4,4X,E12.4,3X,F11.2)
  160    CONTINUE
C
C             PRINT OUT OF INTERPOLATION POINT ARRAY XPTS(J,N,KASE)
C             SUPPRESSED HERE TO REDUCE AMOUNT OF OUTPUT
C
  170 CONTINUE


      TIME=FLOAT(ITIME)/1000
      PRINT *, 'DECmpp ELAPSED TIME:',TIME,' SECONDS'

   50 CONTINUE
      STOP
      END

      SUBROUTINE DOIT(NK,ERRMAX,DECAY)
      INTEGER NK,NPTS
      PARAMETER(NPTS=10)
      REAL ERRMAX(NPTS,2),DECAY(NPTS,2)
      INTEGER NMAX,I,K,KNOTSMAX
      PARAMETER (NMAX=2*NPTS,KNOTSMAX=NPTS+2)
      REAL XPTS(NPTS),COEF(NMAX),H
      REAL TKNOTS(KNOTSMAX)
      REAL LOGERR,OLDERR,INTERP,RATIO,A,B
      DATA A,B / -1.,1. /,PI / 3.141592654 /
      INTEGER ERRTYP1,ERRTYP2
      DATA ERRTYP1,ERRTYP2 /31,32/
      REAL ERR1

      DO KASE = 1, 2
C
C        BEGIN             'NUMBER OF POINTS USED'
C
         DO 80 N = 2, NPTS

C
C           BEGIN                'CASES OF POINT DISTRIBUTION'
C
            IF (KASE.EQ.1) THEN
               H = (B-A)/(N-1)
               DO 10 J = 1, N
                  XPTS(J) = A+(J-1)*H
   10          CONTINUE
            ELSE
               DO 20 J = 1, N
                  XPTS(J) = (A+B)/2+(A-B)/2*COS((2*J-1)*PI/(2*N))
   20          CONTINUE
               RATIO = (B-A)/(XPTS(N)-XPTS(1))
               DO 30 J = 1, N
                  XPTS(J) = (XPTS(J)-(A+B)/2)*
     *                                RATIO+(A+B)/2
   30          CONTINUE
            ENDIF

            DO 40 J = 1, N
               COEF(2*J-1) = F(XPTS(J))
   40       CONTINUE
            DO 50 J = 1, N
               COEF(2*J) = FPRIME(XPTS(J))
   50       CONTINUE

            DO 60 J = 1, N
               TKNOTS(J+1) = XPTS(J)
   60       CONTINUE
C
C           ADD THE DUMMY POINTS AT EACH END OF TKNOTS ARRAY
C
            TKNOTS(1) = A-0.1*(ABS(A)+1.)
            TKNOTS(N+2) = B+0.1*(ABS(B)+1.)
            KNOTS = N+2

            CALL INTERPOLATE(N,NK,A,B,COEF,TKNOTS,ERRMAX(N,KASE))


   80    CONTINUE

         DO 90 N = 3, NPTS
            OLDERR = ALOG(ERRMAX(N-1,KASE))
            LOGERR = ALOG(ERRMAX(N,KASE))
            DECAY(N,KASE) = (LOGERR-OLDERR)/ALOG(FLOAT(N)/FLOAT(N-1))
   90    CONTINUE

      ENDDO


      RETURN
      END

      REAL FUNCTION F (X)
      F = 1./(X*X+25.0)
      RETURN
      END
      REAL FUNCTION FPRIME (X)
      FPRIME = -2.0*X/(X*X+25.)**2
      RETURN
      END

      SUBROUTINE FUN(X,NK,F)
      INTEGER NK
      REAL, ARRAY(NK) :: X,F
      F=1./(X*X+25.0)
      RETURN
      END

      SUBROUTINE INTERPOLATE(N,NK,A,B,COEF,T,ERROR)
      INTEGER NPTS,NMAX,KNOTSMAX
      PARAMETER (NPTS=10)
      PARAMETER (NMAX=2*NPTS,KNOTSMAX=NPTS+2)
      REAL COEF(NMAX),T(KNOTSMAX)
      REAL ERROR,A,B
      INTEGER N,NK,J,IKNOT
      REAL, ARRAY(NK) :: X,F,INTERP,DTDX,DTDX2,DX,COEF2
      REAL, ARRAY(NMAX,NK) :: HCUBIC
CMF$  LAYOUT HCUBIC(:SERIAL,:NEWS)
      INTERFACE
       SUBROUTINE FUN(X,NK,F)
         INTEGER NK
         REAL, ARRAY(NK)    ::X,F
      END INTERFACE

      X=A+[0:NK-1]*(B-A)/(NK-1)
      CALL FUN(X,NK,F)
      COEF2=0

      DO J=1,2*N
       DTDX=0
       DTDX2=0
       DX=0

       IKNOT=(J+3)/2
       IF(MOD(J,2).EQ.1) THEN

         WHERE((X.LT.T(IKNOT+1)) .AND. (X.GT.T(IKNOT)))
             DTDX=(T(IKNOT+1)-X)/(T(IKNOT+1)-T(IKNOT))
         ENDWHERE

         WHERE((X.LE.T(IKNOT)) .AND. (X.GT.T(IKNOT-1)))
             DTDX=(X-T(IKNOT-1))/(T(IKNOT)-T(IKNOT-1))
         ENDWHERE

         WHERE((X.LT.T(IKNOT+1)) .AND. (X.GT.T(IKNOT-1)))
             HCUBIC(J,1:NK)=(3.0-2.0*DTDX)*DTDX**2
         ELSEWHERE
           HCUBIC(J,1:NK)=0
         ENDWHERE

      ELSE

         WHERE((X.LT.T(IKNOT+1)) .AND. (X.GT.T(IKNOT)))
           DX=X-T(IKNOT)
           DTDX2=((X-T(IKNOT+1))/(T(IKNOT)-T(IKNOT+1)))**2
         ENDWHERE

         WHERE((X.LE.T(IKNOT)) .AND. (X.GT.T(IKNOT-1)))
           DX=X-T(IKNOT)
           DTDX2=((X-T(IKNOT-1))/(T(IKNOT)-T(IKNOT-1)))**2

         ENDWHERE

         WHERE((X.LT.T(IKNOT+1)) .AND. (X.GT.T(IKNOT-1)))
           HCUBIC(J,1:NK)=DTDX2*DX
         ELSEWHERE
           HCUBIC(J,1:NK)=0
         ENDWHERE

       ENDIF
       COEF2(1:NK)=COEF(J)*HCUBIC(J,1:NK)+COEF2(1:NK)
      ENDDO

C      INTERP=SUM(COEF2,DIM=1)
      INTERP=COEF2
      ERROR=MAXVAL(ABS(F-INTERP))

      RETURN
      END








