cns02.f90

Go to the documentation of this file.

!> @file cns02.f90
!! @ingroup FORT1THREAD_EXAMPLES
!!
!! This is a CONOPT implementation of the GAMS model:
!! 
!! @verbatim
!! variable x, y;
!! equation f, g;
!! model cns02 / f, g /;
!! 
!! f .. x*x + .001*y =e= 4;
!! g .. x + y        =e= 8;
!! 
!! option limrow = 0, limcol = 0, decimals=8;
!! scalar x1, y1, x2, y2, det;
!! * 1000*f - g yields 1000x^2 - x - 3992 = 0
!! det = sqrt(1 + 4 * 1000 * 3992);
!! x1 = (1 + det)/ 2000;
!! x2 = (1 - det)/ 2000;
!! y1 = 8 - x1;
!! y2 = 8 - x2;
!! display x1, x2, y1, y2;
!! @endverbatim
!! 
!! 
!! *  Case 1: No bounds on the variables. The model should solve fine.
!! @verbatim
!! x.lo = -INF; x.up = INF; x.l = 8; solve cns02 using cns;
!! abort$(cns02.solvestat <> 1 or cns02.modelstat <> 16)  'bad return codes';
!! abort$((abs(x.l-x1) <= tol and abs(y.l-y1) <= tol)
!! eqv (abs(x.l-x2) <= tol and abs(y.l-y2) <= tol)) 'x or y is wrong',x.l,y.l;
!! @endverbatim
!! 
!! *  Case 2: No bounds on the variables. The model should again solve
!!           fine, but the solution can be different because the initial
!!           value of x is different.
!! @verbatim
!! x.lo = -INF; x.up = INF; x.l = -8; solve cns02 using cns;
!! abort$(cns02.solvestat <> 1 or cns02.modelstat <> 16)  'bad return codes';
!! abort$((abs(x.l-x1) <= tol and abs(y.l-y1) <= tol)
!! eqv (abs(x.l-x2) <= tol and abs(y.l-y2) <= tol)) 'x or y is wrong',x.l,y.l;
!! @endverbatim
!! 
!! *  Case 3: The bound on x will make solution unique (feasible).
!! @verbatim
!! x.lo = 1; x.up = 3; x.l = 8; solve cns02 using cns;
!! abort$(cns02.solvestat <> 1 or cns02.modelstat <> 16)  'bad return codes';
!! abort$(abs(x.l-x1) > tol or abs(y.l-y1) > tol) 'x or y is wrong',x.l,y.l;
!! @endverbatim
!! 
!! *  Case 4: The bound on x will make the model infeasible.
!! @verbatim
!! x.lo = 2; x.up = 3; x.l = 8; solve cns02 using cns;
!! abort$(cns02.solvestat <> 1 or cns02.modelstat <> 5)  'bad return codes';
!! abort$(cns02.numinfes < 1)                            'wrong .numinfes';
!! @endverbatim
!! 
!! *  Case 5 and 6: X is fixed and the model is not square any more.
!! @verbatim
!! cns02.holdfixed = 0; x.lo = 2; x.up = 2; x.l = 8; solve cns02 using cns;
!! abort$(execerror=0) 'previous solve should have given exec errors';
!! abort$(cns02.solvestat <> 12 or cns02.modelstat <> 14)  'bad return codes';
!! execerror = 0;
!! 
!! cns02.holdfixed = 1; x.lo = 2; x.up = 2; x.l = 8; solve cns02 using cns;
!! abort$(execerror=0) 'previous solve should have given exec errors';
!! abort$(cns02.solvestat <> 12 or cns02.modelstat <> 14)  'bad return codes';
!! execerror = 0;
!! @endverbatim
!! 
!!
!! 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 cns02
 
   Use proginfo
   Use conopt
   Use casedata_num
   implicit None
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: cns02_readmatrix ! Mandatory Matrix definition routine defined below
   Integer, External :: cns02_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 :: Cns02_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Cns02_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
!
!  Solution info
!
   real*8 x1, x2, y1, y2, det, tol
   logical ok
   det = sqrt(1.d0+4000.d0*3992d0)
   x1 = (1.d0+det)/2000.d0; y1 = 8.d0-x1
   x2 = (1.d0-det)/2000.d0; y2 = 8.d0-x2
   tol = 1.d-8
!
!  Create and initialize a Control Vector
!
   call startup
 
   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, 2 ) )  ! # variables
   coi_error = max( coi_error, coidef_numcon( cntvect, 2 ) )  ! # constraints
   coi_error = max( coi_error, coidef_numnz( cntvect, 4 ) )  ! # nonzeros in the Jacobian
   coi_error = max( coi_error, coidef_numnlnz( cntvect, 1 ) )  ! # of which are nonlinear
   coi_error = max( coi_error, coidef_square( cntvect, 1 ) )  ! Square system
   coi_error = max( coi_error, coidef_optfile( cntvect, 'cns02.opt'    ) )
