ident12.f90

Go to the documentation of this file.

!> @file ident12.f90
!! @ingroup FORT1THREAD_EXAMPLES
!!
!!
!!
!! This is a CONOPT implementation of the GAMS model:
!! @verbatim
!! Set i Rows / 1*m/
!! Set j cols / 1*n/
!! parameter a(j); a(j) = sqrt(ord(j))
!! parameter w(i); w(i) = power(-1,ord(i)+1)*sqrt(ord(i));
!! variable x(j), y(i), obj
!! positive variable x(j);
!! equation e(i), ey(i), objdef;
!! e(i) .. w(i)*sum(j,a(j)*x(j)) - y(i) =L= 0;
!! ey(i) .. y(i) =E= abs(w(i))+4*sign(w(i));
!! objdef .. obj =E= sum(j,x(j));
!! model m / all /; solve m using nlp minimizing obj;
!! @endverbatim
!!
!! The left hand sides of the constraints are all proportional but they are this time
!! inequalities. Two of them should be changed into a range-set and the rest should be
!! made redundant.
!!
!! The model is similar to Ident11 but we have added a constant to the right hand side
!! so all constraints are shifted and zero is no longer a solution. The selected
!! range constraint, 5, is also not binding, but constraint 2 is binding.
!! x(6) = (-sqrt(2)+4) / sqrt(2*6) = obj
!!
!!
!! 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 ident12
 
   Use proginfo
   Use conopt
   Use ident
   implicit None
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: ident_readmatrix ! Mandatory Matrix definition routine defined below
   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 :: Ident_ReadMatrix
!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
!
!  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, ncol+nrow ) )
   coi_error = max( coi_error, coidef_numcon( cntvect, 2*nrow+1 ) )
   coi_error = max( coi_error, coidef_numnz( cntvect, (nrow+1)*ncol+2*nrow ) )
   coi_error = max( coi_error, coidef_numnlnz( cntvect, 0 ) )
   coi_error = max( coi_error, coidef_optdir( cntvect, -1 ) )        ! minimize
   coi_error = max( coi_error, coidef_objcon( cntvect, 2*nrow + 1 ) )  ! Objective is last constraint
   coi_error = max( coi_error, coidef_optfile( cntvect, 'ident12.opt'  ) )
!
!  Tell CONOPT about the callback routines:
!
   coi_error = max( coi_error, coidef_readmatrix( cntvect, ident_readmatrix ) )
   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.
!
!  Start CONOPT:
!
   coi_error = coi_solve( cntvect )
 
   write(*,*)
   write(*,*) 'End of Ident 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 == 1 ) ) then
      call flog( "Solver or Model status was not as expected (1,1)", 1 )
   elseif ( abs( obj - (4.0d0-sqrt(2.0d0))/sqrt(12.d0) ) > 1.d-7 ) then
      call flog( "Incorrect objective returned", 1 )
   Else
      Call checkdual( 'Ident12', minimize )
   endif
 
   if ( coi_free(cntvect) /= 0 ) call flog( "Error while freeing control vector",1)
 
   call flog( "Successful Solve", 0 )
!
!  Free solution memory
!
   call finalize
 
End Program ident12
!
! ============================================================================
!  Define information about the model:
!
 
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function ident_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 :: Ident_ReadMatrix
#endif
   Use ident
   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
 
   real*8, dimension(Nrow) :: w
   real*8, dimension(Ncol) :: a
   Integer :: i, j, k
!
!  Information about Variables:
!  Default: Lower = -Inf, Curr = 0, and Upper = +inf.
!  Default: the status information in Vsta is not used.
!
   do j = 1, ncol  ! Only bounds on the x-variables
      lower(j) = 0.0d0
      a(j)     = sqrt(1.0d0*j)
   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 Nrow:
!  Rhs = 0 and type Less than
!  Constraints Nrow+1 to 2*Nrow:
!  Rhs = abs(w(i)) +/- 4 and type Equality
!
   do i = 1, nrow
      w(i)    = (-1)**(i+1)*sqrt(1.0d0*i);
      type(i) = 2
      if ( w(i) < 0.0 ) then
         rhs(nrow+i)  = abs(w(i)) - 4.0d0
      else
         rhs(nrow+i)  = abs(w(i)) + 4.0d0
      endif
      type(nrow+i) = 0
   enddo
!
!  Constraint Nrow + 1 (Objective)
!  Rhs = 0 and type Non binding
!
   type(2*nrow+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)          y(i)
!  e(i):    L=w(i)*a(j)   -1
!  ey(i):                  1
!  obj:     L=1
!
   k = 1
   do j = 1, ncol
      colsta(j) =  k
      do i = 1, nrow
         rowno(k)  = i
         nlflag(k) = 0
         value(k)  = w(i)*a(j)
         k         = k + 1
      enddo
      rowno(k) = 2*nrow+1
      value(k) = 1.0d0
      k        = k + 1
   enddo
   do i = 1, nrow
      j = ncol+i
      colsta(j) = k
      rowno(k)  = i
      nlflag(k) = 0
      value(k)  = -1.0d0
      k         = k + 1
      rowno(k)  = nrow+i
      nlflag(k) = 0
      value(k)  = +1.0d0
      k         = k + 1
   enddo
   colsta(ncol+nrow+1) = k
 
   ident_readmatrix = 0 ! Return value means OK
 
end Function ident_readmatrix