22 Integer,
Parameter ::
nn = 32
26 Integer,
Parameter ::
nobs = 60
27 Real(8),
dimension(NN,NN) ::
q
28 Real(8),
dimension(nn) ::
mean
31 Integer,
Dimension(Nbin) :: binmap
46 real(8),
dimension(NVar) :: values
47 Integer,
dimension(Nbin) :: fixed
48 Type (
bbnode),
Pointer :: prvnode, nxtnode
68 Real(8),
Parameter :: bintol = 1.d-7
78 Integer :: numintfound
80 Real(8),
dimension(NVar) :: solution
100 Integer,
Dimension(:),
Pointer :: cntvect
114 Type (BBNode),
Pointer :: Node
116 if (
associated(node%prvnode) )
then
117 node%prvnode%nxtnode => node%nxtnode
119 actnode => node%nxtnode
121 if (
associated(node%nxtnode) )
then
122 node%nxtnode%prvnode => node%prvnode
124 actnodes = actnodes - 1
133 Type (BBNode),
Pointer :: Node
135 node%nxtnode => freenode
137 freenodes = freenodes + 1
147 Type (BBNode),
Pointer :: Node
149 if (
associated(freenode) )
then
151 freenode => freenode%Nxtnode
152 freenodes = freenodes - 1
155 Nullify(node%NxtNode);
156 Nullify(node%PrvNode);
157 totnodes = totnodes + 1
158 node%Number = totnodes
173 Type (BBNode),
Pointer :: Node, Next
176 do while (
Associated(node) )
190 Type (BBNode),
pointer :: Node
192 actnodes = actnodes + 1
193 if (
associated(actnode) )
then
194 node%NxtNode => actnode
195 actnode%PrvNode => node;
197 Nullify(node%NxtNode);
199 Nullify(node%PrvNode);
209 Type (BBNode),
pointer :: Node
211 workpool = workpool + 1
213 work_node(workpool)%Bnode => node
220 Type (BBNode),
pointer :: Node
222 If (
Associated(node) )
Then
223 write(*,*)
'List_Node: Number=',node%Number,
' Node_Stat=',node%node_stat
225 write(*,*)
'List_Node: Node not associated.'
232 Type (BBNode),
Pointer :: Node
237 do while (
Associated(node) )
238 write(*,
"('List_nodes: Node number=',i4,' Node_Stat=',I2,' Object=',1pe15.8)") node%number, node%node_stat, node%object
242 write(*,*)
'List_nodes: Total=',n
254 Type (BBNode),
Pointer :: Node
256 Real(8) :: Dist, Maxdist
262 dist = min( node%values(j), 1.d0-node%values(j) )
263 if ( dist > maxdist )
then
278 Type (BBNode),
Pointer :: Node, Curr, Next
284 write(*,*)
'Create_WorkPool'
292 do while ( workpool < maxpool )
300 do while (
associated(curr))
302 if ( curr%Object < cutoff )
then
308 if ( .not.
associated(node) )
then
311 if ( node%Node_Stat == 0 )
then
317 write(*,
"('Unsolved node',i5,' selected. Obj=',e24.17)") node%Number,node%Object
334 next%Values = node%Values
335 next%Fixed = node%Fixed
336 next%Object = node%Object
344 node%Depth = node%Depth + 1
350 next%Depth = node%Depth
357 write(*,
"('New-Down node',i5,' selected. Obj=',e24.17)") node%Number,node%Object
362 if ( workpool < maxpool )
then
364 write(*,
"('New-Up Next',i5,' selected. Obj=',e24.17)") next%Number,next%Object
384 Integer,
Dimension(:),
pointer :: CntVect
386 Type (BBNode),
Pointer :: Node, Curr, Next
393 INTEGER,
External :: COIDEF_UsrMem
400 time0 = omp_get_wtime()
401 if ( workpool .le. 1 .OR. .NOT. parallel )
Then
405 node => work_node(nw)%Bnode
407 write(16,*)
'Starting to solve node',node%number,
' using thread',thread
410 coi_error = coidef_usrmem( cntvect, node )
411 If ( coi_error .ne. 0 )
THEN
413 write(*,*)
'**** Fatal Error while loading CONOPT Callback routines.'
417 node%Sthread = thread
419 If ( coi_error .ne. 0 )
THEN
421 write(*,*)
'**** Fatal Error while solving.'
431 node => work_node(nw)%Bnode
432 thread = 1+omp_get_thread_num()
434 write(16,*)
'Starting to solve node',node%number,
' using thread',thread
437 coi_error = coidef_usrmem( cntvect, node )
438 If ( coi_error .ne. 0 )
THEN
440 write(*,*)
'**** Fatal Error while loading CONOPT Callback routines.'
442 call flog(
"Errors encountered during setup for Solve.", 1 )
444 node%Sthread = thread
446 If ( coi_error .ne. 0 )
THEN
448 write(*,*)
'**** Fatal Error while solving.'
450 call flog(
"Errors encountered during Solve.", 1 )
455 stime = stime + omp_get_wtime() - time0
456 numsolve = numsolve + workpool
464 node => work_node(nw)%Bnode
465 write(*,
'("Node=",I5," Solstat=",I3," Modstat=",I3," Object=",F24.17)') &
466 node%number, node%Solstat, node%modstat, node%Object
467 if ( node%Solstat .ne. 1 .or. node%Modstat .gt. 2 )
then
471 call flog(
"Solve did not return an normal solver status.", 1 )
472 Else if ( node%Modstat .gt. 2 )
then
476 write(*,*)
'Not optimal -- fathom'
478 else if ( intfound )
then
479 if ( node%object >= bestint )
then
483 write(*,*)
'Worse than optimal -- fathom'
486 if ( node%IsInteger )
then
487 write(*,*)
'Integer Solution Stored'
488 numintfound = numintfound + 1
490 bestint = node%object
491 solution = node%Values
502 if ( node%IsInteger )
then
503 write(*,*)
'Integer Solution Stored'
504 numintfound = numintfound + 1
507 bestint = node%object
508 solution = node%Values
522 do while (
associated(curr))
524 if ( curr%Object > bestint )
then
525 write(*,
"('Node',I5,' with Object=',F18.12,' fathomed by new objective')") curr%Number, curr%Object
543 Integer,
save :: seed = 12359
545 seed = mod(seed*1027+25,1048576)
546 rndx = float(seed)/float(1048576)
558 Real,
External :: Rndx
560 Real(8),
dimension(NN,Nobs) :: Xobs
564 mean(i) =
Target + (rndx()*4.d0-2.d0)
578 xobs(i,k) = 100.d0*rndx()
579 sum = sum + xobs(i,k)
583 xobs(i,k) = xobs(i,k) - sum
590 sum = sum + xobs(i,k) * xobs(j,k)
592 q(i,j) = sum / (
nobs-1)
653 Integer,
Dimension(:),
Pointer :: cntvect
655 Real(8) time0, time1, stime1
656 Real(8) time8s, time16s, time32s, timeds
657 Real(8) time8p, time16p, time32p, timedp
658 Real(8) stime8s, stime16s, stime32s, stimeds
659 Real(8) stime8p, stime16p, stime32p, stimedp
661 Integer numsolve8s, numsolve16s, numsolve32s, numsolveds
662 Integer numsolve8p, numsolve16p, numsolve32p, numsolvedp
666 write(11,*)
'Initial MaxThread=',
maxthread
667 write(*,*)
'Initial MaxThread=',
maxthread
671 write(11,*)
'Revised MaxThread=',
maxthread
672 write(*,*)
'Revised MaxThread=',
maxthread
676 Open(12,file=
'threadA1.txt',status=
'unknown')
677 Open(16,file=
'nodesA.txt',status=
'unknown')
688 coi_error = coi_createfort( cntvect )
717#if defined(LICENSE_INT_1) && defined(LICENSE_INT_2) && defined(LICENSE_INT_3) && defined(LICENSE_TEXT)
718 coi_error = max( coi_error,
coidef_license( cntvect, license_int_1, license_int_2, license_int_3, license_text) )
721 If ( coi_error .ne. 0 )
THEN
723 write(*,*)
'**** Fatal Error while loading CONOPT Callback routines.'
725 call flog(
"Skipping First Solve due to setup errors", 1 )
727 write(11,*)
'Finished Setup of model'
728 write(*,*)
'Finished Setup of model'
735 time0 = omp_get_wtime()
737 time1 = omp_get_wtime() - time0
740 write(11,*)
'Finished sequential with Maxpool=',maxpool
741 write(*,*)
'Finished sequential with Maxpool=',maxpool
744 Open(12,file=
'threadB1.txt',status=
'unknown')
745 Open(16,file=
'nodesB.txt',status=
'unknown')
751 time0 = omp_get_wtime()
753 time8s = omp_get_wtime() - time0
755 numsolve8s = numsolve
756 write(11,*)
'Finished sequential with Maxpool=',maxpool
757 write(*,*)
'Finished sequential with Maxpool=',maxpool
762 time0 = omp_get_wtime()
764 time16s = omp_get_wtime() - time0
766 numsolve16s = numsolve
767 write(11,*)
'Finished sequential with Maxpool=',maxpool
768 write(*,*)
'Finished sequential with Maxpool=',maxpool
773 time0 = omp_get_wtime()
775 time32s = omp_get_wtime() - time0
777 numsolve32s = numsolve
778 write(11,*)
'Finished sequential with Maxpool=',maxpool
779 write(*,*)
'Finished sequential with Maxpool=',maxpool
784 time0 = omp_get_wtime()
786 timeds = omp_get_wtime() - time0
788 numsolveds = numsolve
789 write(11,*)
'Finished sequential with Maxpool=',maxpool
790 write(*,*)
'Finished sequential with Maxpool=',maxpool
793 Open(12,file=
'threadC1.txt',status=
'unknown')
794 Open(13,file=
'threadC2.txt',status=
'unknown')
795 Open(14,file=
'threadC3.txt',status=
'unknown')
796 Open(15,file=
'threadC4.txt',status=
'unknown')
797 Open(16,file=
'nodesC.txt',status=
'unknown')
805 time0 = omp_get_wtime()
807 time8p = omp_get_wtime() - time0
809 numsolve8p = numsolve
810 write(11,*)
'Finished parallel with Maxpool=',maxpool
811 write(*,*)
'Finished parallel with Maxpool=',maxpool
813 time0 = omp_get_wtime()
815 time16p = omp_get_wtime() - time0
817 numsolve16p = numsolve
818 write(11,*)
'Finished parallel with Maxpool=',maxpool
819 write(*,*)
'Finished parallel with Maxpool=',maxpool
821 time0 = omp_get_wtime()
823 time32p = omp_get_wtime() - time0
825 numsolve32p = numsolve
827 time0 = omp_get_wtime()
829 timedp = omp_get_wtime() - time0
831 numsolvedp = numsolve
832 write(11,*)
'Finished parallel with Maxpool=',maxpool
833 write(*,*)
'Finished parallel with Maxpool=',maxpool
838 coi_error = max( coi_error, coi_free(
cntvect1(thread)%CntVect ) )
842 write(*,*)
'Number of Binary variables in model=',
nn
843 write(*,*)
'Number of threads in parallel experiments=',
maxthread
845 write(*,
"('Loop type MaxPool Total time % Solve Time % Number of Solves Parallel Speedup')")
847 write(*,1000)
'Sequential',
' 1',time1, 100.d0, stime1, 100.d0, numsolve1
848 write(*,1000)
'Sequential',
' 8',time8s ,time8s /time1*100.d0, stime8s ,stime8s /stime1*100.d0, numsolve8s
849 write(*,1000)
'Sequential',
'16',time16s,time16s/time1*100.d0, stime16s,stime16s/stime1*100.d0, numsolve16s
850 write(*,1000)
'Sequential',
'32',time32s,time32s/time1*100.d0, stime32s,stime32s/stime1*100.d0, numsolve32s
851 write(*,1000)
'Sequential',
' D',timeds ,timeds /time1*100.d0, stimeds ,stimeds /stime1*100.d0, numsolveds
852 write(*,1000)
'Parallel ',
' 8',time8p ,time8p /time1*100.d0, stime8p ,stime8p /stime1*100.d0, numsolve8p , time8s/time8p
853 write(*,1000)
'Parallel ',
'16',time16p,time16p/time1*100.d0, stime16p,stime16p/stime1*100.d0, numsolve16p, time16s/time16p
854 write(*,1000)
'Parallel ',
'32',time32p,time32p/time1*100.d0, stime32p,stime32p/stime1*100.d0, numsolve32p, time32s/time32p
855 write(*,1000)
'Parallel ',
' D',timedp ,timedp /time1*100.d0, stimedp ,stimedp /stime1*100.d0, numsolvedp , timeds/timedp
8571000
Format(a,5x,a2,f12.3,f9.1,f12.3,f9.1,i12,f19.2)
860 write(*,*)
'End of QP Branch & Bound example.'
863 call flog(
"Successful Solve", 0 )
865 if ( numsolve8s .ne. numsolve8p )
then
866 call flog(
"Difference between number of nodes in sequential and parallel solve for MaxPool=8", 1 )
867 else if ( numsolve16s .ne. numsolve16p )
then
868 call flog(
"Difference between number of nodes in sequential and parallel solve for MaxPool=16", 1 )
869 else if ( numsolve32s .ne. numsolve32p )
then
870 call flog(
"Difference between number of nodes in sequential and parallel solve for MaxPool=32", 1 )
871 else if ( numsolveds .ne. numsolvedp )
then
872 call flog(
"Difference between number of nodes in sequential and parallel solve for Dynamic MaxPool", 1 )
874 call flog(
"Successful Solve", 0 )
884Integer Function bb_option( ncall, rval, ival, lval, usrmem, name )
888 integer ncall, ival, lval
889 character(Len=*) :: name
913 Type (BBNode),
Pointer :: Node
930 node%Values(i) = 0.d0
938 do while ( actnodes > 0 )
939 write(*,*)
'Creating WorkPool. Active nodes=',actnodes
941 write(*,*)
'WorkPool with ',workpool,
' nodes selected. Remaining Active nodes=',actnodes
942 if ( workpool > 0 )
then
950 write(*,*)
'An integer solution with objective',bestint,
' has been found.'
952 write(*,*)
'Solution values:'
954 write(*,
"(I5,1pe20.12)") i, solution(i)
957 write(*,*)
'No integer solution has been found'
959 write(*,*)
'Number of Solves: ', numsolve
960 write(*,*)
'Number of Integer Solutions:', numintfound
973 colsta, rowno, value, nlflag, n, m, nz, &
981 integer,
intent (in) :: n
982 integer,
intent (in) :: m
983 integer,
intent (in) :: nz
984 real*8,
intent (in out),
dimension(n) :: lower
985 real*8,
intent (in out),
dimension(n) :: curr
986 real*8,
intent (in out),
dimension(n) :: upper
987 integer,
intent (in out),
dimension(n) :: vsta
989 integer,
intent (out),
dimension(m) ::
type
990 integer,
intent (in out),
dimension(m) :: esta
992 real*8,
intent (in out),
dimension(m) :: rhs
993 integer,
intent (in out),
dimension(n+1) :: colsta
994 integer,
intent (out),
dimension(nz) :: rowno
995 integer,
intent (in out),
dimension(nz) :: nlflag
996 real*8,
intent (in out),
dimension(nz) ::
value
1013 curr(i) = usrmem%Values(i)
1019 if ( usrmem%fixed(i) == 0 )
then
1023 elseif ( usrmem%fixed(i) == 1 )
then
1122Integer Function qp_fdeval( x, g, jac, rowno, jcnm, mode, ignerr, errcnt, &
1123 n, nz, thread, usrmem )
1130 integer,
intent (in) :: n
1131 integer,
intent (in) :: rowno
1132 integer,
intent (in) :: nz
1133 real*8,
intent (in),
dimension(n) :: x
1134 real*8,
intent (in out) :: g
1135 real*8,
intent (in out),
dimension(n) :: jac
1136 integer,
intent (in),
dimension(nz) :: jcnm
1138 integer,
intent (in) :: mode
1140 integer,
intent (in) :: ignerr
1142 integer,
intent (in out) :: errcnt
1144 integer,
intent (in) :: thread
1151 if ( rowno .eq. 1 )
then
1155 if ( mode .eq. 1 .or. mode .eq. 3 )
then
1159 sum = sum + x(i)*
q(i,j)*x(j)
1167 if ( mode .eq. 2 .or. mode .eq. 3 )
then
1171 sum = sum +
q(i,j) * x(j)
1191Integer Function qp_2ddir( X, DX, D2G, ROWNO, JCNM, NODRV, N, NJ, THREAD, USRMEM )
1198 INTEGER,
Intent(IN) :: rowno, n,
nj, thread
1199 INTEGER,
Intent(IN OUT) :: nodrv
1200 INTEGER,
Dimension(NJ),
Intent(IN) :: jcnm
1201 real*8,
Dimension(N),
Intent(IN) :: x
1202 real*8,
Dimension(N),
Intent(IN) :: dx
1203 real*8,
Dimension(N),
Intent(OUT) :: d2g
1211 if ( rowno .eq. 1 )
then
1215 sum = sum +
q(i,j) * dx(j)
1237 INTEGER n, m, nhess, nodrv
1238 INTEGER rowno(nhess), colno(nhess)
1259Integer Function qp_2dlagrval( X, U, ROWNO, COLNO, VALUE, NODRV, N, M, NHESS, USRMEM )
1266 INTEGER n, m, nhess, nodrv
1267 real*8 x(n), u(m), value(nhess)
1268 INTEGER rowno(nhess), colno(nhess)
1276 value(k) =
q(i,j)*u(1)
1289Integer Function bb_status( MODSTA, SOLSTA, ITER, OBJVAL, USRMEM )
1299 INTEGER,
Intent(IN) :: modsta, solsta, iter
1300 real*8,
Intent(IN) :: objval
1303 usrmem%Solstat = solsta
1304 usrmem%Modstat = modsta
1305 usrmem%Object = objval
1306 usrmem%Node_Stat = 1
1313Integer Function bb_solution( XVAL, XMAR, XBAS, XSTA, YVAL, YMAR, YBAS, YSTA, N, M, USRMEM )
1324 INTEGER,
Intent(IN) :: n, m
1325 INTEGER,
Intent(IN),
Dimension(N) :: xbas, xsta
1326 INTEGER,
Intent(IN),
Dimension(M) :: ybas, ysta
1327 real*8,
Intent(IN),
Dimension(N) :: xval, xmar
1328 real*8,
Intent(IN),
Dimension(M) :: yval, ymar
1331 Integer :: numintvar, i, j
1335 usrmem%Values(i) = xval(i)
1342 dist = min( usrmem%values(j), 1.d0-usrmem%values(j) )
1343 if ( dist <= bintol ) numintvar = numintvar + 1
1345 usrmem%IsInteger = (numintvar == nn)
1349Integer Function bb_message( SMSG, DMSG, NMSG, LLEN, USRMEM, MSGV )
1361 INTEGER,
Intent(IN) :: smsg, dmsg, nmsg
1362 INTEGER,
Intent(IN) ,
Dimension(*) :: llen
1363 CHARACTER(Len=133),
Intent(IN),
Dimension(*) :: msgv
1366 Integer :: i, thread, num
1372 thread = usrmem%Sthread
1375 write(11+thread,
"(i4,': ',(A))") num,msgv(i)(1:llen(i))
1383Integer Function bb_errmsg( ROWNO, COLNO, POSNO, MSGLEN, USRMEM, MSG )
1394 INTEGER,
Intent(IN) :: rowno, colno, posno, msglen
1395 CHARACTER(Len=*),
Intent(IN) :: msg
integer function coidef_fdeval(cntvect, coi_fdeval)
define callback routine for performing function and derivative evaluations.
integer function coidef_errmsg(cntvect, coi_errmsg)
define callback routine for returning error messages for row, column or Jacobian elements.
integer function coidef_message(cntvect, coi_message)
define callback routine for handling messages returned during the solution process.
integer function coidef_readmatrix(cntvect, coi_readmatrix)
define callback routine for providing the matrix data to CONOPT.
integer function coidef_status(cntvect, coi_status)
define callback routine for returning the completion status.
integer function coidef_solution(cntvect, coi_solution)
define callback routine for returning the final solution values.
integer function coidef_2dlagrstr(cntvect, coi_2dlagrstr)
define callback routine for providing the structure of the second derivatives of the Lagrangian.
integer function coidef_2ddir(cntvect, coi_2ddir)
define callback routine for computing the second derivative for a constraint in a direction.
integer function coidef_option(cntvect, coi_option)
define callback routine for defining runtime options.
integer function coidef_2dlagrval(cntvect, coi_2dlagrval)
define callback routine for computing the values of the second derivatives of the Lagrangian.
integer function coidef_license(cntvect, licint1, licint2, licint3, licstring)
define the License Information.
integer function coidef_numvar(cntvect, numvar)
defines the number of variables in the model.
integer function coidef_objcon(cntvect, objcon)
defines the Objective Constraint.
integer function coidef_numnz(cntvect, numnz)
defines the number of nonzero elements in the Jacobian.
integer function coidef_optdir(cntvect, optdir)
defines the Optimization Direction.
integer function coidef_numnlnz(cntvect, numnlnz)
defines the Number of Nonlinear Nonzeros.
integer function coidef_numhess(cntvect, numhess)
defines the Number of Hessian Nonzeros.
integer function coidef_numcon(cntvect, numcon)
defines the number of constraints in the model.
integer function coi_solve(cntvect)
method for starting the solving process of CONOPT.
void defdata()
Defines the data for the problem.
float rndx()
Defines a pseudo random number between 0 and 1.
int select_var(struct BBNode *node)
void add_activenode(struct BBNode *node)
void return_freenode(struct BBNode *node)
void add2workpool(struct BBNode *node)
struct BBNode * get_freenode()
void remove_activenode(struct BBNode *node)
void list_node(struct BBNode *node)
integer function bb_message(smsg, dmsg, nmsg, llen, usrmem, msgv)
integer function bb_solution(xval, xmar, xbas, xsta, yval, ymar, ybas, ysta, n, m, usrmem)
program qp_threads
Main program. A simple setup and call of CONOPT.
integer function bb_status(modsta, solsta, iter, objval, usrmem)
integer function bb_option(ncall, rval, ival, lval, usrmem, name)
Sets runtime options.
integer function bb_errmsg(rowno, colno, posno, msglen, usrmem, msg)
type(controlv), dimension(:), pointer cntvect1
subroutine deallocate_nodes
subroutine flog(msg, code)
integer, parameter maxbin
real(8), dimension(nn) mean
real(8), dimension(nn, nn) q
integer function qp_readmatrix(lower, curr, upper, vsta, type, rhs, esta, colsta, rowno, value, nlflag, n, m, nz, usrmem)
Define information about the model.
integer function qp_fdeval(x, g, jac, rowno, jcnm, mode, ignerr, errcnt, n, nz, thread, usrmem)
Compute nonlinear terms and non-constant Jacobian elements.
integer function qp_2ddir(x, dx, d2g, rowno, jcnm, nodrv, n, nj, thread, usrmem)
Computes the second derivative of a constraint in a direction.
integer function qp_2dlagrval(x, u, rowno, colno, value, nodrv, n, m, nhess, usrmem)
Compute the Lagrangian of the Hessian.
integer function qp_2dlagrstr(rowno, colno, nodrv, n, m, nhess, usrmem)
Specify the structure of the Lagrangian of the Hessian.