!
!  Tell CONOPT about the callback routines:
!
   coi_error = max( coi_error, coidef_readmatrix( cntvect, cns02_readmatrix ) )
   coi_error = max( coi_error, coidef_fdeval( cntvect, cns02_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
!
!  Save the solution so we can check the duals:
!
   do_allocate = .true.
   write(10,*) 'Solution set 1: x1=',x1,' y1=',y1
   write(10,*) 'Solution set 2: x2=',x2,' y2=',y2
!
!  Start CONOPT:
!
   casenum = 1
   coi_error = coi_solve( cntvect )
 
   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 /= 16 ) then
      call flog( "Case 1: Solver and Model Status was not as expected (1,16)", 1 )
   else
      ok = ( abs(xprim(1)-x1) < tol .and. abs(xprim(2)-y1) < tol ) .or. &
           ( abs(xprim(1)-x2) < tol .and. abs(xprim(2)-y2) < tol )
      if ( .not. ok ) then
         write(10,*) 'Solution for case 1 was x=',xprim(1),' and y=',xprim(2)
         call flog( "Case 1: Solver values were not correct.", 1 )
      endif
   endif
 
 
   casenum = 2
   coi_error = coi_solve( cntvect )
 
   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 /= 16 ) then
      call flog( "Case 2: Solver and Model Status was not as expected (1,16)", 1 )
   else
      ok = ( abs(xprim(1)-x1) < tol .and. abs(xprim(2)-y1) < tol ) .or. &
           ( abs(xprim(1)-x2) < tol .and. abs(xprim(2)-y2) < tol )
      if ( .not. ok ) then
         write(10,*) 'Solution for case 2 was x=',xprim(1),' and y=',xprim(2)
         call flog( "Case 2: Solver values were not correct.", 1 )
      endif
   endif
 
   casenum = 3
   coi_error = coi_solve( cntvect )
 
   If ( coi_error /= 0 ) then
      call flog( "Case 3: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Case 3: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 16 ) then
      call flog( "Case 3: Solver and Model Status was not as expected (1,16)", 1 )
   else
      ok = ( abs(xprim(1)-x1) < tol .and. abs(xprim(2)-y1) < tol )
      if ( .not. ok ) then
         write(10,*) 'Solution for case 3 was x=',xprim(1),' and y=',xprim(2)
         call flog( "Case 3: Solver values were not correct.", 1 )
      endif
   endif
 
   casenum = 4
   coi_error = coi_solve( cntvect )
 
   If ( coi_error /= 0 ) then
      call flog( "Case 4: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Case 4: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 5 ) then
      call flog( "Case 4: Solver and Model Status was not as expected (1,5)", 1 )
   elseif ( c_infeas == 0 ) then
      call flog( "Case 4: Infeasibility count is zero.", 1 )
   endif
 
   write(*,*)
   write(*,*) 'End of Cns02 example. Return code=',coi_error
 
   if ( coi_free(cntvect) /= 0 ) call flog( "Error while freeing control vector",1)
 
   call flog( "Successful Solve", 0 )
 
End Program cns02
!
! ============================================================================
!  Define information about the model:
!
 
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function cns02_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 :: Cns02_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.
!
   if ( casenum == 1 ) then
      curr(1)  = 8.0d0
   else if ( casenum == 2 ) then
      curr(1)  = -8.0d0
   else if ( casenum == 3 ) then
      lower(1) = 1.0d0
      upper(1) = 3.0d0
      curr(1)  = 8.0d0
   else if ( casenum == 4 ) then
      lower(1) = 2.0d0
      upper(1) = 3.0d0
      curr(1)  = 8.0d0
   endif
!
!  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.
!
   type(1) = 0
   rhs(1)  = 4.0d0
   type(2) = 0
   rhs(2)  = 8.0d0
!
!  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)
!  1:    1     3
!  2:    2     4
!
   colsta(1) =  1
   colsta(2) =  3
   colsta(3) =  5
   rowno(1)  =  1
   rowno(2)  =  2
   rowno(3)  =  1
   rowno(4)  =  2
!
!  Nonlinearity Structure: L = 0 are linear and NL = 1 are nonlinear
!       x(1)  x(2)
!  1:    NL    L
!  2:    L     L
!
   nlflag(1) =  1
   nlflag(2) =  0
   nlflag(3) =  0
   nlflag(4) =  0
!
!  Value (Linear only)
!       x(1)  x(2)
!  1:   NL    0.001
!  2:    1     1
!
   value(2)  =  1.0d0
   value(3)  =  0.001d0
   value(4)  =  1.d0
 
  cns02_readmatrix = 0 ! Return value means OK
 
end Function cns02_readmatrix
!
!==========================================================================
!  Compute nonlinear terms and non-constant Jacobian elements
!
 
!>  Compute nonlinear terms and non-constant Jacobian elements
!!
!! @include{doc} fdeval_params.dox
Integer Function cns02_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 :: Cns02_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: the only nonlinear row: x*x
!
   if ( rowno .eq. 1 ) then
!
!  Mode = 1 or 3. Function value: G = P * Out
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g    = x(1)*x(1)
      endif
!
!  Mode = 2 or 3: Derivative values:
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         jac(1) = x(1)+x(1)
      endif
   endif
   cns02_fdeval = 0
 
end Function cns02_fdeval