mp_qpbandb.f90

Go to the documentation of this file.

!> @file mp_qpbandb.f90
!! @ingroup FORTOPENMP_EXAMPLES
!!
!! 
!! QPBandB: Multi Threaded implementation of QP4 MINLP model
!! The implementation shows how to implement a parallel Branch & Bound
!! procedure that still is deterministic independent of the order in
!! which the sub-models are solved (within the framework defined here
!! with Work-pools).
!! 
!! When testing that the sequential and parallel algorithms uses the same
!! number of nodes we have a test that the solves are identical even
!! when run with different threads. This is some kind of test that
!! CONOPT Solves are reproducable and that there are no un-initialized
!! values used in the algorithm.
!! 
!!
!! 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 qp_data
  Implicit None
  Integer, Parameter        :: NN   = 32   ! Dimension of Q
  Integer, Parameter        :: Nvar = 2*nn ! Number of variables in model
  Integer, Parameter        :: nbin = nn   ! Number of binary variables in model
 
  Integer, Parameter        :: nobs = 60   ! Observations used to estimate Q
  Real(8), dimension(NN,NN) :: q
  Real(8), dimension(nn)    :: mean
  Real(8)                   :: target
  Integer, Parameter        :: maxbin = 15
  Integer, Dimension(Nbin)  :: binmap      ! Mapping from binary to index in all variables
End Module qp_data
 
Module bb_data
!
!  Data used to manage the Branch and Bound Process
!
  Use qp_data
  Use omp_lib
  Implicit None
!
!  Branch and bound node data
!
  Type bbnode
     Integer                  :: Number    ! Sequense number of the node.
     real(8), dimension(NVar) :: Values   ! Solution value or initial values
     Integer, dimension(Nbin) :: fixed    ! Binary: 0=fixed zero, 1=fixed 1, -1 not fixed.
     Type (bbnode), Pointer   :: prvnode, nxtnode ! Pointer forward and backward in
                                                  ! double-linked list of active nodes
     real(8)                  :: object    ! Objective value (if solved) or lower bound if
                                           ! not yet solved
     Integer                  :: node_stat ! Node Status:
                                           ! 0 - created but not solved
                                           ! 1 - solved
     Integer                  :: solstat   ! Solution status. Initialized to -1 and
                                           ! redefined when the model has been solved
     Integer                  :: modstat   ! Model status. Initialized to -1 and
                                           ! redefined when the model has been solved
     Integer                  :: depth     ! Depth of the node in the B&B tree
     Logical                  :: isinteger ! Is true iff the solution is integer
     Integer                  :: sthread   ! The thread being used when solving
  End type bbnode
 
  Type bbnodeptr
     Type (bbnode), Pointer   :: bnode
  End Type bbnodeptr
 
  Real(8), Parameter :: bintol = 1.d-7
  Integer :: totnodes  ! Total number of nodes in use
  Integer :: freenodes ! Number of active nodes
  Integer :: actnodes  ! Number of free nodes in a free list
  Integer :: numsolve  ! Counter for finished Solves.
!
!  Data for best solution
!
  Logical :: intfound  ! Global flag indicating whether we already have an integer
                       ! feasible solution.
  Integer :: numintfound ! Number of integer solutions found
  Real(8) :: bestint   ! Objective value for the best integer solution found so far.
  Real(8), dimension(NVar) :: solution ! The solution corresponding to BestInt
!
!  Data for managing the Branch and Bound Tree
!
  Type (bbnode), Pointer :: actnode   ! Pointer to first node in a list of active untested noded
  Type (bbnode), Pointer :: freenode  ! Pointer to first node in a freelist
!
!  Data for managing a WorkPool with the nodes selected for the next
!  set of parallel Solves
!
  Integer, Parameter :: maxmaxpool = 64 ! Allocation size for WorkPool
  Integer :: maxpool                  ! Largest size of WorkPool selected by user
  Logical :: maxdyn                   ! If .true. then MaxPool above is dynamic
  Integer :: workpool                 ! Size of the actual workpool
  Type (bbnodeptr) :: work_node(maxmaxpool)  ! Pointers to the selected nodes in the WorkPool
!
  integer :: maxthread                ! Maximum number of threads
  Logical :: parallel                 ! If true then run the COI_Solve loop over the nodes
                                      ! in the WorkPool in parallel
  Type controlv
    Integer, Dimension(:), Pointer :: cntvect
  End Type controlv
  Type (controlv), Dimension(:), Pointer :: cntvect1 ! CONOPT Control vectors, one for each thread
  Real(8) :: stime                    ! Accumulated time in Solve Loop
 
  Contains
!
!  Routines used to manage the Branch and Bound Tree
!
Subroutine remove_activenode(Node)
!
!  Remove a node from the list of active nodes.
!
   Implicit None
   Type (BBNode), Pointer :: Node
 
   if ( associated(node%prvnode) ) then
      node%prvnode%nxtnode => node%nxtnode
   else
      actnode => node%nxtnode
   endif
   if ( associated(node%nxtnode) ) then
      node%nxtnode%prvnode => node%prvnode
   endif
   actnodes  = actnodes - 1
 
End Subroutine remove_activenode
 
Subroutine return_freenode(Node)
!
!  Link the node into the free list. The free list does not use prvnode
!
   Implicit None
   Type (BBNode), Pointer :: Node
 
   node%nxtnode => freenode
   freenode     => node
   freenodes    =  freenodes + 1
 
End Subroutine return_freenode
 
