*     E04YBF Example Program Text
*     Mark 14 Revised. NAG Copyright 1989.
*     .. Parameters ..
      INTEGER          MDEC, NDEC, LDFJAC, LB, LIW, LW
      PARAMETER        (MDEC=15,NDEC=3,LDFJAC=MDEC,LB=NDEC*(NDEC+1)/2,
     +                 LIW=1,LW=5*NDEC+MDEC+MDEC*NDEC+NDEC*(NDEC-1)/2)
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
*     .. Arrays in Common ..
      DOUBLE PRECISION T(MDEC,NDEC), Y(MDEC)
*     .. Local Scalars ..
      INTEGER          I, IFAIL, J, K, M, N
*     .. Local Arrays ..
      DOUBLE PRECISION B(LB), FJAC(LDFJAC,NDEC), FVEC(MDEC), W(LW),
     +                 X(NDEC)
      INTEGER          IW(LIW)
*     .. External Subroutines ..
      EXTERNAL         E04YAF, E04YBF, LSQFUN, LSQHES
*     .. Common blocks ..
      COMMON           Y, T
*     .. Executable Statements ..
      WRITE (NOUT,*) 'E04YBF Example Program Results'
*     Skip heading in data file
      READ (NIN,*)
      M = MDEC
      N = NDEC
*     Observations of TJ (J = 1, 2, 3) are held in T(I, J)
*     (I = 1, 2, . . . , 15)
      DO 20 I = 1, M
         READ (NIN,*) Y(I), (T(I,J),J=1,N)
   20 CONTINUE
*     Set up an arbitrary point at which to check the derivatives
      X(1) = 0.19D0
      X(2) = -1.34D0
      X(3) = 0.88D0
*     Check the 1st derivatives
      IFAIL = 1
*
      CALL E04YAF(M,N,LSQFUN,X,FVEC,FJAC,LDFJAC,IW,LIW,W,LW,IFAIL)
      IF (IFAIL.NE.0) THEN
         WRITE (NOUT,99998) IFAIL
         GO TO 60
      END IF
*
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'The test point is'
      WRITE (NOUT,99999) (X(J),J=1,N)
*
*     Check the evaluation of B
      IFAIL = 1
*
      CALL E04YBF(M,N,LSQFUN,LSQHES,X,FVEC,FJAC,LDFJAC,B,LB,IW,LIW,W,LW,
     +            IFAIL)
*
      WRITE (NOUT,*)
      IF (IFAIL.EQ.1) THEN
         WRITE (NOUT,*) 'A parameter is outside its expected range'
      ELSE IF (IFAIL.GE.0) THEN
         IF (IFAIL.EQ.0) THEN
            WRITE (NOUT,*)
     +        'The matrix B is consistent with 1st derivatives'
         ELSE IF (IFAIL.EQ.2) THEN
            WRITE (NOUT,*)
     +        'Probable error in calculation of the matrix B'
         END IF
         WRITE (NOUT,*)
         WRITE (NOUT,*) 'At the test point, LSQFUN gives'
         WRITE (NOUT,*)
         WRITE (NOUT,*)
     +     '      Residuals                   1st derivatives'
         WRITE (NOUT,99997) (FVEC(I),(FJAC(I,J),J=1,N),I=1,M)
         WRITE (NOUT,*)
         WRITE (NOUT,*)
     +     'and LSQHES gives the lower triangle of the matrix B'
         WRITE (NOUT,*)
         K = 1
         DO 40 I = 1, N
            WRITE (NOUT,99997) (B(J),J=K,K+I-1)
            K = K + I
   40    CONTINUE
      ELSE IF (IFAIL.EQ.-399) THEN
         WRITE (NOUT,99996) IFAIL
      ELSE
         WRITE (NOUT,99995) 'IFLAG was set to ', IFAIL,
     +     'in LSQFUN or LSQHES'
      END IF
   60 CONTINUE
*
99999 FORMAT (1X,4F10.5)
99998 FORMAT (1X,/1X,' ** E04YAF returned with IFAIL = ',I5)
99997 FORMAT (1X,1P,4E15.3)
99996 FORMAT (1X,/1X,' ** E04YBF returned with IFAIL = ',I5)
99995 FORMAT (1X,A,I3,A)
      END
*
      SUBROUTINE LSQFUN(IFLAG,M,N,XC,FVEC,FJAC,LDFJAC,IW,LIW,W,LW)
*     Routine to evaluate the residuals and their 1st derivatives
*     .. Parameters ..
      INTEGER          MDEC, NDEC
      PARAMETER        (MDEC=15,NDEC=3)
*     .. Scalar Arguments ..
      INTEGER          IFLAG, LDFJAC, LIW, LW, M, N
*     .. Array Arguments ..
      DOUBLE PRECISION FJAC(LDFJAC,N), FVEC(M), W(LW), XC(N)
      INTEGER          IW(LIW)
*     .. Arrays in Common ..
      DOUBLE PRECISION T(MDEC,NDEC), Y(MDEC)
*     .. Local Scalars ..
      DOUBLE PRECISION DENOM, DUMMY
      INTEGER          I
*     .. Common blocks ..
      COMMON           Y, T
*     .. Executable Statements ..
      DO 20 I = 1, M
         DENOM = XC(2)*T(I,2) + XC(3)*T(I,3)
         FVEC(I) = XC(1) + T(I,1)/DENOM - Y(I)
         FJAC(I,1) = 1.0D0
         DUMMY = -1.0D0/(DENOM*DENOM)
         FJAC(I,2) = T(I,1)*T(I,2)*DUMMY
         FJAC(I,3) = T(I,1)*T(I,3)*DUMMY
   20 CONTINUE
      RETURN
      END
*
      SUBROUTINE LSQHES(IFLAG,M,N,FVEC,XC,B,LB,IW,LIW,W,LW)
*     Routine to compute the lower triangle of the matrix B
*     (stored by rows in the array B)
*     .. Parameters ..
      INTEGER          MDEC, NDEC
      PARAMETER        (MDEC=15,NDEC=3)
*     .. Scalar Arguments ..
      INTEGER          IFLAG, LB, LIW, LW, M, N
*     .. Array Arguments ..
      DOUBLE PRECISION B(LB), FVEC(M), W(LW), XC(N)
      INTEGER          IW(LIW)
*     .. Arrays in Common ..
      DOUBLE PRECISION T(MDEC,NDEC), Y(MDEC)
*     .. Local Scalars ..
      DOUBLE PRECISION DUMMY, SUM22, SUM32, SUM33
      INTEGER          I
*     .. Common blocks ..
      COMMON           Y, T
*     .. Executable Statements ..
      B(1) = 0.0D0
      B(2) = 0.0D0
      SUM22 = 0.0D0
      SUM32 = 0.0D0
      SUM33 = 0.0D0
      DO 20 I = 1, M
         DUMMY = 2.0D0*T(I,1)/(XC(2)*T(I,2)+XC(3)*T(I,3))**3
         SUM22 = SUM22 + FVEC(I)*DUMMY*T(I,2)**2
         SUM32 = SUM32 + FVEC(I)*DUMMY*T(I,2)*T(I,3)
         SUM33 = SUM33 + FVEC(I)*DUMMY*T(I,3)**2
   20 CONTINUE
      B(3) = SUM22
      B(4) = 0.0D0
      B(5) = SUM32
      B(6) = SUM33
      RETURN
      END