mp_lincns.f90

Go to the documentation of this file.

!> @file mp_lincns.f90
!! @ingroup FORTOPENMP_EXAMPLES
!!
!! 
!! Large Linear Dense CNS model.
!! Used to test the inversion routine for a very dense model.
!! 
!!
!! 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 lincns
   Use proginfop
   Use conopt
   Use omp_lib
   Implicit none
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: lincns_readmatrix ! Mandatory Matrix definition routine defined below
   Integer, External :: std_status     ! Standard callback for displaying solution status
   Integer, External :: lincns_solution   ! callback for displaying and testing solution values
   Integer, External :: std_message    ! Standard callback for managing messages
   Integer, External :: std_errmsg     ! Standard callback for managing error messages
   Integer, External :: std_triord     ! Standard callback for triangular order
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lincns_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Status
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lincns_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 :: Std_TriOrd
#endif
!
!  Control vector
!
   INTEGER, Dimension(:), Pointer :: cntvect
   INTEGER :: coi_error
!
!   Locals
!
   Integer :: n
!
!  Create and initialize a Control Vector
!
   real*8 time0, time1, time2
 
   call startup
   coi_error = coi_create( cntvect )
!
!  Define the number of variables and constraints
!
   n = 1000
!
!  Tell CONOPT about the size of the model by populating the Control Vector:
!
!  Number of variables = N
!
   coi_error = max( coi_error, coidef_numvar( cntvect, n ) )
!
!  Number of equations = M
!
   coi_error = max( coi_error, coidef_numcon( cntvect, n ) )
!
!  Number of nonzeros in the Jacobian: N*N
!
   coi_error = max( coi_error, coidef_numnz( cntvect, n*n ) )
!
!  Number of nonlinear nonzeros. 0 -- the model is linear
!
   coi_error = max( coi_error, coidef_numnlnz( cntvect, 0 ) )
!
!  Square system ( No objective is defined)
!
   coi_error = max( coi_error, coidef_square( cntvect, 1 ) )
!
!  Options file is test.opt (if available)
!
   coi_error = max( coi_error, coidef_optfile( cntvect, 'mp_lincns.opt' ) )
!
!  Tell CONOPT about the callback routines:
!
   coi_error = max( coi_error, coidef_readmatrix( cntvect, lincns_readmatrix ) )
   coi_error = max( coi_error, coidef_status( cntvect, std_status     ) )
   coi_error = max( coi_error, coidef_solution( cntvect, lincns_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_triord( cntvect, std_triord     ) )
 
#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 -- first with the default single thread:
!
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time1 = omp_get_wtime() - time0
 
   write(*,*)
   write(*,*) 'Single Thread: End of LinCns. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "Single Thread: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Single Thread: Status or Solution routine was not called", 1 )
   elseif ( mstat /= 15 .or. sstat /= 1 ) then      ! 15 since the model is linear
      call flog( "Single Thread: The model or solver status was not (15,1) as expected", 1 )
   elseif ( miter /= 1 ) then ! we expect one iteration for a linear CNS
      call flog( "Single Thread: The iteration count was not 1 as expected", 1 )
   endif
!
!  Start CONOPT again using multiple threads everywhere:
!
   coi_error = max( coi_error, coidef_threads( cntvect, 0 ) )  ! 0 means use as many threads as you can
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time2 = omp_get_wtime() - time0
 
   write(*,*)
   write(*,*) 'Multiple Threads: End of LinCns. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "Multiple Threads: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Multiple Threads: Status or Solution routine was not called", 1 )
   elseif ( mstat /= 15 .or. sstat /= 1 ) then      ! 15 since the model is linear
      call flog( "Multiple Threads: The model or solver status was not (15,1) as expected", 1 )
   elseif ( miter /= 1 ) then ! we expect one iteration for a linear CNS
      call flog( "Multiple Threads: The iteration count was not 1 as expected", 1 )
   endif
 
   if ( coi_free( cntvect ) /= 0 ) call flog( "Error while freeing control vector", 1 )
 
   write(*,*)
   write(*,"('Time for single thread',f10.3)") time1
   write(*,"('Time for multi  thread',f10.3)") time2
   write(*,"('Speedup               ',f10.3)") time1/time2
   write(*,"('Efficiency            ',f10.3)") time1/time2/omp_get_max_threads()
 
   call flog( "Successful Solve", 0 )
 
end Program lincns
!
! =====================================================================
!  Define information about the model structure
!
 
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function lincns_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 :: Lincns_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 :: i, j, k
!
!  Define the information for the columns. All are unbounded with initial value I
!
   do i = 1, n
      curr(i) = i
   enddo
!
!  All are equalities:
!
   do i = 1, m
      type(i) = 0
   enddo
!
!  Right hand sides: Initialize to 0 and revise so the solution becomes
!  x(i) = 1.0
!
   do i = 1, m
      rhs(i) = 0.d0
   enddo
!
!  Define the structure and content of the Jacobian:
!  We have a pattern with the largest values on the cross-diagonal from
!  (1,N) to (N,1) with gradually smaller value away from this diagonal.
!  The off-diagonal values are bounded away from zero to force a
!  completely dense model even if small values are filtered out.
!
   k = 1
   do i = 1, n
      colsta(i) = k
      do j = 1, n
         rowno(k) = j
         value(k) = (-1)**(i+j)*max(1.d-4, 2.0d0**(-abs(i+j-(n+1))))
         rhs(j)   = rhs(j) + value(k)
         k        = k + 1
      enddo
   enddo
   colsta(n+1) = k
 
   lincns_readmatrix = 0
 
end Function lincns_readmatrix
 
Integer Function lincns_solution( XVAL, XMAR, XBAS, XSTA, YVAL, YMAR, YBAS, YSTA, N, M, USRMEM )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Lincns_Solution
#endif
!
!  We write the solution values to the 'Documentation file' on unit 10
!  and test that the solution values are 1.0
!
   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(*)
 
   INTEGER i
   CHARACTER*5, Parameter, Dimension(4) :: stat = (/ 'Lower','Upper','Basic','Super' /)
 
   Logical error
 
   error = .false.
   lincns_solution = 0
   WRITE(10,"(/' Variable    Solution value    Reduced cost    B-stat')")
   DO i = 1, n
      WRITE(10,"(1P,I7,E20.6,E16.6,4X,A5 )") i, xval(i), xmar(i), stat(1+xbas(i))
      If ( abs(xval(i)-1.d0) > 1.d-7 ) then
         write(10,*) '**** Bad solution value'
         error = .true.
         lincns_solution = 1
      Endif
   ENDDO
 
   WRITE(10,"(/' Constrnt    Activity level    Marginal cost   B-stat')")
   DO i = 1, m
      WRITE(10,"(1P,I7,E20.6,E16.6,4X,A5 )") i, yval(i), ymar(i), stat(1+ybas(i))
   ENDDO
 
   solcalls = solcalls + 1
 
END Function lincns_solution