*     E02DCF Example Program Text
*     Mark 20 Revised. NAG Copyright 2001.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          MXMAX, MYMAX, NXEST, NYEST, NMAX, LWRK, LIWRK
      PARAMETER        (MXMAX=11,MYMAX=9,NXEST=MXMAX+4,NYEST=MYMAX+4,
     +                 NMAX=7,LWRK=4*(MXMAX+MYMAX)+11*(NXEST+NYEST)
     +                 +NXEST*MYMAX+NXEST+54,LIWRK=3+MXMAX+MYMAX+NXEST+
     +                 NYEST)
*     .. Local Scalars ..
      DOUBLE PRECISION DELTA, FP, S, XHI, XLO, YHI, YLO
      INTEGER          I, IFAIL, J, MX, MY, NPX, NPY, NX, NY
      CHARACTER        START
*     .. Local Arrays ..
      DOUBLE PRECISION C((NXEST-4)*(NYEST-4)), F(MXMAX*MYMAX),
     +                 FG(NMAX*NMAX), LAMDA(NXEST), MU(NYEST), PX(NMAX),
     +                 PY(NMAX), WRK(LWRK), X(MXMAX), Y(MYMAX)
      INTEGER          IWRK(LIWRK)
*     .. External Subroutines ..
      EXTERNAL         CPRINT, E02DCF, E02DFF
*     .. Intrinsic Functions ..
      INTRINSIC        MAX, MIN
*     .. Executable Statements ..
      WRITE (NOUT,*) 'E02DCF Example Program Results'
*     Skip heading in data file
      READ (NIN,*)
*     Input the number of X, Y co-ordinates MX, MY.
      READ (NIN,*) MX, MY
      IF (MX.GT.0 .AND. MX.LE.MXMAX .AND. MY.GT.0 .AND. MY.LE.MYMAX)
     +    THEN
*        Input the X co-ordinates followed by the Y co-ordinates.
         READ (NIN,*) (X(I),I=1,MX)
         READ (NIN,*) (Y(I),I=1,MY)
*        Input the MX*MY function values F at the grid points.
         READ (NIN,*) (F(I),I=1,MX*MY)
         START = 'C'
         READ (NIN,*,END=100) S
*        Determine the spline approximation.
         IFAIL = 0
*
         CALL E02DCF(START,MX,X,MY,Y,F,S,NXEST,NYEST,NX,LAMDA,NY,MU,C,
     +               FP,WRK,LWRK,IWRK,LIWRK,IFAIL)
*
         WRITE (NOUT,*)
         WRITE (NOUT,99999) 'Calling with smoothing factor S =', S,
     +     ': NX =', NX, ', NY =', NY, '.'
         WRITE (NOUT,*)
         WRITE (NOUT,*)
     +     '               I   Knot LAMDA(I)      J     Knot MU(J)'
         WRITE (NOUT,*)
         DO 20 J = 4, MAX(NX,NY) - 3
            IF (J.LE.NX-3 .AND. J.LE.NY-3) THEN
               WRITE (NOUT,99997) J, LAMDA(J), J, MU(J)
            ELSE IF (J.LE.NX-3) THEN
               WRITE (NOUT,99997) J, LAMDA(J)
            ELSE IF (J.LE.NY-3) THEN
               WRITE (NOUT,99996) J, MU(J)
            END IF
   20    CONTINUE
         CALL CPRINT(C,NY,NX,NOUT)
         WRITE (NOUT,*)
         WRITE (NOUT,99998) 'Sum of squared residuals FP =', FP
         IF (FP.EQ.0.0D+0) THEN
            WRITE (NOUT,*) '(The spline is an interpolating spline)'
         ELSE IF (NX.EQ.8 .AND. NY.EQ.8) THEN
            WRITE (NOUT,*)
     +        '(The spline is the least-squares bi-cubic polynomial)'
         END IF
*        Evaluate the spline on a rectangular grid at NPX*NPY points
*        over the domain (XLO to XHI) x (YLO to YHI).
         READ (NIN,*) NPX, XLO, XHI
         READ (NIN,*) NPY, YLO, YHI
         IF (NPX.LE.NMAX .AND. NPY.LE.NMAX) THEN
            DELTA = (XHI-XLO)/(NPX-1)
            DO 40 I = 1, NPX
               PX(I) = MIN(XLO+(I-1)*DELTA,XHI)
   40       CONTINUE
            DELTA = (YHI-YLO)/(NPY-1)
            DO 60 I = 1, NPY
               PY(I) = MIN(YLO+(I-1)*DELTA,YHI)
   60       CONTINUE
*
            CALL E02DFF(NPX,NPY,NX,NY,PX,PY,LAMDA,MU,C,FG,WRK,LWRK,IWRK,
     +                  LIWRK,IFAIL)
*
            WRITE (NOUT,*)
            WRITE (NOUT,*) 'Values of computed spline:'
            WRITE (NOUT,*)
            WRITE (NOUT,99995) '          X', (PX(I),I=1,NPX)
            WRITE (NOUT,*) '     Y'
            DO 80 I = NPY, 1, -1
               WRITE (NOUT,99994) PY(I), (FG(NPY*(J-1)+I),J=1,NPX)
   80       CONTINUE
         END IF
      END IF
  100 CONTINUE
      STOP
*
99999 FORMAT (1X,A,1P,D13.4,A,I5,A,I5,A)
99998 FORMAT (1X,A,1P,D13.4)
99997 FORMAT (1X,I16,F12.4,I11,F12.4)
99996 FORMAT (1X,I39,F12.4)
99995 FORMAT (1X,A,7F8.2)
99994 FORMAT (1X,F8.2,3X,7F8.2)
      END
*
      SUBROUTINE CPRINT(C,NY,NX,NOUT)
*     .. Scalar Arguments ..
      INTEGER           NOUT, NX, NY
*     .. Array Arguments ..
      DOUBLE PRECISION  C(NY-4,NX-4)
*     .. Local Scalars ..
      INTEGER           I, J
*     .. Executable Statements ..
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'The B-spline coefficients:'
      WRITE (NOUT,*)
      DO 20 I = 1, NY - 4
         WRITE (NOUT,99999) (C(I,J),J=1,NX-4)
   20 CONTINUE
      RETURN
*
99999 FORMAT (1X,8F9.4)
      END