mp_polygon2.f90

Go to the documentation of this file.

!> @file mp_polygon2.f90
!! @ingroup FORTOPENMP_EXAMPLES
!!
!! 
!! Polygon model from COPS test set.
!! The model is run sequentially and with various number of threads but with ThreadC defined.
!! The models with various number of threads should produce identical solutions and we test
!! this.
!! 
!!
!! 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 pdata2
   Integer, Parameter :: NP = 50 ! Number of polygon points
   Integer, Parameter :: NV = 2*np ! Number of variables
   Integer, Parameter :: no = np-1 ! Number of Ordered constraints
   Integer, Parameter :: nd = np*(np-1)/2 ! Number of distance constraints
   Integer, dimension(ND) :: rowmapi, rowmapj ! Map from distance index to i and j
   Integer, dimension(NP,NP) :: rowmapk
   real*8, dimension(NV) :: xprim, xprim1, xprim2, xprim4
end module pdata2
 
!> Main program. A simple setup and call of CONOPT
!!
Program polygon
 
   Use proginfop
   Use conopt
   Use pdata2
   Use omp_lib
   implicit None
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: poly_readmatrix ! Mandatory Matrix definition routine defined below
   Integer, External :: poly_fdeval     ! Function and Derivative evaluation routine
                                        ! needed a nonlinear model.
   Integer, External :: std_status      ! Standard callback for displaying solution status
   Integer, External :: poly_solution   ! Special 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
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Poly_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Poly_FDEval
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Std_Status
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Poly_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
 
   real*8 time0, time1, time2, time4
   Integer :: iter1, iter2, iter4
   Integer :: i
   Logical :: same12, same14, same24
!
!  The model shown here in GAMS format
!
! set i sides / i1*i%nv% /;
!
! alias (i,j);
!
! scalar pi a famous constant;
!
!
! positive variables
!    r(i)      polar radius (distance to fixed vertex)
!    theta(i)  polar angle (measured from fixed direction)
! variable  polygon_area
!
! equations
!    obj
!    distance(i,j)
!    ordered(i) ;
!
! obj.. polygon_area =e= 0.5*sum(j(i+1), r(i+1)*r(i)*sin(theta(i+1)-theta(i)));
!
! ordered(i+1).. theta(i) =l= theta(i+1);
!
! distance(i,j)$(ord(j) > ord(i))..
!    sqr(r(i)) + sqr(r(j)) - 2*r(i)*r(j)*cos(theta(j)-theta(i)) =l= 1;
!
!
! pi = 2*arctan(inf);
!
! r.up(i)     = 1;
! theta.up(i) = pi;
!
! r.fx('i%nv%')    =  0;
! theta.fx('i%nv%') = pi;
!
! r.l(i)     = 4*ord(i)*(card(i)+ 1- ord(i))/sqr(card(i)+1);
! theta.l(i) = pi*ord(i)/card(i);
!
! model polygon /all/;
!
! $if set workspace polygon.workspace = %workspace%;
!
! solve polygon using nlp maximizing polygon_area;
!
!  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:
!
   coi_error = max( coi_error, coidef_numvar( cntvect, nv ) )           ! # variables
   coi_error = max( coi_error, coidef_numcon( cntvect, no + nd + 1 ) )  ! # constraints
   coi_error = max( coi_error, coidef_numnz( cntvect, 2*no+4*nd+nv ) ) ! # nonzeros in the Jacobian
   coi_error = max( coi_error, coidef_numnlnz( cntvect, 4*nd+nv ) )      ! # of which are nonlinear
   coi_error = max( coi_error, coidef_optdir( cntvect, +1 ) )           ! Maximize
   coi_error = max( coi_error, coidef_objcon( cntvect, no + nd + 1 ) )  ! Objective is last constraint
   coi_error = max( coi_error, coidef_optfile( cntvect, 'mp_polygon2.opt' ) )