Subroutine get_freenode(Node)
!
!  Get a free node to be used in the list of active nodes and initialize
!  the status.
!
   Implicit None
   Type (BBNode), Pointer :: Node
 
   if ( associated(freenode) ) then
      node      => freenode
      freenode  => freenode%Nxtnode
      freenodes = freenodes - 1
   else
      Allocate(node)
      Nullify(node%NxtNode);
      Nullify(node%PrvNode);
      totnodes    = totnodes + 1
      node%Number = totnodes
   endif
   node%Node_Stat = 0
   node%Solstat   =-1
   node%Modstat   =-1
   node%Sthread   =-1
 
End Subroutine get_freenode
 
Subroutine deallocate_nodes
!
!  Deallocates the nodes in the free list between solves so the next
!  solve starts without any nodes.
!
   Implicit None
   Type (BBNode), Pointer :: Node, Next
 
   node => freenode
   do while ( Associated(node) )
      next        => node%Nxtnode
      Deallocate(node);
      node        => next
   enddo
   Nullify(freenode)
 
End Subroutine deallocate_nodes
 
Subroutine add_activenode(Node)
!
!  Add a node to the front of the list of Active nodes
!
   Implicit None
   Type (BBNode), pointer :: Node
 
   actnodes = actnodes + 1
   if ( associated(actnode) ) then
      node%NxtNode    => actnode
      actnode%PrvNode => node;
   else
      Nullify(node%NxtNode);
   endif
   Nullify(node%PrvNode);
   actnode => node
 
End Subroutine add_activenode
 
Subroutine add2workpool( Node )
!
!  Add a node to the Workpool list
!
   Implicit None
   Type (BBNode), pointer :: Node
 
   workpool     =  workpool + 1
!   write(*,*) 'Work_node(',Workpool,') has number',node%number
   work_node(workpool)%Bnode => node
 
end Subroutine add2workpool
 
Subroutine list_node(node)
 
   Implicit None
   Type (BBNode), pointer :: Node
 
   If ( Associated(node) ) Then
      write(*,*) 'List_Node: Number=',node%Number,' Node_Stat=',node%node_stat
   Else
      write(*,*) 'List_Node: Node not associated.'
   Endif
End Subroutine list_node
 
Subroutine list_nodes
 
   Implicit None
   Type (BBNode), Pointer :: Node
   Integer :: N
 
   n = 0
   node => actnode
   do while ( Associated(node) )
      write(*,"('List_nodes: Node number=',i4,' Node_Stat=',I2,' Object=',1pe15.8)") node%number, node%node_stat, node%object
      node => node%nxtnode
      n = n + 1
   enddo
   write(*,*) 'List_nodes: Total=',n
 
End Subroutine list_nodes
 
Subroutine select_var(Node, IB)
!
!  Select variable to branch on, counted in the space of binary variables.
!  Currently there is only one rule, namely to select based on max
!  infeasibility. Other Variable selection rules could be implemented
!  in this routine.
!
   Implicit None
   Type (BBNode), Pointer :: Node
   integer :: IB, i, j
   Real(8) :: Dist, Maxdist
 
   maxdist = 0.d0
   ib = 0
   do i = 1, nn
      j = binmap(i)
      dist = min( node%values(j), 1.d0-node%values(j) )
      if ( dist > maxdist ) then
         maxdist = dist
         ib      = i
      endif
   enddo
 
End Subroutine select_var
 
subroutine create_workpool
!
!  Create a workpool of up to MaxPool nodes to be processed in
!  parallel.
!  The nodes are removed from the active list.
!
   Implicit None
   Type (BBNode), Pointer :: Node, Curr, Next
   Real(8) :: Cutoff
   integer :: IB   ! Branching variable
 
   Nullify(node)
   workpool = 0
   write(*,*) 'Create_WorkPool'
!
!  If the max size of the workpool is dynamic then use the formula
!  MaxThread*min(8,1+ActNodes/(MaxThread*5)). With few nodes
!  (<MaxThread*5) we will select MaxThread nodes. With more nodes we
!  will select gradually more up to a maximum of MaxThread*8.
!
   if ( maxdyn ) maxpool = maxthread*min(8,1+actnodes/(maxthread*5))
   do while ( workpool < maxpool )
      cutoff = 1.d30
      curr   => actnode
!
!  Run through the linked list of active nodes and find the best.
!  We use the objective but other criteria could be used
!
      Nullify(node)
      do while (associated(curr))
         next => curr%NxtNode
         if ( curr%Object < cutoff ) then
            node   => curr
            cutoff =  curr%Object
         endif
         curr => next
      enddo
      if ( .not. associated(node) ) then
         exit
      endif
      if ( node%Node_Stat == 0 ) then
!
!  The node has been created but not solved.
!  Take it immediately: remove it from the active list and add it
!  to the work_node vector
!
         write(*,"('Unsolved node',i5,' selected. Obj=',e24.17)") node%Number,node%Object
         Call remove_activenode(node)
         call add2workpool(node)
      else
!
!  The node has been solved. Split it in two and take one or both
!  depending on the space left in the workpool. The total number
!  of nodes will grow by one so a dynamic MaxPool must be adjusted
!  before the space test below
!
         if ( maxdyn ) maxpool = maxthread*min(8,1+(actnodes+1)/(maxthread*5))
!
!  Get a new node and, select a branching variable, and copy the
!  old node into the new node.
!
         call get_freenode(next)
         call select_var( node, ib )
         next%Values = node%Values
         next%Fixed  = node%Fixed
         next%Object = node%Object
