tutorial2.f90

Go to the documentation of this file.

!> @file tutorial2.f90
!! @ingroup FORT1THREAD_EXAMPLES
!! @brief This model is a revision of Tutorial in which we have added a
!! set of 2nd derivative routines, Tut_2DLagrStr and Tut_2DLagrVal.
!!
!! For more information about the individual callbacks, please have a look at the source code.
 
#if defined(_WIN32) && !defined(_WIN64)
#define dec_directives_win32
#endif
 
!> Main program. A simple setup and call of CONOPT
!!
Program tutorial2
 
   Use proginfo
   Use conopt
   implicit None
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: tut_readmatrix ! Mandatory Matrix definition routine defined below
   Integer, External :: tut_fdeval     ! Function and Derivative evaluation routine
                                       ! needed a nonlinear model.
   Integer, External :: std_status     ! Standard callback for displaying solution status
   Integer, External :: std_solution   ! Standard callback for displaying solution values
   Integer, External :: std_message    ! Standard callback for managing messages
   Integer, External :: std_errmsg     ! Standard callback for managing error messages
   Integer, External :: tut_2dlagrstr  ! Callback for second derivatives structure
   Integer, External :: tut_2dlagrval  ! Callback for second derivatives values
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Tut_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Tut_FDEval
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Status
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Solution
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Message
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_ErrMsg
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Tut_2DLagrStr
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Tut_2DLagrVal
#endif
!
!  Control vector
!
   INTEGER, Dimension(:), Pointer :: cntvect
   INTEGER :: coi_error
 
   call startup
!
!  Create and initialize a Control Vector
!
   coi_error = coi_create( cntvect )
!
!  Tell CONOPT about the size of the model by populating the Control Vector:
!
   coi_error = max( coi_error, coidef_numvar( cntvect, 4 ) )  ! 4 variables
   coi_error = max( coi_error, coidef_numcon( cntvect, 3 ) )  ! 3 constraints
   coi_error = max( coi_error, coidef_numnz( cntvect, 9 ) )  ! 9 nonzeros in the Jacobian
   coi_error = max( coi_error, coidef_numnlnz( cntvect, 4 ) )  ! 4 of which are nonlinear
   coi_error = max( coi_error, coidef_numhess( cntvect, 4 ) )  ! Numher of elements in the Hessian of the Lagrangian
   coi_error = max( coi_error, coidef_optdir( cntvect, 1 ) )  ! Maximize
   coi_error = max( coi_error, coidef_objcon( cntvect, 1 ) )  ! Objective is constraint 1
   coi_error = max( coi_error, coidef_debug2d( cntvect, 0 ) )  ! Debug 2nd derivatives on or off
   coi_error = max( coi_error, coidef_optfile( cntvect, 'tutorial2.opt' ) )
!
!  Tell CONOPT about the callback routines:
!
   coi_error = max( coi_error, coidef_readmatrix( cntvect, tut_readmatrix ) )
   coi_error = max( coi_error, coidef_fdeval( cntvect, tut_fdeval     ) )
   coi_error = max( coi_error, coidef_status( cntvect, std_status     ) )
   coi_error = max( coi_error, coidef_solution( cntvect, std_solution   ) )
   coi_error = max( coi_error, coidef_message( cntvect, std_message    ) )
   coi_error = max( coi_error, coidef_errmsg( cntvect, std_errmsg     ) )
   coi_error = max( coi_error, coidef_2dlagrstr( cntvect, tut_2dlagrstr  ) )
   coi_error = max( coi_error, coidef_2dlagrval( cntvect, tut_2dlagrval  ) )
 
#if defined(CONOPT_LICENSE_INT_1) && defined(CONOPT_LICENSE_INT_2) && defined(CONOPT_LICENSE_INT_3) && defined(CONOPT_LICENSE_TEXT)
   coi_error = max( coi_error, coidef_license( cntvect, conopt_license_int_1, conopt_license_int_2, conopt_license_int_3, conopt_license_text) )
