square6.f90

Go to the documentation of this file.

!> @file square6.f90
!! @ingroup FORT1THREAD_EXAMPLES
!!
!! A square model where we pretend that the last three constraints are nonlinear except
!! for their last entry.
!! 
!! \f{eqnarray*}{
!! &x1 + x2           = 10 \\
!! &x1 - x2           = 0 \\
!! &x1 + x2 + x3      = 9 \\
!! &- x2 + x3 + x4 = 1
!! \f}
!! 
!! In a second case we add a lower bound on x3 and x4 of 0 so the model
!! becomes infeasible in both the last equations.
!! The model is similar to square5, but the last equation is binding
!! upper instead of lower.
!! 
!!
!! 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 square
 
   Use proginfo
   Use conopt
   Use casedata_num
   implicit None
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: sq_readmatrix  ! Mandatory Matrix definition routine defined below
   Integer, External :: sq_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
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Sq_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Sq_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
#endif
!
!  Control vector
!
   INTEGER, Dimension(:), Pointer :: cntvect
   INTEGER :: coi_error
 
   Integer :: i
   Logical :: error
 
   real*8, dimension(4) :: xsol1 = (/ 5.d0, 5.d0, -1.d0, 7.d0 /)
   real*8, dimension(4) :: usol1 = (/ 10.d0, 0.d0, 9.d0, 1.d0 /)
 
   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 ) )  ! # variables
   coi_error = max( coi_error, coidef_numcon( cntvect, 4 ) )  ! # constraints
   coi_error = max( coi_error, coidef_numnz( cntvect, 10 ) ) ! # nonzeros in the Jacobian
   coi_error = max( coi_error, coidef_numnlnz( cntvect, 6 ) )  ! 4 of which are nonlinear
   coi_error = max( coi_error, coidef_square( cntvect, 1 ) )  ! 1 means Square system
   coi_error = max( coi_error, coidef_optfile( cntvect, 'square6.opt' ) )
!
!  Tell CONOPT about the callback routines:
!
   coi_error = max( coi_error, coidef_readmatrix( cntvect, sq_readmatrix ) )
   coi_error = max( coi_error, coidef_fdeval( cntvect, sq_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     ) )
 
#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
!
!  Ask Std_Solution to allocate space for the solution and status
!  vectors and keep this information
!
   do_allocate = .true.
!
!  Start CONOPT:
!
   casenum = 1
   coi_error = coi_solve( cntvect )
 
   write(*,*)
   write(*,*) 'End of Square6 - case 1 example. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "Case 1: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Case 1: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat < 15 .or. mstat > 16 ) then
      call flog( "Case 1: Solver or Model status not as expected (1,15) or (1,16)", 1 )
   elseif ( obj /= 0.d0 ) Then
      call flog( "Case 1: Objective for square model was not as expected 0.0", 1 )
   else
!
!  Check the primal and dual solution itself
!
      error = .false.
      do i = 1, 4
         if ( abs(xprim(i)-xsol1(i)) > 1.d-7 ) error = .true.
         if ( abs(xdual(i)) > 1.d-7 ) error = .true.
         if ( abs(udual(i)) > 1.d-7 ) error = .true.
         if ( abs(uprim(i)- usol1(i)) > 1.e-7 ) error = .true.
      enddo
      if ( error ) call flog( "Case 1: Numerical solution was not as expected.", 1 )
!
!  Check the status information
!
      do i = 1, 4
         if ( xbasc(i) /= 2 ) error = .true.  ! Basic
         if ( xstat(i) /= 0 ) error = .true.  ! Normal
!         if ( Ubasc(i) /= 0 ) error = .true.  ! Lower
         if ( ustat(i) /= 0 ) error = .true.  ! Normal
      enddo
      if ( error ) call flog( "Case 1: Status information was not as expected.", 1 )
   endif
!
!  Start CONOPT with second case where the right hand side in equation 3 is 11:
!
   casenum = 2
   coi_error = coi_solve( cntvect )
 
   write(*,*)
   write(*,*) 'End of Square6 example case 2. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "Case 2: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Case 2: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat < 4 .or. mstat > 5 ) then
      call flog( "Case 2: Solver or Model status not as expected (1,4) or (1,5)", 1 )
   elseif ( obj /= 0.d0 ) Then
      call flog( "Case 2: Objective for square model was not as expected 0.0", 1 )
   else
!
!  Check the primal and dual solution itself
!
      error = .false.
      do i = 1, 4
         if ( abs(xdual(i)) > 1.d-7 ) error = .true.
         if ( abs(udual(i)) > 1.d-7 ) error = .true.
      enddo