!
!  Make the two node unsolved and fix the branching variable down and up
!  respectively.
!
         node%Node_Stat = 0
         node%Solstat   = -1
         node%Modstat   = -1
         node%Depth     = node%Depth + 1
         node%Fixed(ib) = 0  ! Branch down
 
         next%Node_Stat = 0
         next%Solstat   = -1
         next%Modstat   = -1
         next%Depth     = node%Depth
         next%Fixed(ib) = 1  ! Branch up
!
!  Add the first node to the workpool
!
         Call remove_activenode(node)
         call add2workpool(node)
         write(*,"('New-Down node',i5,' selected. Obj=',e24.17)") node%Number,node%Object
!
!  Add the second node to the workpool or to the pool of active
!  nodes, depending on space in the workpool
!
         if ( workpool < maxpool ) then
            call add2workpool(next)
            write(*,"('New-Up   Next',i5,' selected. Obj=',e24.17)") next%Number,next%Object
         else
            call add_activenode(next)
         endif
      endif
   enddo
 
end subroutine create_workpool
 
subroutine solve_workpool
!
!  Solve the models in the WorkPool and integrate the solutions into
!  the active list of nodes. Update the best solution if one is found.
!
   Use proginfop
   Use conopt
   implicit none
 
   integer nw ! index for active node in the WorkPool
   integer Thread
   Integer, Dimension(:), pointer :: CntVect ! CONOPT Control vector for current node
   integer COI_Error
   Type (BBNode), Pointer :: Node, Curr, Next
   Logical :: NewInt
   Real(8) Time0
 
!   write(*,*) 'Calling Solve_WorkPool with',WorkPool,' nodes.'
   time0 = omp_get_wtime()
   if ( workpool .le. 1 .OR. .NOT. parallel ) Then
! This is a serial version
      thread = 1
      do nw = 1, workpool
         node => work_node(nw)%Bnode
         if ( debug ) then
            write(16,*) 'Starting to solve node',node%number,' using thread',thread
         endif
         cntvect => cntvect1(thread)%CntVect
         coi_error = coidef_usrmem( cntvect, node ) ! Pass the current node through Usrmem
         If ( coi_error .ne. 0 ) THEN
            write(*,*)
            write(*,*) '**** Fatal Error while loading CONOPT Callback routines.'
            write(*,*)
            stop
         ENDIF
         node%Sthread = thread
         coi_error = coi_solve( cntvect )
         If ( coi_error .ne. 0 ) THEN
            write(*,*)
            write(*,*) '**** Fatal Error while solving.'
            write(*,*)
            stop
         ENDIF
      enddo
   Else
!  And this is the corresponding Parallel version of the SOLVE loop.
!$OMP PARALLEL DEFAULT(SHARED) Private(COI_Error,Thread, Node)
!$OMP DO PRIVATE(Thread, COI_Error, Node, CntVect) SCHEDULE(Guided)
      do nw = 1, workpool
         node => work_node(nw)%Bnode
         thread = 1+omp_get_thread_num()
         if ( debug ) then
            write(16,*) 'Starting to solve node',node%number,' using thread',thread
         endif
         cntvect => cntvect1(thread)%CntVect
         coi_error = coidef_usrmem( cntvect, node ) ! Pass the current node through Usrmem
         If ( coi_error .ne. 0 ) THEN
            write(*,*)
            write(*,*) '**** Fatal Error while loading CONOPT Callback routines.'
            write(*,*)
            call flog( "Errors encountered during setup for Solve.", 1 )
         ENDIF
         node%Sthread = thread
         coi_error = coi_solve( cntvect )
         If ( coi_error .ne. 0 ) THEN
            write(*,*)
            write(*,*) '**** Fatal Error while solving.'
            write(*,*)
            call flog( "Errors encountered during Solve.", 1 )
         ENDIF
      enddo
!$OMP END PARALLEL
   Endif
   stime = stime + omp_get_wtime() - time0
   numsolve = numsolve + workpool
!
!  Process the solutions. The loop is not parallelized so the order here
!  is fixed, independent of the order in which the nodes finished
!  solving.
!
   newint = .false. ! No new integer solution in this pass
   do nw = 1, workpool
      node => work_node(nw)%Bnode
      write(*,'("Node=",I5," Solstat=",I3," Modstat=",I3," Object=",F24.17)') &
         node%number, node%Solstat, node%modstat, node%Object
      if ( node%Solstat .ne. 1 .or. node%Modstat .gt. 2 ) then
!
! Node was not solved properly
!
         call flog( "Solve did not return an normal solver status.", 1 )
      Else if ( node%Modstat .gt. 2 ) then
!
! Node was not solved to local optimality. Fathom.
!
         write(*,*) 'Not optimal -- fathom'
         call return_freenode(node)
      else if ( intfound ) then
         if ( node%object >= bestint ) then
!
!  A better solution has been found. Fathom.
!
            write(*,*) 'Worse than optimal -- fathom'
            call return_freenode(node)
         else
            if ( node%IsInteger ) then
               write(*,*) 'Integer Solution Stored'
               numintfound = numintfound + 1
               newint   = .true.
               bestint  = node%object
               solution = node%Values
               call return_freenode(node)
            Else
!               write(*,*) 'Non-Integer. Kept'
               call add_activenode(node)
            Endif
         endif
      else
!
!  An Integer solution has not yet been found. We cannot fathom it.
!
         if ( node%IsInteger ) then
            write(*,*) 'Integer Solution Stored'
            numintfound = numintfound + 1
            intfound = .true.
            newint   = .true.
            bestint  = node%object
            solution = node%Values
            call return_freenode(node)
         Else
!            write(*,*) 'Non-Integer. Kept'
            call add_activenode(node)
         Endif
      endif
   enddo