#endif
 
   If ( coi_error .ne. 0 ) THEN
      write(*,*)
      write(*,*) '**** Fatal Error while loading CONOPT Callback routines.'
      write(*,*)
      call flog( "Skipping Solve due to setup errors", 1 )
   ENDIF
!
!  Start CONOPT:
!
   coi_error = coi_solve( cntvect )
 
   write(*,*)
   write(*,*) 'End of Tutorial2 example. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Solver and Model Status was not as expected (1,2)", 1 )
   elseif ( abs( obj-0.572943d0 ) > 0.000001d0 ) then
      call flog( "Incorrect objective returned", 1 )
   endif
 
   if ( coi_free(cntvect) /= 0 ) call flog( "Error while freeing control vector",1)
 
   call flog( "Successful Solve", 0 )
!
!  Free solution memory
!
   call finalize
 
End Program tutorial2
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function tut_readmatrix( lower, curr, upper, vsta, type, rhs, esta,      &
                                 colsta, rowno, value, nlflag, n, m, nz,  &
                                 usrmem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Tut_ReadMatrix
#endif
   implicit none
   integer, intent (in)                      :: n      ! number of variables
   integer, intent (in)                      :: m      ! number of constraints
   integer, intent (in)                      :: nz     ! number of nonzeros
   real*8,  intent (in out), dimension(n)    :: lower  ! vector of lower bounds
   real*8,  intent (in out), dimension(n)    :: curr   ! vector of initial values
   real*8,  intent (in out), dimension(n)    :: upper  ! vector of upper bounds
   integer, intent (in out), dimension(n)    :: vsta   ! vector of initial variable status
                                                       ! (not defined here)
   integer, intent (out),    dimension(m)    :: type   ! vector of equation types
   integer, intent (in out), dimension(m)    :: esta   ! vector of initial equation status
                                                       ! (not defined here)
   real*8,  intent (in out), dimension(m)    :: rhs    ! vector of right hand sides
   integer, intent (in out), dimension(n+1)  :: colsta ! vector with start of column indices
   integer, intent (out),    dimension(nz)   :: rowno  ! vector of row numbers
   integer, intent (in out), dimension(nz)   :: nlflag ! vector of nonlinearity flags
   real*8,  intent (in out), dimension(nz)   :: value  ! vector of matrix values
   real*8   usrmem(*)                                  ! optional user memory
!
!  Information about Variables:
!  Default: Lower = -Inf, Curr = 0, and Upper = +inf.
!  Default: the status information in Vsta is not used.
!
!  Lower bound on L   = X(1) = 0.1 and initial value = 0.5:
!
   lower(1) = 0.1d0
   curr(1)  = 0.5d0
!
!  Lower bound on INP = X(2) = 0.1 and initial value = 0.5:
!
   lower(2) = 0.1d0
   curr(2)  = 0.5d0
!
!  Lower bound on OUT = X(3) and P = X(4) are both 0 and the
!  default initial value of 0 is used:
!
   lower(3) = 0.d0
   lower(4) = 0.d0
!
!  Information about Constraints:
!  Default: Rhs = 0
!  Default: the status information in Esta and the function
!           value in FV are not used.
!  Default: Type: There is no default.
!           0 = Equality,
!           1 = Greater than or equal,
!           2 = Less than or equal,
!           3 = Non binding.
!
!  Constraint 1 (Objective)
!  Rhs = -0.1 and type Non binding
!
   rhs(1)  = -0.1d0
   type(1) = 3
!
!  Constraint 2 (Production Function)
!  Rhs = 0 and type Equality
!
   type(2) = 0
!
!  Constraint 3 (Price equation)
!  Rhs = 4.0 and type Equality
!
   rhs(3)  = 4.d0
   type(3) = 0
!
!  Information about the Jacobian. CONOPT expects a columnwise
!  representation in Rowno, Value, Nlflag and Colsta.
!
!  Colsta = Start of column indices (No Defaults):
!  Rowno  = Row indices
!  Value  = Value of derivative (by default only linear
!                                derivatives are used)
!  Nlflag = 0 for linear and 1 for nonlinear derivative
!           (not needed for completely linear models)
!
!  Indices
!       x(1)  x(2)  x(3)  x(4)
!  1:    1     3     5     8
!  2:    2     4     6
!  3:                7     9
!
   colsta(1) =  1
   colsta(2) =  3
   colsta(3) =  5
   colsta(4) =  8
   colsta(5) = 10
   rowno(1)  =  1
   rowno(2)  =  2
   rowno(3)  =  1
   rowno(4)  =  2
   rowno(5)  =  1
   rowno(6)  =  2
   rowno(7)  =  3
   rowno(8)  =  1
   rowno(9)  =  3
!
!  Nonlinearity Structure: L = 0 are linear and NL = 1 are nonlinear
!       x(1)  x(2)  x(3)  x(4)
!  1:    L     L     NL    NL
!  2:    NL    NL    L
!  3:                L     L
!
   nlflag(1) =  0
   nlflag(2) =  1
   nlflag(3) =  0
   nlflag(4) =  1
   nlflag(5) =  1
   nlflag(6) =  0
   nlflag(7) =  0
   nlflag(8) =  1
   nlflag(9) =  0
!
!  Value (Linear only)
!       x(1)  x(2)  x(3)  x(4)
!  1:   -1    -1     NL    NL
!  2:    NL    NL    -1
!  3:                 1     2
!
   value(1)  = -1.d0
   value(3)  = -1.d0
   value(6)  = -1.d0
   value(7)  =  1.d0
   value(9)  =  2.d0
 
  tut_readmatrix = 0 ! Return value means OK
 
end Function tut_readmatrix
!>  Compute nonlinear terms and non-constant Jacobian elements
!!
!! @include{doc} fdeval_params.dox
Integer Function tut_fdeval( x, g, jac, rowno, jcnm, mode, ignerr, errcnt, &
                             n, nz, thread, usrmem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Tut_FDEval
#endif
   implicit none
   integer, intent (in)                      :: n      ! number of variables
   integer, intent (in)                      :: rowno  ! number of the row to be evaluated
   integer, intent (in)                      :: nz     ! number of nonzeros in this row
   real*8,  intent (in),     dimension(n)    :: x      ! vector of current solution values
   real*8,  intent (in out)                  :: g      ! constraint value
   real*8,  intent (in out), dimension(n)    :: jac    ! vector of derivatives for current constraint
   integer, intent (in),     dimension(nz)   :: jcnm   ! list of variables that appear nonlinearly
                                                       ! in this row. Ffor information only.
   integer, intent (in)                      :: mode   ! evaluation mode: 1 = function value
                                                       ! 2 = derivatives, 3 = both
   integer, intent (in)                      :: ignerr ! if 1 then errors can be ignored as long
                                                       ! as errcnt is incremented
   integer, intent (in out)                  :: errcnt ! error counter to be incremented in case
                                                       ! of function evaluation errors.
   integer, intent (in)                      :: thread
   real*8   usrmem(*)                                  ! optional user memory
!
!  Declare local copies of the optimization variables. This is
!  just for convenience to make the expressions easier to read.
!
   real*8            :: l, inp, out, p
!
!  Declare parameters and their data values.
!
   real*8, parameter :: w    = 1.0d0
   real*8, parameter :: l0   = 0.1d0
   real*8, parameter :: pinp = 1.0d0
   real*8, parameter :: al   = 0.16d0
   real*8, parameter :: ak   = 2.0d0
   real*8, parameter :: ainp = 0.16d0
   real*8, parameter :: rho  = 1.0d0
   real*8, parameter :: k    = 4.0d0
   real*8            :: hold1, hold2, hold3 ! Intermediate results
!
!  Move the optimization variables from the X vector to a set
!  of local variables with the same names as the variables in
!  the model description. This is not necessary, but it should make
!  the equations easier to recognize.
!
   l   = x(1)
   inp = x(2)
   out = x(3)
   p   = x(4)
!
!  Row 1: the objective function is nonlinear
!
   if ( rowno .eq. 1 ) then
!
!  Mode = 1 or 3. Function value: G = P * Out
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g    = p * out
      endif
!
!  Mode = 2 or 3: Derivative values:
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         jac(3) = p     ! derivative w.r.t. Out = X(3)
         jac(4) = out   ! derivative w.r.t. P   = X(4)
      endif
!
!  Row 2: The production function is nonlinear
!
   elseif ( rowno .eq. 2 ) then
!
!  Compute some common terms
!
      hold1 = (al*l**(-rho) + ak*k**(-rho) + ainp*inp**(-rho))
      hold2 = hold1 ** ( -1.d0/rho )
!
!  Mode = 1 or 3: Function value
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g = hold2
      endif
!
!  Mode = 2 or 3: Derivatives
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         hold3  = hold2 / hold1
         jac(1) = hold3 * al   * l   ** (-rho-1.d0) ! derivative w.r.t. L   = X(1)
         jac(2) = hold3 * ainp * inp ** (-rho-1.d0) ! derivative w.r.t. Inp = X(2)
      endif
!
!  Row = 3: The row is linear and will not be called.
!
   endif
   tut_fdeval = 0
 
end Function tut_fdeval
 
!>  Specify the structure of the Lagrangian of the Hessian
!!
!! @include{doc} 2DLagrStr_params.dox
Integer Function tut_2dlagrstr( HSRW, HSCL, NODRV, N, M, NHESS, UsrMem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Tut_2DLagrStr
#endif
   Implicit None
 
   Integer, Intent (IN)                   :: n, m, nhess
   Integer, Intent (IN OUT)               :: nodrv
   Integer, Dimension(Nhess), Intent(Out) :: hsrw, hscl
   real*8,  Intent(IN OUT)                :: usrmem(*)
!
!  There are two nonlinear constraints, P * Out in row 1 and Sdef and
!  Al*L**(-Rho) + Ak*K**(-Rho) + Ainp*Inp**(-Rho)) ** ( -1.d0/Rho ) in
!  row 2 and they have the following Hessians:
!
!  Row1: P * Out
!  Hessian (diagonal and lower triangle): A total of 1 nonzero element
!      Out    P
!       3     4
!  3
!  4    1
!
!  Row 2: Al*L**(-Rho) + Ak*K**(-Rho) + Ainp*Inp**(-Rho)) ** ( -1.d0/Rho )
!         = hold1**(-1.d0/Rh0) = hold2
!  with L and Inp being the variables.
!  Hessian (diagonal and lower triangle): A total of 3 nonzero element
!       L    Inp
!       1     2
!   1   d1
!   2   d2   d3
!
!  where d1, d2, and d3 are derived below.
!
!  The two Hessian matrices do not overlap so the total number of nonzeros
!  is 4. The structure using numerical indices for rows and columns and
!  running indices for the Hessian element is:
!
!       1   2   3   4
!  1    1
!  2    2   3
!  3
!  4            4
!
!  and the same with names of variables and numbers of contributing
!  equations:
!
!       L  Inp Out  P
!  L    2
!  Inp  2   2
!  Out
!  P            1
!
!
!  Define structure of Hessian
!
   hsrw(1) = 1; hscl(1) = 1
   hsrw(2) = 2; hscl(2) = 1
   hsrw(3) = 2; hscl(3) = 2
   hsrw(4) = 4; hscl(4) = 3
 
   tut_2dlagrstr = 0
 
End Function tut_2dlagrstr
 
!>  Compute the Lagrangian of the Hessian
!!
!! @include{doc} 2DLagrVal_params.dox
Integer Function tut_2dlagrval( X, U, HSRW, HSCL, HSVL, NODRV, N, M, NHESS, UsrMem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Tut_2DLagrVal
#endif
   Implicit None
 
   Integer, Intent (IN)                    :: n, m, nhess
   Integer, Intent (IN OUT)                :: nodrv
   real*8,  Dimension(N), Intent (IN)      :: x
   real*8,  Dimension(M), Intent (IN)      :: u
   Integer, Dimension(Nhess), Intent (In)  :: hsrw, hscl
   real*8,  Dimension(Nhess), Intent (Out) :: hsvl
   real*8,  Intent(IN OUT)                 :: usrmem(*)
!
!  Declare local copies of the optimization variables. This is
!  just for convenience to make the expressions easier to read.
!
   real*8            :: l, inp, out, p
!
!  Declare parameters and their data values.
!
   real*8, parameter :: w    = 1.0d0
   real*8, parameter :: l0   = 0.1d0
   real*8, parameter :: pinp = 1.0d0
   real*8, parameter :: al   = 0.16d0
   real*8, parameter :: ak   = 2.0d0
   real*8, parameter :: ainp = 0.16d0
   real*8, parameter :: rho  = 1.0d0
   real*8, parameter :: k    = 4.0d0
   real*8            :: hold1, hold2, hold3, hold4 ! Intermediate results
   integer           :: i
!
!  There are two nonlinear constraints, P * Out in row 1 and Sdef and
!  Al*L**(-Rho) + Ak*K**(-Rho) + Ainp*Inp**(-Rho)) ** ( -1.d0/Rho ) in
!  row 2 and they have the following Hessians:
!
!  Row1: P * Out
!  Hessian (diagonal and lower triangle): A total of 1 nonzero element
!      Out    P
!       3     4
!  3
!  4    1
!
!  Row 2: Al*L**(-Rho) + Ak*K**(-Rho) + Ainp*Inp**(-Rho)) ** ( -1.d0/Rho )
!         = hold1**(-1.d0/Rh0) = hold2
!  with L and Inp being the variables.
!  Hessian (diagonal and lower triangle): A total of 3 nonzero element
!       L    Inp
!       1     2
!   1   d1
!   2   d2   d3
!
!  where d1, d2, and d3 are derived below.
!
!  The two Hessian matrices do not overlap so the total number of nonzeros
!  is 4. The structure using numerical indices for rows and columns and
!  running indices for the Hessian element is:
!
!       1   2   3   4
!  1    1
!  2    2   3
!  3
!  4            4
!
!  and the same with names of variables and numbers of contributing
!  equations:
!
!       L  Inp Out  P
!  L    2
!  Inp  2   2
!  Out
!  P            1
!
!
!  Evaluation mode
!
   hsvl(4) = u(1) ! Second derivative of constraint 1 is 1
   if ( u(2) .ne. 0.d0 ) Then
      l   = x(1)
      inp = x(2)
      out = x(3)
      p   = x(4)
      hold1 = (al*l**(-rho) + ak*k**(-rho) + ainp*inp**(-rho))
      hold2 = hold1 ** ( -1.d0/rho )
      hold3 = hold2 / hold1
