      USE DAESL_INT
      USE MP_TYPES
      IMPLICIT NONE 

!              NEQ = Number of equations
      INTEGER, PARAMETER :: NEQ=100
      REAL T, Y(NEQ), YPRIME(NEQ), TEND, ATOL(NEQ), CON
      INTEGER NOUT, IDO, I, NSTEPS
      REAL, TARGET :: RPARAM(NEQ)
      INTEGER, TARGET :: IPARAM(1)
      TYPE (S_FCN_DATA) USER_DATA
      EXTERNAL GCN
      
!              Pass NEQ and A^T eigenvector V to GCN 
      IPARAM(1) = NEQ 
      RPARAM(1) = 1.0
      DO I = 2, NEQ 
         RPARAM(I) = RPARAM(I-1) / FLOAT(I-1) 
      END DO
      USER_DATA%RDATA=>RPARAM
      USER_DATA%IDATA=>IPARAM
      
!              Initial values for y, y'
      Y = 0.0
      Y(1) = 1.0
      
      YPRIME = 0.0
      YPRIME(1) = -1.0
      YPRIME(2) =  1.0

      ATOL = 1.0E-4
      NSTEPS = 10
      
!              Always set IDO=1 on first call
      IDO = 1
      DO I = 1, NSTEPS
!              Output solution at T=1,2,...,10
        T = I-1
        TEND = I 
!              Set IDO = 3 on last call
        IF (I == NSTEPS) IDO = 3
!              User-defined Jacobian matrix structure (MATSTR=2) 
        CALL DAESL (T, TEND, IDO, Y, YPRIME, GCN, IYPR=0, MATSTR=2, &
                    NCON=1, ATOL=ATOL, FCN_DATA=USER_DATA)
      END DO 
      
!              Check if solution satisfies constraint
      CON = 0.0
      DO I = 1, NEQ
         CON = CON + RPARAM(I) * Y(I)
      END DO
      
      CALL UMACH(2, NOUT)
      WRITE (NOUT, *) ' V dot Y =', CON 
      END
  
  
      SUBROUTINE GCN (T, Y, YPRIME, DELTA, D, LDD, IRES, FCN_DATA)
      USE MP_TYPES
      IMPLICIT NONE
      
      REAL T, Y(*), YPRIME(*), DELTA(*), D(LDD,*), CON, CJ
      INTEGER IRES, LDD, I, NEQ
      SAVE CJ
      TYPE (S_FCN_DATA), OPTIONAL, INTENT(INOUT) :: FCN_DATA

      NEQ = FCN_DATA%IDATA(1)
      SELECT CASE (IRES)
      
!              F_I defined here
      CASE(1)
        DELTA(1) = YPRIME(1) + Y(1)
        DO I = 2, NEQ-1
           DELTA(I) = YPRIME(I) - Y(I-1) + Y(I) / FLOAT(I)   
        END DO
        DELTA(NEQ) = YPRIME(NEQ) - Y(NEQ-1)
        
!              Constraint is V dot Y = 1
      CASE(5)
        CON = -1.0
        DO I = 1, NEQ
           CON = CON + FCN_DATA%RDATA(I) * Y(I)
           D(1,I) = FCN_DATA%RDATA(I)
        END DO 
        DELTA(1) = CON
        
!              Normally, compute matrix A = dF/dY + CJ*dF/dY'  
!              = -B + CJ*I here.  Only CJ needs to be saved 
!              in this case, however, since B is bidiagonal, 
!              so A*x=DELTA can be solved (IRES=7) without 
!              saving or factoring B.
      CASE(6)
        CJ = DELTA(1)
!              If CJ > 0 not close to zero, A is nonsingular,
!              so set IRES = 0.
        IF (CJ >= 1.0E-4) IRES = 0
        
!              Solve A*x=DELTA and return x in DELTA.
      CASE(7)
        DELTA(1) = DELTA(1) / (1.0 + CJ)
        DO I = 2, NEQ-1
           DELTA(I) = (DELTA(I) + DELTA(I-1)) / (1.0 / FLOAT(I) + CJ)  
        END DO
        DELTA(NEQ) = (DELTA(NEQ) + DELTA(NEQ-1)) / CJ
        
      END SELECT
      END