!
!  If we did find a new integer solution then scan the list of active
!  nodes to check if some of them should be fathomed.
!
   if ( newint ) then
      curr => actnode
      do while (associated(curr))
         next => curr%NxtNode
         if ( curr%Object > bestint ) then
            write(*,"('Node',I5,' with Object=',F18.12,' fathomed by new objective')") curr%Number, curr%Object
            call remove_activenode(curr)
            call return_freenode(curr)
         endif
         curr => next
      enddo
   endif
 
end subroutine solve_workpool
 
End Module bb_data
 
REAL FUNCTION rndx( )
!
!  Defines a pseudo random number between 0 and 1
!
   IMPLICIT NONE
 
   Integer, save :: seed = 12359
 
   seed = mod(seed*1027+25,1048576)
   rndx = float(seed)/float(1048576)
 
END FUNCTION rndx
 
subroutine defdata
!
!  Define values for A, B, and Obs
!
   Use qp_data
   IMPLICIT NONE
 
   Integer :: i, j, k
   Real, External :: Rndx
   Real(8) :: Sum
   Real(8), dimension(NN,Nobs) :: Xobs
 
   Target = 10.0d0
   do i = 1, nn
      mean(i) = Target + (rndx()*4.d0-2.d0)
   enddo
!
!  Create a positive definite matrix computed as the sum of some
!  outer products
!
   do i = 1, nn
      do j = 1, nn
         q(i,j) = 0.d0
      enddo
   enddo
   do i = 1, nn
      sum = 0.d0
      do k = 1, nobs
         xobs(i,k) = 100.d0*rndx()
         sum       = sum + xobs(i,k)
      enddo
      sum = sum / nobs
      do k = 1, nobs
         xobs(i,k) = xobs(i,k) - sum
      enddo
   enddo
   do i = 1, nn
      do j = 1, i
         sum = 0.d0
         do k = 1, nobs
            sum = sum + xobs(i,k) * xobs(j,k)
         enddo
         q(i,j) = sum / (nobs-1)
         q(j,i) = q(i,j)
      enddo
   enddo
!
!  Create BinMap
!
   do i = 1, nn
      binmap(i) = nn+i
   enddo
 
end subroutine defdata
!
!  Main program. A simple setup and call of CONOPT for a QP model
!
!> Main program. A simple setup and call of CONOPT
!!
Program qp_threads
 
   Use proginfop
   Use conopt
   use bb_data
   use qp_data
   Use omp_lib
   Implicit None
!
!  Data needed for Branch and bound
!
   Integer :: thread
!
!  Declare the user callback routines as Integer, External:
!
   Integer, External :: qp_readmatrix  ! Mandatory Matrix definition routine defined below
   Integer, External :: qp_fdeval      ! Function and Derivative evaluation routine
                                       ! needed a nonlinear model.
   Integer, External :: qp_2dlagrstr   ! 2nd derivative routine for structure
   Integer, External :: qp_2dlagrval   ! 2nd derivative routine for values
   Integer, External :: qp_2ddir       ! Directional 2nd derivative routine
   Integer, External :: bb_status      ! callback for displaying solution status
   Integer, External :: bb_solution    ! callback for displaying solution values
   Integer, External :: bb_message     ! callback for managing messages
   Integer, External :: bb_errmsg      ! callback for managing error messages
   Integer, External :: bb_option      ! Define options using a callback. Files cannot be used
                                       ! for multithreading since they can be shared and not
                                       ! read in the proper order.
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_ReadMatrix
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_FDEval
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_2DLagrStr
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_2DLagrVal
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_2DDir
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: BB_Status
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: BB_Solution
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: BB_Message
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: BB_ErrMsg
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: BB_Option
#endif
!
!  Control vector
!
   INTEGER :: coi_error
   Integer, Dimension(:), Pointer :: cntvect ! CONOPT Control vector for current node
 
   Real(8) time0, time1, stime1
   Real(8) time8s, time16s, time32s, timeds
   Real(8) time8p, time16p, time32p, timedp
   Real(8) stime8s, stime16s, stime32s, stimeds
   Real(8) stime8p, stime16p, stime32p, stimedp
   Integer numsolve1
   Integer numsolve8s, numsolve16s, numsolve32s, numsolveds
   Integer numsolve8p, numsolve16p, numsolve32p, numsolvedp
 
   call startup
   maxthread = omp_get_max_threads()
   write(11,*) 'Initial MaxThread=',maxthread
   write(*,*) 'Initial MaxThread=',maxthread
   If ( maxthread .lt. 2 ) then
      maxthread = 2
      call omp_set_num_threads(maxthread)
      write(11,*) 'Revised MaxThread=',maxthread
      write(*,*) 'Revised MaxThread=',maxthread
   Endif
!   Debug = .true. ! Turn on if debugging output is needed
   If ( debug ) Then
      Open(12,file='threadA1.txt',status='unknown')
      Open(16,file='nodesA.txt',status='unknown')
   Endif
!
   call defdata
!
!  Create a vector of Control Vectors in CntVect1 and allocate and
!  initialize the individual vectors for each thread.
!
   Allocate( cntvect1( maxthread ) )
   coi_error = 0
   do thread = 1, maxthread
      coi_error = coi_create( cntvect )
      cntvect1(thread)%CntVect => cntvect
