program main

!*****************************************************************************80
!
!! MAIN is the main program for DQED_PRB3.
!
!  Discussion:
!
!    DQED_PRB3 tests DQED.
!
!    This program illustrates the use of the Hanson-Krogh nonlinear 
!    least squares solver DQED 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
!     number of evaluations of parameter model =    14
!     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 ), 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)
  real ( kind = 8 ) fnorm
  integer ( kind = 4 ) i
  integer ( kind = 4 ) igo
  integer ( kind = 4 ) ind(nvars+mcon)
  integer ( kind = 4 ) iopt(24)
  integer ( kind = 4 ) iwork(liwork)
  integer ( kind = 4 ) niters
  real ( kind = 8 ) ropt(1)
  external dqedev
  real ( kind = 8 ), dimension ( mequa ) :: t = (/ &
    0.05D+00, 0.10D+00, 0.40D+00, 0.50D+00, 1.00D+00 /)
  real ( kind = 8 ) work(lwork)
  real ( kind = 8 ) x(nvars)

  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DQED_PRB3'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  A test for DQED'
!
!  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 at all, 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
!
!  Initialize the call counter.
!
  niters = 0
!
!  Use reverse commumication to evaluate the derivatives.
!
  iopt(1) = 16
  iopt(2) = 1
!
!  No more options.
!
  iopt(3) = 99

  do

    call dqed ( dqedev, mequa, nvars, mcon, ind, bl, bu, x, fj, ldfj, fnorm, &
      igo, iopt, ropt, iwork, work )

    if ( 1 < igo ) then
      exit
    end if
!
!  Count function evaluations.
!
    niters = niters + 1
! 
!  Set the value of the constraint, and the residual functions.
!
    fj(1,5) = x(2) - x(4)

    do i = 1, mequa

      fj(mcon+i,5) = 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

  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Model:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '    h(t) = x(1)*exp(-t*x(2)) + x(3)*exp(t*x(4))'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Constrained minimizing X:'
  write ( *, '(a)' ) ' '
  write ( *, '(g14.6)' ) x(1:nvars)
  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) ' Residual of the fitted data, FNORM = ',fnorm
  write ( *, '(a,i6)' ) '  Number of model evaluations NITERS = ',niters
  write ( *, '(a,i6)' ) '  DQED output flag IGO = ',igo
!
!  Terminate.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'DQED_PRB3'
  write ( *, '(a)' ) '  Normal end of execution.'

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

  stop
end