*     E01DAF Example Program Text
*     Mark 14 Release. NAG Copyright 1989.
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          MXMAX, MYMAX
      PARAMETER        (MXMAX=20,MYMAX=MXMAX)
      INTEGER          LIWRK, LWRK
      PARAMETER        (LIWRK=MXMAX+2*(MXMAX-3)*(MYMAX-3),LWRK=(MXMAX+6)
     +                 *(MYMAX+6))
*     .. Local Scalars ..
      DOUBLE PRECISION STEP, XHI, XLO, YHI, YLO
      INTEGER          I, IFAIL, J, MX, MY, NX, NY, PX, PY
*     .. Local Arrays ..
      DOUBLE PRECISION C(MXMAX*MYMAX), F(MXMAX*MYMAX), FG(MXMAX*MYMAX),
     +                 LAMDA(MXMAX+4), MU(MYMAX+4), TX(MXMAX),
     +                 TY(MYMAX), WRK(LWRK), X(MXMAX), Y(MYMAX)
      INTEGER          IWRK(LIWRK)
      CHARACTER*10     CLABS(MYMAX), RLABS(MXMAX)
*     .. External Subroutines ..
      EXTERNAL         E01DAF, E02DFF, X04CBF
*     .. Intrinsic Functions ..
      INTRINSIC        MAX, MIN
*     .. Executable Statements ..
      WRITE (NOUT,*) 'E01DAF Example Program Results'
*     Skip heading in data file
      READ (NIN,*)
*     Read the number of X points, MX, and the values of the
*     X co-ordinates.
      READ (NIN,*) MX
      READ (NIN,*) (X(I),I=1,MX)
*     Read the number of Y points, MY, and the values of the
*     Y co-ordinates.
      READ (NIN,*) MY
      READ (NIN,*) (Y(I),I=1,MY)
*     Read the function values at the grid points.
      DO 20 J = 1, MY
         READ (NIN,*) (F(MY*(I-1)+J),I=1,MX)
   20 CONTINUE
      IFAIL = 1
*
*     Generate the (X,Y,F) interpolating bicubic B-spline.
      CALL E01DAF(MX,MY,X,Y,F,PX,PY,LAMDA,MU,C,WRK,IFAIL)
*
      IF (IFAIL.EQ.0) THEN
*
*        Print the knot sets, LAMDA and MU.
         WRITE (NOUT,*)
         WRITE (NOUT,*)
     +     '               I   Knot LAMDA(I)      J     Knot MU(J)'
         DO 40 J = 4, MAX(PX,PY) - 3
            IF (J.LE.PX-3 .AND. J.LE.PY-3) THEN
               WRITE (NOUT,99997) J, LAMDA(J), J, MU(J)
            ELSE IF (J.LE.PX-3) THEN
               WRITE (NOUT,99997) J, LAMDA(J)
            ELSE IF (J.LE.PY-3) THEN
               WRITE (NOUT,99996) J, MU(J)
            END IF
   40    CONTINUE
*        Print the spline coefficients.
         WRITE (NOUT,*)
         WRITE (NOUT,*) 'The B-Spline coefficients:'
         WRITE (NOUT,99999) (C(I),I=1,MX*MY)
         WRITE (NOUT,*)
*        Evaluate the spline on a regular rectangular grid at NX*NY
*        points over the domain (XLO to XHI) x (YLO to YHI).
         READ (NIN,*) NX, XLO, XHI
         READ (NIN,*) NY, YLO, YHI
         IF (NX.LE.MXMAX .AND. NY.LE.MYMAX) THEN
            STEP = (XHI-XLO)/(NX-1)
            DO 60 I = 1, NX
*              Generate NX equispaced X co-ordinates.
               TX(I) = MIN(XLO+(I-1)*STEP,XHI)
*              Generate X axis labels for printing results.
               WRITE (CLABS(I),99998) TX(I)
   60       CONTINUE
            STEP = (YHI-YLO)/(NY-1)
            DO 80 I = 1, NY
               TY(I) = MIN(YLO+(I-1)*STEP,YHI)
               WRITE (RLABS(I),99998) TY(I)
   80       CONTINUE
*
*           Evaluate the spline.
            IFAIL = 0
            CALL E02DFF(NX,NY,PX,PY,TX,TY,LAMDA,MU,C,FG,WRK,LWRK,IWRK,
     +                  LIWRK,IFAIL)
*
*           Print the results.
            IFAIL = 0
            CALL X04CBF('General','X',NY,NX,FG,NY,'F8.3',
     +          'Spline evaluated on a regular mesh (X across, Y down):'
     +                  ,'Character',RLABS,'Character',CLABS,80,0,IFAIL)
*
         END IF
      ELSE
         WRITE (NOUT,99995) IFAIL
      END IF
*
99999 FORMAT (1X,8F9.4)
99998 FORMAT (F5.2)
99997 FORMAT (1X,I16,F12.4,I11,F12.4)
99996 FORMAT (1X,I39,F12.4)
99995 FORMAT (1X,' ** E01DAF returned with IFAIL = ',I5)
      END