!
!  Tell CONOPT about the size of the model by populating the Control Vector:
!  We will create a QP model with one constraint, sum(i, x(i) ) = 1
!  and NN variables. NN and the Q matrix are defined in Module QP_Data.
!
      coi_error = max( coi_error, coidef_numvar( cntvect, 2*nn ) ) ! NN variables
      coi_error = max( coi_error, coidef_numcon( cntvect, 4+nn ) ) ! 3+NN constraint + 1 objective
      coi_error = max( coi_error, coidef_numnz( cntvect, 6*nn ) ) ! 6*NN nonzeros in the Jacobian
      coi_error = max( coi_error, coidef_numnlnz( cntvect, nn ) )   ! NN of which are nonlinear
      coi_error = max( coi_error, coidef_numhess( cntvect, nn*(nn+1)/2 ) )  ! Hessian is dense in the NN nonlinear variables
      coi_error = max( coi_error, coidef_optdir( cntvect, -1 ) )   ! Minimize
      coi_error = max( coi_error, coidef_objcon( cntvect, 1 ) )    ! Objective is constraint 1
!
!  Tell CONOPT about the callback routines:
!
      coi_error = max( coi_error, coidef_readmatrix( cntvect, qp_readmatrix ) )
      coi_error = max( coi_error, coidef_fdeval( cntvect, qp_fdeval     ) )
      coi_error = max( coi_error, coidef_2dlagrstr( cntvect, qp_2dlagrstr  ) )
      coi_error = max( coi_error, coidef_2dlagrval( cntvect, qp_2dlagrval  ) )
      coi_error = max( coi_error, coidef_2ddir( cntvect, qp_2ddir      ) )
      coi_error = max( coi_error, coidef_status( cntvect, bb_status     ) )
      coi_error = max( coi_error, coidef_solution( cntvect, bb_solution   ) )
      coi_error = max( coi_error, coidef_message( cntvect, bb_message    ) )
      coi_error = max( coi_error, coidef_errmsg( cntvect, bb_errmsg     ) )
      If ( debug ) &
      coi_error = max( coi_error, coidef_option( cntvect, bb_option   ) )
 
#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
   enddo
   If ( coi_error .ne. 0 ) THEN
      write(*,*)
      write(*,*) '**** Fatal Error while loading CONOPT Callback routines.'
      write(*,*)
      call flog( "Skipping First Solve due to setup errors", 1 )
   ENDIF
   write(11,*) 'Finished Setup of model'
   write(*,*) 'Finished Setup of model'
!
!  Solve the B and B process with MaxPool = 1
!
   parallel = .false.
   maxdyn   = .false.
   maxpool  = 1
   time0 = omp_get_wtime()
   Call bandb_solve
   time1 = omp_get_wtime() - time0
   stime1 = stime
   numsolve1 = numsolve
   write(11,*) 'Finished sequential with Maxpool=',maxpool
   write(*,*) 'Finished sequential with Maxpool=',maxpool
   If ( debug ) Then
      close(12); close(16)
      Open(12,file='threadB1.txt',status='unknown')
      Open(16,file='nodesB.txt',status='unknown')
   Endif
!
!  Redo the B and B process with MaxPool = 8
!
   maxpool  = 8
   time0 = omp_get_wtime()
   Call bandb_solve
   time8s = omp_get_wtime() - time0
   stime8s = stime
   numsolve8s = numsolve
   write(11,*) 'Finished sequential with Maxpool=',maxpool
   write(*,*) 'Finished sequential with Maxpool=',maxpool
!
!  Redo the B and B process with MaxPool = 16
!
   maxpool  = 16
   time0 = omp_get_wtime()
   Call bandb_solve
   time16s = omp_get_wtime() - time0
   stime16s = stime
   numsolve16s = numsolve
   write(11,*) 'Finished sequential with Maxpool=',maxpool
   write(*,*) 'Finished sequential with Maxpool=',maxpool
!
!  Redo the B and B process with MaxPool = 32
!
   maxpool  = 32
   time0 = omp_get_wtime()
   Call bandb_solve
   time32s = omp_get_wtime() - time0
   stime32s = stime
   numsolve32s = numsolve
   write(11,*) 'Finished sequential with Maxpool=',maxpool
   write(*,*) 'Finished sequential with Maxpool=',maxpool
!
!  Redo the B and B process with Dynamic MaxPool
!
   maxdyn   = .true.
   time0 = omp_get_wtime()
   Call bandb_solve
   timeds = omp_get_wtime() - time0
   stimeds = stime
   numsolveds = numsolve
   write(11,*) 'Finished sequential with Maxpool=',maxpool
   write(*,*) 'Finished sequential with Maxpool=',maxpool
   If ( debug ) Then
      close(12); close(16)
      Open(12,file='threadC1.txt',status='unknown')
      Open(13,file='threadC2.txt',status='unknown')
      Open(14,file='threadC3.txt',status='unknown')
      Open(15,file='threadC4.txt',status='unknown')
      Open(16,file='nodesC.txt',status='unknown')
   Endif
!
!  Redo again with same three MaxPool sizes but with a parallel solve loop.
!
   parallel = .true.
   maxdyn   = .false.
   maxpool  = 8
   time0 = omp_get_wtime()
   Call bandb_solve
   time8p = omp_get_wtime() - time0
   stime8p = stime
   numsolve8p = numsolve
   write(11,*) 'Finished parallel with Maxpool=',maxpool
   write(*,*) 'Finished parallel with Maxpool=',maxpool
   maxpool     = 16
   time0 = omp_get_wtime()
   Call bandb_solve
   time16p = omp_get_wtime() - time0
   stime16p = stime
   numsolve16p = numsolve
   write(11,*) 'Finished parallel with Maxpool=',maxpool
   write(*,*) 'Finished parallel with Maxpool=',maxpool
   maxpool     = 32
   time0 = omp_get_wtime()
   Call bandb_solve
   time32p = omp_get_wtime() - time0
   stime32p = stime
   numsolve32p = numsolve
   maxdyn     = .true.
   time0 = omp_get_wtime()
   Call bandb_solve
   timedp = omp_get_wtime() - time0
   stimedp = stime
   numsolvedp = numsolve
   write(11,*) 'Finished parallel with Maxpool=',maxpool
   write(*,*) 'Finished parallel with Maxpool=',maxpool
