minimax07.f90

Go to the documentation of this file.

!> @file minimax07.f90
!! @ingroup FORT1THREAD_EXAMPLES
!!
!! 
!! Minimax07 is a minimax model with multiple minimax variables
!! and some side constraints.
!! 
!! Compared to Minimax06 we have added higher order terms on z(k)
!! so the limits on res are
!! \f{eqnarray*}{
!! k = 1: && z_1 \\
!! k = 2: && z_1^2 + z_2 \\
!! k = 3: && z_2^2 + z_3
!! \f}
!! 
!! We solve the following nonlinear minimax model:
!! 
!! \f[
!! \min \sum_k \max_{i} |res_{ik}| \\
!! \sum_{j} (a_{ijk}x_{jk} + b_{ijk}x_{jk}^2) + res_{ik} = obs_{ik} \quad \forall i, k \\
!! \sum_{j} x_{jk} \leq 1 \quad \forall k
!! \f]
!! 
!! where \f$a\f$, \f$b\f$, and \f$obs\f$ are known data, and \f$res\f$ and \f$x\f$ are the
!! variables of the model.
!! 
!! To get a model with continuous derivatives we introduce an
!! objective variable, \f$z\f$, and the constraints
!! \f[
!! res_{ik} - z_k \leq 0 \\
!! -res_{ik} - z_k \leq 0
!! \f]
!! Note opposite sign from Minimax01
!! and minimize  \f$\sum_k z_k\f$
!! 
!!
!! 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
 
REAL FUNCTION rndx( )
!
!  Defines a pseudo random number between 0 and 1
!
   IMPLICIT NONE
 
   Integer, save :: seed = 12359
 
   seed = mod(seed*1027+25,1048576)
   rndx = float(seed)/float(1048576)
 
END FUNCTION rndx
 
subroutine defdata
!
!  Define values for A, B, and Obs
!
   Use lsq_7k
   IMPLICIT NONE
 
   Integer :: i, j, k
   real*8, Parameter ::  xtarg = -1.0
   real*8, Parameter ::  noise = 1.0
   Real, External :: Rndx
   real*8 :: o
 
   do i = 1, nobs
      o = 0.d0
      do k = 1, dimk
         do j = 1, dimx
            a(i,k,j) = rndx()
            b(i,k,j) = rndx()
            o = o + a(i,k,j) * xtarg + b(i,k,j) * xtarg**2
         enddo
         obs(i,k) = o + noise * rndx()
      enddo
   enddo
 
end subroutine defdata
!
!  Main program.
!
!> Main program. A simple setup and call of CONOPT
!!
Program minimax07
 
   Use proginfo
   Use conopt
   Use lsq_7k
   implicit None
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: mm_readmatrix ! Mandatory Matrix definition routine defined below
   Integer, External :: mm_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 :: MM_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: MM_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
 
   call startup
!
!  Define data
!
   Call defdata
!
!  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, nobs*dimk + dimx*dimk + dimk ) )
   coi_error = max( coi_error, coidef_numcon( cntvect, 3*nobs*dimk+dimk+1  ) )
   coi_error = max( coi_error, coidef_numnz( cntvect, nobs * dimx * dimk + 5* nobs*dimk + dimk + dimk*dimx + 4*nobs ) )
   coi_error = max( coi_error, coidef_numnlnz( cntvect, nobs * dimx * dimk + 4*nobs ) )
   coi_error = max( coi_error, coidef_optdir( cntvect, -1 ) ) !  Minimize
   coi_error = max( coi_error, coidef_objcon( cntvect, 3*nobs*dimk + dimk+ 1 ) )  ! Objective is last constraint
   coi_error = max( coi_error, coidef_optfile( cntvect, 'Minimax07.opt'  ) )