!
!  The code for the first derivative was:
!        jac(1) = hold3 * Al   * L   ** (-Rho-1.d0)
!        jac(2) = hold3 * Ainp * Inp ** (-Rho-1.d0)
!  Define Hold4 as the derivative of Hold3 with respect to Hold1
!
     hold4 = hold3 / hold1 * (-1.d0/rho-1.d0)
!
!  The second derivatives are computed as:
!  (1) The derivative of Hold3 multiplied by the other terms PLUS
!  (2) Hold3 multiplied by the derivative of the other terms
!  where the derivative of Hold3 is Hold4 * the derivative of Hold1
!  with respect to the variable in question.
!
      hsvl(1) = hold4 * (-rho) * (al*l**(-rho-1.d0))**2 +     &
                hold3 * al   * (-rho-1.d0)*l**(-rho-2.d0)
      hsvl(2) = hold4 * (-rho) * (al*l**(-rho-1.d0)) *         &
                                 (ainp*inp**(-rho-1.d0))
      hsvl(3) = hold4 * (-rho) * (ainp*inp**(-rho-1.d0))**2 + &
                hold3 * ainp * (-rho-1.d0)*inp**(-rho-2.d0)
!
!  Scale the derivatives with the Lagrange multiplier, U(2):
!
      do i = 1, 3
         hsvl(i) = hsvl(i) * u(2)
      enddo
   endif
 
   tut_2dlagrval = 0
 
End Function tut_2dlagrval