mp_rosex.f90

Go to the documentation of this file.

!> @file mp_rosex.f90
!! @ingroup FORTOPENMP_EXAMPLES
!!
!! @copydoc rosex.f90
 
#if defined(_WIN32) && !defined(_WIN64)
#define dec_directives_win32
#endif
 
Module tlim
   Integer :: Thread2d
End module tlim
 
!> Main program. A simple setup and call of CONOPT
!!
Program rosex
 
   Use proginfop
   Use conopt
   Use omp_lib
   Use tlim
   IMPLICIT NONE
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: ros_readmatrix ! Mandatory Matrix definition routine defined below
   Integer, External :: ros_fdeval     ! Function and Derivative evaluation routine
                                       ! needed a nonlinear model.
   Integer, External :: std_status     ! Standard callback for displaying solution status
   Integer, External :: ros_solution   ! Special callback for displaying solution values that does not write the solution
   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 :: Ros_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Ros_FDEval
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Status
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Ros_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
!
!  The model is an unconstrained problem in which the objective is the
!  sum of N identical copies of the Rosenbrock function placed after
!  each other:
!
!  min sum(i, (1-x(2*i-1))**2+100*(x(2*i)-x(2*i-1)**2)**2)
!
!  The number of Rosenborck functions
!
   Integer :: nr
!
!  Thread Info
!
   Integer :: maxthread
   real*8 time0, time1, time4
!
!  Create and initialize a Control Vector
!
   call startup
 
   coi_error = coi_create( cntvect )
!
!  Tell CONOPT about the size of the model by populating the Control Vector:
!
   nr = 16000
   coi_error = max( coi_error, coidef_numvar( cntvect, 2*nr ) ) ! # variables
   coi_error = max( coi_error, coidef_numcon( cntvect, 1 ) )    ! # constraints
   coi_error = max( coi_error, coidef_numnz( cntvect, 2*nr ) ) ! # nonzeros in the Jacobian
   coi_error = max( coi_error, coidef_numnlnz( cntvect, 2*nr ) ) ! # of which are nonlinear
   coi_error = max( coi_error, coidef_optdir( cntvect, -1 ) )   ! Minimize
   coi_error = max( coi_error, coidef_objcon( cntvect, 1 ) )    ! Objective is constraint 1
   coi_error = max( coi_error, coidef_optfile( cntvect, 'mp_rosex.opt' ) )
