      program main

c*********************************************************************72
c
cc toms660_test() tests toms660().
c
c  Discussion:
c
C      ALGORITHM 660, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 14, NO. 2, P.149.
C
C                     QS2TEST
C
C   THIS PROGRAM TESTS THE SCATTERED DATA INTERPOLATION
C PACKAGE QSHEP2D BY PRINTING THE MAXIMUM ERRORS ASSOCIATED
C WITH INTERPOLATED VALUES AND GRADIENTS ON A 10 BY 10
C UNIFORM GRID IN THE UNIT SQUARE.  THE DATA SET CONSISTS
C OF 36 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 QS2VAL AND QS2GRD ARE COMPARED FOR AGREEMENT.
C
      INTEGER LCELL(3,3), LNEXT(36)
      REAL    X(36), Y(36), F(36), RSQ(36), A(5,36), P(10)
C
C QSHEP2 PARAMETERS AND LOGICAL UNIT FOR OUTPUT
C
      DATA N/36/,  NQ/13/,  NW/19/,  NR/3/,  LOUT/6/
C
C QUADRATIC TEST FUNCTION AND PARTIAL DERIVATIVES
C
      FQ(XX,YY) = ((XX + 2.*YY)/3.)**2
      FX(XX,YY) = 2.*(XX + 2.*YY)/9.
      FY(XX,YY) = 4.*(XX + 2.*YY)/9.

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'toms660_test():'
      write ( *, '(a)' ) '  Fortran77 version.'
      write ( *, '(a)' ) '  Test toms660().'
C
C GENERATE A 6 BY 6 GRID OF NODES IN THE UNIT SQUARE WITH
C   THE NATURAL ORDERING.
C
      K = 0
      DO J = 1,6
      YK = FLOAT(6-J)/5.
      DO I = 1,6
        K = K + 1
        X(K) = FLOAT(I-1)/5.
             Y(K) = YK
        end do
      end do
C
C COMPUTE THE DATA VALUES.
C
      DO K = 1,N
           F(K) = FQ(X(K),Y(K))
      end do
C
C COMPUTE PARAMETERS DEFINING THE INTERPOLANT Q.
C
      CALL QSHEP2 (N,X,Y,F,NQ,NW,NR, LCELL,LNEXT,XMIN,YMIN,
     &  DX,DY,RMAX,RSQ,A,IER)
      IF (IER .NE. 0) GO TO 6
C
C GENERATE A 10 BY 10 UNIFORM GRID OF INTERPOLATION POINTS
C   (P(I),P(J)) IN THE UNIT SQUARE.  THE FOUR CORNERS COIN-
C   CIDE WITH NODES.
C
      DO I = 1,10
           P(I) = FLOAT(I-1)/9.
      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 QS2VAL AND QS2GRD.
C
      EQ = 0.
      EQX = 0.
      EQY = 0.
      DO J = 1,10
      PY = P(J)
      DO I = 1,10
        PX = P(I)
        Q1 = QS2VAL (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN,
     &      YMIN,DX,DY,RMAX,RSQ,A)
        CALL QS2GRD (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN,
     &      YMIN,DX,DY,RMAX,RSQ,A, Q,QX,QY,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)-Q))
        EQX = AMAX1(EQX,ABS(FX(PX,PY)-QX))
             EQY = AMAX1(EQY,ABS(FY(PX,PY)-QY))
        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
      STOP
  100 FORMAT (///1H ,31HMAXIMUM ABSOLUTE ERRORS IN THE ,
     &  25HINTERPOLANT Q AND PARTIAL/
     &  1H ,31HDERIVATIVES QX AND QY 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)
C
C ERROR IN QSHEP2
C
    6 WRITE (LOUT,200) IER
      STOP
  200 FORMAT (///1H ,28H*** ERROR IN QSHEP2 -- IER =,I2,
     &  4H ***)
C
C ERROR IN QS2GRD
C
    7 WRITE (LOUT,210) IER
      STOP
  210 FORMAT (///1H ,28H*** ERROR IN QS2GRD -- IER =,I2,
     &  4H ***)
C
C VALUES RETURNED BY QS2VAL AND QS2GRD DIFFER BY A RELATIVE
C   AMOUNT GREATER THAN 3*EPS.
C
    8 WRITE (LOUT,220) Q1, Q

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'toms660_test():'
      write ( *, '(a)' ) '  Normal end of execution.'

      STOP
  220 FORMAT (///1H ,33H*** ERROR -- INTERPOLATED VALUES ,
     &  37HQ1 (QS2VAL) AND Q2 (QS2GRD) DIFFER --//
     &  1H ,5X,5HQ1 = ,E21.14,5X,5HQ2 = ,E21.14)
      END