!
!  Tell CONOPT about the callback routines:
!
   coi_error = max( coi_error, coidef_readmatrix( cntvect, mm_readmatrix ) )
   coi_error = max( coi_error, coidef_fdeval( cntvect, mm_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_debugfv( cntvect, 0 ) )  ! Debug Fdeval on or off
 
#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.
!
!  Start CONOPT:
!
   coi_error = coi_solve( cntvect )
 
   write(*,*)
   write(*,*) 'End of Minimax07 example 1. Return code=',coi_error
 
   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 ( .not. ( sstat == 1 .and. mstat == 2 ) ) then
      call flog( "Solver or Model status was not as expected (1,2)", 1 )
!   elseif ( abs( OBJ - 19.44434311d0 ) > 1.d-7 ) then
!      call flog( "Incorrect objective returned", 1 )
   Else
      Call checkdual( 'Minimax07', minimize )
   endif
 
   if ( coi_free(cntvect) /= 0 ) call flog( "Error while freeing control vector",1)
 
   call flog( "Successful Solve", 0 )
 
End Program minimax07
!
! ============================================================================
!  Define information about the model:
!
 
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function mm_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 :: MM_ReadMatrix
#endif
   Use lsq_7k
   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, k, l, nzc
!
!  Information about Variables:
!  Default: Lower = -Inf, Curr = 0, and Upper = +inf.
!  Default: the status information in Vsta is not used.
!
   l = 0
   do i = 1, dimx
      do k = 1, dimk
         l = l + 1
         curr(l) =  +0.9d0 ! Make the initial point infeasible to the sum(j,(x(j,k)) =L= 1;
         lower(l) = -1.0d0
      enddo
   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.
!
!  Constraints 1 to Nobs*Dimk:
!  Rhs = Obs(i,k) and type Equality
!
   l = 0
   do i = 1, nobs
      do k = 1, dimk
         l = l + 1
         rhs(l)  = obs(i,k)
         type(l) = 0
      enddo
   enddo
!
!  Constraint Nobs*Dimk+1 to 2*Nobs*Dimk
!  Rhs = 0 (default) and type Less than or equal
!
   do i = nobs*dimk+1, 2*nobs*dimk
      type(i) = 2
   enddo
!
!  Constraint 2*Nobs*Dimk+1 to 3*Nobs*Dimk
!  Rhs = 0 (default) and type Greater than or equal
!
   do i = 2*nobs*dimk+1, 3*nobs*dimk
      type(i) = 1
   enddo
!
!  Constraint 3*Nobs*Dimk+1 to 3*Nobs*Dimk+Dimk
!  Rhs = 1 and type Less than or equal
!
   do i = 3*nobs*dimk+1, 3*nobs*dimk+dimk
      rhs(i)  = 1.0d0
      type(i) = 2
   enddo
!
!  Objective, type Nonbinding
!
   type(3*nobs*dimk+dimk+1) = 3
!
!  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(j,k)  res(i,k)   z(k)
!  i,k:    NL    L=1
!  i,k:          L=1       nl and L=-1
!  i,k:          L=1       nl and L=+1
!  k       L=1
!  obj                     L=+1
!
!  We map (i,k) -> k+DimK*(i-1)
!     and (j,k) -> k+DimK*(j-1)
!
   nzc = 1
   do j = 1, dimx  ! x(j,k)
      do k = 1, dimk
         l = k+dimk*(j-1)
         colsta(l) =  nzc
         do i = 1, nobs
            rowno(nzc)  = k+dimk*(i-1)
            nlflag(nzc) = 1
            nzc         = nzc + 1
         enddo
         rowno(nzc) = 3*nobs*dimk+k
         nlflag(nzc) = 0
         value(nzc)  = 1.d0
         nzc         = nzc + 1
      enddo
   enddo
   do i = 1, nobs  ! res(i,k)
      do k = 1, dimk
         colsta(dimx*dimk+k+dimk*(i-1)) = nzc
         rowno(nzc)  = k+dimk*(i-1)
         nlflag(nzc) = 0
         value(nzc)  = 1.d0
         nzc         = nzc + 1
         rowno(nzc)  = nobs*dimk+k+dimk*(i-1)
         nlflag(nzc) = 0
         value(nzc)  = 1.d0
         nzc         = nzc + 1
         rowno(nzc)  = 2*nobs*dimk+k+dimk*(i-1)
         nlflag(nzc) = 0
         value(nzc)  = 1.d0
         nzc         = nzc + 1
      enddo
   enddo
   do k = 1, dimk  ! z(k)
      colsta(dimx*dimk+nobs*dimk+k) = nzc
      do i = 1, nobs
         rowno(nzc)  = nobs*dimk+k+dimk*(i-1)
         nlflag(nzc) = 0
         value(nzc)  = -1.d0
         nzc         = nzc + 1
         rowno(nzc)  = 2*nobs*dimk+k+dimk*(i-1)
         nlflag(nzc) = 0
         value(nzc)  = +1.d0
         nzc         = nzc + 1
      enddo
      if ( k < dimk ) then ! also nonlinear terms for z(k-1)
         do i = 1, nobs
            rowno(nzc)  = nobs*dimk+k+1+dimk*(i-1)
            nlflag(nzc) = 1
            nzc         = nzc + 1
            rowno(nzc)  = 2*nobs*dimk+k+1+dimk*(i-1)
            nlflag(nzc) = 1
            nzc         = nzc + 1
         enddo
      endif
      rowno(nzc) = 3*nobs*dimk+dimk+1 ! Objective
      nlflag(nzc) = 0
      value(nzc)  = +1.d0
      nzc         = nzc + 1
   enddo
   colsta(dimx*dimk+nobs*dimk+dimk+1) = nzc
 
   mm_readmatrix = 0 ! Return value means OK
 
end Function mm_readmatrix
!
!==========================================================================
!  Compute nonlinear terms and non-constant Jacobian elements
!
 
!>  Compute nonlinear terms and non-constant Jacobian elements
!!
!! @include{doc} fdeval_params.dox
Integer Function mm_fdeval( x, g, jac, rowno, jcnm, mode, ignerr, errcnt, &
                             n, nzc, thread, usrmem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: MM_FDEval
#endif
   use lsq_7k
   implicit none
   integer, intent (in)                      :: n      ! number of variables
   integer, intent (in)                      :: rowno  ! number of the row to be evaluated
   integer, intent (in)                      :: nzc     ! number of nonzceros 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(nzc)   :: 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,l
   real*8  :: s
 
   mm_fdeval = 0
   if ( rowno .le. nobs*dimk ) then
!
!  We map (i,k): rowno -> k+DimK*(i-1) so
!  i = (Dimk+rowno-1))/Dimk
!  k = rowno-Dimk*(i-1)
!     and (j,k) -> k+DimK*(j-1)
!  Mode = 1 or 3: Function value - x-part only
!
      i = (dimk+rowno-1)/dimk
      k = rowno-dimk*(i-1)
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         s = 0.d0
         do j = 1, dimx
            l = k+dimk*(j-1)
            s = s + a(i,k,j)*x(l) + b(i,k,j)*x(l)**2
         enddo
         g = s
      endif
!
!  Mode = 2 or 3: Derivatives
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         do j = 1, dimx
            l = k+dimk*(j-1)
            jac(l) = a(i,k,j) + 2.d0*b(i,k,j)*x(l)
         enddo
      endif
   Else if ( rowno .le. 2*nobs*dimk ) then
      i = (dimk+rowno-nobs*dimk-1)/dimk
      k = rowno-nobs*dimk-dimk*(i-1)-1
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g = x(nobs*dimk + dimx*dimk+k )**2
      endif
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         jac(nobs*dimk + dimx*dimk+k) = 2.0d0*x(nobs*dimk + dimx*dimk+k )
      endif
   Else if ( rowno .le. 3*nobs*dimk ) then
      i = (dimk+rowno-2*nobs*dimk-1)/dimk
      k = rowno-2*nobs*dimk-dimk*(i-1)-1
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g = x(nobs*dimk + dimx*dimk+k )**2
      endif
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         jac(nobs*dimk + dimx*dimk+k) = 2.0d0*x(nobs*dimk + dimx*dimk+k )
      endif
   Else
      mm_fdeval = 1 ! Illegal row number
   endif
 
end Function mm_fdeval