      PROGRAM NODE15
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 MIMD/iPSC860-1.00-NODE   (F77+mp)
C       ==================================================
C
      INCLUDE 'picl_def.h'
      DOUBLE PRECISION HOST_SEC, NODE_SEC, STARTTIME
      INTEGER*4 TIMESIZE, TIMETYPE
      DATA TIMESIZE,TIMETYPE /8,10/
      INTEGER*4 ME,HOST,ALLNODES
      INTEGER*4 NPROCS,TOP,ORD,DIR
      INTEGER*4 INITTYPE
      DATA INITTYPE /1/
      DATA TOP,ORD,DIR,ALLNODES /1,1,1,-1/
C -----------------------------------------------------------
C           parameters of the initial broadcast
C -----------------------------------------------------------
      INTEGER*4 MSG(1),MSGSIZE
      DATA MSGSIZE /4/
      EQUIVALENCE (MSG(1),NK)
C -----------------------------------------------------------
C      parameters of messages with partial results
C -----------------------------------------------------------
      INTEGER*4 PARTTYP1,PARTTYP2
      DATA PARTTYP1,PARTTYP2 /20,21/
      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)
      COMMON /POINTS/ N,KASE
      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
      CALL open0(NPROCS,ME,HOST)
      CALL setarc0(NPROCS,TOP,ORD,DIR)

  100 CONTINUE
      STARTTIME=clock0()
      CALL bcast0(MSG,MSGSIZE,INITTYPE,HOST)
      IF(NK.GT.0) THEN
        IF (ME.LT.NPROCS) THEN


      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

            ERR1 = -10E20
            DO 70 I=ME+1, NK, NPROCS
               X = A+(I-1)*(B-A)/(NK-1)
               ERROR = ABS(F(X)-INTERP(X,COEF,TKNOTS,KNOTS,NPTS,NMAX,
     *            KNOTSMAX))
               IF (ERROR.GT.ERR1) THEN
                  ERR1 = ERROR
               ENDIF
   70       CONTINUE
            CALL gmax0(ERR1,1,4,ERRTYP1,0)
            IF(ME.EQ.0) ERRMAX(N,KASE)=ERR1
            CALL bcast0(ERRMAX(N,KASE),4,ERRTYP2,0)


C
C           END               'NUMBER OF POINTS USED'
C
   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

      IF(ME.EQ.0) CALL send0(ERRMAX,NPTS*8,PARTTYP1,HOST)
      IF(ME.EQ.0) CALL send0(DECAY,NPTS*8,PARTTYP2,HOST)
      NODE_SEC=clock0()-STARTTIME

      IF(ME.EQ.0) CALL send0(NODE_SEC,TIMESIZE,TIMETYPE,HOST)
      GOTO 100
      ENDIF

      ELSE

      CALL close0(1)
      WRITE(6,*) 'end of the job on node', ME
      STOP

      ENDIF
      END


      REAL FUNCTION INTERP (X,COEF,T,KNOTS,NPTS,NMAX,KNOTSMAX)
      REAL COEF(NMAX),T(KNOTSMAX)
      COMMON /POINTS/ N,KASE
C
C          PUT THE APPROPRIATE ARRAY OF KNOTS INTO THE T ARRAY

      INTERP = 0
      DO 10 J = 1, 2*N
         INTERP = INTERP+COEF(J)*HCUBIC(J,X,T,KNOTS)
   10 CONTINUE
      IF (INTERP.EQ.0) THEN
         INTERP = F(X)
      ENDIF
      RETURN
      END

      REAL FUNCTION HCUBIC (J,X,T,KNOTS)
C
C         THE HERMITE CUBIC KNOTS ARE IN THE ARRAY T OF LENGTH KNOTS
C           THERE ARE 2 SUCH ARRAYS, 1 FOR EACH CASE OF POINT DIST.
C         J = INDEX OF THE BASIS FUNCTION
C         X = POINT OF EVALUATION OF THE BASIS FUNCTION
C
      REAL T(KNOTS)
      IKNOT = (J+3)/2
      IF (MOD(J,2).EQ.1) THEN
         HCUBIC = HERMC(T,KNOTS,X,IKNOT)
      ELSE
         HCUBIC = HERMC1(T,KNOTS,X,IKNOT)
      ENDIF
      RETURN
      END

      REAL FUNCTION HERMC (T,KNOTS,X,IKNOT)
      REAL T(KNOTS)
      IF ((X.LT.T(IKNOT+1)).AND.(X.GT.T(IKNOT-1))) THEN
         HERMC = 0.0
         IF (X.GT.T(IKNOT)) THEN
            DT = T(IKNOT+1)-T(IKNOT)
            DX = T(IKNOT+1)-X
         ELSE
            DT = T(IKNOT)-T(IKNOT-1)
            DX = X-T(IKNOT-1)
         ENDIF
         HERMC = (3.-2.*DX/DT)*(DX**2)/DT**2
      ELSE
         HERMC = 0.0
      ENDIF
      RETURN
      END

      REAL FUNCTION HERMC1 (T,KNOTS,X,IKNOT)
      REAL T(KNOTS)
      IF ((X.LT.T(IKNOT+1)).AND.(X.GT.T(IKNOT-1))) THEN
         DX = X-T(IKNOT)
         IF (X.GT.T(IKNOT)) THEN
            DT2 = (T(IKNOT)-T(IKNOT+1))**2
            DX2 = (X-T(IKNOT+1))**2
         ELSE
            DT2 = (T(IKNOT)-T(IKNOT-1))**2
            DX2 = (X-T(IKNOT-1))**2
         ENDIF
         HERMC1 = DX2*DX/DT2
      ELSE
         HERMC1 = 0.0
      ENDIF
      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









