stress1.f90

Go to the documentation of this file.

!> @file stress1.f90
!! @ingroup FORT1THREAD_EXAMPLES
!!
!! 
!! This is a CONOPT implementation of the gams model
!! 
!! @verbatim
!! set i / i1*i10000/
!! variable x(i), obj;
!! equation objdef; objdef .. obj =E= prod(i, x(i)*x(i) );
!! x.lo(i) = 1;
!! model m / all /;
!! solve m using nlp minimizing obj;
!! @endverbatim
!! 
!!
!! 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
 
Module stress1data
   integer, Parameter :: N = 15000
End Module stress1data
 
!> Main program. A simple setup and call of CONOPT
!!
Program stress1
   Use proginfo
   Use conopt
   Use stress1data
   Implicit none
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: stress_readmatrix ! Mandatory Matrix definition routine defined below
   Integer, External :: stress_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
   Integer, External :: stress_2dlagrstr  ! Hessian structure evaluation routine
   Integer, External :: stress_2dlagrval  ! Hessian value evaluation routine
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Stress_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Stress_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
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Stress_2DLagrStr
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Stress_2DLagrVal
#endif
!
!  Control vector
!
   INTEGER, Dimension(:), Pointer :: cntvect
   INTEGER :: coi_error
 
   call startup
!
!  Create and initialize a Control Vector
!
   coi_error = coi_create( cntvect )
!
!  Number of variables (excl. slacks): N
!
   coi_error = max( coi_error, coidef_numvar( cntvect, n ) )
!
!  Number of equations: 1 objective
!
   coi_error = max( coi_error, coidef_numcon( cntvect, 1 ) )
!
!  Number of nonzeros in the Jacobian: N
!
   coi_error = max( coi_error, coidef_numnz( cntvect, n ) )
!
!  Number of nonlinear nonzeros. N
!
   coi_error = max( coi_error, coidef_numnlnz( cntvect, n ) )
!
!  Number of nonzeros in the Hessian. N*(N+1)/2
!
   coi_error = max( coi_error, coidef_numhess( cntvect, n*(n+1)/2 ) )
!
!  Direction: -1 = minimization.
!
   coi_error = max( coi_error, coidef_optdir( cntvect, -1 ) )
!
!  Objective: Constraint no 1
!
   coi_error = max( coi_error, coidef_objcon( cntvect, 1 ) )
!
!  Turn Hessian debugging off
!
   coi_error = max( coi_error, coidef_debug2d( cntvect, 0 ) )
!
!  Options file = stress1.opt
!
   coi_error = max( coi_error, coidef_optfile( cntvect, 'stress1.opt' ) )
