*     D02AGF Example Program Text
*     Mark 14 Revised. NAG Copyright 1989.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
*     .. External Subroutines ..
      EXTERNAL         EX1, EX2
*     .. Executable Statements ..
      WRITE (NOUT,*) 'D02AGF Example Program Results'
      CALL EX1
      CALL EX2
      STOP
      END
*
      SUBROUTINE EX1
*     .. Parameters ..
      INTEGER          N, M1
      PARAMETER        (N=2,M1=6)
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
*     .. Scalars in Common ..
      INTEGER          IPRINT
*     .. Local Scalars ..
      DOUBLE PRECISION DUM, H, R, X, X1
      INTEGER          I, IFAIL, J, N1
*     .. Local Arrays ..
      DOUBLE PRECISION C(M1,N), COPY(N,N), ERROR(N), G(N), G1(N),
     +                 MAT(N,N), PARAM(N), PARERR(N), WSPACE(N,9)
*     .. External Subroutines ..
      EXTERNAL         AUX1, BCAUX1, D02AGF, PRSOL, RNAUX1
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Common blocks ..
      COMMON           /BLOCK1/IPRINT
*     .. Executable Statements ..
      WRITE (NOUT,*)
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Case 1'
      WRITE (NOUT,*)
*     * Set IPRINT to 1 to obtain output from PRSOL at each iteration *
      IPRINT = 0
      PARAM(1) = 0.2D0
      PARAM(2) = 0.0D0
      N1 = 2
      H = 0.1D0
      PARERR(1) = 1.0D-5
      PARERR(2) = 1.0D-3
      ERROR(1) = 1.0D-4
      ERROR(2) = 1.0D-4
      IFAIL = 1
*
      CALL D02AGF(H,ERROR,PARERR,PARAM,C,N,N1,M1,AUX1,BCAUX1,RNAUX1,
     +            PRSOL,MAT,COPY,WSPACE,G,G1,IFAIL)
*
      IF (IFAIL.EQ.0) THEN
         WRITE (NOUT,*) 'Final parameters'
         WRITE (NOUT,99998) (PARAM(I),I=1,N1)
         WRITE (NOUT,*)
         WRITE (NOUT,*) 'Final solution'
         WRITE (NOUT,*) 'X-value     Components of solution'
         CALL RNAUX1(X,X1,R,PARAM)
         H = (X1-X)/5.0D0
         DO 20 I = 1, 6
            DUM = X + DBLE(I-1)*H
            WRITE (NOUT,99997) DUM, (C(I,J),J=1,N)
   20    CONTINUE
      ELSE
         WRITE (NOUT,99999) 'IFAIL = ', IFAIL
      END IF
      RETURN
*
99999 FORMAT (1X,A,I3)
99998 FORMAT (1X,3D16.6)
99997 FORMAT (1X,F7.2,3D13.4)
      END
*
      SUBROUTINE AUX1(F,Y,X,PARAM)
*     .. Scalar Arguments ..
      DOUBLE PRECISION X
*     .. Array Arguments ..
      DOUBLE PRECISION F(2), PARAM(2), Y(2)
*     .. Executable Statements ..
      F(1) = Y(2)
      F(2) = (Y(1)**3-Y(2))/(2.0D0*X)
      RETURN
      END
*
      SUBROUTINE RNAUX1(X,X1,R,PARAM)
*     .. Scalar Arguments ..
      DOUBLE PRECISION  R, X, X1
*     .. Array Arguments ..
      DOUBLE PRECISION  PARAM(2)
*     .. Executable Statements ..
      X = 0.1D0
      X1 = 16.0D0
      R = 16.0D0
      RETURN
      END
*
      SUBROUTINE BCAUX1(G,G1,PARAM)
*     .. Array Arguments ..
      DOUBLE PRECISION  G(2), G1(2), PARAM(2)
*     .. Local Scalars ..
      DOUBLE PRECISION  Z
*     .. Intrinsic Functions ..
      INTRINSIC         SQRT
*     .. Executable Statements ..
      Z = 0.1D0
      G(1) = 0.1D0 + PARAM(1)*SQRT(Z)*0.1D0 + 0.01D0*Z
      G(2) = PARAM(1)*0.05D0/SQRT(Z) + 0.01D0
      G1(1) = 1.0D0/6.0D0
      G1(2) = PARAM(2)
      RETURN
      END
*
      SUBROUTINE EX2
*     .. Parameters ..
      INTEGER          N, M1
      PARAMETER        (N=3,M1=6)
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
*     .. Scalars in Common ..
      INTEGER          IPRINT
*     .. Local Scalars ..
      DOUBLE PRECISION DUM, H, R, X, X1
      INTEGER          I, IFAIL, J
*     .. Local Arrays ..
      DOUBLE PRECISION C(M1,N), COPY(N,N), ERROR(N), G(N), G1(N),
     +                 MAT(N,N), PARAM(N), PARERR(N), WSPACE(N,9)
