tria04.f90

Go to the documentation of this file.

!> @file tria04.f90
!! @ingroup FORT1THREAD_EXAMPLES
!!
!! 
!! Triangular Demo model 04
!! 
!! This is a CONOPT implementation of the GAMS model:
!! 
!! @verbatim
!! variable x1, x2, x3  ;
!! equation e1, e2, e3  ;
!! 
!! e1 .. x1 + sqr(x2) =e= 5;
!! 
!! e2 .. 5*x2         =e= 4;
!! 
!! e3 .. x3           =e= x2 + x1;
!! @endverbatim
!! 
!! This model is triangular and the solution is known to be unique
!! even though it is nonlinear. x2 is determined uniquely from e2 and it
!! is no longer variable when it is used nonlinearly in e1.
!! CONOPT will therefore return model status = 1 or "Optimal" and
!! not the usual model status = 2 or "Locally Optimal".
!! 
!! @verbatim
!! model tria04 / all /;
!! 
!! solve tria04 using nlp minimizing x3;
!! @endverbatim
!! 
!! 
!! Add a lower bound on x1 that is too high and solve again.
!! The model will now be infeasible due to the bound.
!! The infeasibility is not local since x2 is unique from e2 and
!! x1 appears linearly in e1.
!! 
!! @verbatim
!! x1.lo = 5;
!! 
!! solve tria04 using nlp minimizing x3;
!! @endverbatim
!! 
!!
!! For more information about the individual callbacks, please have a look at the source code.
 
 
!> Main program. A simple setup and call of CONOPT
!!
Program tria04
 
   Use proginfo
   Use coidef
   use casedata_num
   implicit None
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: tria_readmatrix ! Mandatory Matrix definition routine defined below
   Integer, External :: tria_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 :: std_triord     ! callback for triangular order
#if defined(itl)
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Tria_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Tria_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 :: Std_TriOrd
#endif
!
!  Control vector
!
   INTEGER :: numcallback
   INTEGER, Dimension(:), Pointer :: cntvect
   INTEGER :: coi_error
 
   call startup
!
!  Create and initialize a Control Vector
!
   numcallback = coidef_size()
   Allocate( cntvect(numcallback) )
   coi_error = coidef_inifort( cntvect )
!
!  Tell CONOPT about the size of the model by populating the Control Vector:
!
   coi_error = max( coi_error, coidef_numvar( cntvect, 3 ) )  ! 3 variables
   coi_error = max( coi_error, coidef_numcon( cntvect, 3 ) )  ! 3 constraints
   coi_error = max( coi_error, coidef_numnz( cntvect, 6 ) )  ! 6 nonzeros in the Jacobian
   coi_error = max( coi_error, coidef_numnlnz( cntvect, 1 ) )  ! 1 of which are nonlinear
   coi_error = max( coi_error, coidef_optdir( cntvect, -1 ) ) ! Minimize
   coi_error = max( coi_error, coidef_objvar( cntvect, 3 ) )  ! Objective is variable 3
   coi_error = max( coi_error, coidef_optfile( cntvect, 'tria04.opt' ) )
