mp_lincns2.f90

Go to the documentation of this file.

!> @file mp_lincns2.f90
!! @ingroup FORTOPENMP_EXAMPLES
!!
!! 
!! Large Linear Dense CNS model.
!! Used to test the inversion routine for a very dense model.
!! The model is similar to Lincns.f90, but many nonzeros are so small
!! that they are removed by the factorization routine. The
!! factorization will therefore have to handle many fairly large
!! blocks with many pivots and streamlining the updates is a
!! challenge.
!! 
!!
!! 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   ! 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 :: 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
 
   real*8 time0, time1, time2
!
!  Create and initialize a Control Vector
!
   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
!
   coi_error = max( coi_error, coidef_optfile( cntvect, 'mp_lincns2.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 with a single thread:
!
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time1 = omp_get_wtime() - time0
 
   write(*,*)
   write(*,*) 'Single Thread: End of LinCns2. 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 LinCns2. 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.
!
   k = 1
   do i = 1, n
      colsta(i) = k
      do j = 1, n
         rowno(k) = j
         value(k) = 2.0d0**(-abs(i+j-(n+1))) ! many value are really small
         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