mp_polygon.f90

Go to the documentation of this file.

!> @file mp_polygon.f90
!! @ingroup FORTOPENMP_EXAMPLES
!!
!! 
!! Polygon model from COPS test set.
!! 
!!
!! 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 pdata
   Integer, Parameter :: NP = 150 ! 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
end module pdata
 
!> Main program. A simple setup and call of CONOPT
!!
Program polygon
 
   Use proginfop
   Use conopt
   Use pdata
   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 :: 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
#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 :: Std_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
   Integer :: iter1, iter2
!
!  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_polygon.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, std_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
   coi_error = max( coi_error, coidef_itlim( cntvect, 50  ) ) ! only 50 iterations to get reasonable time
!
!  Start CONOPT -- First time using a single thread:
!
   time0 = omp_get_wtime()
   coi_error = coi_solve( cntvect )
   time1 = omp_get_wtime() - time0
 
   write(*,*)
   write(*,*) 'End of Polygon example with one thread. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "One thread: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "One thread: Status or Solution routine was not called", 1 )
!   elseif ( sstat /= 1 .or. mstat /= 2 ) then
!      call flog( "One thread: Solver and Model Status was not as expected (1,2)", 1 )
   endif
   iter1 = miter
!
!  Start CONOPT again -- this time using multiple threads:
!
   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(*,*) 'End of Polygon example with multiple threads. Return code=',coi_error
 
   If ( coi_error /= 0 ) then
      call flog( "Multi thread: Errors encountered during solution", 1 )
   elseif ( stacalls == 0 .or. solcalls == 0 ) then
      call flog( "Multi thread: Status or Solution routine was not called", 1 )
!   elseif ( sstat /= 1 .or. mstat /= 2 ) then
!      call flog( "Multi thread: Solver and Model Status was not as expected (1,2)", 1 )
   endif
   iter2 = miter
 
   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()
   write(*,"('Speedup/iteration     ',f10.3)") time1/time2/iter1*iter2
   write(*,"('Efficiency/iteration  ',f10.3)") time1/time2/iter1*iter2/omp_get_max_threads()
 
   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 pdata
   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 pdata
   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