!
!  Tell CONOPT about the callback routines:
!
   coi_error = max( coi_error, coidef_readmatrix( cntvect, ros_readmatrix ) )
   coi_error = max( coi_error, coidef_fdeval( cntvect, ros_fdeval     ) )
   coi_error = max( coi_error, coidef_status( cntvect, std_status     ) )
   coi_error = max( coi_error, coidef_solution( cntvect, ros_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
!
!  Start CONOPT with a single thread:
!
   thread2d = 1
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time1 = omp_get_wtime() - time0
   If ( coi_error /= 0 ) then
      call flog( "Solve 1: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Solve 1: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Solve 1: Solver and Model Status was not as expected (1,2)", 1 )
   elseif ( abs( obj-0.0d0 ) > 0.000001d0 ) then
      call flog( "Solve 1: Incorrect objective returned", 1 )
   endif
!
   maxthread = omp_get_max_threads() ! Define ThreadC as the largest number of threads we can get
   coi_error = max( coi_error, coidef_threadc( cntvect, maxthread ) )
   coi_error = max( coi_error, coidef_threadf( cntvect, 1 ) )  ! Functions must run with one thread
#if defined (complete)
!
!  Start CONOPT -- Multi-thread mode with a single thread:
!
   thread2d = 1
   coi_error = max( coi_error, coidef_threads( cntvect, 1 ) )  ! 1 means use 1 thread
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time1a = omp_get_wtime() - time0
   If ( coi_error /= 0 ) then
      call flog( "Solve 1a: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Solve 1a: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Solve 1a: Solver and Model Status was not as expected (1,2)", 1 )
   elseif ( abs( obj-0.0d0 ) > 0.000001d0 ) then
      call flog( "Solve 1a: Incorrect objective returned", 1 )
   endif
!
!  Start CONOPT -- Multi-thread mode with two threads:
!
   thread2d = 2
   coi_error = max( coi_error, coidef_threads( cntvect, 2 ) )  ! 2 threads
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time2 = omp_get_wtime() - time0
   If ( coi_error /= 0 ) then
      call flog( "Solve 2: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Solve 2: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Solve 2: Solver and Model Status was not as expected (1,2)", 1 )
   elseif ( abs( obj-0.0d0 ) > 0.000001d0 ) then
      call flog( "Solve 2: Incorrect objective returned", 1 )
   endif
#endif
   if ( maxthread >= 4 ) then
!
!  Start CONOPT -- Multi-thread mode with four threads:
!
   thread2d = 4
   coi_error = max( coi_error, coidef_threads( cntvect, 4 ) )  ! 4 threads
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time4 = omp_get_wtime() - time0
   If ( coi_error /= 0 ) then
      call flog( "Solve 4: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Solve 4: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Solve 4: Solver and Model Status was not as expected (1,2)", 1 )
   elseif ( abs( obj-0.0d0 ) > 0.000001d0 ) then
      call flog( "Solve 4: Incorrect objective returned", 1 )
   endif
   endif
 
   if ( coi_free( cntvect ) /= 0 ) call flog( "Error while freeing control vector", 1 )
 
   write(*,*)
   write(*,"('Time for single thread ',f10.3)") time1
#if defined (complete)
   write(*,"('Time for single thread ',f10.3)") time1a
   write(*,"('Time for dual  threads ',f10.3)") time2
#endif
   if ( maxthread >= 4 ) then
   write(*,"('Time for quad  threads ',f10.3)") time4
   endif
   write(*,*)
   write(*,*) 'End of Extended Rosenbrock Function example. Return code=',coi_error
 
   call flog( "Successful Solve", 0 )
 
End Program rosex
!
! ============================================================================
!  Define information about the model:
!
 
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function ros_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 :: Ros_ReadMatrix
#endif
      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 :: nr
   Integer :: i
 
   nr = n / 2
!
!  Information about Variables:
!  Default: Lower = -Inf, Curr = 0, and Upper = +inf.
!  Default: the status information in Vsta is not used.
!
   DO i = 1, nr
      curr(2*i-1) = -1.d0
      curr(2*i  ) = +1.d0
   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.
!
!  Constraint 1 (Objective)
!  Rhs = 0.0 and type Non binding
!
   type(1) = 3
!
!  Information about the Jacobian. We have to define 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(1)  x(2)  x(3)  x(4)
!  1:    1     2     3     4    etc
!
!  Nonlinearity Structure: L = 0 are linear and NL = 1 are nonlinear
!       x(1)  x(2)  x(3)  x(4)
!  1:    NL    NL    NL    NL  etc
!
!  Value (Linear only)
!
   DO i = 1, nr*2+1
      colsta(i) = i
   enddo
   DO i = 1, nr*2
      rowno(i)  = 1
      nlflag(i) = 1
   enddo
 
   ros_readmatrix = 0 ! Return value means OK
 
end Function ros_readmatrix
!
!==========================================================================
!  Compute nonlinear terms and non-constant Jacobian elements
!
 
!>  Compute nonlinear terms and non-constant Jacobian elements
!!
!! @include{doc} fdeval_params.dox
Integer Function ros_fdeval( x, g, jac, rowno, jcnm, mode, ignerr, errcnt, &
                             n, nz, thread, usrmem )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Ros_FDEval
#endif
   Use tlim
      IMPLICIT NONE
   integer, intent (in)                      :: n      ! 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(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(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 :: nr
   Integer :: i
 
   nr = n / 2
!
!  Row 1: the objective function is nonlinear
!
   if ( rowno .eq. 1 ) then
!
!  Mode = 1 or 3. Function value:
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g = 0.d0
         do i = 1, nr
            g = g + (x(2*i-1)-1.d0)**2 + 1.d2*(x(2*i)-x(2*i-1)**2)**2
         enddo
      endif
!
!  Mode = 2 or 3: Derivative values:
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
!$OMP PARALLEL DEFAULT(Shared) IF( Thread2D > 1 ) NUM_THREADS(Thread2D)
!$OMP DO SCHEDULE(Static)
         do i = 1, nr
            jac(2*i-1) = 2.d0*(x(2*i-1)-1.d0) - 4.d2*(x(2*i)-x(2*i-1)**2)*x(2*i-1)
            jac(2*i  ) = 2.d2*(x(2*i)-x(2*i-1)**2)
         enddo
!$OMP END DO
!$OMP END PARALLEL
      endif
   endif
   ros_fdeval = 0
 
end Function ros_fdeval
 
Integer Function ros_solution( XVAL, XMAR, XBAS, XSTA, YVAL, YMAR, YBAS, YSTA, N, M, USRMEM )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Ros_Solution
#endif
!
!  Simple implementation in which we write the solution values to
!  the 'Documentation file' on unit 10.
!
   Use proginfop
   IMPLICIT NONE
   INTEGER, Intent(IN)                  :: n, m
   INTEGER, Intent(IN),    Dimension(N) :: xbas, xsta
   INTEGER, Intent(IN),    Dimension(M) :: ybas, ysta
   real*8,  Intent(IN),    Dimension(N) :: xval, xmar
   real*8,  Intent(IN),    Dimension(M) :: yval, ymar
   real*8,  Intent(IN OUT)              :: usrmem(*)
 
 
   solcalls = solcalls + 1
   ros_solution = 0
 
END Function ros_solution