!
!  Clean up
!
   DO thread = 1, maxthread
      coi_error = max( coi_error, coi_free( cntvect1(thread)%CntVect ) )
   Enddo
   Deallocate( cntvect1 )
 
   write(*,*) 'Number of Binary variables in model=',nn
   write(*,*) 'Number of threads in parallel experiments=',maxthread
   write(*,*)
   write(*,"('Loop type  MaxPool  Total time    %      Solve Time    %    Number of Solves    Parallel Speedup')")
   write(*,*)
   write(*,1000) 'Sequential',' 1',time1,                100.d0, stime1,                  100.d0, numsolve1
   write(*,1000) 'Sequential',' 8',time8s ,time8s /time1*100.d0, stime8s ,stime8s /stime1*100.d0, numsolve8s
   write(*,1000) 'Sequential','16',time16s,time16s/time1*100.d0, stime16s,stime16s/stime1*100.d0, numsolve16s
   write(*,1000) 'Sequential','32',time32s,time32s/time1*100.d0, stime32s,stime32s/stime1*100.d0, numsolve32s
   write(*,1000) 'Sequential',' D',timeds ,timeds /time1*100.d0, stimeds ,stimeds /stime1*100.d0, numsolveds
   write(*,1000) 'Parallel  ',' 8',time8p ,time8p /time1*100.d0, stime8p ,stime8p /stime1*100.d0, numsolve8p , time8s/time8p
   write(*,1000) 'Parallel  ','16',time16p,time16p/time1*100.d0, stime16p,stime16p/stime1*100.d0, numsolve16p, time16s/time16p
   write(*,1000) 'Parallel  ','32',time32p,time32p/time1*100.d0, stime32p,stime32p/stime1*100.d0, numsolve32p, time32s/time32p
   write(*,1000) 'Parallel  ',' D',timedp ,timedp /time1*100.d0, stimedp ,stimedp /stime1*100.d0, numsolvedp , timeds/timedp
 
1000 Format(a,5x,a2,f12.3,f9.1,f12.3,f9.1,i12,f19.2)
 
   write(*,*)
   write(*,*) 'End of QP Branch & Bound example.'
 
#if defined (nocomp)
   call flog( "Successful Solve", 0 )
#else
   if ( numsolve8s .ne. numsolve8p ) then
      call flog( "Difference between number of nodes in sequential and parallel solve for MaxPool=8", 1 )
   else if ( numsolve16s .ne. numsolve16p ) then
      call flog( "Difference between number of nodes in sequential and parallel solve for MaxPool=16", 1 )
   else if ( numsolve32s .ne. numsolve32p ) then
      call flog( "Difference between number of nodes in sequential and parallel solve for MaxPool=32", 1 )
   else if ( numsolveds .ne. numsolvedp ) then
      call flog( "Difference between number of nodes in sequential and parallel solve for Dynamic MaxPool", 1 )
   else
      call flog( "Successful Solve", 0 )
   endif
#endif
 
End Program qp_threads
 
 
!>  Sets runtime options
!!
!! @include{doc} option_params.dox
Integer Function bb_option( ncall, rval, ival, lval, usrmem, name )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: BB_Option
#endif
   integer  ncall, ival, lval
   character(Len=*) :: name
   real*8   rval
   real*8   usrmem(*)                                  ! optional user memory
 
   Select Case (ncall)
   case (1)
      name = 'locon3'
      ival = 0
   case (2)
      name = 'lotria'
      ival = 0
   case default
      name = ' '
   end Select
   bb_option = 0
 
end Function bb_option
 
Subroutine bandb_solve
 
   Use bb_data
   Use qp_data
   implicit none
   Integer i
   Type (BBNode), Pointer :: Node
!
!  Initialize the first node for B&B tree
!
   numsolve    = 0
   numintfound = 0
   Nullify(freenode)
   Nullify(actnode)
   intfound  = .false.
   totnodes  = 0 ! Total number of nodes in use
   actnodes  = 0 ! Number of active nodes
   freenodes = 0 ! Number of free nodes in a free list
   stime     = 0 ! Reset Solve loop timer
!
   Call get_freenode(node) ! Create a root node with no binaries fixed
   node%Depth     = 0;
   do i = 1, nvar
      node%Values(i) = 0.d0
   enddo
   do i = 1, nbin
      node%Fixed(i) = -1
   enddo
   node%Object = 1.e29 ! We need a value for the first selection before
                       ! the node has been solved
   Call add_activenode(node)
   do while ( actnodes > 0 )
      write(*,*) 'Creating WorkPool. Active nodes=',actnodes
      Call create_workpool
      write(*,*) 'WorkPool with ',workpool,' nodes selected. Remaining Active nodes=',actnodes
      if ( workpool > 0 ) then
         Call solve_workpool
      endif
   Enddo
   Call deallocate_nodes ! Clean up again in the node list
 
   if ( intfound ) then
      write(*,*)
      write(*,*) 'An integer solution with objective',bestint,' has been found.'
      write(*,*)
      write(*,*) 'Solution values:'
      do i = 1, 2*nn
         write(*,"(I5,1pe20.12)") i, solution(i)
      enddo
   else
      write(*,*) 'No integer solution has been found'
   endif
   write(*,*) 'Number of Solves:           ', numsolve
   write(*,*) 'Number of Integer Solutions:', numintfound
   write(*,*)
 