!  Model  x1 + x2           = 10  Check the activity levels based on primal values.
!         x1 - x2           = 0
!         x1 + x2 + x3      = 9
!            - x2 + x3 + x4 = 1
      if ( abs(xprim(1)+xprim(2)-uprim(1)) > 1.d-7 ) error = .true.
      if ( abs(xprim(1)-xprim(2)-uprim(2)) > 1.d-7 ) error = .true.
      if ( abs(xprim(1)+xprim(2)+xprim(3)-uprim(3)) > 1.d-7 ) error = .true.
      if ( abs(-xprim(2)+xprim(3)+xprim(3)-uprim(4)) > 1.d-7 ) error = .true.
      if ( error ) call flog( "Case 2: Numerical solution was not as expected.", 1 )
!
!  Check the status information
!
      do i = 1, 4
         if ( xbasc(i) /= 2 ) error = .true.  ! Basic
         if ( xstat(i) /= 0 ) error = .true.  ! Normal
!         if ( Ubasc(i) /= 0 ) error = .true.  ! Lower
      enddo
      if ( ustat(1) /= 2 .and. ustat(3) /= 2 .and. ustat(4) /= 2 ) error = .true. ! One of them must be infeasible
!      if ( Ubasc(3) /= 1 ) error = .true.  ! Upper
!      if ( Ubasc(4) /= 0 ) error = .true.  ! Lower
      if ( error ) call flog( "Case 2: Status information was not as expected.", 1 )
   endif
 
   if ( coi_free(cntvect) /= 0 ) call flog( "Error while freeing control vector",1)
 
   call flog( "Successful Solve", 0 )
 
End Program square
!
! ============================================================================
!  Define information about the model:
!
 
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function sq_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 :: Sq_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.
!
!  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.
   if ( casenum == 2 ) then
      lower(3) = 0.d0
      upper(4) = 0.d0
   endif
!
!  Constraint 1
!  Rhs = 10 and type Equal
!
   rhs(1)  = 10.d0
   type(1) = 0
!
!  Constraint 2
!  Rhs = 0 and type Equality
!
   type(2) = 0
!
!  Constraint 3
!  Rhs = 9 and type Equality
!
   rhs(3)  = 9.d0
   type(3) = 0
!
!  Constraint 4
!  Rhs = 1 and type Equality
!
   rhs(4)  = 1.d0
   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
!       x(1)  x(2)   x(3)  x(4)
!  1:    1     4
!  2:    2     5
!  3:    3     6      8
!  4:          7      9     10
!
   colsta(1) =  1
   colsta(2) =  4
   colsta(3) =  8
   colsta(4) = 10
   colsta(5) = 11
   rowno(1)  =  1
   rowno(2)  =  2
   rowno(3)  =  3
   rowno(4)  =  1
   rowno(5)  =  2
   rowno(6)  =  3
   rowno(7)  =  4
   rowno(8)  =  3
   rowno(9)  =  4
   rowno(10) =  4
!
!  Nonlinearity Structure: L = 0 are linear and NL = 1 are nonlinear
!       x(1)  x(2)  x(3)  x(4)
!  1:    L     L
!  2:    NL    NL
!  3:    NL    NL    L
!  4:          NL    NL    L
!
   nlflag(1) =  0
   nlflag(2) =  1
   nlflag(3) =  1
   nlflag(4) =  0
   nlflag(5) =  1
   nlflag(6) =  1
   nlflag(7) =  1
   nlflag(8) =  0
   nlflag(9) =  1
   nlflag(10)=  0
!
!  Value (Linear only)
!       x(1)  x(2)  x(3)  x(4)
!  1:    1     1
!  2:    NL    NL
!  3:    NL    NL    1
!  4:          NL    NL    1
!
   value(1)  = 1.d0
   value(4)  = 1.d0
   value(8)  = 1.d0
   value(10) = 1.d0
 
  sq_readmatrix = 0 ! Return value means OK
 
end Function sq_readmatrix
!
!==========================================================================
!  Compute nonlinear terms and non-constant Jacobian elements
!
 
!>  Compute nonlinear terms and non-constant Jacobian elements
!!
!! @include{doc} fdeval_params.dox
Integer Function sq_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 :: Sq_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: Is declared as linear and should not be called.
!
   if ( rowno .eq. 1 ) then
      sq_fdeval = 1
      return
!
!  Row 2: x1 + x2 assumed to be nonlinear
!
   elseif ( rowno .eq. 2 ) then
!
!  Mode = 1 or 3: Function value
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g = x(1) - x(2)
      endif
!
!  Mode = 2 or 3: Derivatives
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         jac(1) =  1.d0
         jac(2) = -1.d0
      endif
   elseif ( rowno .eq. 3 ) then
!
!  Mode = 1 or 3: Function value
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g = x(1) + x(2)
      endif
!
!  Mode = 2 or 3: Derivatives
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         jac(1) =  1.d0
         jac(2) =  1.d0
      endif
   elseif ( rowno .eq. 4 ) then
!
!  Mode = 1 or 3: Function value
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g = -x(2) + x(3)
      endif
!
!  Mode = 2 or 3: Derivatives
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         jac(2) =  -1.d0
         jac(3) =  1.d0
      endif
   endif
   sq_fdeval = 0
 
end Function sq_fdeval