!
!  Tell CONOPT about the callback routines:
!
   coi_error = max( coi_error, coidef_readmatrix( cntvect, stress_readmatrix ) )
   coi_error = max( coi_error, coidef_fdeval( cntvect, stress_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_2dlagrstr( cntvect, stress_2dlagrstr  ) )
   coi_error = max( coi_error, coidef_2dlagrval( cntvect, stress_2dlagrval  ) )
 
#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 First Solve due to setup errors", 1 )
   ENDIF
!
!  Start CONOPT:
!
   coi_error = coi_solve( cntvect )
   If ( coi_error .ne. 0 ) THEN
      write(*,*)
      write(*,*) '**** Fatal Error while Solving model. Return Code=',coi_error
      write(*,*)
      call flog( "Errors encountered during Solve.", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Status or Solution routine was not called during Solve", 1 )
   else if ( mstat /= 2 .or. sstat /= 1 ) then
      call flog( "Incorrect Model or Solver Status during solve", 1 )
   ENDIF
 
   write(*,*)
   write(*,*) 'End of Stress1 Model. Return code=',coi_error
 
   if ( coi_free(cntvect) /= 0 ) call flog( "Error while freeing control vector",1)
 
   call flog( "Successful Solve", 0 )
 
end Program stress1
!
! =====================================================================
!  Define information about the model structure
!
 
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function stress_readmatrix( lower, curr, upper, vsta, type, rhs, esta,      &
                                 colsta, rowno, value, nlflag, nc, mc, nz, usrmem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Stress_ReadMatrix
#endif
   Use stress1data
   implicit none
   integer, intent (in)                      :: nc     ! number of variables
   integer, intent (in)                      :: mc     ! number of constraints
   integer, intent (in)                      :: nz     ! number of nonzeros
   real*8,  intent (in out), dimension(nc)   :: lower  ! vector of lower bounds
   real*8,  intent (in out), dimension(nc)   :: curr   ! vector of initial values
   real*8,  intent (in out), dimension(nc)   :: upper  ! vector of upper bounds
   integer, intent (in out), dimension(nc)   :: vsta   ! vector of initial variable status
                                                       ! Defined during restarts
   integer, intent (out),    dimension(mc)   :: type   ! vector of equation types
   integer, intent (in out), dimension(mc)   :: esta   ! vector of initial equation status
                                                       ! Defined during restarts
   real*8,  intent (in out), dimension(mc)   :: rhs    ! vector of right hand sides
   integer, intent (in out), dimension(nc+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
!
!  Define the information for the columns.
!
!  We should not supply status information, vsta.
!
   do i = 1, n
      lower(i) = 1.d0
      curr(i)  = 1.d0
   enddo
!
!  Define the information for the rows
!
!  We order the constraints as follows: The objective is first,
!  followed by tddef, sdef, csdef, ddef, rdef, and revdef for
!  the first period, the same for the second period, etc.
!
!  The objective is a nonbinding constraint:
!
   type(1) = 3
!
!  Define the structure and content of the Jacobian:
!
   do i = 1, n
      colsta(i) = i
      rowno(i)  = 1
      nlflag(i) = 1
   enddo
   colsta(n+1) = n+1
 
   stress_readmatrix = 0
 
end Function stress_readmatrix
!
! =====================================================================
!  Compute nonlinear terms and non-constant Jacobian elements
!
 
!>  Compute nonlinear terms and non-constant Jacobian elements
!!
!! @include{doc} fdeval_params.dox
Integer Function stress_fdeval( x, g, jac, rowno, jcnm, mode, ignerr, errcnt, &
                             nc, nz, thread, usrmem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Stress_FDEval
#endif
   Use stress1data
   implicit none
   integer, intent (in)                      :: nc     ! 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(nc)   :: x      ! vector of current solution values
   real*8,  intent (in out)                  :: g      ! constraint value
   real*8,  intent (in out), dimension(nc)   :: 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
   real*8  :: obj
!
   if ( rowno == 1 ) then
      obj = 1.0d0
      do i = 1, n
         obj = obj * x(i)*x(i)
      enddo
      if ( mode == 1 .or. mode == 3 ) then
         g = obj
      endif
      if ( mode == 2 .or. mode == 3 ) then
         do i = 1, n
            jac(i) = 2.d0*obj/x(i)
         enddo
      endif
   else
!
!  Error - this equation is not nonlinear
!
      write(*,*)
      write(*,*) 'Error. FDEval called with rowno=',rowno
      write(*,*)
      stress_fdeval = 1
      Return
   endif
   stress_fdeval = 0
 
end Function stress_fdeval
 
 
!>  Specify the structure of the Lagrangian of the Hessian
!!
!! @include{doc} 2DLagrStr_params.dox
Integer Function stress_2dlagrstr( HSRW, HSCL, NODRV, Nc, Mc, NHESS, UsrMem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Stress_2DLagrStr
#endif
   USE stress1data
   Implicit None
 
   Integer, Intent (IN)                    :: nc, mc, nhess
   Integer, Intent (IN OUT)                :: nodrv
   Integer, Dimension(NHess), Intent (Out) :: hsrw, hscl
   real*8,  Intent(IN OUT)                 :: usrmem(*)
 
   Integer :: i, j, iz
!
!  Define structure of Hessian
!
   if ( n /= nc ) then
      write(*,*) 'Struct: Number of variables should be=',n,' but CONOPT uses',nc
      stress_2dlagrstr = 1; Return
      stop
   endif
   if ( mc /= 1 ) then
      write(*,*) 'Struct: Number of equations should be=',+1,' but CONOPT uses',mc
      stress_2dlagrstr = 1; Return
   endif
   if ( nhess /= n*(n+1)/2 ) then
      write(*,*) 'Struct: Number of Hessian elements should be=',n*(n+1)/2,' but CONOPT uses',nhess
      stress_2dlagrstr = 1; Return
   endif
   iz = 0
   Do i = 1, n
      Do j = i, n
         iz = iz + 1
         hscl(iz) = i
         hsrw(iz) = j
      Enddo
   Enddo
 
   stress_2dlagrstr = 0
 
End Function stress_2dlagrstr
 
 
!>  Compute the Lagrangian of the Hessian
!!
!! @include{doc} 2DLagrVal_params.dox
Integer Function stress_2dlagrval( X, U, HSRW, HSCL, HSVL, NODRV, Nc, Mc, NHESS, UsrMem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Stress_2DLagrVal
#endif
   USE stress1data
   Implicit None
 
   Integer, Intent (IN)                    :: nc, mc, nhess
   Integer, Intent (IN OUT)                :: nodrv
   real*8,  Dimension(Nc), Intent (IN)     :: x
   real*8,  Dimension(Mc), Intent (IN)     :: u
   Integer, Dimension(Nhess), Intent (IN)  :: hsrw, hscl
   real*8,  Dimension(NHess), Intent (Out) :: hsvl
   real*8,  Intent(IN OUT)                 :: usrmem(*)
 
   Integer :: i, j, iz
!
!  Normal Evaluation mode -- we are testing memory only so we return a unit matrix.
!
   iz = 0
   Do i = 1, n
      iz = iz + 1
      hsvl(iz) = 1.0d0
      Do j = i+1, n
         iz = iz + 1
         hsvl(iz) = 0.0d0
      Enddo
   Enddo
 
   stress_2dlagrval = 0
 
End Function stress_2dlagrval