* E02DDF Example Program Text * Mark 20 Revised. NAG Copyright 2001. * .. Parameters .. INTEGER NIN, NOUT PARAMETER (NIN=5,NOUT=6) INTEGER MMAX, NXEST, NYEST, NMAX, LWRK, LIWRK PARAMETER (MMAX=30,NXEST=14,NYEST=14,NMAX=7, + LWRK=(7*(NXEST-4)*(NYEST-4)+25*(NXEST-4)) + *(NXEST-3)+2*((NXEST-4)+(NYEST-4)+4*MMAX) + +23*(NXEST-4)+56,LIWRK=MMAX+2*(NXEST-7)*(NYEST-7) + ) * .. Local Scalars .. DOUBLE PRECISION DELTA, FP, S, XHI, XLO, YHI, YLO INTEGER I, IFAIL, J, M, NPX, NPY, NX, NY, RANK CHARACTER START * .. Local Arrays .. DOUBLE PRECISION C((NXEST-4)*(NYEST-4)), F(MMAX), FG(NMAX*NMAX), + LAMDA(NXEST), MU(NYEST), PX(NMAX), PY(NMAX), + W(MMAX), WRK(LWRK), X(MMAX), Y(MMAX) INTEGER IWRK(LIWRK) * .. External Subroutines .. EXTERNAL CPRINT, E02DDF, E02DFF * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. Executable Statements .. WRITE (NOUT,*) 'E02DDF Example Program Results' * Skip heading in data file READ (NIN,*) * Input the number of data-points M. READ (NIN,*) M IF (M.GT.0 .AND. M.LE.MMAX) THEN * Input the data-points and the weights. DO 20 I = 1, M READ (NIN,*) X(I), Y(I), F(I), W(I) 20 CONTINUE START = 'C' READ (NIN,*,END=120) S * Determine the spline approximation. IFAIL = 0 * CALL E02DDF(START,M,X,Y,F,W,S,NXEST,NYEST,NX,LAMDA,NY,MU,C,FP, + RANK,WRK,LWRK,IWRK,LIWRK,IFAIL) * WRITE (NOUT,*) WRITE (NOUT,99999) 'Calling with smoothing factor S =', S, + ': NX =', NX, ', NY =', NY, ',' WRITE (NOUT,99998) 'rank deficiency =', (NX-4)*(NY-4) - RANK WRITE (NOUT,*) WRITE (NOUT,*) + ' I Knot LAMDA(I) J Knot MU(J)' WRITE (NOUT,*) DO 40 J = 4, MAX(NX,NY) - 3 IF (J.LE.NX-3 .AND. J.LE.NY-3) THEN WRITE (NOUT,99996) J, LAMDA(J), J, MU(J) ELSE IF (J.LE.NX-3) THEN WRITE (NOUT,99996) J, LAMDA(J) ELSE IF (J.LE.NY-3) THEN WRITE (NOUT,99995) J, MU(J) END IF 40 CONTINUE CALL CPRINT(C,NY,NX,NOUT) WRITE (NOUT,*) WRITE (NOUT,99997) ' Sum of squared residuals FP =', FP 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 60 I = 1, NPX PX(I) = MIN(XLO+(I-1)*DELTA,XHI) 60 CONTINUE DO 80 I = 1, NPY PY(I) = MIN(YLO+(I-1)*DELTA,YHI) 80 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,99994) ' X', (PX(I),I=1,NPX) WRITE (NOUT,*) ' Y' DO 100 I = NPY, 1, -1 WRITE (NOUT,99993) PY(I), (FG(NPY*(J-1)+I),J=1,NPX) 100 CONTINUE END IF END IF 120 CONTINUE STOP * 99999 FORMAT (1X,A,1P,D13.4,A,I5,A,I5,A) 99998 FORMAT (1X,A,I5) 99997 FORMAT (1X,A,1P,D13.4,A) 99996 FORMAT (1X,I16,F12.4,I11,F12.4) 99995 FORMAT (1X,I39,F12.4) 99994 FORMAT (1X,A,7F8.2) 99993 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,7F9.2) END