sqp3.f90

Go to the documentation of this file.

!> @file sqp3.f90
!! @ingroup FORT1THREAD_EXAMPLES
!!
!! 
!! SQP3: Sparse QP model with `2DLagrStr` and `2DLagrVal` routines
!! 
!!
!! 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 for a QP model
!!
Program qp
 
   Use proginfo
   Use conopt
   Use qp_data_sp
   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 :: qp_2dlagrstr   ! 2nd derivative routine for structure
   Integer, External :: qp_2dlagrval   ! 2nd derivative routine for values
   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 :: QP_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_FDEval
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_2DLagrStr
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_2DLagrVal
!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
!
!  Local variables used to define Q
!
   Integer :: i,j,k
 
   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:
!  We will create a QP model with one constraint, sum(i, x(i) ) = 1
!  and NN variables. NN and the Q matrix are defined in Module QP_Data_Sp.
!
   j = 0
   do i = 1, nn
      j       = j + 1
      q(j)    = 1.d0
      qrow(j) = i
      qcol(j) = i
      do k = i+1, nn, 3
         j       = j + 1
         q(j)    = q(j-1)*0.25 ! Make sure the matrix is diagonally dominant.
         qrow(j) = k
         qcol(j) = i
      enddo
   enddo
   nq = j
   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_numhess( cntvect, nq ) )   ! NQ elements in the Hessian
   coi_error = max( coi_error, coidef_optdir( cntvect, -1 ) )   ! Minimize
   coi_error = max( coi_error, coidef_objcon( cntvect, 1 ) )    ! Objective is constraint 1
   coi_error = max( coi_error, coidef_hessfac( cntvect, 1.d0 ) ) ! Only allow a very sparse Hessian in 2DLagr
   coi_error = max( coi_error, coidef_optfile( cntvect, 'sqp3.opt' ) )
!
!  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_2dlagrstr( cntvect, qp_2dlagrstr  ) )
   coi_error = max( coi_error, coidef_2dlagrval( cntvect, qp_2dlagrval  ) )
   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
!
!  Start CONOPT:
!
   coi_error = coi_solve( cntvect )
 
   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 )
   endif
!
!  Repeat but this time allow a dense Hessian
!
   coi_error = max( coi_error, coidef_hessfac( cntvect, 0.d0 ) ) ! Allow a dense Hessian in 2DLagr. 0.0 means no limit
   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
 
   coi_error = coi_solve( cntvect )
 
   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 )
   endif
!
!  Repeat once again but this time with a MaxHeap limit that should prevent using the Hessian
!
   coi_error = max( coi_error, coidef_maxheap( cntvect, 2.5d0 ) ) ! 2.5 MBytes is not enough for the Hessian
   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
 
   coi_error = coi_solve( cntvect )
 
   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 )
   endif
 
   write(*,*)
   write(*,*) 'End of SQP example 3. Return code=',coi_error
 
   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:
!
 
!>  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 )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_ReadMatrix
#endif
   Use qp_data_sp
   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. 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)
!
   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
!
 
!>  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 )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_FDEval
#endif
   Use qp_data_sp
   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)) / 2.d0
            else
               sum = sum + (x(i)-target(i))*q(k)*(x(j)-target(j))
            endif
         enddo
         g    = sum
      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.0d0
         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
 
 
!>  Specify the structure of the Lagrangian of the Hessian
!!
!! @include{doc} 2DLagrStr_params.dox
Integer Function qp_2dlagrstr( ROWNO, COLNO, NODRV, N, M, NHESS, USRMEM )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_2DLagrStr
#endif
   Use qp_data_sp
   implicit none
   INTEGER n, m, nhess, nodrv
   real*8  usrmem(*)
   INTEGER rowno(nhess), colno(nhess)
 
   Integer :: j
 
   Do j = 1, nq
      rowno(j) = qrow(j)
      colno(j) = qcol(j)
   Enddo
   qp_2dlagrstr = 0
 
End function qp_2dlagrstr
 
 
!>  Compute the Lagrangian of the Hessian
!!
!! @include{doc} 2DLagrVal_params.dox
Integer Function qp_2dlagrval( X, U, ROWNO, COLNO, VALUE, NODRV, N, M, NHESS, USRMEM )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_2DLagrVal
#endif
   Use qp_data_sp
   implicit none
   INTEGER n, m, nhess, nodrv
   real*8  x(n), u(m), value(nhess), usrmem(*)
   INTEGER rowno(nhess), colno(nhess)
 
   Integer :: j
 
   Do j = 1, nq
      value(j) = q(j)*u(1)
   Enddo
   qp_2dlagrval = 0
 
End function qp_2dlagrval