qp1.f90

Go to the documentation of this file.

!> @file qp1.f90
!! @ingroup FORT1THREAD_EXAMPLES
!! @brief The current model is a simple QP model with a sparse Q matrix,
!! bounded variables, and one constraint. 
!!
!! The number of superbasic variables is larger than the default limit of 500
!! and the model does not solve nicely. To solve it faster, there are these
!! possibilities:
!! 1. Increase the limit on the number of superbasics.
!!    This can be done with a call to coidef_maxsup() or with an option.
!!    Both approaches are shown in extra solves in this file.
!! 2. Use directional 2nd derivatives, see qp2
!! 3. Use 2nd derivatives as a matrix, see qp3
!! 4. 2 and 3 combined, see qp4
!! 5. Use 2nd derivatives computed using perturbations, see qp5
!! 6. Use 2nd derivatives computed using perturbations defined using
!!    an options routine, see qp6
!! 7. As 5 but with the objective defined as a positive variable, i.e.
!!    it cannot be removed in the post-triangle.
!!
!! For more information about the individual callbacks, please have a look at the source code.
 
 
!> Main program. A simple setup and call of CONOPT for a QP model
!!
Program qp
 
   Use proginfo
   Use coidef
   Use qp_data
   implicit None
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: qp_readmatrix  ! Mandatory Matrix definition routine defined below
   Integer, External :: qp_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
#if defined(itl)
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_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 :: numcallback
   INTEGER, Dimension(:), Pointer :: cntvect
   INTEGER :: coi_error
!
!  Local variables used to define Q
!
   Integer :: i,j
   real*8  :: obj_1 ! Objective from first solve
 
   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:
!  We will create a QP model with one constraint, sum(i, x(i) ) = 1
!  and NN variables. NN and the Q matrix are declared in Module QP_Data.
!
   j = 1
   do i = 1, nn
      q(j)    = 1.d0
      qrow(j) = i
      qcol(j) = i
      j       = j + 1
      if ( i < nn ) then
         q(j)    = 0.1d0
         qrow(j) = i+1
         qcol(j) = i
         j       = j + 1
      endif
   enddo
   do i = 1, nn
      target(i) = 10.d0
   enddo
 
   coi_error = max( coi_error, coidef_numvar( cntvect, nn ) )   ! NN variables
   coi_error = max( coi_error, coidef_numcon( cntvect, 2 ) )    ! 1 constraint + 1 objective
   coi_error = max( coi_error, coidef_numnz( cntvect, 2*nn ) ) ! 2*NN nonzeros in the Jacobian
   coi_error = max( coi_error, coidef_numnlnz( cntvect, nn ) )   ! NN of which are nonlinear
   coi_error = max( coi_error, coidef_optdir( cntvect, -1 ) )   ! Minimize
   coi_error = max( coi_error, coidef_objcon( cntvect, 1 ) )    ! Objective is constraint 1
!
!  Tell CONOPT about the callback routines:
!
   coi_error = max( coi_error, coidef_readmatrix( cntvect, qp_readmatrix ) )
   coi_error = max( coi_error, coidef_fdeval( cntvect, qp_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     ) )
!
!  Default options file for debugging. Will be overwritten later.
!
   coi_error = max( coi_error, coidef_optfile( cntvect, 'qp1.opt'  ) )
 
#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 first Solve due to setup errors", 1 )
   ENDIF
!
!  Save the solution so we can check the duals:
!
   do_allocate = .true.
!
!  Start CONOPT:
!
   coi_error = coi_solve( cntvect )
 
   write(*,*)
   write(*,*) 'End of first solve in QP example 1. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "Errors encountered during first solve", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Status or Solution routine was not called during first solve", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Solver and Model Status in first call was not as expected (1,2)", 1 )
   Else
      Call checkdual( 'Solve 1', minimize )
   endif
   obj_1 = obj ! Remember the objective to test if next calls give the same
   stacalls = 0; solcalls = 0;
 
!
!  Define a limit on the number of superbasics of NN
!
   coi_error = max( coi_error, coidef_maxsup( cntvect, nn           ) )
   If ( coi_error .ne. 0 ) THEN
      call flog( "Skipping second Solve due to setup errors", 1 )
   ENDIF
   coi_error = coi_solve( cntvect )
   write(*,*)
   write(*,*) 'End of second solve in QP example 1. Return code=',coi_error
   If ( coi_error /= 0 ) then
      call flog( "Errors encountered during second solve", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Status or Solution routine was not called during second solve", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Solver and Model Status in second call was not as expected (1,2)", 1 )
   elseif ( abs( obj-obj_1 ) > 1.d-8*obj_1 ) then
      call flog( "Objective returned from second call was not the same as from first call", 1 )
   Else
      Call checkdual( 'Solve 2', minimize )
   endif
   stacalls = 0; solcalls = 0;
