      program main

c*********************************************************************72
c
cc toms661_test() tests toms661().
c
c  Discussion:
c
C      ALGORITHM 661, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 14, NO. 2, P.151.
C
C                     QS3TEST
C
C   THIS PROGRAM TESTS THE SCATTERED DATA INTERPOLATION
C PACKAGE QSHEP3D BY PRINTING THE MAXIMUM ERRORS ASSOCIATED
C WITH INTERPOLATED VALUES AND GRADIENTS ON A 5 BY 5 BY 5
C UNIFORM GRID IN THE UNIT CUBE.  THE DATA SET CONSISTS
C OF 64 NODES WITH DATA VALUES TAKEN FROM A QUADRATIC FUNC-
C TION FOR WHICH THE METHOD IS EXACT.  THE RATIO OF MAXIMUM
C INTERPOLATION ERROR RELATIVE TO THE MACHINE PRECISION IS
C ALSO PRINTED.  THIS SHOULD BE O(1).  THE INTERPOLATED
C VALUES FROM QS3VAL AND QS3GRD ARE COMPARED FOR AGREEMENT.
C
      INTEGER LCELL(3,3,3), LNEXT(64)
      REAL X(64), Y(64), Z(64), F(64), RSQ(64), A(9,64),
     .         XYZMIN(3), XYZDEL(3)
      REAL P(5)
C
C QSHEP3 PARAMETERS AND LOGICAL UNIT FOR OUTPUT
C
      DATA N/64/,  NQ/17/,  NW/32/,  NR/3/,  LOUT/6/
C
C QUADRATIC TEST FUNCTION AND PARTIAL DERIVATIVES
C
      FQ(XX,YY,ZZ) = ((XX + 2.*YY + 3.*ZZ)/6.)**2
      FX(XX,YY,ZZ) = (XX + 2.*YY + 3.*ZZ)/18.
      FY(XX,YY,ZZ) = (XX + 2.*YY + 3.*ZZ)/9.
      FZ(XX,YY,ZZ) = (XX + 2.*YY + 3.*ZZ)/6.

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'toms661_test():'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) '  Test toms661().'
C
C GENERATE A 4 BY 4 BY 4 GRID OF NODES IN THE UNIT CUBE.
C
      L = 0
      DO K = 1,4
        ZL = FLOAT(K-1)/3.
        DO J = 1,4
          YL = FLOAT(J-1)/3.
          DO I = 1,4
            L = L + 1
            X(L) = FLOAT(I-1)/3.
            Y(L) = YL
            Z(L) = ZL
          end do
        end do
      end do
C
C COMPUTE THE DATA VALUES.
C
      DO L = 1,N
        F(L) = FQ(X(L),Y(L),Z(L))
      end do
C
C COMPUTE PARAMETERS DEFINING THE INTERPOLANT Q.
C
      CALL QSHEP3 (N,X,Y,Z,F,NQ,NW,NR, LCELL,LNEXT,XYZMIN,
     .               XYZDEL,RMAX,RSQ,A,IER)
      IF (IER .NE. 0) GO TO 6
C
C GENERATE A 5 BY 5 BY 5 UNIFORM GRID OF INTERPOLATION
C   POINTS (P(I),P(J),P(K)) IN THE UNIT CUBE.  THE EIGHT
C   CORNERS COINCIDE WITH NODES.
C
      DO I = 1,5
        P(I) = FLOAT(I-1)/4.
      end do
C
C COMPUTE THE MACHINE PRECISION EPS.
C
      EPS = 1.
    4 EPS = EPS/2.
      EP1 = EPS + 1.
      IF (STORE(EP1) .GT. 1.) GO TO 4
      EPS = EPS*2.
C
C COMPUTE INTERPOLATION ERRORS AND TEST FOR AGREEMENT IN THE
C   Q VALUES RETURNED BY QS3VAL AND QS3GRD.
C
      EQ = 0.
      EQX = 0.
      EQY = 0.
      EQZ = 0.
      DO K = 1,5
        PZ = P(K)
        DO J = 1,5
          PY = P(J)
          DO I = 1,5
            PX = P(I)
            Q1 = QS3VAL (PX,PY,PZ,N,X,Y,Z,F,NR,LCELL,LNEXT,
     .                   XYZMIN,XYZDEL,RMAX,RSQ,A)
            CALL QS3GRD (PX,PY,PZ,N,X,Y,Z,F,NR,LCELL,LNEXT,
     .                   XYZMIN,XYZDEL,RMAX,RSQ,A, Q,QX,QY,
     .                   QZ,IER)
            IF (IER .NE. 0) GO TO 7
            IF (ABS(Q1-Q) .GT. 3.*ABS(Q)*EPS) GO TO 8
            EQ = AMAX1(EQ,ABS(FQ(PX,PY,PZ)-Q))
            EQX = AMAX1(EQX,ABS(FX(PX,PY,PZ)-QX))
            EQY = AMAX1(EQY,ABS(FY(PX,PY,PZ)-QY))
            EQZ = AMAX1(EQZ,ABS(FZ(PX,PY,PZ)-QZ))
          end do
        end do
      end do
C
C PRINT ERRORS AND THE RATIO EQ/EPS.
C
      RQ = EQ/EPS
      WRITE (LOUT,100)
      WRITE (LOUT,110) EQ, RQ
      WRITE (LOUT,120) EQX
      WRITE (LOUT,130) EQY
      WRITE (LOUT,140) EQZ
      STOP
  100 FORMAT (///1H ,31HMAXIMUM ABSOLUTE ERRORS IN THE ,
     .            25HINTERPOLANT Q AND PARTIAL/
     .            1H ,32HDERIVATIVES (QX,QY,QZ) RELATIVE ,
     .            24HTO MACHINE PRECISION EPS//
     .            1H ,10X,8HFUNCTION,3X,9HMAX ERROR,3X,
     .            13HMAX ERROR/EPS/)
  110 FORMAT (1H ,13X,1HQ,7X,E9.3,7X,F4.2)
  120 FORMAT (1H ,13X,2HQX,6X,E9.3)
  130 FORMAT (1H ,13X,2HQY,6X,E9.3)
  140 FORMAT (1H ,13X,2HQZ,6X,E9.3)
C
C ERROR IN QSHEP3
C
    6 WRITE (LOUT,200) IER
      STOP
  200 FORMAT (///1H ,28H*** ERROR IN QSHEP3 -- IER =,I2,
     .            4H ***)
C
C ERROR IN QS3GRD
C
    7 WRITE (LOUT,210) IER
      STOP
  210 FORMAT (///1H ,28H*** ERROR IN QS3GRD -- IER =,I2,
     .            4H ***)
C
C VALUES RETURNED BY QS3VAL AND QS3GRD DIFFER BY A RELATIVE
C   AMOUNT GREATER THAN 3*EPS.
C
    8 WRITE (LOUT,220) Q1, Q
      STOP
  220 FORMAT (///1H ,33H*** ERROR -- INTERPOLATED VALUES ,
     .            37HQ1 (QS3VAL) AND Q2 (QS3GRD) DIFFER --//
     .            1H ,5X,5HQ1 = ,E21.14,5X,5HQ2 = ,E21.14)
      END