!
!  Tell CONOPT about the callback routines:
!
   coi_error = max( coi_error, coidef_readmatrix( cntvect, poly_readmatrix ) )
   coi_error = max( coi_error, coidef_fdeval( cntvect, poly_fdeval     ) )
   coi_error = max( coi_error, coidef_status( cntvect, std_status     ) )
   coi_error = max( coi_error, coidef_solution( cntvect, poly_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 -- First time use ThreadC < ThreadS. We should get error 115.
!
   coi_error = max( coi_error, coidef_threads( cntvect, 2 ) )  !
   coi_error = max( coi_error, coidef_threadc( cntvect, 1 ) )  !
   coi_error = coi_solve( cntvect )
 
   write(*,*)
   write(*,*) 'End of Polygon example with one thread. Return code=',coi_error
 
   If ( coi_error /= 115 ) then
      call flog( "ThreadC = 1 and ThreadS = 2 did not give error return 115.", 1 )
   endif
   coi_error = max( coi_error, coidef_threadc( cntvect, 8 ) )  ! Cut work in 8 when testing threads
!
!  Start CONOPT again -- this time exactly 1 thread
!
   coi_error = max( coi_error, coidef_threads( cntvect, 1 ) )  ! Use 1 thread
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time1 = omp_get_wtime() - time0
 
   write(*,*)
   write(*,*) 'End of Polygon example with 1 threads. 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 ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Single thread: Solver and Model Status was not as expected (1,2)", 1 )
   endif
   iter1 = miter
   xprim1 = xprim
!
!  Start CONOPT again -- this time exactly 2 threads -- we keep ThreadC = 8
!
   coi_error = max( coi_error, coidef_threads( cntvect, 2 ) )   ! Use 2 thread
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time2 = omp_get_wtime() - time0
 
   write(*,*)
   write(*,*) 'End of Polygon example with 2 threads. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "Double threads: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Double threads: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Double threads: Solver and Model Status was not as expected (1,2)", 1 )
   endif
   iter2 = miter
   xprim2 = xprim
!
!  Start CONOPT again -- this time exactly 4 threads -- we keep ThreadC = 16
!
   coi_error = max( coi_error, coidef_threads( cntvect, 4 ) )   ! Use 4 thread
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time4 = omp_get_wtime() - time0
 
   write(*,*)
   write(*,*) 'End of Polygon example with 4 threads. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "Quad threads: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Quad threads: Status or Solution routine was not called", 1 )
   elseif ( sstat /= 1 .or. mstat /= 2 ) then
      call flog( "Quad threads: Solver and Model Status was not as expected (1,2)", 1 )
   endif
   iter4 = miter
   xprim4 = xprim
 
   same12 = .true.
   same14 = .true.
   same24 = .true.
   Do i = 1, nv
      If ( xprim2(i) /= xprim1(i) ) then
         if ( same12 ) then
            write(*,*) 'Solution 2 is different from Solution 1. First diff for index',i
         endif
         same12 = .false.
      endif
      If ( xprim4(i) /= xprim1(i) ) then
         if ( same14 ) then
            write(*,*) 'Solution 4 is different from Solution 1. First diff for index',i
         endif
         same14 = .false.
      endif
      If ( xprim4(i) /= xprim2(i) ) then
         if ( same24 ) then
            write(*,*) 'Solution 4 is different from Solution 2. First diff for index',i
         endif
         same24 = .false.
      endif
   Enddo
 
   write(*,*)
   write(*,"('Time for single thread',f10.3)") time1
   write(*,"('Time for double thread',f10.3)") time2
   write(*,"('Time for quad   thread',f10.3)") time4
   write(*,"('Speedup  double       ',f10.3)") time1/time2
   write(*,"('Speedup  quad         ',f10.3)") time1/time4
   if ( same12 .and. same14 .and. same24 ) then
      write(*,*) 'All solutions are the same'
   else
      write(*,*) 'The solutions are NOT the same'
!      call flog( "Solutions are not the same", 1 )
   Endif
 
   if ( coi_free( cntvect ) /= 0 ) call flog( "Error while freeing control vector", 1 )
 
   call flog( "Successful Solve", 0 )
 
End Program polygon
!
! ============================================================================
!  Define information about the model:
!
 
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function poly_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 :: Poly_ReadMatrix
#endif
   Use pdata2
   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
   real*8, parameter :: pi = 3.141592
 
!
!  Information about Variables:
!  Default: Lower = -Inf, Curr = 0, and Upper = +inf.
!  Default: the status information in Vsta is not used.
!
! r.up(i)     = 1;
! theta.up(i) = pi;
!
! r.fx('i%nv%')    =  0;
! theta.fx('i%nv%') = pi;
!
! r.l(i)     = 4*ord(i)*(card(i)+ 1- ord(i))/sqr(card(i)+1);
! theta.l(i) = pi*ord(i)/card(i);
   DO i = 1, np-1
      lower(i)    = 0.d0
      lower(i+np) = 0.d0
      upper(i)    = 1.d0
      upper(i+np) = pi
      curr(i)     = 4.d0*i*(np+1-i)/(1+np)**2
      curr(i+np)  = pi*i/np
   enddo
   lower(np)   = 0.d0
   upper(np)   = 0.d0
   curr(np)    = 0.d0
   lower(2*np) = pi
   upper(2*np) = pi
   curr(2*np)  = pi
!
!  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.
!
! ordered(i+1).. theta(i) =l= theta(i+1);
!
! distance(i,j)$(ord(j) > ord(i))..
!    sqr(r(i)) + sqr(r(j)) - 2*r(i)*r(j)*cos(theta(j)-theta(i)) =l= 1;
!
! obj.. polygon_area =e= 0.5*sum(j(i+1), r(i+1)*r(i)*sin(theta(i+1)-theta(i)));
!
   k = 0
   DO i = 1, no
      k       = k + 1
      Type(k) = 2 ! Ordered constraint Less than and RHS = 0
   Enddo
   rowmapk = 0 ! Initialize so we can recognize (i,j) combinations that are not equations
   DO i = 1, np-1
      DO j = i+1, np
         k       = k + 1
         Type(k) = 2 ! Dist(i,j) Less than
         rhs(k)  = 1.d0
         rowmapi(k-no) = i ! To be used for easier lookup
         rowmapj(k-no) = j
         rowmapk(i,j)  = k
      Enddo
   Enddo
   k = k + 1
   type(k) = 3 ! Objective of type Non binding
!
!  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
!
!       r(1)  r(2)  r(3)  r(4)  .. Th(1)  th(2)  th(3)  th(4)
!  1:     1    -1
!  2:           1    -1
!  3:                 1    -1
!  ..
!  NP:   NL    NL                   NL      NL
!  NP+1: NL          NL             NL             NL
!  NP+2  NL                NL       NL                    NL
!  ..
!              NL    NL                     NL     NL
!              NL          NL               NL            NL
!  ..
!  Obj:  NL    NL   NL      NL      NL    NL
!
!  Nonlinearity Structure: L = 0 are linear and NL = 1 are nonlinear
!  All nonzeros are nonlinear
!
   k = 0
   DO i = 1, np        ! r variables
      colsta(i) = k+1
      DO j = i+1, np
         k = k + 1
         rowno(k)  = rowmapk(i,j) ! Distance - i index
         nlflag(k) = 1
      enddo
      DO j = 1,np
         IF ( rowmapk(j,i) > 0 ) Then
            k = k + 1
            rowno(k)  = rowmapk(j,i) ! Distance - j index
            nlflag(k) = 1
         endif
      enddo
      k = k + 1
      rowno(k)  = no+nd+1  ! Objective constraint
      nlflag(k) = 1
   enddo
   DO i = 1, np         ! theta variables
      colsta(np+i) = k+1
      if ( i < np ) then
         k = k + 1
         rowno(k)  = i    ! ordered constraint index i
         nlflag(k) = 0
         value(k)  = 1.d0
      endif
      if ( i > 1 ) then
         k = k + 1
         rowno(k)  = i-1  ! ordered constraint index i+1
         nlflag(k) = 0
         value(k)  = -1.d0
      endif
      DO j = i+1, np
         k = k + 1
         rowno(k)  = rowmapk(i,j) ! Distance - i index
         nlflag(k) = 1
      enddo
      DO j = 1,np
         IF ( rowmapk(j,i) > 0 ) Then
            k = k + 1
            rowno(k)  = rowmapk(j,i) ! Distance - j index
            nlflag(k) = 1
         endif
      enddo
      k = k + 1
      rowno(k)  = no+nd+1  ! Objective constraint
      nlflag(k) = 1
   enddo
   colsta(nv+1) = k+1
 
   poly_readmatrix = 0 ! Return value means OK
 
end Function poly_readmatrix
!
!==========================================================================
!  Compute nonlinear terms and non-constant Jacobian elements
!
 
!>  Compute nonlinear terms and non-constant Jacobian elements
!!
!! @include{doc} fdeval_params.dox
Integer Function poly_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 :: Poly_FDEval
#endif
   Use pdata2
   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 :: ne
   Integer :: i, j
   real*8 :: s, c
 
   ne = n / 3
!
!  Row NO+ND+1: the objective function is nonlinear
!
   if ( rowno == no+nd+1 ) then
!
!  Mode = 1 or 3. Function value:
!  obj.. polygon_area =e= 0.5*sum(j(i+1), r(i+1)*r(i)*sin(theta(i+1)-theta(i)));
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g = 0.d0
         do i = 1, np-1
            g = g + x(i+1)*x(i)*sin(x(np+i+1)-x(np+i))
         enddo
         g = g * 0.5d0
      endif
!
!  Mode = 2 or 3: Derivative values: w.r.t. x: 2*(x(i)-x(j))*(-1/2)*dist(-3/2)
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         do i = 1, nv
            jac(i) = 0.d0
         enddo
         do i = 1, np-1
!            g = g + 0.5* x(i+1)*x(i)*Sin(x(NP+i+1)-x(NP+i))
            s           = 0.5d0*sin(x(np+i+1)-x(np+i))
            jac(i+1)    = jac(i+1) + x(i)*s
            jac(i)      = jac(i)   + x(i+1)*s
            c           = 0.5d0*x(i+1)*x(i)*cos(x(np+i+1)-x(np+i))
            jac(np+i+1) = jac(np+i+1) + c
            jac(np+i)   = jac(np+i)   - c
         enddo
      endif
!
   Else if ( rowno > no ) then ! this is a distance constraint
!
!  Mode = 1 or 3. Function value:
!    sqr(r(i)) + sqr(r(j)) - 2*r(i)*r(j)*cos(theta(j)-theta(i)) =l= 1;
!
      i = rowmapi(rowno-no)
      j = rowmapj(rowno-no)
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         g = x(i)**2 + x(j)**2 -2.d0*x(i)*x(j)*cos(x(j+np)-x(i+np))
      endif
!
!  Mode = 2 or 3: Derivative values:
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         c         = cos(x(j+np)-x(i+np))
         jac(i)    = 2.d0*(x(i)-x(j)*c)
         jac(j)    = 2.d0*(x(j)-x(i)*c)
         s         = 2.d0*sin(x(j+np)-x(i+np))*x(i)*x(j)
         jac(i+np) = -s
         jac(j+np) = s
      endif
   Else
      write(*,*) 'Fdeval called with rowno=',rowno
      poly_fdeval = 1; return ! error in row number
   endif
   poly_fdeval = 0
 
end Function poly_fdeval
 
Integer Function poly_solution( XVAL, XMAR, XBAS, XSTA, YVAL, YMAR, YBAS, YSTA, N, M, USRMEM )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: Poly_Solution
#endif
!
!  Simple implementation in which we write the solution values to
!  the 'Documentation file' on unit 10.
!
   Use proginfop
   Use pdata2
   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' /)
 
   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))
   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
   xprim(1:n) = xval(1:n)
   Call flush(10)
 
   solcalls = solcalls + 1
   poly_solution = 0
 
END Function poly_solution