End Subroutine bandb_solve
!
! ============================================================================
!  Define information about the model:
!
 
!>  Define information about the model
!!
!! @include{doc} readMatrix_params.dox
Integer Function qp_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 :: QP_ReadMatrix
#endif
   Use bb_data
   Use qp_data
   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
   Type (bbnode) :: usrmem
   integer :: i, j
!
!  Information about Variables:
!  Default: Lower = -Inf, Curr = 0, and Upper = +inf.
!  Default: the status information in Vsta is not used.
!
!  Nondefault: Lower bound = 0 for all variables
!
   do i = 1, nvar
     lower(i) = 0.0d0
   enddo
   do i = 1, nn
      upper(nn+i) = 1.d0
   enddo
   do i = 1, nvar
      curr(i) = usrmem%Values(i)
   enddo
!
!  Fix the binary variables if indicated by Fixed
!
   do i = 1, nbin
      if ( usrmem%fixed(i) == 0 ) then
         j        = binmap(i)
         lower(j) = 0.d0
         upper(j) = 0.d0
      elseif ( usrmem%fixed(i) == 1 ) then
         j        = binmap(i)
         lower(j) = 1.d0
         upper(j) = 1.d0
      endif
   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 (default) and type Non binding
!
   type(1) = 3
!
!  Constraint 2 (Sum of all X-vars = 1)
!  Rhs = 1 and type Equality
!
   rhs(2)  = 1.d0
   type(2) = 0
!
!  Constraint 3 (average of all mean = target)
!  Rhs = Target and type Eguality
!
   rhs(3)  = target
   type(3) = 0
!
!  Constraint 4 (sum of all Bin-vars =L= Maxbin)
!  Rhs = Maxbin and type L
!
   rhs(4)  = maxbin
   type(4) = 2
!
!  Next N constraints ( x(i) - bin(i) =L= 0 )
!  Rhs = 0 (default) and type L
!
   do i = 1,nn
      type(4+i) = 2
   enddo
!
!  Information about the Jacobian.
!
!  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)
!
   j = 1 ! counter for current nonzero
   do i = 1, nn
      colsta(i) = j
      rowno(j)  = 1       ! Objective - nonlinear
      nlflag(j) = 1
      j         = j + 1
      rowno(j)  = 2       ! sum(x) = 1
      value(j)  = 1.d0
      nlflag(j) = 0
      j         = j + 1
      rowno(j)  = 3       ! sum(mean*x) = target
      value(j)  = mean(i)
      nlflag(j) = 0
      j         = j + 1
      rowno(j)  = 4+i     ! x(i) - bin(i) =l= 0
      value(j)  = 1.d0
      nlflag(j) = 0
      j         = j + 1
   enddo
   do i = 1, nn
      colsta(nn+i) = j
      rowno(j)  = 4       ! sum(bin) =l= Maxbin
      value(j)  = 1.d0
      nlflag(j) = 0
      j         = j + 1
      rowno(j)  = 4+i     ! x(i) - bin(i) =l= 0
      value(j)  = -1.d0
      nlflag(j) = 0
      j         = j + 1
   enddo
   colsta(2*nn+1) = j
 
  qp_readmatrix = 0 ! Return value means OK
 
end Function qp_readmatrix
!
!==========================================================================
!  Compute nonlinear terms and non-constant Jacobian elements
!
 
!>  Compute nonlinear terms and non-constant Jacobian elements
!!
!! @include{doc} fdeval_params.dox
Integer Function qp_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 :: QP_FDEval
#endif
   Use bb_data
   Use qp_data
   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 !
   Type (bbnode) :: usrmem
   Integer :: i, j
   real*8  :: sum
!
!  Row 1: The objective function is nonlinear
!
   if ( rowno .eq. 1 ) then
!
!  Mode = 1 or 3. Function value: G = x * Q * x / 2
!
      if ( mode .eq. 1 .or. mode .eq. 3 ) then
         sum = 0.d0
         do i = 1, nn
            do j = 1, nn
               sum = sum + x(i)*q(i,j)*x(j)
            enddo
         enddo
         g    = sum / 2.d0
      endif
!
!  Mode = 2 or 3: Derivative values: Q*x
!
      if ( mode .eq. 2 .or. mode .eq. 3 ) then
         do i = 1, nn
            sum = 0
            do j = 1, nn
               sum = sum + q(i,j) * x(j)
            enddo
            jac(i) = sum
         enddo
      endif
      qp_fdeval = 0
   Else
!
!  Other rows are linear and will not be called.
!
      qp_fdeval = 1
   endif
   qp_fdeval = 0
 
end Function qp_fdeval
 
 
!>  Computes the second derivative of a constraint in a direction
!!
!! @include{doc} 2DDir_params.dox
Integer Function qp_2ddir( X, DX, D2G, ROWNO, JCNM, NODRV, N, NJ, THREAD, USRMEM )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_2DDir
#endif
   Use qp_data
   Use bb_data
   implicit none
   INTEGER, Intent(IN)                       :: rowno, n, nj, thread
   INTEGER, Intent(IN OUT)                   :: nodrv
   INTEGER, Dimension(NJ), Intent(IN)        :: jcnm
   real*8,  Dimension(N),  Intent(IN)        :: x
   real*8,  Dimension(N),  Intent(IN)        :: dx
   real*8,  Dimension(N),  Intent(OUT)       :: d2g
   Type (bbnode) :: usrmem
   integer :: i,j
   real(8) :: sum
