program main

!*****************************************************************************80
!
!! MAIN is the main program for DQED_PRB2.
!
!  Discussion:
!
!    DQED_PRB2 tests DQED.
!
!    The program illustrates the use of the Hanson-Krogh nonlinear least 
!    squares solver QED for fitting two exponentials to data.
! 
!    The problem is to find values for the four variables x(1),...,x(4),
!    which specify the model function
! 
!      h(t) = x(1)*exp(x(2)*t) + x(3)*exp(x(4)*t)
!
!    which mininize the sum of the squares of the error between h(t) and
!    recorded values of h at five values of t:
!
!      t = 0.05, 0.1, 0.4, 0.5, and 1.0.
!
!    We also have problem constraints that 
!
!          0 <= x(1)
!      -25.0 <= x(2) <= 0
!          0 <= x(3)
!      -25.0 <= x(4) <= 0. 
!
!    and a minimal separation requirement between x(2) and x(4) that
!    can be expressed as
!
!       0.05 <= x(2) - x(4)
!
!
!  Output:
! 
!    model is h(t) = x(1)*exp(-t*x(2)) + x(3)*exp(t*x(4))
!    x(1),x(2),x(3),x(4)  = 
!       1.999475    -.999801     .500057   -9.953988
!     residual after the fit =   4.2408d-04
!     output flag from solver =                      4
! 
  implicit none

  integer ( kind = 4 ), parameter :: liwork = 84
  integer ( kind = 4 ), parameter :: lwork = 640
  integer ( kind = 4 ), parameter :: mcon = 1
  integer ( kind = 4 ), parameter :: mequa = 5
  integer ( kind = 4 ), parameter :: nvars = 4

  integer ( kind = 4 ), parameter :: ldfj = mcon + mequa

  real ( kind = 8 ) bl(nvars+mcon)
  real ( kind = 8 ) bu(nvars+mcon)
  real ( kind = 8 ) fj(ldfj,nvars+1)
  real ( kind = 8 ) fnorm
  integer ( kind = 4 ) igo
  integer ( kind = 4 ) ind(nvars+mcon)
  integer ( kind = 4 ) iopt(24)
  integer ( kind = 4 ) iwork(liwork)
  real ( kind = 8 ) ropt(1)
  external dqedhd
  real ( kind = 8 ) work(lwork)
  real ( kind = 8 ) x(nvars)
  real ( kind = 8 ), dimension ( nvars ) :: x_expected = &
    (/ 1.999475D+00, -0.999801D+00, 0.500057D+00, -9.953988D+00 /)

  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DQED_PRB2'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  A sample calling program for DQED.'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  The problem is to find values for the four variables'
  write ( *, '(a)' ) '  x(1), x(2), x(3), x(4), which specify the model function:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    h(t) = x(1) * exp(x(2)*t) + x(3) * exp(x(4)*t)'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  which mininize the sum of the squares of the error '
  write ( *, '(a)' ) '  between h(t) and recorded values at five values of t:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    t = 0.05, 0.1, 0.4, 0.5, and 1.0.'
!
!  Define the bounding constraints on the variables X(1) through X(4),
!  and then the linear constraints on the separation condition.
!
!  1: 0.0 <= X(1)
!
  ind(1) = 1
  bl(1) = 0.0D+00
  bu(1) = 0.0D+00
!
!  2: -25.0 <= X(2) <= 0.0
!
  ind(2) = 3
  bl(2) = -25.0D+00
  bu(2) = 0.0D+00
!
!  3: 0.0 <= X(3)
!
  ind(3) = 1
  bl(3) = 0.0D+00
  bu(3) = 0.0D+00
!
!  4: -25.0 <= X(4) <= 0.0
!
  ind(4) = 3
  bl(4) = -25.0D+00
  bu(4) = 0.0D+00
!
!  5: 0.05 <= X(2) - X(4)
!
  ind(5) = 1
  bl(5) = 0.05D+00
  bu(5) = 0.0D+00
!
!  Define the initial values of the variables.
!  We don't know anything more, so all variables are set zero.
!
  x(1:nvars) = 0.0D+00
