program main

c*********************************************************************72
c
cc MAIN is the main program for TOMS660_PRB.
c
c  Discussion:
c
c    TOMS660_PRB tests the TOMS660 library.
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.
C
C GENERATE A 6 BY 6 GRID OF NODES IN THE UNIT SQUARE WITH
C   THE NATURAL ORDERING.
C
      K = 0
      DO 1 J = 1,6
	YK = FLOAT(6-J)/5.
	DO 1 I = 1,6
	  K = K + 1
	  X(K) = FLOAT(I-1)/5.
    1	  Y(K) = YK
C
C COMPUTE THE DATA VALUES.
C
      DO 2 K = 1,N
    2	F(K) = FQ(X(K),Y(K))
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 3 I = 1,10
    3	P(I) = FLOAT(I-1)/9.
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 5 J = 1,10
	PY = P(J)
	DO 5 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))
    5	  EQY = AMAX1(EQY,ABS(FY(PX,PY)-QY))
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
      STOP
  220 FORMAT (///1H ,33H*** ERROR -- INTERPOLATED VALUES ,
     .	      37HQ1 (QS2VAL) AND Q2 (QS2GRD) DIFFER --//
     .	      1H ,5X,5HQ1 = ,E21.14,5X,5HQ2 = ,E21.14)
      END