subroutine graph_arc_euler_circ ( nnode, nedge, inode, jnode, loop ) ! !******************************************************************************* ! !! GRAPH_ARC_EULER_CIRC finds an Euler circuit in an Eulerian graph. ! ! ! Discussion: ! ! An Euler circuit of a graph is a path that uses each edge exactly once. ! ! A graph is Eulerian if it has an Euler circuit. ! ! An Eulerian graph may have many circuits; this routine only finds one. ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 21 July 2000 ! ! Parameters: ! ! Input, integer NNODE, the number of nodes in the graph. ! ! Input, integer NEDGE, the number of edges in the graph. ! ! Input, integer INODE(NEDGE), JNODE(NEDGE), the two end nodes of each edge. ! ! Output, integer LOOP(NEDGE), the Euler circuit, as a series of nodes. ! implicit none ! integer nedge integer nnode ! logical copyon logical found integer i integer ibase integer iforwd integer inode(nedge) integer insert integer ipivot integer iwork1(nedge) integer iwork2(nedge) integer iwork3(nnode) integer iwork4(nnode) integer iwork5(nnode) integer iwork6(nnode) integer j integer jnode(nedge) integer k integer l integer locbas integer loop(nedge) integer nbreak integer ncopy integer numarc integer numnode ! ! The number of times each node has been visited begins at 0. ! iwork3(1:nnode) = 0 loop(1:nedge) = 0 iwork1(1:nedge) = 0 iwork2(1:nedge) = 0 ! ! Begin the Euler circuit with the edge INODE(1), JNODE(1). ! numarc = 1 iwork2(1) = 1 numnode = 1 i = inode(1) iwork1(numnode) = i iwork3(i) = 1 numnode = numnode + 1 j = jnode(1) iwork1(numnode) = j iwork3(j) = 1 ibase = j nbreak = 0 ! ! Look for the next arc. ! 30 continue do i = 2, nedge if ( iwork2(i) == 0 ) then if ( inode(i) == ibase ) then found = .true. ibase = jnode(i) else if ( jnode(i) == ibase ) then found = .true. ibase = inode(i) else found = .false. end if if ( found ) then iwork2(i) = 1 numarc = numarc + 1 numnode = numnode + 1 if ( numnode <= nedge ) then iwork1(numnode) = ibase end if iwork3(ibase) = 1 go to 30 end if end if end do ! ! A cycle has been found. ! if ( nbreak > 0 ) then numnode = numnode - 1 iwork5(nbreak) = numnode end if if ( numarc < nedge ) then iwork1(numnode) = ibase ! ! Find a node in the current Euler circuit. ! do i = 2, nedge if ( iwork2(i) == 0 ) then if ( iwork3(inode(i)) /= 0 ) then found = .true. j = inode(i) k = jnode(i) else if ( iwork3(jnode(i)) /= 0 ) then found = .true. j = jnode(i) k = inode(i) else found = .false. end if ! ! Identify a path which will be added to the circuit. ! if ( found ) then nbreak = nbreak + 1 iwork6(nbreak) = j ibase = k iwork3(k) = 1 numnode = numnode + 1 iwork4(nbreak) = numnode iwork1(numnode) = ibase iwork2(i) = 1 numarc = numarc + 1 go to 30 end if end if end do end if ! ! Form the Euler circuit. ! if ( nbreak == 0 ) then numnode = numnode - 1 loop(1:numnode) = iwork1(1:numnode) return end if insert = 1 ipivot = iwork6(insert) iforwd = 0 70 continue ncopy = 1 ibase = iwork1(1) locbas = 1 loop(ncopy) = ibase ! ! A path identified before is added to the circuit. ! 80 continue if ( ibase == ipivot ) then j = iwork4(insert) + iforwd k = iwork5(insert) + iforwd do l = j, k ncopy = ncopy + 1 loop(ncopy) = iwork1(l) iwork1(l) = 0 end do ncopy = ncopy + 1 ! ! Add the intersecting node to the circuit. ! loop(ncopy) = ibase iforwd = iforwd + 1 if ( ncopy < numnode ) then do if ( ncopy >= nedge ) then exit end if locbas = locbas + 1 if ( locbas >= nedge ) then exit end if ibase = iwork1(locbas) if ( ibase /= 0 ) then ncopy = ncopy + 1 loop(ncopy) = ibase end if end do end if else ncopy = ncopy + 1 if ( ncopy <= numnode ) then locbas = locbas + 1 ibase = iwork1(locbas) loop(ncopy) = ibase go to 80 end if end if ! ! Check if more paths are to be added to the circuit. ! copyon = .false. insert = insert + 1 if ( insert <= nbreak ) then copyon = .true. ipivot = iwork6(insert) end if if ( copyon ) then iwork1(1:nedge) = loop(1:nedge) go to 70 end if return end subroutine graph_arc_min_path ( nnode, nedge, inode, jnode, arcost, & istart, last, num_path, ispath, xlen ) ! !******************************************************************************* ! !! GRAPH_ARC_MIN_PATH finds the shortest path between two nodes. ! ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! ! Input, integer NNODE, the number of nodes in the graph. ! ! Input, integer NEDGE, the number of edges in the graph. ! ! Input, integer INODE(NEDGE), JNODE(NEDGE), the edges of the graph, ! describe by pairs of nodes. ! ! Input, real ARCOST(NEDGE), the distance or cost of each edge. ! ! Input, integer ISTART, LAST, are the two nodes between which a ! shortest path is desired. ! ! Output, integer NUM_PATH, the number of nodes in the shortest path. ! NUM_PATH is zero if no path could be found. ! ! Output, integer ISPATH(NNODE), lists the nodes in the shortest path, ! from ISPATH(1) to ISPATH(NUM_PATH). ! ! Output, real XLEN, the length of the shortest path from ISTART to LAST. ! implicit none ! integer nedge integer nnode ! real arcost(nedge) real d integer i integer ic integer ient logical ifin integer ij integer inode(nedge) integer ispath(nnode) integer istart logical iwork1(nnode) integer iwork2(nnode) integer iwork3(nedge) integer j integer jnode(nedge) integer k integer l integer last integer num_path real wk4(nnode) real xlen ! wk4(1:nnode) = huge ( xlen ) iwork1(1:nnode) = .true. iwork2(1:nnode) = 0 wk4(istart) = 0.0E+00 i = istart iwork1(istart) = .false. xlen = 0 ! ! For each forward arc originating at node I calculate ! the length of the path to node I. ! 10 continue ic = 0 do k = 1, nedge if ( inode(k) == i ) then ic = ic + 1 iwork3(ic) = k ispath(ic) = jnode(k) end if if ( jnode(k) == i ) then ic = ic + 1 iwork3(ic) = k ispath(ic) = inode(k) end if end do if ( ic > 0 ) then do l = 1, ic k = iwork3(l) j = ispath(l) if ( iwork1(j) ) then d = wk4(i) + arcost(k) if ( d < wk4(j) ) then wk4(j) = d iwork2(j) = k end if end if end do end if ! ! Find the minimum potential. ! d = huge ( d ) ient = 0 ifin = .false. do i = 1, nnode if ( iwork1(i) ) then ifin = .true. if ( wk4(i) < d ) then d = wk4(i) ient = i end if end if end do ! ! Include the node in the current path. ! if ( d < huge ( d ) ) then iwork1(ient) = .false. if ( ient /= last ) then i = ient go to 10 end if else if ( ifin ) then num_path = 0 return end if end if ij = last num_path = 1 ispath(1) = last do k = iwork2(ij) if ( inode(k) == ij ) then ij = jnode(k) else ij = inode(k) end if num_path = num_path + 1 ispath(num_path) = ij if ( ij == istart ) then exit end if end do l = num_path / 2 j = num_path do i = 1, l call i_swap ( ispath(i), ispath(j) ) j = j - 1 end do xlen = wk4(last) return end subroutine graph_arc_min_span_tree ( nnode, nedge, inode, jnode, cost, & itree, jtree, tree_cost ) ! !******************************************************************************* ! !! GRAPH_ARC_MIN_SPAN_TREE finds a minimum spanning tree of a graph. ! ! ! Discussion: ! ! The input graph is represented by a list of edges. ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 21 July 2000 ! ! Parameters: ! ! Input, integer NNODE, the number of nodes in the graph. ! ! Input, integer NEDGE, the number of edges in the graph. ! ! Input, integer INODE(NEDGE), JNODE(NEDGE), the start and end nodes ! of the edges. ! ! Input, real COST(NEDGE), the cost or length of each edge. ! ! Output, integer ITREE(NNODE-1), JTREE(NNODE-1), the pairs of nodes that ! form the edges of the spanning tree. ! ! Output, real TREE_COST, the total cost or length of the spanning tree. ! implicit none ! integer nedge integer nnode ! integer best real cost(nedge) real d logical free(nnode) integer i integer ic integer ij integer inode(nedge) integer itr integer itree(nnode-1) integer iwork1(nnode) integer iwork2(nnode) integer iwork4(nedge) integer iwork5(nedge) integer j integer jnode(nedge) integer jtree(nnode-1) integer jj integer k integer kk integer l real tree_cost real wk6(nnode) ! wk6(1:nnode) = huge ( tree_cost ) free(1:nnode) = .true. iwork1(1:nnode) = 0 iwork2(1:nnode) = 0 itree(1:nnode-1) = 0 jtree(1:nnode-1) = 0 ! ! Find the first non-zero arc. ! do ij = 1, nedge if ( inode(ij) /= 0 ) then i = inode(ij) exit end if end do wk6(i) = 0.0E+00 free(i) = .false. tree_cost = 0.0E+00 do jj = 1, nnode - 1 wk6(1:nnode) = huge ( tree_cost ) do i = 1, nnode ! ! For each forward arc originating at node I ! calculate the length of the path to node I. ! if ( .not. free(i) ) then ic = 0 do k = 1, nedge if ( inode(k) == i ) then ic = ic + 1 iwork5(ic) = k iwork4(ic) = jnode(k) end if if ( jnode(k) == i ) then ic = ic + 1 iwork5(ic) = k iwork4(ic) = inode(k) end if end do if ( ic > 0 ) then do l = 1, ic k = iwork5(l) j = iwork4(l) if ( free(j) ) then d = tree_cost + cost(k) if ( d < wk6(j) ) then wk6(j) = d iwork1(j) = i iwork2(j) = k end if end if end do end if end if end do ! ! Identify the free node of minimum potential. ! d = huge ( d ) best = 0 do i = 1, nnode if ( free(i) ) then if ( wk6(i) < d ) then d = wk6(i) best = i itr = iwork1(i) kk = iwork2(i) end if end if end do ! ! Add that node to the tree. ! if ( d < huge ( d ) ) then free(best) = .false. tree_cost = tree_cost + cost(kk) itree(jj) = itr jtree(jj) = best end if end do return end subroutine graph_arc_print ( nedge, inode, jnode, title ) ! !******************************************************************************* ! !! GRAPH_ARC_PRINT prints out a graph from an edge list. ! ! ! Modified: ! ! 04 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NEDGE, the number of edges. ! ! Input, integer INODE(NEDGE), JNODE(NEDGE), the beginning and end ! nodes of the edges. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none ! integer nedge ! integer i integer inode(nedge) integer jnode(nedge) character ( len = * ) title ! if ( len_trim ( title ) /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if write ( *, '(a)' ) ' ' do i = 1, nedge write ( *, '(i6,4x,2i6)' ) i, inode(i), jnode(i) end do return end subroutine graph_arc_weight_print ( nedge, inode, jnode, wnode, title ) ! !******************************************************************************* ! !! GRAPH_ARC_WEIGHT_PRINT prints out a weighted graph from an edge list. ! ! ! Modified: ! ! 04 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NEDGE, the number of edges. ! ! Input, integer INODE(NEDGE), JNODE(NEDGE), the beginning and end ! nodes of the edges. ! ! Input, real WNODE(NEDGE), the weights of the edges. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none ! integer nedge ! integer i integer inode(nedge) integer jnode(nedge) character ( len = * ) title real wnode(nedge) ! if ( len_trim ( title ) /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if write ( *, '(a)' ) ' ' do i = 1, nedge write ( *, '(i6,4x,2i6,g14.6)' ) i, inode(i), jnode(i), wnode(i) end do return end subroutine graph_dist_min_span_tree3 ( lda, nnode, dist, inode, jnode ) ! !******************************************************************************* ! !! GRAPH_DIST_MIN_SPAN_TREE3 finds a minimum spanning tree. ! ! ! Discussion: ! ! The input graph is represented by a distance matrix. ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 05 April 2001 ! ! Parameters: ! ! Input, integer LDA, the leading dimension of DIST, which should be ! at least NNODE. ! ! Input, integer NNODE, the number of nodes. ! ! Input, real DIST(LDA,NNODE), an NNODE by NNODE distance matrix. ! DIST(I,J) is the distance from node I to node J. The matrix ! should be symmetric. If there is no arc from node I to node J, ! set DIST(I,J) = HUGE(1.0). ! ! Output, integer INODE(NNODE), JNODE(NNODE); entries 1 through NNODE-1 ! describe the edges of the spanning tree as pairs of nodes. ! implicit none ! integer lda integer nnode ! real d real dist(lda,nnode) integer i integer ient integer ij integer inode(nnode) integer itr integer iwork1(nnode) integer iwork2(nnode) integer j integer jj integer jnode(nnode) integer k integer kj integer kk4 real tree_length real work(nnode) ! work(1:nnode) = huge ( d ) iwork1(1:nnode) = 0 iwork2(1:nnode) = 0 ! ! Find the first non-zero arc. ! i = 0 do ij = 1, nnode do kj = 1, nnode if ( dist(ij,kj) < huge ( d ) ) then i = ij exit end if end do if ( i /= 0 ) then exit end if end do if ( i == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRAPH_DIST_MIN_SPAN_TREE - Fatal error!' write ( *, '(a)' ) ' There are no nonzero arcs in this graph.' stop end if work(i) = 0.0E+00 iwork1(i) = 1 tree_length = 0.0E+00 kk4 = nnode - 1 do jj = 1, kk4 work(1:nnode) = huge ( d ) do i = 1, nnode ! ! For each forward arc originating at node I calculate ! the length of the path to node I ! if ( iwork1(i) == 1 ) then do j = 1, nnode if ( dist(i,j) < huge ( 1.0E+00 ) .and. iwork1(j) == 0 ) then d = tree_length + dist(i,j) if ( d < work(j) ) then work(j) = d iwork2(j) = i end if end if end do end if end do ! ! Find the minimum potential. ! d = huge ( d ) ient = 0 do i = 1, nnode if ( iwork1(i) == 0 .and. work(i) < d ) then d = work(i) ient = i itr = iwork2(i) end if end do ! ! Include the node in the current path. ! if ( d < huge ( d ) ) then iwork1(ient) = 1 tree_length = tree_length + dist(itr,ient) inode(jj) = itr jnode(jj) = ient end if end do return end subroutine graph_dist_print ( dist, lda, nnode, title ) ! !******************************************************************************* ! !! GRAPH_DIST_PRINT prints a distance matrix. ! ! ! Modified: ! ! 04 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real DIST(LDA,NNODE), the distance matrix. ! DIST(I,J) is the distance from node I to node J. ! ! Input, integer LDA, the leading dimension of DIST, which must be at ! least NNODE. ! ! Input, integer NNODE, the number of nodes. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none ! integer lda integer nnode ! real dist(lda,nnode) integer ihi integer ilo integer jhi integer jlo integer ncol integer nrow character ( len = * ) title ! if ( len_trim ( title ) /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if write ( *, '(a)' ) ' ' ilo = 1 ihi = nnode jlo = 1 jhi = nnode ncol = nnode nrow = nnode call rmat_print ( dist, ihi, ilo, jhi, jlo, lda, ncol, nrow ) return end subroutine i_swap ( i, j ) ! !******************************************************************************* ! !! I_SWAP switches two integer values. ! ! ! Modified: ! ! 30 November 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer I, J. On output, the values of I and ! J have been interchanged. ! implicit none ! integer i integer j integer k ! k = i i = j j = k return end subroutine int_lp ( m, n, a, b, c, lda, x, infs ) ! !******************************************************************************* ! !! INT_LP is a heuristic algorithm for the integer linear programming problem. ! ! ! Discussion: ! ! The integer linear programming problem has the form: ! ! Maximize the objective function: ! ! Sum ( 1 <= J <= N ) C(J) * X(J) ! ! Subject to the constraint equations: ! ! Sum ( 1 <= J <= N ) A(I,J) * X(J) <= B(I), for 1 <= I <= M, ! ! and positivity and integer conditions: ! ! X(J) >= 0, for 1 <= J <= N, ! X(J) integer, for 1 <= J <= N. ! ! The heuristic algorithm described here will usually find a vector ! X which is feasible, and gives a "reasonably good" objective ! function value. ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! ! Input, integer M, the number of constraints. ! ! Input, integer N, the number of variables. ! ! Input, real A(LDA,N+M+1), contains the coefficients of the constraints ! in the first M rows and N columns. There is also an M+1 row used ! for working storage. The last M+1 columns are also used as working ! storage. ! ! Input, real B(M), the right hand side of the constraints. ! ! Input, real C(N), the coefficients of the objective function. ! ! Input, integer LDA, the leading dimension of the A array. ! LDA must be at least M+1. ! ! Output, real X(max ( M+1, N ) ), contains the solution. ! ! Output, integer INFS, error indicator. ! 0, no error, a solution was found. ! 1, the linear program is infeasible. No solution was found. ! implicit none ! integer lda integer m integer n ! real a(lda,n+m+1) real b(m) real c(n) real eps integer i integer ii integer infs integer init integer isub integer iwk18(n) integer iwk19(m+1) integer iwk20(n+m) integer j integer jj integer jval integer k integer kval integer ll1 integer ll2 integer ll3 integer more integer moveon integer n1 real row_norm real sum2 real temp real tol real wk1(n,n) real wk2(m+1,m+1) real wk3(n) real wk4(n) real wk5(n) real wk6(n) real wk7(n) real wk8(n) real wk9(m) real wk10(m) real wk11(m) real wk12(m) real wk13(m+1) real wk14(m+1) real wk15(m+1) real wk16(n) logical wk17(n) real x(n) ! eps = sqrt ( epsilon ( eps ) ) ! ! Normalize the coefficients. ! do i = 1, m row_norm = sqrt ( sum ( a(i,1:n)**2 ) ) if ( row_norm /= 0.0E+00 ) then a(i,1:n) = a(i,1:n) / row_norm b(i) = b(i) / row_norm end if end do ! ! Set up the data for the simplex subroutine. ! do i = m, 1, -1 a(i+1,1:n) = a(i,1:n) end do a(1,1:n) = - c(1:n) a(1,n+1:n+m) = 0.0E+00 do i = 2, m + 1 a(i,n+1:n+m) = 0.0E+00 end do do i = 2, m + 1 j = n + i - 1 a(i,j) = 1.0E+00 end do ! ! Set the right hand side. ! wk16(1) = 0.0E+00 wk16(2:m+1) = b(1:m) ! ! Use the simplex method to find an optimal noninteger solution. ! call simplex ( n, m, a, x, infs, wk2, wk13, wk14, wk16, iwk19, iwk20 ) if ( infs /= 0 ) then return end if ! ! Restore the matrix. ! do i = 1, m a(i,1:n) = a(i+1,1:n) end do wk15(1:m+1) = x(1:m+1) do j = 1, n isub = iwk20(j) if ( isub > 0 ) then wk7(j) = x(isub) else wk7(j) = 0.0E+00 end if end do jj = 1 ii = 0 do j = 1, m + 1 do i = 1, m + 1 ii = ii + 1 if ( ii > m + 1 ) then ii = 1 jj = jj + 1 end if a(i,n+j) = wk2(ii,jj) end do end do do i = 1, m do j = 1, m wk2(i,j) = a(i+1,n+j+1) end do end do ! ! Select a second solution corresponding to the noninteger ! optimal solution by tightening the binding constraints ! sufficiently so that a rounded solution still satisfies ! the original constraints. ! do i = 1, m isub = iwk20(n+i) if ( isub > 0 ) then if ( x(isub) > 0.0E+00 ) then wk12(i) = b(i) go to 220 end if end if temp = 0.0E+00 do j = 1, n if ( iwk20(j) > 0 ) then temp = temp + abs(a(i,j)) end if end do wk12(i) = b(i) - 0.5E+00 * temp 220 continue end do ! ! The second solution is stored in WK8. ! wk8(1:n) = 0.0E+00 do i = 1, m isub = iwk19(i+1) if ( isub <= n ) then temp = 0.0E+00 do k = 1, m temp = temp + wk2(i,k) * wk12(k) end do wk8(isub) = temp end if end do ! ! Perform the linear search along the line joined by two points ! WK7 and WK8 found above. ! call parse1 ( n, m, lda, a, b, x, wk4, wk6, wk7, wk8, wk9, wk16, wk17 ) ! ! Search for better feasible solutions. ! call parse2 ( n, m, lda, a, b, c, x, wk3, wk5, wk10, iwk18, init ) ! ! Prepare for the interchange of two variables. ! if ( init > 1 ) then call parse4 ( n, c, wk1, iwk18, init ) end if ! ! Increase or decrease one variable by one as long as the ! resulting better solution is feasible. ! 260 continue call parse3 ( n, m, lda, a, c, x, wk10, kval ) if ( init <= 1 ) then return end if wk6(1:n) = x(1:n) jval = 1 kval = n ! ! Check whether the changed value of a variable is negative. ! 280 continue call parse5 ( n, m, lda, a, x, wk5, wk10, wk11, iwk18, moveon, jval ) if ( moveon /= 1 ) then kval = jval + 1 go to 320 end if 290 continue isub = iwk18(kval) ll1 = 0 ll2 = 0 do i = 1, m if ( wk11(i) < 0.0E+00 ) then if ( a(i,isub) == 0.0E+00 ) then go to 310 else if ( a(i,isub) < 0.0E+00 ) then ll1 = 1 else ll2 = 2 end if end if end do ll3 = ll1 + ll2 ! ! Investigate the possibility of changing X(JVAL) in one favorable ! direction and seek for a change in another variable X(KVAL). ! ! Try to decrease X(KVAL) to get a better solution. ! if ( ll3 == 2 ) then call parse6 ( n, m, lda, a, c, x, wk1, wk3, wk5, wk10, wk11, & iwk18, more, jval, kval ) if ( more == 1 ) then go to 280 end if ! ! Try to increase X(KVAL) to get a better solution. ! else if ( ll3 < 2 ) then call parse7 ( n, m, lda, a, c, x, wk1, wk3, wk5, wk10, wk11, & iwk18, more, jval, kval ) if ( more == 1 ) then go to 280 end if end if 310 continue if ( jval /= kval-1 ) then kval = kval - 1 go to 290 end if ! ! Consider changing X(JVAL) by one in the other direction and ! make a small integer change of X(KVAL). ! 320 continue call parse8 ( n, m, lda, a, x, wk1, wk5, wk10, wk11, iwk18, init, & jval, kval ) n1 = min ( n - 1, init ) ! ! If one or more improved solution is found then repeat the search. ! if ( jval /= n1 ) then jval = jval + 1 kval = n go to 280 end if do i = 1, n if ( x(i) /= wk6(i) ) then go to 260 end if end do return end subroutine k_center ( nnode, cost, kmax, lda, knum, kset ) ! !******************************************************************************* ! !! K_CENTER is a heuristic algorithm for the K-center problem. ! ! ! Discussion: ! ! Let V be the set of nodes of a complete undirected graph of ! NNODE nodes, with the edge from node I to J having nonnegative ! cost COST(I,J), and COST(I,I) = 0. ! ! Given an integer K between 1 and N, the K-center location problem ! is to find a subset S of the nodes V, of size at most K, such ! that ! ! Z(S) = max ( i in V ) min ( j in S ) COST(I,J) ! ! is minimized. ! ! ! In other words, for any subset S of the nodes, we measure the ! distance from every node in V to a node in S; the maximum of ! these values is the "cost" of the subset S. We are seeking the ! set S with minimum cost. This is a sort of discrete minimization ! problem in the L-infinity norm. ! ! ! In general, this is a hard problem. However, in the particular case ! where the COST function satisfies the triangle inequality, that is, ! ! COST(I,J) + COST(J,K) >= COST(I,K) for all I, J, K, ! ! then solutions close to the optimum can be found with this algorithm. ! In fact, the computed solution is guaranteed to have a value of Z ! no more than twice that of the optimal solution. ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! ! Input, integer NNODE, the number of nodes. ! ! Input, real COST(LDA,NNODE), the cost matrix, which contains the ! "cost" for each pair of nodes, or the distance between them. ! ! Input, integer KMAX, the maximum size of the subset of nodes to be found. ! ! Input, integer LDA, the leading dimension of COST. ! ! Output, integer KNUM, the size of the subset of nodes found. ! KNUM will be no greater than KMAX. ! ! Output, integer KSET(KMAX), contains in entries 1 through KNUM ! the list of nodes forming the subset. ! implicit none ! integer lda integer m integer nnode ! real big real cost(lda,nnode) real great integer high integer i integer ifirst integer inode integer iwk1(nnode,nnode) integer iwk2(nnode) integer iwk3(nnode) integer iwk4((nnode*(nnode-1))/2) integer iwk5((nnode*(nnode-1))/2) integer j integer k integer kmax integer knum integer kset(nnode) integer l integer low integer mid integer ncheck integer num1 integer num2 integer numk real small real work6((nnode*(nnode-1))/2) ! m = ( nnode * ( nnode - 1 ) ) / 2 great = 1.0E+00 do i = 1, nnode - 1 do j = i + 1, nnode great = great + cost(i,j) end do end do small = great do i = 1, nnode big = - great do j = 1, nnode if ( j /= i ) then if ( cost(i,j) > big ) then big = cost(i,j) numk = i end if end if end do if ( big < small ) then small = big knum = numk end if end do ifirst = knum ! ! Return the optimal solution if KMAX = 1. ! if ( kmax == 1 ) then knum = 1 kset(1) = ifirst return end if ! ! Sort the edges in order of increasing cost. ! k = 0 do i = 1, nnode - 1 do j = i + 1, nnode k = k + 1 work6(k) = cost(i,j) end do end do ! ! Sort the edges in order of increasing cost. ! call rvec_sort_heap_index_a ( m, work6, iwk4 ) do i = 1, m j = iwk4(i) iwk5(j) = i end do k = 0 do i = 1, nnode - 1 do j = i + 1, nnode k = k + 1 iwk1(i,j) = iwk5(k) iwk1(j,i) = iwk5(k) end do end do ! ! Binary search. ! low = 1 high = m 100 continue if ( high /= low + 1 ) then mid = ( high + low ) / 2 ! ! Restrict to the subgraph with the original N nodes ! but only having the first MID number of edges. ! numk = 0 ncheck = nnode iwk2(1:nnode) = 0 inode = ifirst 120 continue numk = numk + 1 ! ! Include node INODE into the solution set. ! iwk3(numk) = inode ! ! Consider all nodes adjacent to node INODE. ! do k = 1, nnode if ( k /= inode ) then if ( iwk1(k,inode) <= mid ) then ! ! Node K is adjacent to node INODE; ! delete node K from the subgraph. ! if ( iwk2(k) == 0 ) then ncheck = ncheck - 1 iwk2(k) = 2 end if ! ! Delete all nodes adjacent to node K from the subgraph. ! do l = 1, nnode if ( iwk2(l) == 0 ) then if ( l /= k ) then if ( iwk1(k,l) <= mid ) then ! ! Node L is adjacent to node K; ! delete node L from the subgraph. ! ncheck = ncheck - 1 iwk2(l) = 2 end if end if end if end do end if end if end do ! ! Mark node INODE as being already selected. ! if ( iwk2(inode) == 0 ) then ncheck = ncheck - 1 end if iwk2(inode) = 1 ! ! Continue the binary search if the subgraph is nonempty. ! if ( ncheck > 0 ) then ! ! Pick the next center by the greedy heuristic. ! if ( ncheck <= 2 ) then do i = 1, nnode if ( iwk2(i) == 0 ) then inode = i go to 120 end if end do end if small = great do i = 1, nnode if ( iwk2(i) == 0 ) then big = - great do j = 1, nnode if ( iwk2(j) == 0 ) then if ( j /= i ) then if ( cost(i,j) > big ) then big = cost(i,j) num1 = i end if end if end if end do if ( big < small ) then small = big num2 = num1 end if end if end do inode = num2 go to 120 end if ! ! Store up the temporary solution set. ! if ( numk <= kmax ) then high = mid knum = numk kset(1:knum) = iwk3(1:knum) else low = mid end if go to 100 end if return end subroutine k_median ( m, n, c, k, lda, isol ) ! !******************************************************************************* ! !! K_MEDIAN is a heuristic algorithm for the K-median problem. ! ! ! Discussion: ! ! If C is an M by N matrix, and K is an integer between 1 and M, ! the K median location problem is to select K rows of C which ! form a subset S that maximizes the value: ! ! F(S) = sum ( 1 <= J <= N ) max ( I in S ) C(I,J) ! ! The heuristic algorithm given here finds a near-optimum solution in ! time proportional to the number of entries of C. ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns of the cost matrix. ! ! Input, real C(LDA,N), the M by N cost matrix. ! ! Input, integer K, the number of rows to be used in the maximizing set. ! ! Input, integer LDA, the leading dimension of the cost matrix. ! ! Output, integer ISOL(M), contains in entries 1 through K the elements ! of the maximizing set. ! implicit none ! integer lda integer m integer n ! real big real c(lda,n) real eps integer i integer iflag integer ii integer ij integer index integer ir integer isol(m) integer iw integer iwk5(n) integer iwk6(n) integer j integer k integer l real optval real row_sum real temp real tmin real tol real tot1 real w real work1(n) real work2(n) real work3(m) real work4(m) real z ! eps = sqrt ( epsilon ( eps ) ) big = huge ( big ) work3(1:m) = 0.0E+00 ! ! Obtain the optimal solution when K = 1. ! if ( k == 1 ) then optval = - big do i = 1, m row_sum = sum ( c(i,1:n) ) if ( row_sum > optval ) then optval = row_sum index = i end if end do isol(1) = index return end if ! ! Find the largest and second-largest elements among the first ! M elements of column J; store them in WORK1(J) and WORK2(J). ! do i = 1, n work1(i) = - big do j = 1, k if ( c(j,i) > work1(i) ) then work1(i) = c(j,i) iwk5(i) = j end if end do work2(i) = - big do j = 1, k if ( c(j,i) > work2(i) .and. j /= iwk5(i) ) then work2(i) = c(j,i) iwk6(i) = j end if end do end do do i = 1, k tot1 = 0.0E+00 do j = 1, n if ( iwk5(j) == i ) then tot1 = tot1 + work1(j) - work2(j) end if end do work3(i) = tot1 end do do i = 1, m isol(i) = i end do iflag = 0 ir = k do if ( ir + 1 > m ) then ir = k end if ir = ir + 1 do i = 1, k work4(isol(i)) = 0.0E+00 end do work4(isol(ir)) = 0.0E+00 ! ! Find out which of the K+l rows will make the objective ! value decrease the least if the row is removed. ! do j = 1, n z = c(isol(ir),j) if ( z > work2(j) .and. z <= work1(j) ) then work4(iwk5(j)) = work4(iwk5(j)) - z + work2(j) else if ( z > work1(j) ) then work4(isol(ir)) = work4(isol(ir)) + z - work1(j) work4(iwk5(j)) = work4(iwk5(j)) - work1(j) + work2(j) end if end if end do ! ! Determine which row to remove. ! tmin = big l = 0 do ij = 1, k i = isol(ij) temp = work3(i) + work4(i) if ( temp < tmin ) then tmin = temp l = ij end if end do i = isol(ir) temp = work3(i) + work4(i) if ( temp <= tmin .or. abs ( temp - tmin ) <= eps ) then iflag = iflag + 1 if ( iflag == m - k ) then return end if cycle end if iflag = 0 tot1 = 0.0E+00 do i = 1, n if ( c(isol(ir),i) > work1(i) ) then tot1 = tot1 + c(isol(ir),i) - work1(i) end if end do work3(isol(ir)) = tot1 ! ! Update all arrays after exchanging two rows. ! do j = 1, n z = c(isol(ir),j) if ( z <= work2(j) ) then if ( isol(l) == iwk5(j) ) then ! ! Replacing the largest element by something no better than ! the third largest. ! work1(j) = work2(j) iwk5(j) = iwk6(j) w = - big do ij = 1, k ii = isol(ij) if ( c(ii,j) > w .and. ii /= iwk5(j) ) then w = c(ii,j) iw = ii end if end do ii = isol(ir) if ( c(ii,j) > w .and. ii /= iwk5(j) ) then w = c(ii,j) iw = ii end if work3(iwk6(j)) = work3(iwk6(j)) - w + work2(j) work2(j) = w iwk6(j) = iw end if if ( isol(l) == iwk6(j) ) then ! ! Replacing the second largest element by something no ! better than the third largest. ! w = - big do ij = 1, k ii = isol(ij) if ( c(ii,j) > w .and. ii /= iwk5(j) .and. ii /= iwk6(j) ) then w = c(ii,j) iw = ii end if end do ii = isol(ir) if ( c(ii,j) > w .and. ii /= iwk5(j) .and. ii /= iwk6(j) ) then w = c(ii,j) iw = ii end if work3(iwk5(j)) = work3(iwk5(j)) - w + work2(j) work2(j) = w iwk6(j) = iw end if else if ( z > work2(j) .and. z <= work1(j) ) then if ( isol(l) == iwk5(j) ) then ! ! Replacing the largest element by a new and smaller largest element. ! work1(j) = z iwk5(j) = isol(ir) work4(isol(ir)) = work4(isol(ir)) + z - work2(j) else ! ! Z becomes the new second largest element. ! work3(iwk5(j)) = work3(iwk5(j)) - z + work2(j) work2(j) = z iwk6(j) = isol(ir) end if else if ( z > work1(j) ) then if ( isol(l) == iwk5(j) ) then ! ! Replacing the largest element by a new largest element. ! work4(isol(ir)) = work4(isol(ir)) + work1(j) - work2(j) iwk5(j) = isol(ir) work1(j) = z else ! ! Z becomes the largest element. ! work3(iwk5(j)) = work3(iwk5(j)) - work1(j) + work2(j) work2(j) = work1(j) work1(j) = z iwk6(j) = iwk5(j) iwk5(j) = isol(ir) end if end if end if end if end do ! ! Iterate until no improvement by exchange can be found. ! call i_swap ( isol(l), isol(ir) ) work3(isol(l)) = work4(isol(l)) do i = k + 1, m work3(isol(i)) = 0.0E+00 end do end do return end subroutine knapsack ( n, a, b, c, eps, m, objval, numsol, isol ) ! !******************************************************************************* ! !! KNAPSACK is a heuristic algorithm for the zero-one knapsack problem. ! ! ! Discussion: ! ! The zero-one knapsack problem seeks to maximize ! ! Sum ( 1 <= J <= N ) C(J) * X(J) ! ! subject to ! ! Sum ( 1 <= J <= N ) A(J) * X(J) <= B ! ! where ! ! A(J) and C(J) are nonnegative for each J; ! B is nonnegative; ! X(J) = 0 or 1 for each J. ! ! Intuitively, C(J) is the value of X(J), A(J) is the "size" or "cost" ! of X(J), and B is an overall size or spending limit. Thus, ! the ratio C(J)/A(J) is a heuristic indication of the relative ! value of X(J). ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! ! Input, integer N, the number of variables. ! ! Input, real A(N), the coefficients of the constraints. ! ! Input, real B, the right hand side of the constraints. ! ! Input, real C(N), the coefficients of the objective function. ! ! Input, real EPS, a positive real number prescribing the required ! degree of accuracy in the solution. ! ! Input, integer M, the smallest integer greater than ( 3 / EPS )**2. ! ! Output, real OBJVAL, the value of the objective function for the ! suggested solution. ! ! Output, integer NUMSOL, the number of nonzero variables in the ! suggested solution. ! ! Output, integer ISOL(N), contains the indices of the nonzero variables ! in entries 1 through NUMSOL. ! implicit none ! integer m integer n ! real a(n) real b real big real c(n) real eps real est real gvalue integer i integer icol1 integer icol2 integer igp1 integer igp2 integer ii integer index integer irow integer isol(n) integer isub integer iwk1(n) integer iwk2(n) integer iwk3(n) integer iwk4(n) logical iwk5(m) integer iwk6(m) integer iwk7(n,m) integer j integer k integer numsol real objval real parm1 real parm2 real rhs logical switch real value real work1(2,m) real work2(2,m) real work3(n) real work4(n) real work5(n) real xtemp real ytemp ! ! Reorder the variables. ! work3(1:n) = c(1:n) / a(1:n) big = 1.0E+00 + sum ( c ) call rvec_sort_heap_index_d ( n, work3, iwk1 ) ! ! Calculate the initial parameters. ! numsol = 0 objval = 0.0E+00 rhs = b j = 1 20 continue index = iwk1(j) if ( rhs >= a(index) ) then rhs = rhs - a(index) j = j + 1 if ( j <= n ) then go to 20 end if end if if ( j > n ) then numsol = n do i = 1, n objval = objval + c(j) end do do i = 1, n isol(i) = i end do return end if est = 0.0E+00 do i = 1, j est = est + c(iwk1(i)) end do xtemp = eps / 3.0E+00 parm2 = xtemp * est parm1 = xtemp * parm2 ! ! Split the variables into two groups. ! igp1 = 0 igp2 = 0 do i = 1, n index = iwk1(i) if ( c(index) >= parm2 ) then igp1= igp1 + 1 iwk4(igp1) = c(index) / parm1 work3(igp1) = a(index) iwk2(igp1) = index else igp2 = igp2 + 1 work5(igp2) = a(index) work4(igp2) = c(index) iwk3(igp2) = index end if end do if ( igp1 == 0 .or. m <= 1 ) then rhs = b j = 1 do while ( j <= n ) index = iwk1(j) if ( rhs >= a(index) ) then numsol = numsol + 1 isol(numsol) = index objval = objval + c(index) rhs = rhs - a(index) end if j = j + 1 end do return end if ! ! Solve the sequence of problems in group one. ! work1(1,1) = 0.0E+00 work1(2,1) = 0.0E+00 work2(1,1) = 0.0E+00 work2(2,1) = 0.0E+00 work1(1,2:m) = big work2(1,2:m) = 0.0E+00 i = iwk4(1) + 1 work1(1,i) = work3(1) index = iwk2(1) work2(1,i) = c(index) switch = .true. icol2 = 1 k = 2 do while ( k <= igp1 ) if ( switch ) then icol1 = 1 icol2 = 2 else icol1 = 2 icol2 = 1 end if switch = .not. switch do i = 2, m irow = i - 1 - iwk4(k) if ( irow < 0 ) then xtemp = big else irow = irow + 1 xtemp = work1(icol1,irow) if ( xtemp < big ) then xtemp = xtemp + work3(k) end if end if ytemp = work1(icol1,i) if ( ytemp <= xtemp ) then work1(icol2,i) = ytemp work2(icol2,i) = work2(icol1,i) else work1(icol2,i) = xtemp index = iwk2(k) work2(icol2,i) = work2(icol1,irow) + c(index) end if end do k = k + 1 end do ! ! Find the maximum objective value. ! do i = 1, m value = work1(icol2,i) if ( value <= b ) then gvalue = 0.0E+00 if ( igp2 > 0 ) then rhs = b - value if ( rhs < 0 ) then work2(icol2,i) = 0.0E+00 else j = 1 do while ( j <= igp2 ) if ( rhs >= work5(j) ) then gvalue = gvalue + work4(j) rhs = rhs - work5(j) end if j = j + 1 end do end if end if gvalue = gvalue + work2(icol2,i) if ( gvalue > objval ) then objval = gvalue isub = i end if end if end do ! ! Obtain the best solution. ! if ( isub == 1 ) then if ( igp2 > 0 ) then rhs = b j = 1 do while ( j <= igp2 ) index = iwk3(j) if ( rhs >= work5(j) ) then numsol = numsol + i isol(numsol) = index rhs = rhs - work5(j) end if j = j + 1 end do end if return end if iwk5(1:m) = .false. iwk6(1:m) = 0 work1(1,1) = 0.0E+00 work1(2,1) = 0.0E+00 work2(1,1) = 0.0E+00 work1(1,2:isub) = big i = iwk4(1) + 1 work1(1,i) = work3(1) iwk6(i) = 1 iwk7(1,i) = 1 switch = .true. icol2 = 1 k = 2 do while ( k <= igp1 ) if ( switch ) then icol1 = 1 icol2 = 2 else icol1 = 2 icol2 = 1 end if switch = .not. switch do i = 2, isub irow = i - 1 - iwk4(k) if ( irow < 0 ) then xtemp = big else irow = irow + 1 xtemp = work1(icol1,irow) if ( xtemp < big ) then xtemp = xtemp + work3(k) end if end if ytemp = work1(icol1,i) if ( ytemp <= xtemp ) then work1(icol2,i) = ytemp else work1(icol2,i) = xtemp end if if ( work1(icol1,i) /= work1(icol2,i) ) then iwk5(i) = .true. end if end do do ii = 2, isub i = isub - ii + 2 if ( iwk5(i) ) then iwk5(i) = .false. irow = i - iwk4(k) j = iwk6(irow) do while ( j > 0 ) iwk7(j,i) = iwk7(j,irow) j = j - 1 end do index = iwk6(irow) + 1 iwk6(i) = index iwk7(index,i) = k end if end do k = k + 1 end do j = iwk6(isub) do while ( j > 0 ) numsol = numsol + 1 index = iwk7(j,isub) isol(numsol) = iwk2(index) j = j - 1 end do if ( igp2 > 0 ) then rhs = b - work1(icol2,isub) j = 1 do while ( j <= igp2 ) if ( rhs >= work5(j)) then numsol = numsol + 1 isol(numsol) = iwk3(j) rhs = rhs - work5(j) end if j = j + 1 end do end if return end subroutine multi_knap ( m, n, a, b, c, lda, objval, numsol, isol ) ! !******************************************************************************* ! !! MULTI_KNAP is a heuristic algorithm for the multi-dimensional 0/1 knapsack problem. ! ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! ! Input, integer M, the number of constraints. ! ! Input, integer N, the number of variables. ! ! Input, real A(LDA,N), the M by N matrix of coefficients of the constraints. ! ! Input, real B(M), the right hand sides of the constraints. ! ! Input, real C(N), the coefficients of the objective function. ! ! Input, integer LDA, the leading dimension of the matrix A. ! ! Output, real OBJVAL, the value of the objective function for the ! suggested solution. ! ! Output, integer NUMSOL, the number of nonzero variables in the ! suggested solution. ! ! Output, integer ISOL(N), contains the indices of the nonzero variables ! in entries 1 through NUMSOL. ! implicit none ! integer lda integer m integer n ! real a(lda,n) real b(m) real big real c(n) logical check real cmax real col_sum real eps real grad real grdmax integer i integer isol(n) integer iwk3(n) integer j integer jgdmax integer jgm integer jj integer ncan integer numcan integer numsol real objval real orimov real rootm real small real sum2 real sumsq real tol real wk1(m) real wk2(m) ! eps = sqrt ( epsilon ( eps ) ) ! ! Compute the machine infinity. ! big = huge ( big ) small = tiny ( small ) ! ! Initialize. ! do i = 1, m a(i,1:n) = a(i,1:n) / b(i) end do wk1(1:m) = 0.0E+00 objval = 0.0E+00 numsol = 0 do j = 1, n iwk3(j) = j end do rootm = sqrt ( real ( m ) ) numcan = n 60 continue ncan = numcan numcan = 0 grdmax = - big ! ! When the resource requirement vector is zero. ! do j = 1, ncan ! ! Select variables. ! jj = iwk3(j) do i = 1, m if ( wk1(i) + a(i,jj) - 1.0E+00 > eps ) then go to 90 end if end do numcan = numcan + 1 iwk3(numcan) = jj ! ! Compute the effective gradients. ! col_sum = sum ( a(1:m,jj) ) if ( col_sum <= small ) then grad = big else grad = c(jj) * rootm / col_sum end if ! ! Find the variable whose effective gradient is the largest. ! if ( grad > grdmax ) then grdmax = grad jgdmax = jj jgm = numcan end if 90 continue end do ! ! Accept the variable whose effective gradient is the largest. ! if ( numcan <= 0 ) then return end if if ( numcan == 1 ) then objval = objval + c(jgdmax) wk1(1:m) = wk1(1:m) + a(1:m,jgdmax) numsol = numsol + 1 isol(numsol) = jgdmax return end if objval = objval + c(jgdmax) numsol = numsol + 1 isol(numsol) = jgdmax iwk3(jgm) = iwk3(numcan) numcan = numcan - 1 check = .true. do i = 1, m wk1(i) = wk1(i) + a(i,jgdmax) if ( wk1(i) > eps ) then check = .false. end if end do if ( check ) then go to 60 end if do cmax = wk1(1) do i = 2, m cmax = max ( cmax, wk1(i) ) end do orimov = cmax * cmax sumsq = 0.0E+00 do i = 1, m wk2(i) = wk1(i) - orimov if ( wk2(i) <= eps ) then wk2(i) = 0.0E+00 end if sumsq = sumsq + wk2(i) * wk2(i) end do sumsq = sqrt ( sumsq ) ncan = numcan numcan = 0 grdmax = - big ! ! When the resource requirement vector is non-zero. ! do j = 1, ncan ! ! Select variables. ! jj = iwk3(j) do i = 1, m if ( wk1(i) + a(i,jj) - 1.0E+00 > eps ) then go to 170 end if end do numcan = numcan + 1 iwk3(numcan) = jj ! ! Compute the effective gradients. ! sum2 = dot_product ( wk2(1:m), a(1:m,jj) ) if ( sum2 <= small ) then grad = big else grad = c(jj) * sumsq / sum2 end if ! ! Find the variable whose effective gradient is the largest. ! if ( grad > grdmax ) then grdmax = grad jgdmax = jj jgm = numcan end if 170 continue end do ! ! Accept the variable whose effective gradient is the largest. ! if ( numcan <= 0 ) then exit end if objval = objval + c(jgdmax) numsol = numsol + 1 isol(numsol) = jgdmax wk1(1:m) = wk1(1:m) + a(1:m,jgdmax) if ( numcan == 1 ) then exit end if iwk3(jgm) = iwk3(numcan) numcan = numcan - 1 end do return end subroutine parse1 ( n, m, lda, a, b, x, wk4, wk6, wk7, wk8, wk9, wk16, wk17 ) ! !******************************************************************************* ! !! PARSE1 is used by INT_LP. ! ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! implicit none ! integer lda integer m integer n ! real a(lda,n+m+1) real b(m) real delps real delta real eps integer i integer incr integer j integer jhigh integer k integer kk integer ll real temp logical wk17(n) real wk4(n) real wk6(n) real wk7(n) real wk8(n) real wk9(m) real wk16(n) real x(n) real xwork real xx1 real xx2 real yy1 real yy2 ! ! Linear search. ! eps = - huge ( eps ) delps = 0.05E+00 delta = 0.0E+00 do temp = 1.0E+00 - delta ! ! Form a linear combination of WK7 and WK8. ! wk6(1:n) = wk7(1:n) + delta * ( wk8(1:n) - wk7(1:n) ) do j = 1, n if ( wk6(j) <= 0.0E+00 ) then x(j) = 0.0E+00 else x(j) = aint ( wk6(j) + 0.5E+00 ) end if end do ! ! Compute WK9, the degree of infeasibility. ! do i = 1, m wk9(i) = dot_product ( a(i,1:n), x(1:n) ) - b(i) end do xwork = 0.0E+00 do i = 1, m if ( wk9(i) > 0.0E+00 ) then xwork = xwork + wk9(i) end if end do 70 continue if ( xwork <= 0.0E+00 ) then return end if ll = m ! ! Compute the change in solution X that improves the objective value. ! This is indicated by WK4. ! incr = 0 xx1 = eps kk = 0 do j = 1, n yy1 = 0.0E+00 yy2 = 0.0E+00 do i = 1, ll if ( wk9(i) > 0.0E+00 ) then yy1 = yy1 + a(i,j) end if end do xx2 = abs ( yy1 ) if ( yy1 == 0.0E+00 ) then wk4(j) = 0.0E+00 wk17(j) = .false. else if ( yy1 < 0.0E+00 ) then wk4(j) = 1.0E+00 do i = 1, ll yy1 = wk9(i) + a(i,j) if ( yy1 > 0.0E+00 ) then yy2 = yy2 + yy1 end if end do else if ( yy1 > 0.0E+00 ) then if ( x(j) > 0.0E+00 ) then wk4(j) = - 1.0E+00 do i = 1, ll yy1 = wk9(i) - a(i,j) if ( yy1 > 0.0E+00 ) then yy2 = yy2 + yy1 end if end do end if end if wk16(j) = yy2 if ( wk16(j) <= 0.0E+00 ) then if ( xx2 > xx1 ) then xx1 = xx2 kk = j end if else if ( yy2 < xwork ) then wk17(j) = .true. k = j incr = incr + 1 else wk17(j) = .false. end if end if end if end do if ( kk /= 0 ) then x(kk) = x(kk) + wk4(kk) return end if ! ! Compute the improvement of the objective function. ! if ( incr > 0 ) then if ( incr > 1 ) then temp = eps do j = 1, n if ( wk17(j) ) then if ( xwork - wk16(j) > temp ) then jhigh = j temp = xwork - wk16(j) end if end if end do k = jhigh end if x(k) = x(k) + wk4(k) wk9(1:ll) = wk9(1:ll) + a(1:ll,k) * wk4(k) xwork = wk16(k) go to 70 end if ! ! Continue with the linear search. ! delta = delta + delps if ( delta > 1.0E+00 ) then exit end if end do return end subroutine parse2 ( n, m, lda, a, b, c, x, wk3, wk5, wk10, iwk18, init ) ! !******************************************************************************* ! !! PARSE2 is used by INT_LP. ! ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! implicit none ! integer lda integer m integer n ! real a(lda,n+m+1) real b(m) real c(n) integer i integer init integer iwk18(n) integer j integer ll real wk3(n) real wk5(n) real wk10(m) real temp real x(n) ! ! Improve the feasible solution by changing some variable by one. ! This is the initialization part. ! do j = 1, n iwk18(j) = j end do wk3(1:n) = c(1:n) do i = 1, n - 1 temp = abs ( wk3(i) ) do j = i + 1, n if ( abs ( wk3(j) ) > temp ) then temp = abs ( wk3(j) ) call i_swap ( iwk18(i), iwk18(j) ) call r_swap ( wk3(i), wk3(j) ) end if end do end do do j = 1, n if ( abs ( wk3(j) ) > 0.0E+00 ) then init = j end if end do ! ! Set WK5 to indicate whether the objective coefficient is positive. ! do j = 1, init ll = iwk18(j) if ( wk3(j) < 0.0E+00 ) then wk5(ll) = - 1.0E+00 else if ( wk3(j) > 0.0E+00 ) then wk5(ll) = 1.0E+00 end if end do ! ! Compute WK10, the feasibility test slacks. ! do i = 1, m wk10(i) = b(i) - dot_product ( a(i,1:n), x(1:n) ) end do return end subroutine parse3 ( n, m, lda, a, c, x, wk10, kval ) ! !******************************************************************************* ! !! PARSE3 is used by INT_LP. ! ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 21 July 2000 ! ! Parameters: ! implicit none ! integer lda integer m integer n ! real a(lda,n+m+1) real c(n) integer i integer j integer kval real temp real wk10(m) real x(n) real xx1 real xx2 ! ! Iteratively change the value of some variable by 1 to improve ! the objective value. ! do kval = 1 xx1 = 0.0E+00 ! ! The variable that gives a large increase in the objective value ! is chosen to be changed. ! do j = 1, n if ( c(j) /= 0.0E+00 ) then if ( c(j) > 0.0E+00 .or. x(j) > 0.0 ) then temp = huge ( 1.0E+00 ) do i = 1, m if ( c(j) * a(i,j) > 0.0E+00 ) then xx2 = wk10(i) / abs ( a(i,j) ) if ( xx2 < 0.0E+00 ) then if ( xx2 /= aint ( xx2 ) ) then xx2 = aint ( xx2 ) - 1.0E+00 end if else if ( xx2 > 0.0E+00 ) then xx2 = aint ( xx2 ) end if temp = min ( temp, xx2 ) end if end do if ( xx1 <= abs ( c(j) ) * temp ) then xx1 = abs ( c(j) ) * temp kval = j end if end if end if end do ! ! Increase or decrease X(KVAL) by one. ! if ( xx1 == 0.0E+00 ) then exit end if if ( c(kval) <= 0.0E+00 ) then x(kval) = x(kval) - 1.0E+00 wk10(1:m) = wk10(1:m) + a(1:m,kval) else x(kval) = x(kval) + 1.0E+00 wk10(1:m) = wk10(1:m) - a(1:m,kval) end if end do return end subroutine parse4 ( n, c, wk1, iwk18, init ) ! !******************************************************************************* ! !! PARSE4 is used by INT_LP. ! ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! implicit none ! integer n ! real c(n) integer init integer iwk18(n) integer j integer j2 integer j3 integer k real temp real wk1(n,n) ! ! Prepare for the interchange of two variables. ! ! Determine the bounds on the size of changes. ! ! The ratios of the objective coefficients are computed and stored in WK1. ! do j = 1, init - 1 j2 = iwk18(j) do k = j + 1, init j3 = iwk18(k) temp = abs ( c(j2) / c(j3) ) if ( temp <= 1.0E+00 ) then wk1(j,k) = aint ( temp - 1.0E+00 ) else if ( aint ( temp ) == temp ) then wk1(j,k) = aint ( temp - 1.0E+00 ) else wk1(j,k) = aint ( temp - 1.0E+00 ) + 1.0E+00 end if end if if ( temp >= -1.0E+00 ) then wk1(k,j) = aint ( temp + 1.0E+00 ) else if ( aint ( temp ) == temp ) then wk1(k,j) = aint ( temp + 1.0E+00 ) else wk1(k,j) = aint ( temp + 1.0E+00 ) - 1.0E+00 end if end if end do end do return end subroutine parse5 ( n, m, lda, a, x, wk5, wk10, wk11, iwk18, moveon, jval ) ! !******************************************************************************* ! !! PARSE5 is used by INT_LP. ! ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! implicit none ! integer lda integer m integer n ! real a(lda,n+m+1) integer i integer icon integer isub integer iwk18(n) integer jval integer moveon real wk5(n) real wk10(m) real wk11(m) real x(n) ! ! Check if the changed value of X(JVAL) is negative. ! isub = iwk18(jval) do if ( x(isub) < - wk5(isub) ) then moveon = 0 exit end if ! ! If the change is negative, then check whether this change is ! feasible without changing another variable. ! icon = 0 do i = 1, m wk11(i) = wk10(i) - wk5(isub) * a(i,isub) if ( wk11(i) < 0.0E+00 ) then icon = 1 end if end do if ( icon == 1 ) then moveon = 1 exit end if x(isub) = x(isub) + wk5(isub) wk10(1:m) = wk11(1:m) end do return end subroutine parse6 ( n, m, lda, a, c, x, wk1, wk3, wk5, wk10, wk11, & iwk18, more, jval, kval ) ! !******************************************************************************* ! !! PARSE6 is used by INT_LP. ! ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! implicit none ! integer lda integer m integer n ! real a(lda,n+m+1) real c(n) integer i integer isub integer iwk18(n) integer jval integer jsub integer kval integer more real temp real wk1(n,n) real wk3(n) real wk5(n) real wk10(m) real wk11(m) real x(n) real xwk1 real xwk2 real xwk3 ! ! Check for the change of variables X(JVAL) and X(KVAL) and identify a ! better solution. ! isub = iwk18(kval) if ( wk3(kval) > 0.0E+00 ) then xwk2 = - min ( x(isub), wk1(jval,kval) ) else xwk2 = - x(isub) end if xwk3 = huge ( 1.0E+00 ) isub = iwk18(kval) do i = 1, m if ( wk11(i) < 0.0E+00 ) then temp = wk11(i) / a(i,isub) if ( temp /= aint ( temp ) ) then temp = aint ( temp ) - 1.0E+00 end if xwk3 = min ( xwk3, temp ) end if end do ! ! Check whether feasibility can be restored. ! if ( xwk2 > xwk3 ) then more = 0 return end if isub = iwk18(kval) xwk1 = - huge ( 1.0E+00 ) do i = 1, m if ( a(i,isub) < 0.0E+00 ) then temp = wk11(i) / a(i,isub) if ( temp < 0.0E+00 ) then temp = aint ( temp ) else if ( temp > 0.0E+00 ) then if ( temp /= aint ( temp ) ) then temp = aint ( temp ) + 1.0E+00 end if end if end if xwk1 = max ( xwk1, temp ) end if end do ! ! Check if an improved solution can be obtained. ! xwk2 = max ( xwk2, xwk1 ) if ( xwk2 > xwk3 ) then more = 0 return end if isub = iwk18(kval) if ( c(isub) <= 0.0E+00 ) then jsub = iwk18(jval) x(jsub) = x(jsub) + wk5(jsub) x(isub) = x(isub) + xwk2 wk10(1:m) = wk11(1:m) - xwk2 * a(1:m,isub) else x(isub) = x(isub) + xwk3 wk10(1:m) = wk11(1:m) - xwk3 * a(1:m,isub) isub = iwk18(jval) x(isub) = x(isub) + wk5(isub) end if more = 1 return end subroutine parse7 ( n, m, lda, a, c, x, wk1, wk3, wk5, wk10, wk11, & iwk18, more, jval, kval ) ! !******************************************************************************* ! !! PARSE7 is used by INT_LP. ! ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! implicit none ! integer lda integer m integer n ! real a(lda,n+m+1) real c(n) integer i integer isub integer iwk18(n) integer jval integer jsub integer kval integer more real temp real wk1(n,n) real wk3(n) real wk5(n) real wk10(m) real wk11(m) real x(n) real xwk1 real xwk2 real xwork ! ! Try to increase X(KVAL) to get a better solution. ! if ( wk3(kval) < 0.0E+00 ) then xwk1 = wk1(jval,kval) else xwk1 = huge ( 1.0E+00 ) end if isub = iwk18(kval) temp = - huge ( 1.0E+00 ) do i = 1, m if ( wk11(i) < 0.0E+00 ) then xwork = wk11(i) / a(i,isub) if ( aint ( xwork ) /= xwork ) then xwork = aint ( xwork ) + 1.0E+00 end if temp = max ( temp, xwork ) end if end do ! ! Check whether feasibility can be restored. ! xwk2 = temp if ( xwk2 > xwk1 ) then more = 0 return end if isub = iwk18(kval) temp = huge ( 1.0E+00 ) do i = 1, m if ( a(i,isub) > 0.0E+00 ) then xwork = wk11(i) / a(i,isub) if ( xwork > 0.0E+00 ) then xwork = aint ( xwork ) else if ( xwork < 0.0E+00 ) then if ( xwork /= aint ( xwork ) ) then xwork = aint ( xwork ) - 1.0E+00 end if end if temp = min ( temp, xwork ) end if end do ! ! Check if an improved solution can be obtained. ! xwk1 = min ( xwk1, temp ) if ( xwk2 > xwk1 ) then more = 0 return end if isub = iwk18(kval) if ( c(isub) <= 0.0E+00 ) then jsub = iwk18(jval) x(jsub) = x(jsub) + wk5(jsub) x(isub) = x(isub) + xwk2 wk10(1:m) = wk11(1:m) - xwk2 * a(1:m,isub) ! ! A better solution is found. ! else x(isub) = x(isub) + xwk1 wk10(1:m) = wk11(1:m) - xwk1 * a(1:m,isub) isub = iwk18(jval) x(isub) = x(isub) + wk5(isub) end if more = 1 return end subroutine parse8 ( n, m, lda, a, x, wk1, wk5, wk10, wk11, iwk18, init, & jval, kval ) ! !******************************************************************************* ! !! PARSE8 is used by INT_LP. ! ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! implicit none ! integer lda integer m integer n ! real a(lda,n+m+1) integer i integer init integer isub integer iwk18(n) integer jval integer kval real wk1(n,n) real wk5(n) real wk10(m) real wk11(m) real x(n) ! ! Consider changing X(JVAL) and a small integer change of X(KVAL). ! do isub = iwk18(jval) if ( x(isub) < wk5(isub) ) then exit end if wk11(1:m) = wk10(1:m) + a(1:m,isub) * wk5(isub) 30 continue if ( kval > init ) then return end if ! ! Check whether the simultaneous change of X(JVAL) and X(KVAL) is feasible. ! isub = iwk18(kval) if ( x(isub) + wk5(isub) * wk1(kval,jval) < 0.0E+00 ) then kval = kval + 1 go to 30 end if isub = iwk18(kval) do i = 1, m if ( wk11(i) < wk5(isub) * wk1(kval,jval) * a(i,isub) ) then kval = kval + 1 go to 30 end if end do ! ! Make the simultaneous change of X(JVAL) and X(KVAL). ! isub = iwk18(jval) x(isub) = x(isub) - wk5(isub) isub = iwk18(kval) x(isub) = x(isub) + wk5(isub) * wk1(kval,jval) wk10(1:m) = wk11(1:m) - a(1:m,isub) * wk5(isub) * wk1(kval,jval) end do return end subroutine partition ( n, cost, lda, init, ip, iq, kp, kq, tcost ) ! !******************************************************************************* ! !! PARTITION is a heuristic algorithm for the graph partitioning problem. ! ! ! Discussion: ! ! Let V be the set of 2 * N nodes of a complete graph with an ! associated 2*N by 2*N symmetric cost matrix C on its edges. ! The graph partitioning problem is to partion the nodes into ! two sets P and Q, each containing N nodes, such that we minimize ! the total cost of the edges in the cut set: ! ! Sum ( I in P, J in Q ) C(I,J) ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! ! Input, integer N, the number of nodes in each partition. ! ! Input, real COST(LDA,2*N), the 2*N by 2*N edge cost matrix. ! ! Input, integer LDA, the leading dimension of COST, which must ! be at least 2*N. ! ! Input, integer INIT: ! If TRUE, then the initial partition will be generated by this routine; ! If FALSE, then the initial partition will be supplied on input. ! ! Input, integer IP(N), IQ(N). ! If INIT is TRUE, then IP and IQ are not input quantities, but simply ! workspace controlled by the routine. ! If INIT is FALSE, then the user must supply in IP and IQ and initial ! partition of the nodes. ! ! Output, integer KP(N), KQ(N), the final partition of the nodes. ! ! Output, real TCOST, the total cost of the edges in the cut set. ! implicit none ! integer lda integer n ! real cost(lda,2*n) real gain integer i integer ia integer ib integer ind1 integer ind2 logical init integer ip(n) integer iq(n) logical iwk4(n) logical iwk5(n) integer j integer k integer kp(n) integer kq(n) real small real tcost real tmax real tot1 real tot2 real wk1(n) real wk2(n) real wk3(n) ! ! Initial partitioning. ! if ( init ) then do i = 1, n ip(i) = i iq(i) = i + n end do end if ! ! Set the flags. ! do iwk4(1:n) = .true. iwk5(1:n) = .true. tcost = 0.0E+00 do i = 1, n do j = 1, n tcost = tcost + cost(ip(i),iq(j)) end do end do small = - 2.0E+00 * tcost ! ! Calculate the external cost of each element in the first partition. ! do i = 1, n tot1 = 0.0E+00 do j = 1, n tot1 = tot1 + cost(ip(i),iq(j)) end do ! ! Calculate the internal cost of each element in the first partition. ! tot2 = 0.0E+00 do k = 1, n tot2 = tot2 + cost(ip(i),ip(k)) end do ! ! Calculate the difference between external and internal costs. ! wk1(i) = tot1 - tot2 end do do i = 1, n ! ! Calculate the external cost of each element in the second partition. ! tot1 = 0.0E+00 do j = 1, n tot1 = tot1 + cost(iq(i),ip(j)) end do ! ! Calculate the internal cost of each element in the second partition. ! tot2 = 0.0E+00 do k = 1, n tot2 = tot2 + cost(iq(i),iq(k)) end do ! ! Calculate the difference between external and internal costs. ! wk2(i) = tot1 - tot2 end do do i = 1, n ! ! Choose IA from the first partition and IB from the second ! partition such that the gain is maximum. ! tmax = small do j = 1, n if ( iwk4(j) ) then do k = 1, n if ( iwk5(k) ) then gain = wk1(j) + wk2(k) - 2.0E+00 * cost(ip(j),iq(k)) if ( gain > tmax ) then tmax = gain ia = ip(j) ib = iq(k) ind1 = j ind2 = k end if end if end do end if end do wk3(i) = tmax kp(i) = ia kq(i) = ib iwk4(ind1) = .false. iwk5(ind2) = .false. ! ! Recalculate the cost differences. ! do j = 1, n if ( iwk4(j) ) then wk1(j) = wk1(j) + 2.0E+00 * cost(ip(j),ia) - 2.0*cost(ip(j),ib) end if if ( iwk5(j) ) then wk2(i) = wk2(j) + 2.0E+00 * cost(iq(j),ib) - 2.0 * cost(iq(j),ia) end if end do end do ! ! Choose K such that WK3(K) is maximal. ! tmax = small do i = 1, n tot1 = 0.0E+00 do j = 1, i tot1 = tot1 + wk3(j) end do if ( tot1 > tmax ) then tmax = tot1 k = i end if end do ! ! Exchange the two elements found above. ! Iterate until no reduction in cost can be obtained. ! if ( tmax <= 0.0E+00 ) then exit end if do i = 1, n if ( i <= k ) then ip(i) = kq(i) iq(i) = kp(i) else ip(i) = kp(i) iq(i) = kq(i) end if end do end do return end subroutine pmatch ( n, cost, ipair ) ! !******************************************************************************* ! !! PMATCH finds a minimum weight perfect matching in a graph. ! ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 15 September 1999 ! ! Parameters: ! ! Input, integer N, the number of nodes in the complete graph. ! N is assumed to be even. ! ! Input, real COST(N*(N-1)/2), the strict upper triangle of the cost ! matrix, stored by rows. In other words, the first elements of ! COST are the costs C(1,2), C(1,3), ..., C(1,N), C(2,3), C(2,4), ! ..., C(2,N). ! ! Output, integer IPAIR(N), contains the minimum weight perfect matching. ! Node I is connected to node IPAIR(I), for I = 1 to N. ! implicit none ! integer n ! real cost((n*(n-1))/2) real cst real cstlow real cswk real cwk2 integer i integer ihead integer index integer ipair(n) integer isub integer j integer jwk1(n) integer jwk2(n) integer jwk3(n) integer jwk4(n) integer jwk5(n) integer jwk6(n) integer jwk7(n) integer jwk8(n) integer jwk9(n) integer kk1 integer kk2 integer kk3 integer kk4 integer kk5 integer kk6 integer ll1 integer ll2 integer ll3 integer ll4 integer ll5 integer max integer min integer mm1 integer mm2 integer mm3 integer mm4 integer mm5 integer nn integer nn2 real xcst real value real work1(n) real work2(n) real work3(n) real work4(n) real xwk2 real xwk3 real xwork ! ! Initialization ! nn2 = ( n * ( n - 1 ) ) / 2 jwk1(2) = 0 do i = 3, n jwk1(i) = jwk1(i-1) + i - 2 end do ihead = n + 2 do i = 1, n jwk2(i) = i jwk3(i) = i jwk5(i) = i end do jwk4(1:n) = 0 jwk6(1:n) = ihead jwk7(1:n) = ihead jwk8(1:n) = ihead ipair(1:n) = ihead work1(1:n) = huge ( 1.0E+00 ) work2(1:n) = 0.0E+00 work3(1:n) = 0.0E+00 work4(1:n) = huge ( 1.0E+00 ) ! ! Start procedure. ! do i = 1, n if ( ipair(i) == ihead ) then nn = 0 cwk2 = huge ( 1.0E+00 ) do j = 1, n min = i max = j if ( i /= j ) then if ( i > j ) then max = i min = j end if isub = jwk1(max) + min xcst = cost(isub) cswk = cost(isub) - work2(j) if ( cswk <= cwk2 ) then if ( cswk == cwk2 ) then if ( nn == 0 ) then go to 30 else go to 40 end if end if cwk2 = cswk nn = 0 30 continue if ( ipair(j) == ihead ) then nn = j end if end if end if 40 continue end do if ( nn /= 0 ) then work2(i) = cwk2 ipair(i) = nn ipair(nn) = i end if end if end do ! ! Initial labeling. ! nn = 0 do i = 1, n if ( ipair(i) == ihead ) then nn = nn + 1 jwk6(i) = 0 work4(i) = 0.0E+00 xwk2 = work2(i) do j = 1, n min = i max = j if ( i /= j ) then if ( i > j ) then max = i min = j end if isub = jwk1(max) + min xcst = cost(isub) cswk = cost(isub) - xwk2 - work2(j) if ( cswk < work1(j) ) then work1(j) = cswk jwk4(j) = i end if end if end do end if end do if ( nn <= 1 ) then go to 340 end if ! ! Examine the labeling and prepare for the next step. ! 80 continue cstlow = huge ( 1.0E+00 ) do i = 1, n if ( jwk2(i) == i ) then cst = work1(i) if ( jwk6(i) < ihead ) then cst = 0.5E+00 * ( cst + work4(i) ) if ( cst <= cstlow ) then index = i cstlow = cst end if else if ( jwk7(i) < ihead ) then if ( jwk3(i) /= i ) then cst = cst + work2(i) if ( cst < cstlow ) then index = i cstlow = cst end if end if else if ( cst < cstlow ) then index = i cstlow = cst end if end if end if end if end do if ( jwk7(index) < ihead ) then go to 190 end if if ( jwk6(index) < ihead ) then ll4 = jwk4(index) ll5 = jwk5(index) kk4 = index kk1 = kk4 kk5 = jwk2(ll4) kk2 = kk5 do jwk7(kk1) = kk2 mm5 = jwk6(kk1) if ( mm5 == 0 ) then exit end if kk2 = jwk2(mm5) kk1 = jwk7(kk2) kk1 = jwk2(kk1) end do ll2 = kk1 kk1 = kk5 kk2 = kk4 110 continue if ( jwk7(kk1) >= ihead ) then jwk7(kk1) = kk2 mm5 = jwk6(kk1) if ( mm5 == 0 ) then go to 280 end if kk2 = jwk2(mm5) kk1 = jwk7(kk2) kk1 = jwk2(kk1) go to 110 end if 120 continue if ( kk1 == ll2 ) then go to 130 end if mm5 = jwk7(ll2) jwk7(ll2) = ihead ll1 = ipair(mm5) ll2 = jwk2(ll1) go to 120 end if ! ! Growing an alternating tree, add two edges. ! jwk7(index) = jwk4(index) jwk8(index) = jwk5(index) ll1 = ipair(index) ll3 = jwk2(ll1) work4(ll3) = cstlow jwk6(ll3) = ipair(ll3) call pmatch_sub_b ( ll3, n, nn2, cost, jwk1, jwk2, jwk3, jwk4, & jwk5, jwk7, jwk9, work1, work2, work3, work4 ) go to 80 ! ! Shrink a blossom. ! 130 continue xwork = work2(ll2) + cstlow - work4(ll2) work2(ll2) = 0.0E+00 mm1 = ll2 do work3(mm1) = work3(mm1) + xwork mm1 = jwk3(mm1) if ( mm1 == ll2 ) then exit end if end do mm5 = jwk3(ll2) if ( ll2 /= kk5 ) go to 160 150 continue kk5 = kk4 kk2 = jwk7(ll2) 160 continue jwk3(mm1) = kk2 ll1 = ipair(kk2) jwk6(kk2) = ll1 xwk2 = work2(kk2) + work1(kk2) - cstlow mm1 = kk2 do mm2 = mm1 work3(mm2) = work3(mm2) + xwk2 jwk2(mm2) = ll2 mm1 = jwk3(mm2) if ( mm1 == kk2 ) then exit end if end do jwk5(kk2) = mm2 work2(kk2) = xwk2 kk1 = jwk2(ll1) jwk3(mm2) = kk1 xwk2 = work2(kk1) + cstlow - work4(kk1) mm2 = kk1 do mm1 = mm2 work3(mm1) = work3(mm1) + xwk2 jwk2(mm1) = ll2 mm2 = jwk3(mm1) if ( mm2 == kk1 ) then exit end if end do jwk5(kk1) = mm1 work2(kk1) = xwk2 if ( kk5 /= kk1 ) then kk2 = jwk7(kk1) jwk7(kk1) = jwk8(kk2) jwk8(kk1) = jwk7(kk2) go to 160 end if if ( kk5 /= index ) then jwk7(kk5) = ll5 jwk8(kk5) = ll4 if ( ll2 /= index ) then go to 150 end if else jwk7(index) = ll4 jwk8(index) = ll5 end if jwk3(mm1) = mm5 kk4 = jwk3(ll2) jwk4(kk4) = mm5 work4(kk4) = xwork jwk7(ll2) = ihead work4(ll2) = cstlow call pmatch_sub_b ( ll2, n, nn2, cost, jwk1, jwk2, jwk3, jwk4, & jwk5, jwk7, jwk9, work1, work2, work3, work4 ) go to 80 ! ! Expand a T-labeled blossom. ! 190 continue kk4 = jwk3(index) kk3 = kk4 ll4 = jwk4(kk4) mm2 = kk4 do mm1 = mm2 ll5 = jwk5(mm1) xwk2 = work2(mm1) do jwk2(mm2) = mm1 work3(mm2)= work3(mm2) - xwk2 if ( mm2 == ll5 ) then exit end if mm2 = jwk3(mm2) end do mm2 = jwk3(ll5) jwk3(ll5) = mm1 if ( mm2 == ll4 ) then exit end if end do xwk2 = work4(kk4) work2(index) = xwk2 jwk3(index) = ll4 mm2 = ll4 do work3(mm2) = work3(mm2) - xwk2 if ( mm2 == index ) then exit end if mm2 = jwk3(mm2) end do mm1 = ipair(index) kk1 = jwk2(mm1) mm2 = jwk6(kk1) ll2 = jwk2(mm2) if ( ll2 /= index ) then kk2 = ll2 do mm5 = jwk7(kk2) kk1 = jwk2(mm5) if ( kk1 == index ) then exit end if kk2 = jwk6(kk1) kk2 = jwk2(kk2) end do jwk7(ll2) = jwk7(index) jwk7(index) = jwk8(kk2) jwk8(ll2) = jwk8(index) jwk8(index) = mm5 mm3 = jwk6(ll2) kk3 = jwk2(mm3) mm4 = jwk6(kk3) jwk6(ll2) = ihead ipair(ll2) = mm1 kk1 = kk3 do mm1 = jwk7(kk1) mm2 = jwk8(kk1) jwk7(kk1) = mm4 jwk8(kk1) = mm3 jwk6(kk1) = mm1 ipair(kk1) = mm1 kk2 = jwk2(mm1) ipair(kk2) = mm2 mm3 = jwk6(kk2) jwk6(kk2) = mm2 if ( kk2 == index ) then exit end if kk1 = jwk2(mm3) mm4 = jwk6(kk1) jwk7(kk2) = mm3 jwk8(kk2) = mm4 end do end if mm2 = jwk8(ll2) kk1 = jwk2(mm2) work1(kk1) = cstlow kk4 = 0 if ( kk1 /= ll2 ) then mm1 = jwk7(kk1) kk3 = jwk2(mm1) jwk7(kk1) = jwk7(ll2) jwk8(kk1) = mm2 do mm5 = jwk6(kk1) jwk6(kk1) = ihead kk2 = jwk2(mm5) mm5 = jwk7(kk2) jwk7(kk2) = ihead kk5 = jwk8(kk2) jwk8(kk2) = kk4 kk4 = kk2 work4(kk2) = cstlow kk1 = jwk2(mm5) work1(kk1) = cstlow if ( kk1 == ll2 ) then exit end if end do jwk7(ll2) = kk5 jwk8(ll2) = mm5 jwk6(ll2) = ihead if ( kk3 == ll2 ) go to 270 end if kk1 = 0 kk2 = kk3 do mm5 = jwk6(kk2) jwk6(kk2) = ihead jwk7(kk2) = ihead jwk8(kk2) = kk1 kk1 = jwk2(mm5) mm5 = jwk7(kk1) jwk6(kk1) = ihead jwk7(kk1) = ihead jwk8(kk1) = kk2 kk2 = jwk2(mm5) if ( kk2 == ll2 ) then exit end if end do call pmatch_sub_a ( kk1, n, nn2, cost, jwk1, jwk2, jwk3, jwk4, & jwk5, jwk6, jwk8, work1, work2, work3, work4 ) 270 continue if ( kk4 == 0 ) then go to 80 end if ll2 = kk4 call pmatch_sub_b ( ll2, n, nn2, cost, jwk1, jwk2, jwk3, jwk4, & jwk5, jwk7, jwk9, work1, work2, work3, work4 ) kk4 = jwk8(ll2) jwk8(ll2) = ihead go to 270 ! ! Augmentation of the matching. ! ! Exchange the matching and non-matching edges along the augmenting path. ! 280 continue ll2 = kk4 mm5 = ll4 do kk1 = ll2 do ipair(kk1) = mm5 mm5 = jwk6(kk1) jwk7(kk1) = ihead if ( mm5 == 0 ) then exit end if kk2 = jwk2(mm5) mm1 = jwk7(kk2) mm5 = jwk8(kk2) kk1 = jwk2(mm1) ipair(kk2) = mm1 end do if ( ll2 /= kk4 ) then exit end if ll2 = kk5 mm5 = ll5 end do ! ! Remove all labels of non-exposed base nodes. ! do i = 1, n if ( jwk2(i) == i ) then if ( jwk6(i) < ihead ) then cst = cstlow - work4(i) work2(i) = work2(i) + cst jwk6(i) = ihead if ( ipair(i) /= ihead ) then work4(i) = huge ( 1.0E+00 ) else jwk6(i) = 0 work4(i) = 0.0E+00 end if else if ( jwk7(i) < ihead ) then cst = work1(i) - cstlow work2(i) = work2(i) + cst jwk7(i) = ihead jwk8(i) = ihead end if work4(i) = huge ( 1.0E+00 ) end if work1(i) = huge ( 1.0E+00 ) end if end do nn = nn - 2 ! ! Determine the new WORK1 values. ! if ( nn <= 1 ) then go to 340 end if do i = 1, n kk1 = jwk2(i) if ( jwk6(kk1) == 0 ) then xwk2 = work2(kk1) xwk3 = work3(i) do j = 1, n kk2 = jwk2(j) if ( kk1 /= kk2 ) then min = i max = j if ( i /= j ) then if ( i > j ) then max = i min = j end if isub = jwk1(max) + min xcst = cost(isub) cswk = cost(isub) - xwk2 - xwk3 cswk = cswk - work2(kk2) - work3(j) if ( cswk < work1(kk2) ) then jwk4(kk2) = i jwk5(kk2) = j work1(kk2) = cswk end if end if end if end do end if end do go to 80 ! ! Generate the original graph by expanding all shrunken blossoms. ! 340 continue value = 0.0E+00 do i = 1, n if ( jwk2(i) == i ) then if ( jwk6(i) >= 0 ) then kk5 = ipair(i) kk2 = jwk2(kk5) kk4 = ipair(kk2) jwk6(i) = -1 jwk6(kk2) = -1 min = kk4 max = kk5 if ( kk4 /= kk5 ) then if ( kk4 > kk5 ) then max = kk4 min = kk5 end if isub = jwk1(max) + min xcst = cost(isub) value = value + xcst end if end if end if end do do i = 1, n 360 continue ll2 = jwk2(i) if ( ll2 == i ) then cycle end if mm2 = jwk3(ll2) ll4 = jwk4(mm2) kk3 = mm2 xwork = work4(mm2) do mm1 = mm2 ll5 = jwk5(mm1) xwk2 = work2(mm1) do jwk2(mm2) = mm1 work3(mm2) = work3(mm2) - xwk2 if ( mm2 == ll5 ) then exit end if mm2 = jwk3(mm2) end do mm2 = jwk3(ll5) jwk3(ll5) = mm1 if ( mm2 == ll4 ) then exit end if end do work2(ll2) = xwork jwk3(ll2) = ll4 mm2 = ll4 do work3(mm2) = work3(mm2) - xwork if ( mm2 == ll2 ) then exit end if mm2 = jwk3(mm2) end do mm5 = ipair(ll2) mm1 = jwk2(mm5) mm1 = ipair(mm1) kk1 = jwk2(mm1) if ( ll2 /= kk1 ) then ipair(kk1) = mm5 kk3 = jwk7(kk1) kk3 = jwk2(kk3) do mm3 = jwk6(kk1) kk2 = jwk2(mm3) mm1 = jwk7(kk2) mm2 = jwk8(kk2) kk1 = jwk2(mm1) ipair(kk1) = mm2 ipair(kk2) = mm1 min = mm1 max = mm2 if ( mm1 == mm2 ) go to 360 if ( mm1 > mm2 ) then max = mm1 min = mm2 end if isub = jwk1(max) + min xcst = cost(isub) value = value + xcst if ( kk1 == ll2 ) then exit end if end do if ( kk3 == ll2 ) go to 360 end if do kk5 = jwk6(kk3) kk2 = jwk2(kk5) kk6 = jwk6(kk2) min = kk5 max = kk6 if ( kk5 == kk6 ) go to 360 if ( kk5 > kk6 ) then max = kk5 min = kk6 end if isub = jwk1(max) + min xcst = cost(isub) value = value + xcst kk6 = jwk7(kk2) kk3 = jwk2(kk6) if ( kk3 == ll2 ) then go to 360 end if end do end do return end subroutine r_swap ( x, y ) ! !******************************************************************************* ! !! R_SWAP switches two real values. ! ! ! Modified: ! ! 30 November 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real X, Y. On output, the values of X and ! Y have been interchanged. ! implicit none ! real x real y real z ! z = x x = y y = z return end subroutine rmat_print ( a, ihi, ilo, jhi, jlo, lda, ncol, nrow ) ! !******************************************************************************* ! !! RMAT_PRINT prints out a portion of a dense matrix. ! ! ! Modified: ! ! 01 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A(LDA,NCOL), an NROW by NCOL matrix to be printed. ! ! Input, integer IHI, ILO, the first and last rows to print. ! ! Input, integer JHI, JLO, the first and last columns to print. ! ! Input, integer LDA, the leading dimension of A. ! ! Input, integer NCOL, NROW, the number of rows and columns ! in the matrix A. ! implicit none ! integer, parameter :: incx = 5 integer lda integer ncol ! real a(lda,ncol) character ( len = 14 ) ctemp(incx) integer i integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2 integer j2hi integer j2lo integer jhi integer jlo integer nrow ! write ( *, '(a)' ) ' ' do j2lo = jlo, jhi, incx j2hi = j2lo + incx - 1 j2hi = min ( j2hi, ncol ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo write ( *, '(a)' ) ' ' do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(i7,7x)' ) j end do write ( *, '(''Columns '',5a14)' ) ( ctemp(j2), j2 = 1, inc ) write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ' i2lo = max ( ilo, 1 ) i2hi = min ( ihi, nrow ) do i = i2lo, i2hi do j2 = 1, inc j = j2lo - 1 + j2 if ( a(i,j) == 0.0E+00 ) then ctemp(j2) = ' 0.0' else write ( ctemp(j2), '(g14.6)' ) a(i,j) end if end do write ( *, '(i5,1x,5a14)' ) i, ( ctemp(j), j = 1, inc ) end do end do write ( *, '(a)' ) ' ' return end subroutine rvec_sort_heap_index_a ( n, a, indx ) ! !******************************************************************************* ! !! RVEC_SORT_HEAP_INDEX_A does an indexed heap ascending sort of a real vector. ! ! ! Discussion: ! ! The sorting is not actually carried out. Rather an index array is ! created which defines the sorting. This array may be used to sort ! or index the array, or to sort or index related arrays keyed on the ! original array. ! ! Once the index array is computed, the sorting can be carried out ! "implicitly: ! ! A(INDX(I)), I = 1 to N is sorted, ! ! or explicitly, by the call ! ! call RVEC_PERMUTE ( N, A, INDX ) ! ! after which A(I), I = 1 to N is sorted. ! ! Modified: ! ! 17 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of entries in the array. ! ! Input, real A(N), an array to be index-sorted. ! ! Output, integer INDX(N), contains the sort index. The ! I-th element of the sorted array is A(INDX(I)). ! implicit none ! integer n ! real a(n) real aval integer i integer indx(n) integer indxt integer ir integer j integer l ! do i = 1, n indx(i) = i end do l = n / 2 + 1 ir = n do if ( l > 1 ) then l = l - 1 indxt = indx(l) aval = a(indxt) else indxt = indx(ir) aval = a(indxt) indx(ir) = indx(1) ir = ir - 1 if ( ir == 1 ) then indx(1) = indxt exit end if end if i = l j = l + l do while ( j <= ir ) if ( j < ir ) then if ( a(indx(j)) < a(indx(j+1)) ) then j = j + 1 end if end if if ( aval < a(indx(j)) ) then indx(i) = indx(j) i = j j = j + j else j = ir + 1 end if end do indx(i) = indxt end do return end subroutine rvec_sort_heap_index_d ( n, a, indx ) ! !******************************************************************************* ! !! RVEC_SORT_HEAP_INDEX_D does an indexed heap descending sort of a real vector. ! ! ! Discussion: ! ! The sorting is not actually carried out. Rather an index array is ! created which defines the sorting. This array may be used to sort ! or index the array, or to sort or index related arrays keyed on the ! original array. ! ! Once the index array is computed, the sorting can be carried out ! "implicitly: ! ! A(INDX(I)), I = 1 to N is sorted, ! ! or explicitly, by the call ! ! call RVEC_PERMUTE ( N, A, INDX ) ! ! after which A(I), I = 1 to N is sorted. ! ! Modified: ! ! 07 November 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of entries in the array. ! ! Input, real A(N), an array to be index-sorted. ! ! Output, integer INDX(N), contains the sort index. The ! I-th element of the sorted array is A(INDX(I)). ! implicit none ! integer n ! real a(n) real aval integer i integer indx(n) integer indxt integer ir integer j integer l ! do i = 1, n indx(i) = i end do l = n / 2 + 1 ir = n do if ( l > 1 ) then l = l - 1 indxt = indx(l) aval = a(indxt) else indxt = indx(ir) aval = a(indxt) indx(ir) = indx(1) ir = ir - 1 if ( ir == 1 ) then indx(1) = indxt exit end if end if i = l j = l + l do while ( j <= ir ) if ( j < ir ) then if ( a(indx(j)) > a(indx(j+1)) ) then j = j + 1 end if end if if ( aval > a(indx(j)) ) then indx(i) = indx(j) i = j j = j + j else j = ir + 1 end if end do indx(i) = indxt end do return end subroutine simplex ( n, m, a, x, infs, wk2, wk13, wk14, wk16, iwk19, iwk20 ) ! !******************************************************************************* ! !! SIMPLEX is used by INT_LP. ! ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 13 September 2000 ! ! Parameters: ! ! Input, integer N, the number of variables. ! ! Input, integer M, the number of constraints. ! ! Input, real A(LDA,N+M+1), contains the coefficients of the constraints ! in the first M rows and N columns. There is also an M+1 row used ! for working storage. The last M+1 columns are also used as working ! storage. ! ! Output, real X(max ( M+1, N ) ), contains the solution. ! ! Output, integer INFS, error indicator. ! 0, no error, a solution was found. ! 1, the linear program is infeasible. No solution was found. ! ! Workspace, real WK2((M+1)*(M+1)). ! ! Workspace, real WK13(M+1). ! ! Workspace, real WK14(M+1). ! ! Workspace, real WK16(max(M+1,N)). ! ! Workspace, integer IWK19(M+1). ! ! Output, integer IWK20(N+M). ! implicit none ! integer m integer n ! real a(*) real beta real delta1 real delta2 real delta3 real delta4 real delta5 real delta6 real delta7 real delta8 real delta9 real eps integer i integer ic1 integer ic2 integer ic3 integer ic4 integer ic5 integer ic6 integer ic7 integer ic8 integer ic9 integer inc1 integer inc2 integer inc3 integer inc4 integer infs integer ispec integer it integer iter integer iwk19(m+1) integer iwk20(n+m) integer j integer jc1 integer jc2 integer jc3 integer jc4 integer jc5 integer k integer l integer ll integer loop1 integer loop2 integer loop3 integer, parameter :: maxit = 900 integer m2 integer m99 integer n99 real wk2(*) real wk13(m+1) real wk14(m+1) real wk16(n) real x(n) ! ! Compute the machine epsilon. ! eps = sqrt ( epsilon ( 1.0E+00 ) ) ! ! Simplex method for linear programming. ! n99 = n + m ic5 = m + 1 m99 = m + 1 ic6 = 2 ic7 = 1 inc1 = 0 k = 0 iter = 0 inc2 = 0 ic1 = 0 jc5 = 0 delta1 = eps delta2 = eps delta3 = - eps * 2.0E+00 jc1 = 0 delta4 = 0.0E+00 delta6 = 0.0E+00 m2 = m99 + m99 infs = 1 ispec = 7777 ! ! Start phase one. ! iwk19(1:m99) = 0 ic9 = 0 do j = 1, n99 iwk20(j) = 0 ic8 = ic9 + ic6 ll = ic9 + m99 jc2 = 0 do l = ic8, ll if ( a(l) /= 0.0E+00 ) then jc2 = jc1 + 1 jc1 = l end if end do if ( jc2 == 1 ) then jc3 = jc1 - ic9 if ( iwk19(jc3) == 0 ) then if ( a(jc1) * wk16(jc3) >= 0.0E+00 ) then iwk19(jc3) = j iwk20(j) = jc3 ic9 = ic9 + ic5 end if end if end if end do ! ! Form the inverse from IWK20. ! 40 continue loop3 = 1 loop2 = 1 if ( jc5 <= 0 ) then inc2 = 0 end if ! ! Old code: ! ! wk2(1:m2) ! ! New code: ! wk2(1:(m+1)*(m+1)) = 0.0E+00 ic8 = 1 do i = 1, m99 wk2(ic8) = 1.0E+00 ic8 = ic8 + m99 + 1 end do x(1:m99) = wk16(1:m99) do i = ic6, m99 if ( iwk19(i) /= 0 ) then iwk19(i) = ispec end if end do infs = 1 ic1 = 1 ! ! Form the inverse. ! 80 continue if ( iwk20(ic1) == 0 ) then go to 110 else go to 250 end if 90 continue ! ! Reset the artificials. ! delta5 = 0.0E+00 do i = ic6, m99 if ( iwk19(i) == ispec ) then if ( abs ( wk14(i) ) > delta5 ) then ic3 = i delta5 = abs ( wk14(i) ) end if end if end do if ( delta5 < delta1 ) then iwk20(ic1) = 0 else iwk19(ic3) = ic1 iwk20(ic1) = ic3 go to 350 end if 110 continue ic1 = ic1 + 1 if ( ic1 <= n99 ) then go to 80 end if 120 continue do i = 1, m99 if ( iwk19(i) == ispec ) then iwk19(i) = 0 end if end do loop1 = 1 loop2 = 2 loop3 = 2 ! ! Perform one iteration. ! 140 continue inc3 = 0 inc4 = 0 do i = ic6, m99 if ( abs ( x(i) ) < delta2 ) then x(i) = 0.0E+00 else if ( x(i) /= 0.0E+00 ) then if ( x(i) > 0.0E+00 ) then if ( iwk19(i) == 0 ) then inc3 = 1 end if else inc4 = 1 inc3 = 1 end if end if end if end do ! ! If infeasible, then invert again. ! if ( infs < inc3 ) then go to 40 end if if ( infs > inc3 ) then infs = 0 delta4 = 0.0E+00 end if ! ! Obtain prices. ! 160 continue ic8 = ic7 do j = 1, m99 wk13(j) = wk2(ic8) ic8 = ic8 + m99 end do ! ! Pricing. ! if ( infs /= 0 ) then wk13(1:m99) = wk13(1:m99) * delta4 do i = ic6, m99 ic8 = i if ( x(i) < 0.0E+00 ) then do j = 1, m99 wk13(j) = wk13(j) + wk2(ic8) ic8 = ic8 + m99 end do else if ( iwk19(i) == 0 ) then do j = 1, m99 write ( *, '(a,i6)' ) 'J=', j write ( *, '(a,g14.6)' ) 'WK13(j)=', wk13(j) write ( *, '(a,i6)' ) 'IC8=', ic8 write ( *, '(a,g14.6)' ) 'WK2(IC8)=', wk2(ic8) wk13(j) = wk13(j) - wk2(ic8) ic8 = ic8 + m99 end do end if end do end if ! ! Select entering column. ! ic1 = 0 delta7 = delta3 ic2 = 1 220 continue if ( iwk20(ic2) == 0 ) then go to 460 else go to 240 end if 230 continue if ( delta6 < delta7 ) then delta7 = delta6 ic1 = ic2 end if 240 continue ic2 = ic2 + 1 ! ! All costs are non-negative. ! if ( ic2 <= n99 ) then go to 220 end if if ( ic1 <= 0 ) then k = 3 + infs go to 340 end if ! ! Multiply the column into the basis inverse. ! 250 continue wk14(1:m99) = 0.0E+00 ic4 = ( ic1 - 1 ) * ic5 ll = 0 do i = 1, m99 ic4 = ic4 + 1 do j = 1, m99 ll = ll + 1 wk14(j) = wk14(j) + a(ic4) * wk2(ll) end do end do if ( loop2 == 1 ) then go to 90 else if ( loop2 == 3 ) then return end if ! ! Get the maximum value from the pivot column. ! 290 continue ic3 = 0 delta8 = 0.0E+00 jc3 = 0 do i = ic6, m99 if ( x(i) == 0.0E+00 ) then beta = abs ( wk14(i) ) if ( beta > delta1 ) then if ( iwk19(i) == 0 ) then go to 300 end if if ( jc3 == 0 ) then if ( wk14(i) > 0.0E+00 ) then 300 continue if ( jc3 /= 0 ) then if ( beta > delta8 ) then delta8 = beta ic3 = i end if else jc3 = 1 delta8 = beta ic3 = i end if end if end if end if end if end do ! ! Find the maximum pivot from the positive equations. ! if ( ic3 == 0 ) then delta8 = huge ( 1.0E+00 ) do it = ic6, m99 if ( wk14(it) > delta1 ) then if ( x(it) > 0.0E+00 ) then delta9 = x(it) / wk14(it) if ( delta9 < delta8 ) then delta8 = delta9 ic3 = it else if ( delta9 == delta8 ) then if ( iwk19(it) == 0 ) then delta8 = delta9 ic3 = it end if end if end if end if end if end do ! ! Find pivot among negative equations. ! if ( inc4 /= 0 ) then delta7 = - delta1 do i = ic6, m99 if ( x(i) < 0.0E+00 ) then if ( wk14(i) < delta7 ) then if ( wk14(i) * delta8 <= x(i) ) then delta7 = wk14(i) ic3 = i end if end if end if end do end if end if ! ! Test pivot. ! if ( ic3 <= 0 ) then k = 5 340 continue if ( delta4 /= 0.0E+00 ) then delta4 = 0.0E+00 go to 160 end if go to 400 end if ! ! Terminate if too many iterations. ! if ( iter < maxit ) then 350 continue beta = - wk14(ic3) wk14(ic3) = - 1.0E+00 ll = 0 ! ! Transform inverse. ! do l = ic3, m2, m99 if ( wk2(l) == 0.0E+00 ) then ll = ll + m99 else delta9 = wk2(l) / beta wk2(l) = 0.0E+00 do i = 1, m99 ll = ll + 1 wk2(ll) = wk2(ll) + delta9 * wk14(i) end do end if end do ! ! Transform X. ! delta9 = x(ic3) / beta x(ic3) = 0.0E+00 x(1:m99) = x(1:m99) + delta9 * wk14(1:m99) ! ! Restore the pivot. ! wk14(ic3) = - beta if ( loop3 == 1 ) then go to 120 end if 390 continue jc3 = iwk19(ic3) if ( jc3 > 0 ) then iwk20(jc3) = 0 end if iwk20(ic1) = ic3 iwk19(ic3) = ic1 jc5 = 0 iter = iter + 1 inc2 = inc2 + 1 if ( inc1 == inc2 ) then go to 40 else go to 140 end if end if k = 6 ! ! Error checking. ! 400 continue loop1 = 2 wk14(1:m99) = - wk16(1:m99) do i = 1, m99 jc4 = iwk19(i) if ( jc4 /= 0 ) then jc3 = ic5 * ( jc4 - 1 ) wk14(1:m99) = wk14(1:m99) + x(i) * a(jc3+1:jc3+m99) end if end do ! ! Set error flags. ! i = 1 440 continue ic2 = iwk19(i) if ( ic2 == 0 ) then 450 continue i = i + 1 if ( i <= m99 ) then go to 440 end if if ( jc5 == 0 ) then jc5 = 4 end if if ( k /= 5 ) then return end if loop2 = 3 go to 250 end if ! ! Price out a column. ! 460 continue ll = ( ic2 - 1 ) * ic5 delta6 = 0.0E+00 do ic8 = 1, m99 ll = ll + 1 if ( a(ll) /= 0.0E+00 ) then delta6 = delta6 + wk13(ic8) * a(ll) end if end do if ( loop1 == 1 ) then go to 230 else go to 450 end if end subroutine steiner ( n, m, inode, jnode, arcost, ns, spoint, nsp, & istree, jstree, xlen ) ! !******************************************************************************* ! !! STEINER is a heuristic algorithm for the minimal Steiner tree problem. ! ! ! Discussion: ! ! Consider an undirected graph G whose edges have specified lengths. ! Let S be a subset of the nodes. The Steiner tree problem is to find ! a tree of G that includes all the nodes S, and has minimal total ! distance. The nodes in the set S are called Steiner points. ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! ! Input, integer N, the number of nodes in the graph. ! ! Input, integer M, the number of edges in the graph. ! ! Input, integer INODE(M), JNODE(M), the edges of the graph, ! described by pairs of nodes. ! ! Input, real ARCOST(M), the length of each edge. ! ! Input, integer NS, the number of nodes which are Steiner points. ! ! Input, logical SPOINT(N), is TRUE for each node that is a Steiner point. ! ! Output, integer NSP, the number of edges in the Steiner tree. ! ! Output, integer ISTREE(N), JSTREE(N), the edges of the Steiner tree, ! stored as pairs of nodes, in entries 1 through NSP. ! ! Output, real XLEN, the total length of the edges in the Steiner tree. ! implicit none ! integer m integer n integer ns ! real arcost(m) integer i integer ic1 integer ic2 integer iflag integer ii integer ij integer il integer in integer inode(m) integer istree(n) integer it integer iv1 integer iv2 integer iwk1(n) integer iwk2(n) logical iwk3(n) integer iwk4(n) integer iwk5(n) integer iwk6(n) integer iwk7(n) integer iwk8(m) integer iwk9(m) integer iwk10(m) integer iwk11(m) integer iwk12(ns) integer j integer jj integer jk integer jnode(m) integer jstree(n) integer k integer kk integer kk1 integer l integer ll integer mx integer np integer nsp integer nump logical spoint(n) real work13(n) real work14(m) real xlen ! ! Identify the Steiner points. ! it = 0 do i = 1, n if ( spoint(i) ) then it = it + 1 iwk12(it) = i end if end do ! ! Construct the complete graph for the Steiner points. ! il = 0 ii = 2 do i = 1, ns - 1 do j = ii, ns if ( i /= j ) then il = il + 1 iwk10(il) = iwk12(i) iwk11(il) = iwk12(j) call graph_arc_min_path ( n, m, inode, jnode, arcost, iwk12(i), & iwk12(j), nump, iwk4, xlen ) work14(il) = xlen endif end do ii = ii + 1 end do ll = ( ns * ( ns - 1 ) ) / 2 mx = 0 do i = 1, it mx = max ( mx, iwk12(i) ) end do ! ! Construct the complete graph for the Steiner points. ! ! Find a minimum spanning tree of the complete graph. ! call graph_arc_min_span_tree ( mx, ll, iwk10, iwk11, work14, & iwk1, iwk2, xlen ) iwk10(1:m) = 0 ! ! Construct the subgraph by replacing each edge of the ! minimum spanning tree by its corresponding shortest path. ! do i = 1, ns - 1 ii = iwk1(i) jj = iwk2(i) call graph_arc_min_path ( n, m, inode, jnode, arcost, ii, & jj, nump, iwk4, xlen ) do ij = 1, nump - 1 iv1 = iwk4(ij) iv2 = iwk4(ij+1) do jk = 1, m if ( ( inode(jk) == iv1 .and. jnode(jk) == iv2 ) .or. & ( inode(jk) == iv2 .and. jnode(jk) == iv1 ) ) then kk1 = jk end if end do iwk10(kk1) = - 1 end do end do do i = 1, m if ( iwk10(i) == 0 ) then inode(i) = 0 jnode(i) = 0 end if end do ! ! Find a minimum spanning tree of the subgraph. ! call graph_arc_min_span_tree ( n, m, inode, jnode, arcost, & iwk1, iwk2, xlen ) in = 0 do i = 1, n if ( iwk1(i) /= 0 ) then in = in + 1 end if end do ! ! Construct a Steiner tree by deleting edges, if ! such that all leaves are Steiner points. ! do j = 1, in iflag = 0 do i = 1, in if ( iwk1(i) /= 0 .and. iwk2(i) /= 0 ) then l = iwk2(i) ic1 = 0 do k = 1, in if ( iwk1(k) == l ) then ic1 = ic1 + 1 iwk9(ic1) = k iwk8(ic1) = iwk2(k) endif if ( iwk2(k) == l ) then ic1 = ic1 + 1 iwk9(ic1) = k iwk8(ic1) = iwk1(k) endif end do l = iwk1(i) ic2 = 0 do k = 1, in if ( iwk1(k) == l ) then ic2 = ic2 + 1 iwk9(ic2) = k iwk8(ic2) = iwk2(k) endif if ( iwk2(k) == l ) then ic2 = ic2 + 1 iwk9(ic2) = k iwk8(ic2) = iwk1(k) endif end do ii = iwk1(i) jj = iwk2(i) if ( ( ic1 == 1 .and. .not. spoint(jj) ) .or. & ( ic2 == 1 .and. .not. spoint(ii) ) ) then do k = 1, m if ( inode(k) == iwk1(i) .and. jnode(k) == iwk2(i) ) then kk = k exit end if end do iwk1(i) = 0 iwk2(i) = 0 iflag = 1 xlen = xlen - arcost(kk) end if end if end do if ( iflag == 0 ) then exit end if end do ! ! Store the solution. ! i = 0 do j = 1, in if ( iwk1(j) /= 0 .and. iwk2(j) /= 0 ) then i = i + 1 istree(i) = iwk1(j) jstree(i) = iwk2(j) endif end do nsp = i return end subroutine pmatch_sub_a ( kk, n, nn2, cost, jwk1, jwk2, jwk3, jwk4, & jwk5, jwk6, jwk8, work1, work2, work3, work4 ) ! !******************************************************************************* ! !! PMATCH_SUB_A is used by PMATCH. ! ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! implicit none ! integer n integer nn2 ! real cost(nn2) real cstwk real cswk integer i integer ihead integer isub integer j integer jj1 integer jj2 integer jj3 integer jj4 integer jwk1(n) integer jwk2(n) integer jwk3(n) integer jwk4(n) integer jwk5(n) integer jwk6(n) integer jwk8(n) integer kk integer max integer min real work1(n) real work2(n) real work3(n) real work4(n) real xwk2 real xwk3 ! ihead = n + 2 do jj1 = kk kk = jwk8(jj1) jwk8(jj1) = ihead cstwk = huge ( 1.0E+00 ) jj3 = 0 jj4 = 0 j = jj1 xwk2 = work2(jj1) do xwk3 = work3(j) do i = 1, n jj2 = jwk2(i) if ( jwk6(jj2) < ihead ) then min = j max = i if ( j /= i ) then if ( j > i ) then max = j min = i end if isub = jwk1(max) + min cswk = cost(isub) - xwk2 - xwk3 - work2(jj2) - work3(i) + work4(jj2) if ( cswk < cstwk ) then jj3 = i jj4 = j cstwk = cswk end if end if end if end do j = jwk3(j) if ( j == jj1 ) then exit end if end do jwk4(jj1) = jj3 jwk5(jj1) = jj4 work1(jj1) = cstwk if ( kk == 0 ) then exit end if end do return end subroutine pmatch_sub_b ( kk, n, nn2, cost, jwk1, jwk2, jwk3, jwk4, & jwk5, jwk7, jwk9, work1, work2, work3, work4 ) ! !******************************************************************************* ! !! PMATCH_SUB_B is used by PMATCH. ! ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 11 September 1999 ! ! Parameters: ! implicit none ! integer n integer nn2 ! real cost(nn2) real cswk integer i integer ihead integer ii integer isub integer jj1 integer jj2 integer jj3 integer jwk1(n) integer jwk2(n) integer jwk3(n) integer jwk4(n) integer jwk5(n) integer jwk7(n) integer jwk9(n) integer kk integer max integer min real work1(n) real work2(n) real work3(n) real work4(n) real xwk1 real xwk2 ! ihead = n + 2 xwk1 = work4(kk) - work2(kk) work1(kk) = huge ( 1.0E+00 ) xwk2 = xwk1 - work3(kk) jwk7(kk) = 0 ii = 0 do i = 1, n jj3 = jwk2(i) if ( jwk7(jj3) >= ihead ) then ii = ii + 1 jwk9(ii) = i min = kk max = i if ( kk /= i ) then if ( kk > i ) then max = kk min = i end if isub = jwk1(max) + min cswk = cost(isub) + xwk2 - work2(jj3) - work3(i) if ( cswk < work1(jj3) ) then jwk4(jj3) = kk jwk5(jj3) = i work1(jj3) = cswk end if end if end if end do jwk7(kk) = ihead jj1 = kk jj1 = jwk3(jj1) if ( jj1 == kk ) then return end if do xwk2 = xwk1 - work3(jj1) do i = 1, ii jj2 = jwk9(i) jj3 = jwk2(jj2) min = jj1 max = jj2 if ( jj1 /= jj2 ) then if ( jj1 > jj2 ) then max = jj1 min = jj2 end if isub = jwk1(max) + min cswk = cost(isub) + xwk2 - work2(jj3) - work3(jj2) if ( cswk < work1(jj3) ) then jwk4(jj3) = jj1 jwk5(jj3) = jj2 work1(jj3) = cswk end if end if end do jj1 = jwk3(jj1) if ( jj1 == kk ) then exit end if end do return end subroutine timestamp ( ) ! !******************************************************************************* ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! ! Example: ! ! May 31 2001 9:45:54.872 AM ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none ! character ( len = 8 ) ampm integer d character ( len = 8 ) date integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s character ( len = 10 ) time integer values(8) integer y character ( len = 5 ) zone ! call date_and_time ( date, time, zone, values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end subroutine tsp ( nnode, dist, lda, isol ) ! !******************************************************************************* ! !! TSP is a heuristic algorithm for the traveling salesman problem. ! ! ! Discussion: ! ! Let G be a complete graph with an associated distance matrix ! DIST(I,J) on its edges. The traveling salesman problem is to ! start from a node in G, visit all the other nodes exactly once, ! and return to the starting node in such a way that the total ! traveled distance is a minimum. ! ! In general, it is very hard to develop efficient algorithms ! that yield good approximate solutions to this problem. However, ! in the special case where the distance matrix is symmetric: ! ! DIST(I,J) = DIST(J,I) ! ! and satisfies the triangle inequality: ! ! DIST(I,J) + DIST(J,K) >= DIST(I,K), ! ! solutions close to the optimum value can be found in a relatively ! short time. The algorithm here will find a circuit of no worse ! than 3/2 the optimum length, in polynomial time. ! ! Reference: ! ! H T Lau, ! Combinatorial Heuristic Algorithms in FORTRAN, ! Springer Verlag, 1986. ! ! Modified: ! ! 21 July 2000 ! ! Parameters: ! ! Input, integer NNODE, the number of nodes. ! ! Input, real DIST(LDA,NNODE), the distance between each pair of nodes. ! ! Input, integer LDA, the leading dimension of the DIST array. ! LDA must be at least NNODE. ! ! Output, integer ISOL(NNODE), contains the nodes in the order in which ! they should be visited. ! implicit none ! integer lda integer nnode ! real dist(lda,nnode) integer i integer ict integer ii integer inode((3*nnode)/2) integer isol(nnode) integer iwk10((3*nnode)/2) integer iwk11(nnode) integer iwk18(nnode) integer j integer jj integer jnode((3*nnode)/2) integer k integer nedge real wk1( ( nnode * ( nnode - 1 ) ) / 2) real wk2(nnode) real wk3(nnode) real wk4(nnode) real wk5(nnode) ! ! From the graph described by DIST, construct a minimum length spanning ! tree, described by INODE, JNODE. ! call graph_dist_min_span_tree3 ( lda, nnode, dist, inode, jnode ) ! ! Compute the degree (in the spanning tree) of each node. ! iwk10(1:nnode) = 0 do i = 1, nnode - 1 ii = inode(i) jj = jnode(i) iwk10(ii) = iwk10(ii) + 1 iwk10(jj) = iwk10(jj) + 1 end do ! ! Construct a list of the nodes of odd degree. ! ict = 0 do i = 1, nnode if ( mod ( iwk10(i), 2 ) /= 0 ) then ict = ict + 1 iwk11(ict) = i end if end do ! ! Determine a minimum weight perfect matching for the set of odd-degree nodes. ! k = 0 do i = 2, ict do j = 1, i - 1 k = k + 1 wk1(k) = dist(iwk11(j),iwk11(i)) end do end do call pmatch ( ict, wk1, iwk18 ) ! ! Store up the edges in the perfect matching. ! do i = 1, ict iwk18(i) = iwk11(iwk18(i)) end do iwk10(1:nnode) = 0 k = nnode - 1 do i = 1, ict if ( iwk10(iwk11(i)) == 0 ) then iwk10(iwk11(i)) = 1 iwk10(iwk18(i)) = 1 k = k + 1 inode(k) = iwk11(i) jnode(k) = iwk18(i) end if end do ! ! Find an Euler circuit. ! nedge = nnode - 1 + ( ict / 2 ) call graph_arc_euler_circ ( nnode, nedge, inode, jnode, iwk10 ) ! ! Form the Hamiltonian circuit. ! isol(1) = iwk10(1) j = 2 isol(j) = iwk10(j) k = 2 do 90 continue j = j + 1 do i = 1, j - 1 if ( iwk10(j) == iwk10(i) ) then go to 90 end if end do k = k + 1 isol(k) = iwk10(j) if ( k >= nnode ) then exit end if end do return end