!
!  Tell how much storage we gave the solver.
!
  iwork(1) = lwork
  iwork(2) = liwork
!
!  No additional options are in use.
!
  iopt(1) = 99
!
!  Call the program.
!
  call dqed ( dqedhd, mequa, nvars, mcon, ind, bl, bu, x, fj, ldfj, &
    fnorm, igo, iopt, ropt, iwork, work )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'h(t) = x(1)*exp(t*x(2)*t) + x(3)*exp(t*x(4))'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Output flag from DQED, IGO = ', igo
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Computed X:'
  write ( *, '(a)' ) ' '
  write ( *, '(g14.6)' ) x(1:nvars)
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  L2 norm of the residual, FNORM = ', fnorm
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Expected X:'
  write ( *, '(a)' ) ' '
  write ( *, '(g14.6)' ) x_expected(1:nvars)
!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DQED_PRB2'
  write ( *, '(a)' ) '  Normal end of execution.'

  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine dqedhd ( x, fj, ldfj, igo, iopt, ropt )

!*****************************************************************************80
!
!! DQEDHD evaluates the functions and derivatives for DQED.
!
!  Discussion:
!
!    The user problem has MCON constraint functions,
!    MEQUA least squares equations, and involves NVARS
!    unknown variables.
! 
!    When this subprogram is entered, the general (near)
!    linear constraint partial derivatives, the derivatives
!    for the least squares equations, and the associated
!    function values are placed into the array FJ(*,*).
!
!    All partials and functions are evaluated at the point
!    in X(*).  Then the subprogram returns to the calling
!    program unit.  Typically one could do the following
!    steps:
! 
!      step 1. Place the partials of the i-th constraint
!      function with respect to variable j in the
!      array FJ(i,j), i = 1,...,MCON, j=1,...,NVARS.
!
!      step 2. Place the values of the i-th constraint
!      equation into FJ(i,NVARS+1).
!
!      step 3. Place the partials of the i-th least squares
!      equation with respect to variable j in the
!      array FJ(MCON+i,j), i = 1,...,MEQUA,
!      j = 1,...,NVARS.
!
!      step 4. Place the value of the i-th least squares
!      equation into FJ(MCON+i,NVARS+1).
!
!      step 5. Return to the calling program unit.
!
  implicit none

  integer ( kind = 4 ) ldfj
  integer ( kind = 4 ), parameter :: mcon = 1
  integer ( kind = 4 ), parameter :: mequa = 5
  integer ( kind = 4 ), parameter :: nvars = 4

  real ( kind = 8 ), save, dimension ( mequa ) :: f = &
    (/ 2.206D+00, 1.994D+00, 1.350D+00, 1.216D+00, 0.7358D+00 /)
  real ( kind = 8 ) fj(ldfj,nvars+1)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) igo
  integer ( kind = 4 ) iopt(*)
  real ( kind = 8 ) ropt(*)
  real ( kind = 8 ), save, dimension ( mequa ) :: t = &
    (/ 0.05D+00, 0.10D+00, 0.40D+00, 0.50D+00, 1.00D+00 /)
  real ( kind = 8 ) x(nvars)
! 
!  Set the value of the constraint.
!
  fj(1,nvars+1) = x(2) - x(4)
! 
!  Set the value of the residual functions.
!
  do i = 1, mequa

    fj(mcon+i,nvars+1) = &
        x(1) * exp ( x(2) * t(i) ) &
      + x(3) * exp ( x(4) * t(i) ) &
      - f(i)

  end do
!
!  If IGO is nonzero, compute the derivatives.
!
  if ( igo /= 0 ) then

    fj(1,1) = 0.0D+00
    fj(1,2) = 1.0D+00
    fj(1,3) = 0.0D+00
    fj(1,4) = -1.0D+00

    do i = 1, mequa

      fj(mcon+i,1) = exp ( x(2) * t(i) )
      fj(mcon+i,2) = x(1) * t(i) * exp ( x(2) * t(i) )
      fj(mcon+i,3) = exp ( x(4) * t(i) )
      fj(mcon+i,4) = x(3) * t(i) * exp ( x(4) * t(i) )

    end do

  end if

  return
end