!
!  Is only called for Irow = 1, the nonlinear objective
!  Return H*Dx = Q*Dx in D2G
!
   if ( rowno .eq. 1 ) then
      do i = 1, nn
         sum = 0.d0
         do j = 1, nn
            sum = sum + q(i,j) * dx(j)
         enddo
         d2g(i) = sum
      enddo
      qp_2ddir = 0
   else
      qp_2ddir = 1
   endif
 
End Function qp_2ddir
 
 
!>  Specify the structure of the Lagrangian of the Hessian
!!
!! @include{doc} 2DLagrStr_params.dox
Integer Function qp_2dlagrstr( ROWNO, COLNO, NODRV, N, M, NHESS, USRMEM )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_2DLagrStr
#endif
   Use qp_data
   Use bb_data
   implicit none
   INTEGER n, m, nhess, nodrv
   INTEGER rowno(nhess), colno(nhess)
   Type (bbnode) :: usrmem
 
   Integer :: i,j,k
 
   k = 1
   Do i = 1, nn
      do j = i,nn
         rowno(k) = j
         colno(k) = i
         k        = k + 1
      Enddo
   Enddo
   qp_2dlagrstr = 0
 
End function qp_2dlagrstr
 
 
!>  Compute the Lagrangian of the Hessian
!!
!! @include{doc} 2DLagrVal_params.dox
Integer Function qp_2dlagrval( X, U, ROWNO, COLNO, VALUE, NODRV, N, M, NHESS, USRMEM )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: QP_2DLagrVal
#endif
   Use qp_data
   Use bb_data
   implicit none
   INTEGER n, m, nhess, nodrv
   real*8  x(n), u(m), value(nhess)
   INTEGER rowno(nhess), colno(nhess)
   Type (bbnode) :: usrmem
 
   Integer :: i,j,k
 
   k = 1
   Do i = 1, nn
      do j = i,nn
         value(k) = q(i,j)*u(1)
         k        = k + 1
      Enddo
   Enddo
   qp_2dlagrval = 0
 
End function qp_2dlagrval
!
!  Toutines for Message, Status, Solution, and ErrMsg called
!  BB_Message, BB_Status, BB_Solution, and BB_ErrMsg.
!  They are very simple and we avoid written output since overlapping
!  threads can give unreadable output.
!
Integer Function bb_status( MODSTA, SOLSTA, ITER, OBJVAL, USRMEM )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: BB_Status
#endif
!
!  Simple implementation in which we just store the status and
!  objective value in the current node.
!
   Use bb_data
   IMPLICIT NONE
   INTEGER, Intent(IN)        :: modsta, solsta, iter
   real*8,  Intent(IN)        :: objval
   Type (bbnode) :: usrmem
 
   usrmem%Solstat = solsta
   usrmem%Modstat = modsta
   usrmem%Object  = objval
   usrmem%Node_Stat = 1
 
   bb_status = 0
   return
 
END Function bb_status
 
Integer Function bb_solution( XVAL, XMAR, XBAS, XSTA, YVAL, YMAR, YBAS, YSTA, N, M, USRMEM )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: BB_Solution
#endif
!
!  Store the solution in the Values field of the solved node.
!  Also check for integrality. If done here it is inside the
!  parallel loop.
!
   Use bb_data
   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
   Type (bbnode) :: usrmem
 
   Integer :: numintvar, i, j
   Real(8) :: dist
 
   DO i = 1, n
      usrmem%Values(i) = xval(i)
   ENDDO
   bb_solution = 0
 
   numintvar = 0
   do i = 1, nn
      j = binmap(i)
      dist = min( usrmem%values(j), 1.d0-usrmem%values(j) )
      if ( dist <= bintol ) numintvar = numintvar + 1
   enddo
   usrmem%IsInteger = (numintvar == nn)
 
END Function bb_solution
 
Integer Function bb_message( SMSG, DMSG, NMSG, LLEN, USRMEM, MSGV )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: BB_Message
#endif
!
!  Simple implementation in which does nothing. Messages from different
!  threads can give a confusing output.
!
   Use bb_data
   Use proginfop
   IMPLICIT NONE
 
   INTEGER,            Intent(IN)                 :: smsg, dmsg, nmsg
   INTEGER,            Intent(IN) ,  Dimension(*) :: llen
   CHARACTER(Len=133), Intent(IN),   Dimension(*) :: msgv
   Type (bbnode)                                  :: usrmem
 
   Integer :: i, thread, num
 
   if ( debug ) then
!
!  write to document file
!
      thread = usrmem%Sthread
      num    = usrmem%Number
      do i = 1, dmsg
         write(11+thread,"(i4,': ',(A))") num,msgv(i)(1:llen(i))
      enddo
   endif
   bb_message = 0
   return
 
END Function bb_message
 
Integer Function bb_errmsg( ROWNO, COLNO, POSNO, MSGLEN, USRMEM, MSG )
#ifdef dec_directives_win32
!DEC$ ATTRIBUTES STDCALL, REFERENCE, NOMIXED_STR_LEN_ARG :: BB_ErrMsg
#endif
!
!  Simple implementation in which does nothing. Messages from different
!  threads can give a confusing output.
!
   Use bb_data
   IMPLICIT NONE
 
   INTEGER,          Intent(IN)                    :: rowno, colno, posno, msglen
   CHARACTER(Len=*), Intent(IN)                    :: msg
   Type (bbnode)                                   :: usrmem
   bb_errmsg = 0
   return
 
END Function bb_errmsg