minimax04.f90

Go to the documentation of this file.

!> @file minimax04.f90
!! @ingroup FORT1THREAD_EXAMPLES
!!
!! 
!! Minimax04 is a minimax model with multiple minimax variables.
!! 
!! 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
!! \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 minimax04
 
   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+1  ) )
   coi_error = max( coi_error, coidef_numnz( cntvect, nobs * dimx * dimk + 5* nobs*dimk + dimk ) )
   coi_error = max( coi_error, coidef_numnlnz( cntvect, nobs * dimx * dimk ) )
   coi_error = max( coi_error, coidef_optdir( cntvect, -1 ) ) !  Minimize
   coi_error = max( coi_error, coidef_objcon( cntvect, 3*nobs*dimk + 1 ) )  ! Objective is last constraint
   coi_error = max( coi_error, coidef_optfile( cntvect, 'Minimax04.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 Minimax04 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( 'Minimax04', minimize )
   endif
 
   if ( coi_free(cntvect) /= 0 ) call flog( "Error while freeing control vector",1)
 
   call flog( "Successful Solve", 0 )
 
End Program minimax04
!
! ============================================================================
!  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.8d0
      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
!
!  Objective, type Nonbinding
!
   type(3*nobs*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       L=-1
!  i,k:          L=1       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
      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
      rowno(nzc) = 3*nobs*dimk+1
      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
      mm_fdeval = 1 ! Illegal row number
   endif
 
end Function mm_fdeval