!
!  Tell CONOPT about the callback routines:
!
   coi_error = max( coi_error, coidef_readmatrix( cntvect, tria_readmatrix ) )
   coi_error = max( coi_error, coidef_fdeval( cntvect, tria_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_triord( cntvect, std_triord     ) )
 
#if defined(LICENSE_INT_1) && defined(LICENSE_INT_2) && defined(LICENSE_INT_3) && defined(LICENSE_TEXT)
   coi_error = max( coi_error, coidef_license( cntvect, license_int_1, license_int_2, license_int_3, 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
!
!  Save the solution so we can check the duals:
!
   do_allocate = .true.
!
!  Start CONOPT:
!
   casenum = 1
   coi_error = coi_solve( cntvect )
 
   If ( coi_error /= 0 ) then
      call flog( "Solve 1: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Solve 1: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 1 ) then
      call flog( "Solve 1: Solver and Model Status was not as expected (1,1)", 1 )
   elseif ( abs( obj-5.16d0 ) > 0.000001d0 ) then
      call flog( "Solve 1: Incorrect objective returned", 1 )
   Else
      Call checkdual( 'Tria04-Opt', minimize )
   endif
 
   casenum = 2
   coi_error = coi_solve( cntvect )
 
   If ( coi_error /= 0 ) then
      call flog( "Solve 2: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Solve 2: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 4 ) then
      call flog( "Solve 2: Solver and Model Status was not as expected (1,4)", 1 )
! No objective test for infeasible model
   Elseif ( abs( xprim(1)+xprim(2)**2-uprim(1) ) > 1.d-7 ) then
      call flog( "Solve 2: Activity level in constraint 1 is not corrrect.", 1 )
   Elseif ( abs( 5*xprim(2)-uprim(2) ) > 1.d-7 ) then
      call flog( "Solve 2: Activity level in constraint 2 is not corrrect.", 1 )
   Elseif ( abs( xprim(3) - xprim(2) - xprim(1) - uprim(3) ) > 1.d-7 ) then
      call flog( "Solve 2: Activity level in constraint 3 is not corrrect.", 1 )
   Else
      Call checkdual( 'Tria04-Inf', infeasible )
   endif
 
   if ( coi_free(cntvect) /= 0 ) call flog( "Error while freeing control vector",1)
 
   write(*,*)
   write(*,*) 'End of Tria04 example. Return code=',coi_error
   call flog( "Successful Solve", 0 )
 
End Program tria04
!
! ============================================================================
!  Define information about the model:
!
 
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function tria_readmatrix( lower, curr, upper, vsta, type, rhs, esta,      &
                                 colsta, rowno, value, nlflag, n, m, nz,  &
                                 usrmem )
#if defined(itl)
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Tria_ReadMatrix
#endif
   use casedata_num
   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.
!
!  The model uses defaults, except for the second solve where we have
!  a lower bound on x1
!
   If ( casenum == 2 ) lower(1) = 5.0d0
!
!  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: e1
!  Rhs = 5.0 and type Equality
!
   rhs(1)  = 5.0d0
   type(1) = 0
!
!  Constraint 2: e2
!  Rhs = 4.0 and type Equality
!
   rhs(2)  = 4.0d0
   type(2) = 0
!
!  Constraint 3: e3
!  Rhs = 0.0 and type Equality
!
   type(3) = 0
!
!  Information about the Jacobian. We use the standard method with
!  Rowno, Value, Nlflag and Colsta and we do not use Colno.
!
!  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)
!  1:    1     3
!  2:          4
!  3:    2     5     6
!
   colsta(1) =  1
   colsta(2) =  3
   colsta(3) =  6
   colsta(4) =  7
   rowno(1)  =  1
   rowno(2)  =  3
   rowno(3)  =  1
   rowno(4)  =  2
   rowno(5)  =  3
   rowno(6)  =  3
!
!  Nonlinearity Structure: L = 0 are linear and NL = 1 are nonlinear
!       x(1)  x(2)  x(3)
!  1:    L     NL
!  2:          L
!  3:    L     L     L
!
   nlflag(1) =  0
   nlflag(2) =  0
   nlflag(3) =  1
   nlflag(4) =  0
   nlflag(5) =  0
   nlflag(6) =  0
!
!  Value (Linear only)
!       x(1)  x(2)  x(3)
!  1:   1.0    NL
!  2:         5.0
!  3:  -1.0  -1.0   1.0
!
   value(1)  =  1.d0
   value(2)  = -1.d0
   value(4)  =  5.d0
   value(5)  = -1.d0
   value(6)  =  1.d0
 
  tria_readmatrix = 0 ! Return value means OK
 
end Function tria_readmatrix
!
!==========================================================================
!  Compute nonlinear terms and non-constant Jacobian elements
!
 
!>  Compute nonlinear terms and non-constant Jacobian elements
!!
!! @include{doc} fdeval_params.dox
Integer Function tria_fdeval( x, g, jac, rowno, jcnm, mode, ignerr, errcnt, &
                              n, nz, thread, usrmem )
#if defined(itl)
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Tria_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
!
!  Row 1: e1
!
   if ( rowno .eq. 1 ) then
!
!  Mode = 1 or 3. G = sqr(x2)
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g    = x(2)*x(2)
      endif
!
!  Mode = 2 or 3: Derivative values:
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         jac(2) = 2.d0*x(2)
      endif
      tria_fdeval = 0
!
!  Row 2 and 3: The rows are linear and will not be called
!
   else
      tria_fdeval = 1
   endif
 
end Function tria_fdeval