!
!  Add an options file with extra superbasics and run again after first
!  resetting MaxSup to a low number (10).
!
   coi_error = max( coi_error, coidef_optfile( cntvect, 'qp.opt'  ) )
!
!  Make sure the options file has the right content by writing it before it is used.
!
   Open(20, file='qp.opt', action='write');
   write(20,*) 'Lfnsup := 1000'
   Close(20)
   coi_error = max( coi_error, coidef_maxsup( cntvect, 10           ) )
   If ( coi_error .ne. 0 ) THEN
      call flog( "Skipping third Solve due to setup errors", 1 )
   ENDIF
   coi_error = coi_solve( cntvect )
   write(*,*)
   write(*,*) 'End of third solve in QP example 1. Return code=',coi_error
   If ( coi_error /= 0 ) then
      call flog( "Errors encountered during third solve", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Status or Solution routine was not called during third solve", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Solver and Model Status in the third call was not as expected (1,2)", 1 )
   elseif ( abs( obj-obj_1 ) > 1.d-8*obj_1 ) then
      call flog( "Objective returned from third call was not the same as from first call", 1 )
   Else
      Call checkdual( 'Solve 3', minimize )
   endif
 
   if ( coi_free(cntvect) /= 0 ) call flog( "Error while freeing control vector",1)
 
   call flog( "Successful Solve", 0 )
 
End Program qp
 
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function qp_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 :: QP_ReadMatrix
#endif
   Use qp_data
   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
   integer :: i, j
!
!  Information about Variables:
!  Default: Lower = -Inf, Curr = 0, and Upper = +inf.
!  Default: the status information in Vsta is not used.
!
!  Nondefault: Lower bound = 0 for all variables
!
   do i = 1, nn
     lower(i) = 0.0d0
   enddo
!
!  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.0 and type Non binding
!
   type(1) = 3
!
!  Constraint 2 (Sum of all vars = 1)
!  Rhs = 1 and type Equality
!
   rhs(2)  = 1.d0
   type(2) = 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)
!
   j = 1 ! counter for current nonzero
   do i = 1, nn
      colsta(i) = j
      rowno(j)  = 1
      nlflag(j) = 1
      j         = j + 1
      rowno(j)  = 2
      value(j)  = 1.d0
      nlflag(j) = 0
      j         = j + 1
   enddo
   colsta(nn+1) = j
 
  qp_readmatrix = 0 ! Return value means OK
 
end Function qp_readmatrix
 
!>  Compute nonlinear terms and non-constant Jacobian elements
!!
!! @include{doc} fdeval_params.dox
Integer Function qp_fdeval( x, g, jac, rowno, jcnm, mode, ignerr, errcnt, &
                            n, nz, thread, usrmem )
#if defined(itl)
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_FDEval
#endif
   Use qp_data
   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
   Integer :: i, j, k
   real*8  :: sum
!
!  Row 1: The objective function is nonlinear
!
   if ( rowno .eq. 1 ) then
!
!  Mode = 1 or 3. Function value: G = x * Q * x / 2
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         sum = 0.d0
         do k = 1, nq
            i = qrow(k)
            j = qcol(k)
            if ( i == j ) then
               sum = sum + (x(i)-target(i))*q(k)*(x(j)-target(j))
            else
               sum = sum + 2*(x(i)-target(i))*q(k)*(x(j)-target(j))
            endif
         enddo
         g    = sum / 2.d0
      endif
!
!  Mode = 2 or 3: Derivative values: Q*x
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         do i = 1, nn
            jac(i) = 0
         enddo
         do k = 1, nq
            i = qrow(k)
            j = qcol(k)
            if ( i == j ) then
               jac(i) = jac(i) + q(k) * (x(i)-target(i))
            else
               jac(i) = jac(i) + q(k) * (x(j)-target(j))
               jac(j) = jac(j) + q(k) * (x(i)-target(i))
            endif
         enddo
      endif
!
!  Row = 2: The row is linear and will not be called.
!
   endif
   qp_fdeval = 0
 
end Function qp_fdeval