tutorial2r.f90

Go to the documentation of this file.

!> @file tutorial2r.f90
!! @ingroup FORT1THREAD_EXAMPLES
!! @brief This is a revised version of Tutorial2 in which the production
!! function is split into two equations using an intermediate variable.
!!
!! The purpose is to make the expressions simpler, in particular the
!! second order terms.
!! The equation
!! Row 2: `Al*L**(-Rho) + Ak*K**(-Rho) + Ainp*Inp**(-Rho)) ** ( -1.d0/Rho ) - Out = 0`
!! is replaced by
!! Row 2 (new): `Al*L**(-Rho) + Ak*K**(-Rho) + Ainp*Inp**(-Rho)) = Int`
!! Row 4 (new): `Int**( -1.d0/Rho ) = Out`
!! where `Int` is a new variable.
!!
!! 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 tutorial2r
 
   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, 5 ) )  ! 5 variables
   coi_error = max( coi_error, coidef_numcon( cntvect, 4 ) )  ! 4 constraints
   coi_error = max( coi_error, coidef_numnz( cntvect, 11 ) ) ! 11 nonzeros in the Jacobian
   coi_error = max( coi_error, coidef_numnlnz( cntvect, 5 ) )  ! 5 of which are nonlinear
   coi_error = max( coi_error, coidef_numhess( cntvect, 4 ) )  ! 4 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_debugfv( cntvect, 0 ) )  ! Debug Fdeval on or off
   coi_error = max( coi_error, coidef_debug2d( cntvect, 0 ) )  ! Debug 2nd derivatives on or off
   coi_error = max( coi_error, coidef_optfile( cntvect, 'tutorial.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 Tutorial2r 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 tutorial2r
 
!>  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
!
!  Parameters needed for new initial value -- could be moved to a
!  Module
!
   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
!
!  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
!
!  The added variable, in the following called Int, is given a lower
!  bound of 0.01 and an initial value consistent with the defining
!  equation, i.e. int = (Al*L**(-Rho) + Ak*K**(-Rho) + Ainp*Inp**(-Rho))
!
    lower(5) = 0.01d0
    curr(5)  = (al*curr(1)**(-rho) + ak*k**(-rho) + ainp*curr(2)**(-rho))
!
!  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
!
!  Constraint 4 (additional equation)
!  Rhs = 0.0 and type Equality
!
   type(4) = 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
!        L     INP  OUT    P    Int
!       x(1)  x(2)  x(3)  x(4)  x(5)
!  1:    1     3     5     8
!  2:    2     4                10
!  3:                6     9
!  4                 7          11
!
   colsta(1) =  1
   colsta(2) =  3
   colsta(3) =  5
   colsta(4) =  8
   colsta(5) = 10
   colsta(6) = 12
   rowno(1)  =  1
   rowno(2)  =  2
   rowno(3)  =  1
   rowno(4)  =  2
   rowno(5)  =  1
   rowno(6)  =  3
   rowno(7)  =  4
   rowno(8)  =  1
   rowno(9)  =  3
   rowno(10) =  2
   rowno(11) =  4
!
!  Nonlinearity Structure: L = 0 are linear and NL = 1 are nonlinear
!       x(1)  x(2)  x(3)  x(4)  x(5)
!  1:    L     L     NL    NL
!  2:    NL    NL                L
!  3:                L     L
!  4:                L           NL
!
   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
   nlflag(10)=  0
   nlflag(11)=  1
!
!  Value (Linear only)
!       x(1)  x(2)  x(3)  x(4)  x(5)
!  1:   -1    -1     NL    NL
!  2:    NL    NL                -1
!  3:                 1     2
!  4:                -1          NL
!
   value(1)  = -1.d0
   value(3)  = -1.d0
   value(6)  =  1.d0
   value(7)  = -1.d0
   value(9)  =  2.d0
   value(10) = -1.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, int
!
!  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
!
!  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)
   int = x(5)
!
!  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:
!         Al*L**(-Rho) + Ak*K**(-Rho) + Ainp*Inp**(-Rho))
!
   elseif ( rowno .eq. 2 ) then
!
!  Mode = 1 or 3: Function value
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g = al*l**(-rho) + ak*k**(-rho) + ainp*inp**(-rho)
      endif
!
!  Mode = 2 or 3: Derivatives
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         jac(1) = -rho * al   * l   ** (-rho-1.d0) ! derivative w.r.t. L   = X(1)
         jac(2) = -rho * ainp * inp ** (-rho-1.d0) ! derivative w.r.t. Inp = X(2)
      endif
!
!  Row = 4: Out = Int**(-1/Rho)
!
   elseif ( rowno .eq. 4 ) then
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g = int**(-1.d0/rho)
      endif
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         jac(5) = (-1.d0/rho)*int**(-1.d0/rho-1.d0) ! derivative w.r.t. Int = X(5)
      endif
 
   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))
!  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
!
!  where d1 = (-Rho)*(-Rho-1)*Al*L**(-Rho-2)
!  and   d2 = (-Rho)*(-Rho-1)*Ainp*Inp**(-Rho-2)
!
!  Row 4: Int**(-1/Rho)
!  Hessian (diagonal and lower triangle): A total of 1 nonzero element
!           Int
!             5
!  5   Int   d3
!
!  where d3 = (-1/Rho)*(-1/Rho-1)*Int**(-1/Rho-2)
!
!  The three 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   5
!  1    1
!  2        2
!  3
!  4            3
!  5                    4
!
!  and the same with names of variables and numbers of contributing
!  equations:
!
!       L  Inp Out  P  Int
!  L    2
!  Inp      2
!  Out
!  P            1
!  Int                  4
!
!
!  Define structure of Hessian
!
   hsrw(1) = 1; hscl(1) = 1
   hsrw(2) = 2; hscl(2) = 2
   hsrw(3) = 4; hscl(3) = 3
   hsrw(4) = 5; hscl(4) = 5
 
   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
!
!  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
!
!  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))
!  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
!
!  where d1 = (-Rho)*(-Rho-1)*Al*L**(-Rho-2)
!  and   d2 = (-Rho)*(-Rho-1)*Ainp*Inp**(-Rho-2)
!
!  Row 4: Int**(-1/Rho)
!  Hessian (diagonal and lower triangle): A total of 1 nonzero element
!           Int
!             5
!  5   Int   d3
!
!  where d3 = (-1/Rho)*(-1/Rho-1)*Int**(-1/Rho-2)
!
!  The three 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   5
!  1    1
!  2        2
!  3
!  4            3
!  5                    4
!
!  and the same with names of variables and numbers of contributing
!  equations:
!
!       L  Inp Out  P  Int
!  L    2
!  Inp      2
!  Out
!  P            1
!  Int                  4
!
!
!  Normal Evaluation mode
!
   hsvl(3) = u(1) ! Second derivative of constraint 1 is 1
   if ( u(2) .ne. 0.d0 ) Then
      l   = x(1)
      inp = x(2)
      hsvl(1) = (-rho) * (-rho-1.d0) * al * l**(-rho-2.d0) * u(2)
      hsvl(2) = (-rho) * (-rho-1.d0) * ainp * inp**(-rho-2.d0) * u(2)
   endif
   if ( u(4) .ne. 0.d0 ) then
      hsvl(4) = (-1.d0/rho)*(-1.d0/rho-1.d0)*x(5)**(-1.d0/rho-2.d0)
   endif
 
   tut_2dlagrval = 0
 
End Function tut_2dlagrval