ident02.f90

Go to the documentation of this file.

!> @file ident02.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), obj
!! positive variable x(j);
!! equation e(i), objdef;
!! e(i) .. w(i)*sum(j,a(j)*x(j)) =L= w(i);
!! objdef .. obj =E= sum(j,x(j));
!! model m / all /; solve m using nlp maximizing obj;
!! @endverbatim
!! 
!! The left hand sides of the constraints are all proportional but they are this time
!! inequalities. The last two (with the largest size) should be singled out and the last
!! changed into an equality and the rest should be made redundant.
!! 
!!
!! 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 ident02
 
   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 ) )
   coi_error = max( coi_error, coidef_numcon( cntvect, nrow+1 ) )
   coi_error = max( coi_error, coidef_numnz( cntvect, (nrow+1)*ncol ) )
   coi_error = max( coi_error, coidef_numnlnz( cntvect, 0 ) )
   coi_error = max( coi_error, coidef_optdir( cntvect, +1 ) )        ! maximize
   coi_error = max( coi_error, coidef_objcon( cntvect, nrow + 1 ) )  ! Objective is last constraint
   coi_error = max( coi_error, coidef_optfile( cntvect, 'ident02.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 - 1.0d0 ) > 1.d-7 ) then
      call flog( "Incorrect objective returned", 1 )
   Else
      Call checkdual( 'Ident02', maximize )
   endif
 
   if ( coi_free(cntvect) /= 0 ) call flog( "Error while freeing control vector",1)
 
   call flog( "Successful Solve", 0 )
 
End Program ident02
!
! ============================================================================
!  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
      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 = Obs(i) and type Equality
!
   do i = 1, nrow
      w(i)    = (-1)**(i+1)*sqrt(1.0d0*i);
      rhs(i)  = w(i)
      type(i) = 2
   enddo
!
!  Constraint Nrow + 1 (Objective)
!  Rhs = 0 and type Non binding
!
   type(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)
!  i:    L=w(i)*a(j)
!  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) = nrow+1
      value(k) = 1.0d0
      k         = k + 1
   enddo
   colsta(ncol+1) = k
 
   ident_readmatrix = 0 ! Return value means OK
 
end Function ident_readmatrix