*     .. External Subroutines ..
      EXTERNAL         AUX2, BCAUX2, D02AGF, PRSOL, RNAUX2
*     .. Intrinsic Functions ..
      INTRINSIC        DBLE
*     .. Common blocks ..
      COMMON           /BLOCK1/IPRINT
*     .. Executable Statements ..
      WRITE (NOUT,*)
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Case 2'
      WRITE (NOUT,*)
*     * Set IPRINT to 1 to obtain output from PRSOL at each iteration *
      IPRINT = 0
      H = 10.0D0
      PARAM(1) = 32.0D0
      PARAM(2) = 6000.0D0
      PARAM(3) = 0.54D0
      PARERR(1) = 1.0D-5
      PARERR(2) = 1.0D-4
      PARERR(3) = 1.0D-4
      ERROR(1) = 1.0D-2
      ERROR(2) = 1.0D-2
      ERROR(3) = 1.0D-2
      IFAIL = 1
*
      CALL D02AGF(H,ERROR,PARERR,PARAM,C,N,N,M1,AUX2,BCAUX2,RNAUX2,
     +            PRSOL,MAT,COPY,WSPACE,G,G1,IFAIL)
*
      IF (IFAIL.EQ.0) THEN
         WRITE (NOUT,*) 'Final parameters'
         WRITE (NOUT,99998) (PARAM(I),I=1,N)
         WRITE (NOUT,*)
         WRITE (NOUT,*) 'Final solution'
         WRITE (NOUT,*) 'X-value     Components of solution'
         CALL RNAUX2(X,X1,R,PARAM)
         H = (X1-X)/5.0D0
         DO 20 I = 1, 6
            DUM = X + DBLE(I-1)*H
            WRITE (NOUT,99997) DUM, (C(I,J),J=1,N)
   20    CONTINUE
      ELSE
         WRITE (NOUT,99999) 'IFAIL = ', IFAIL
      END IF
      RETURN
*
99999 FORMAT (1X,A,I3)
99998 FORMAT (1X,3D16.6)
99997 FORMAT (1X,F7.0,3D13.4)
      END
*
      SUBROUTINE AUX2(F,Y,X,PARAM)
*     .. Scalar Arguments ..
      DOUBLE PRECISION X
*     .. Array Arguments ..
      DOUBLE PRECISION F(3), PARAM(3), Y(3)
*     .. Local Scalars ..
      DOUBLE PRECISION C, S
*     .. Intrinsic Functions ..
      INTRINSIC       COS, SIN
*     .. Executable Statements ..
      C = COS(Y(3))
      S = SIN(Y(3))
      F(1) = S/C
      F(2) = -(PARAM(1)*S+0.00002D0*Y(2)*Y(2))/(Y(2)*C)
      F(3) = -PARAM(1)/(Y(2)*Y(2))
      RETURN
      END
*
      SUBROUTINE RNAUX2(X,X1,R,PARAM)
*     .. Scalar Arguments ..
      DOUBLE PRECISION  R, X, X1
*     .. Array Arguments ..
      DOUBLE PRECISION  PARAM(3)
*     .. Executable Statements ..
      X = 0.0D0
      X1 = PARAM(2)
      R = PARAM(2)
      RETURN
      END
*
      SUBROUTINE BCAUX2(G,G1,PARAM)
*     .. Array Arguments ..
      DOUBLE PRECISION  G(3), G1(3), PARAM(3)
*     .. Executable Statements ..
      G(1) = 0.0D0
      G(2) = 500.0D0
      G(3) = 0.5D0
      G1(1) = 0.0D0
      G1(2) = 450.0D0
      G1(3) = PARAM(3)
      RETURN
      END
*
      SUBROUTINE PRSOL(PARAM,RESID,N1,ERR)
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
*     .. Scalar Arguments ..
      DOUBLE PRECISION RESID
      INTEGER          N1
*     .. Array Arguments ..
      DOUBLE PRECISION ERR(N1), PARAM(N1)
*     .. Scalars in Common ..
      INTEGER          IPRINT
*     .. Local Scalars ..
      INTEGER          I
*     .. Common blocks ..
      COMMON           /BLOCK1/IPRINT
*     .. Executable Statements ..
      IF (IPRINT.NE.0) THEN
         WRITE (NOUT,99999) 'Current parameters   ', (PARAM(I),I=1,N1)
         WRITE (NOUT,99998) 'Residuals   ', (ERR(I),I=1,N1)
         WRITE (NOUT,99998) 'Sum of residuals squared  ', RESID
         WRITE (NOUT,*)
      END IF
      RETURN
*
99999 FORMAT (1X,A,6(D14.6,2X))
99998 FORMAT (1X,A,6(D12.4,1X))
      END
