/* relax4.c (relaxation method of Bertsekas and Tseng) */ /*********************************************************************** * This code is part of GLPK (GNU Linear Programming Kit). * * THIS CODE IS THE RESULT OF TRANSLATION OF THE FORTRAN CODE RELAX4. * * THE TRANSLATION HAS BEEN DONE WITH THE PERMISSION OF THE AUTHOR OF * THE ORIGINAL FORTRAN CODE PROF. DIMITRI P. BERTSEKAS, MASSACHUSETTS * INSTITUTE OF TECHNOLOGY, CAMBRIDGE, MASSACHUSETTS, USA. * * The translation was made by Andrew Makhorin . * * GLPK is free software: you can redistribute it and/or modify it * under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * GLPK is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * License for more details. * * You should have received a copy of the GNU General Public License * along with GLPK. If not, see . ***********************************************************************/ #include "relax4.h" /*********************************************************************** * WARNING * * A serious bug was *tentatively* fixed in this code (see #if/#endif * marked by 'mao'). * * This bug is inherited from the original Fortran version of the * RELAX-IV code. Unfortunately, the code is very intricate, so this * bug is still under investigation. Thanks to Sylvain Fournier for bug * report. * * RELAX-IV bug details * -------------------- * In the original RELAX-IV code there are four similar fragments in * subroutines ascnt1 and ascnt2 like this: * * C * C DECREASE THE PRICES OF THE SCANNED NODES BY DELPRC. * C ADJUST FLOW TO MAINTAIN COMPLEMENTARY SLACKNESS WITH * C THE PRICES. * C * NB = 0 * DO 6 I=1,NSAVE * . . . * IF (RC(ARC).EQ.0) THEN * DELX=DELX+U(ARC) * NB = NB + 1 * PRDCSR(NB) = ARC * END IF * . . . * * On some instances the variable NB becomes greater than N (the number * of nodes) that leads to indexing error, because the array PRDCSR is * declared as array of N elements (more precisely, as array of MAXNN * elements, however, NB becomes even much greater than MAXNN). ***********************************************************************/ #define false 0 #define true 1 /*********************************************************************** * NAME * * RELAX-IV (version of October 1994) * * PURPOSE * * This routine implements the relaxation method of Bertsekas and Tseng * (see [1], [2]) for linear cost ordinary network flow problems. * * [1] Bertsekas, D. P., "A Unified Framework for Primal-Dual Methods" * Mathematical Programming, Vol. 32, 1985, pp. 125-145. * [2] Bertsekas, D. P., and Tseng, P., "Relaxation Methods for * Minimum Cost" Operations Research, Vol. 26, 1988, pp. 93-114. * * The relaxation method is also described in the books: * * [3] Bertsekas, D. P., "Linear Network Optimization: Algorithms and * Codes" MIT Press, 1991. * [4] Bertsekas, D. P. and Tsitsiklis, J. N., "Parallel and Distributed * Computation: Numerical Methods", Prentice-Hall, 1989. * [5] Bertsekas, D. P., "Network Optimization: Continuous and Discrete * Models", Athena Scientific, 1998. * * RELEASE NOTE * * This version of relaxation code has option for a special crash * procedure for the initial price-flow pair. This is recommended for * difficult problems where the default initialization results in long * running times. crash = 1 corresponds to an auction/shortest path * method * * These initializations are recommended in the absence of any prior * information on a favorable initial flow-price vector pair that * satisfies complementary slackness. * * The relaxation portion of the code differs from the code RELAXT-III * and other earlier relaxation codes in that it maintains the set of * nodes with nonzero deficit in a fifo queue. Like its predecessor * RELAXT-III, this code maintains a linked list of balanced (i.e., of * zero reduced cost) arcs so to reduce the work in labeling and * scanning. Unlike RELAXT-III, it does not use selectively shortest * path iterations for initialization. * * SOURCE * * The original Fortran code was written by Dimitri P. Bertsekas and * Paul Tseng, with a contribution by Jonathan Eckstein in the phase II * initialization. The original Fortran routine AUCTION was written by * Dimitri P. Bertsekas and is based on the method described in the * paper: * * [6] Bertsekas, D. P., "An Auction/Sequential Shortest Path Algorithm * for the Minimum Cost Flow Problem", LIDS Report P-2146, MIT, * Nov. 1992. * * For inquiries about the original Fortran code, please contact: * * Dimitri P. Bertsekas * Laboratory for information and decision systems * Massachusetts Institute of Technology * Cambridge, MA 02139 * (617) 253-7267, dimitrib@mit.edu * * This code is the result of translation of the original Fortran code. * The translation was made by Andrew Makhorin . * * USER GUIDELINES * * This routine is in the public domain to be used only for research * purposes. It cannot be used as part of a commercial product, or to * satisfy in any part commercial delivery requirements to government * or industry, without prior agreement with the authors. Users are * requested to acknowledge the authorship of the code, and the * relaxation method. * * No modification should be made to this code other than the minimal * necessary to make it compatible with specific platforms. * * INPUT PARAMETERS (see notes 1, 2, 4) * * n = number of nodes * na = number of arcs * large = a very large integer to represent infinity * (see note 3) * repeat = true if initialization is to be skipped * (false otherwise) * crash = 0 if default initialization is used * 1 if auction initialization is used * startn[j] = starting node for arc j, j = 1,...,na * endn[j] = ending node for arc j, j = 1,...,na * fou[i] = first arc out of node i, i = 1,...,n * nxtou[j] = next arc out of the starting node of arc j, j = 1,...,na * fin[i] = first arc into node i, i = 1,...,n * nxtin[j] = next arc into the ending node of arc j, j = 1,...,na * * UPDATED PARAMETERS (see notes 1, 3, 4) * * rc[j] = reduced cost of arc j, j = 1,...,na * u[j] = capacity of arc j on input * and (capacity of arc j) - x[j] on output, j = 1,...,na * dfct[i] = demand at node i on input * and zero on output, i = 1,...,n * * OUTPUT PARAMETERS (see notes 1, 3, 4) * * x[j] = flow on arc j, j = 1,...,na * nmultinode = number of multinode relaxation iterations in RELAX4 * iter = number of relaxation iterations in RELAX4 * num_augm = number of flow augmentation steps in RELAX4 * num_ascnt = number of multinode ascent steps in RELAX4 * nsp = number of auction/shortest path iterations * * WORKING PARAMETERS (see notes 1, 4, 5) * * label[1+n], prdcsr[1+n], save[1+na], tfstou[1+n], tnxtou[1+na], * tfstin[1+n], tnxtin[1+na], nxtqueue[1+n], scan[1+n], mark[1+n], * extend_arc[1+n], sb_level[1+n], sb_arc[1+n] * * RETURNS * * 0 = normal return * 1,...,8 = problem is found to be infeasible * * NOTE 1 * * To run in limited memory systems, declare the arrays startn, endn, * nxtin, nxtou, fin, fou, label, prdcsr, save, tfstou, tnxtou, tfstin, * tnxtin, ddpos, ddneg, nxtqueue as short instead. * * NOTE 2 * * This routine makes no effort to initialize with a favorable x from * amongst those flow vectors that satisfy complementary slackness with * the initial reduced cost vector rc. If a favorable x is known, then * it can be passed, together with the corresponding arrays u and dfct, * to this routine directly. This, however, requires that the capacity * tightening portion and the flow initialization portion of this * routine (up to line labeled 90) be skipped. * * NOTE 3 * * All problem data should be less than large in magnitude, and large * should be less than, say, 1/4 the largest int of the machine used. * This will guard primarily against overflow in uncapacitated problems * where the arc capacities are taken finite but very large. Note, * however, that as in all codes operating with integers, overflow may * occur if some of the problem data takes very large values. * * NOTE 4 * * [This note being specific to Fortran was removed.-A.M.] * * NOTE 5 * * ddpos and ddneg are arrays that give the directional derivatives for * all positive and negative single-node price changes. These are used * only in phase II of the initialization procedure, before the linked * list of balanced arcs comes to play. Therefore, to reduce storage, * they are equivalence to tfstou and tfstin, which are of the same size * (number of nodes) and are used only after the tree comes into use. */ static void ascnt1(struct relax4_csa *csa, int dm, int *delx, int *nlabel, int *feasbl, int *svitch, int nscan, int curnode, int *prevnode); static void ascnt2(struct relax4_csa *csa, int dm, int *delx, int *nlabel, int *feasbl, int *svitch, int nscan, int curnode, int *prevnode); static int auction(struct relax4_csa *csa); int relax4(struct relax4_csa *csa) { /* input parameters */ int n = csa->n; int na = csa->na; int large = csa->large; int repeat = csa->repeat; int crash = csa->crash; int *startn = csa->startn; int *endn = csa->endn; int *fou = csa->fou; int *nxtou = csa->nxtou; int *fin = csa->fin; int *nxtin = csa->nxtin; /* updated parameters */ int *rc = csa->rc; int *u = csa->u; int *dfct = csa->dfct; /* output parameters */ int *x = csa->x; # define nmultinode (csa->nmultinode) # define iter (csa->iter) # define num_augm (csa->num_augm) # define num_ascnt (csa->num_ascnt) # define nsp (csa->nsp) /* working parameters */ int *label = csa->label; int *prdcsr = csa->prdcsr; int *save = csa->save; int *tfstou = csa->tfstou; int *tnxtou = csa->tnxtou; int *tfstin = csa->tfstin; int *tnxtin = csa->tnxtin; int *nxtqueue = csa->nxtqueue; char *scan = csa->scan; char *mark = csa->mark; int *ddpos = tfstou; int *ddneg = tfstin; /* local variables */ int arc, augnod, capin, capout, defcit, delprc, delx, dm, dp, dx, feasbl, i, ib, indef, j, lastqueue, maxcap, narc, nb, nlabel, node, node2, node_def, naugnod, nscan, num_passes, numnz, numnz_new, numpasses, nxtarc, nxtbrk, nxtnode, passes, pchange, posit, prevnode, prvarc, quit, rdcost, scapin, scapou, svitch, t, t1, t2, tmparc, tp, trc, ts; /*--------------------------------------------------------------*/ /* Initialization phase I */ /* In this phase, we reduce the arc capacities by as much as * possible without changing the problem; then we set the initial * flow array x, together with the corresponding arrays u and * dfct. */ /* This phase and phase II (from here up to line labeled 90) can * be skipped (by setting repeat to true) if the calling program * places in common user-chosen values for the arc flows, the * residual arc capacities, and the nodal deficits. When this is * done, it is critical that the flow and the reduced cost for * each arc satisfy complementary slackness and the dfct array * properly correspond to the initial arc/flows. */ if (repeat) goto L90; for (node = 1; node <= n; node++) { node_def = dfct[node]; ddpos[node] = node_def; ddneg[node] = -node_def; maxcap = 0; scapou = 0; for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { if (scapou <= large - u[arc]) scapou += u[arc]; else goto L10; } if (scapou <= large - node_def) capout = scapou + node_def; else goto L10; if (capout < 0) { /* problem is infeasible */ /* exogenous flow into node exceeds out capacity */ return 1; } scapin = 0; for (arc = fin[node]; arc > 0; arc = nxtin[arc]) { if (u[arc] > capout) u[arc] = capout; if (maxcap < u[arc]) maxcap = u[arc]; if (scapin <= large - u[arc]) scapin += u[arc]; else goto L10; } if (scapin <= large + node_def) capin = scapin - node_def; else goto L10; if (capin < 0) { /* problem is infeasible */ /* exogenous flow out of node exceeds in capacity */ return 2; } for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { if (u[arc] > capin) u[arc] = capin; } L10: ; } /*--------------------------------------------------------------*/ /* Initialization phase II */ /* In this phase, we initialize the prices and flows by either * calling the routine auction or by performing only single node * (coordinate) relaxation iterations. */ if (crash == 1) { nsp = 0; if (auction(csa) != 0) { /* problem is found to be infeasible */ return 3; } goto L70; } /* Initialize the arc flows to satisfy complementary slackness * with the prices. u[arc] is the residual capacity of arc, and * x[arc] is the flow. These two always add up to the total * capacity for arc. Also compute the directional derivatives for * each coordinate and compute the actual deficits. */ for (arc = 1; arc <= na; arc++) { x[arc] = 0; if (rc[arc] <= 0) { t = u[arc]; t1 = startn[arc]; t2 = endn[arc]; ddpos[t1] += t; ddneg[t2] += t; if (rc[arc] < 0) { x[arc] = t; u[arc] = 0; dfct[t1] += t; dfct[t2] -= t; ddneg[t1] -= t; ddpos[t2] -= t; } } } /* Make 2 or 3 passes through all nodes, performing only single * node relaxation iterations. The number of passes depends on the * density of the network. */ if (na > n * 10) numpasses = 2; else numpasses = 3; for (passes = 1; passes <= numpasses; passes++) for (node = 1; node <= n; node++) { if (dfct[node] == 0) continue; if (ddpos[node] <= 0) { /* Compute delprc, the stepsize to the next breakpoint in * the dual cost as the price of node is increased. * [Since the reduced cost of all outgoing (resp., incoming) * arcs will decrease (resp., increase) as the price of node * is increased, the next breakpoint is the minimum of the * positive reduced cost on outgoing arcs and of the * negative reduced cost on incoming arcs.] */ delprc = large; for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { trc = rc[arc]; if ((trc > 0) && (trc < delprc)) delprc = trc; } for (arc = fin[node]; arc > 0; arc = nxtin[arc]) { trc = rc[arc]; if ((trc < 0) && (trc > -delprc)) delprc = -trc; } /* If no breakpoint is left and dual ascent is still * possible, the problem is infeasible. */ if (delprc >= large) { if (ddpos[node] == 0) continue; return 4; } /* delprc is the stepsize to next breakpoint. Increase * price of node by delprc and compute the stepsize to the * next breakpoint in the dual cost. */ L53: nxtbrk = large; /* Look at all arcs out of node. */ for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { trc = rc[arc]; if (trc == 0) { t1 = endn[arc]; t = u[arc]; if (t > 0) { dfct[node] += t; dfct[t1] -= t; x[arc] = t; u[arc] = 0; } else t = x[arc]; ddneg[node] -= t; ddpos[t1] -= t; } /* Decrease the reduced costs on all outgoing arcs. */ trc -= delprc; if ((trc > 0) && (trc < nxtbrk)) nxtbrk = trc; else if (trc == 0) { /* Arc goes from inactive to balanced. Update the rate * of dual ascent at node and at its neighbor. */ ddpos[node] += u[arc]; ddneg[endn[arc]] += u[arc]; } rc[arc] = trc; } /* Look at all arcs into node. */ for (arc = fin[node]; arc > 0; arc = nxtin[arc]) { trc = rc[arc]; if (trc == 0) { t1 = startn[arc]; t = x[arc]; if (t > 0) { dfct[node] += t; dfct[t1] -= t; u[arc] = t; x[arc] = 0; } else t = u[arc]; ddpos[t1] -= t; ddneg[node] -= t; } /* Increase the reduced cost on all incoming arcs. */ trc += delprc; if ((trc < 0) && (trc > -nxtbrk)) nxtbrk = -trc; else if (trc == 0) { /* Arc goes from active to balanced. Update the rate * of dual ascent at node and at its neighbor. */ ddneg[startn[arc]] += x[arc]; ddpos[node] += x[arc]; } rc[arc] = trc; } /* If price of node can be increased further without * decreasing the dual cost (even the dual cost doesn't * increase), return to increase the price further. */ if ((ddpos[node] <= 0) && (nxtbrk < large)) { delprc = nxtbrk; goto L53; } } else if (ddneg[node] <= 0) { /* Compute delprc, the stepsize to the next breakpoint in * the dual cost as the price of node is decreased. * [Since the reduced cost of all outgoing (resp., incoming) * arcs will increase (resp., decrease) as the price of node * is decreased, the next breakpoint is the minimum of the * negative reduced cost on outgoing arcs and of the * positive reduced cost on incoming arcs.] */ delprc = large; for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { trc = rc[arc]; if ((trc < 0) && (trc > -delprc)) delprc = -trc; } for (arc = fin[node]; arc > 0; arc = nxtin[arc]) { trc = rc[arc]; if ((trc > 0) && (trc < delprc)) delprc = trc; } /* If no breakpoint is left and dual ascent is still * possible, the problem is infeasible. */ if (delprc == large) { if (ddneg[node] == 0) continue; return 5; } /* delprc is the stepsize to next breakpoint. Decrease * price of node by delprc and compute the stepsize to the * next breakpoint in the dual cost. */ L63: nxtbrk = large; /* Look at all arcs out of node. */ for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { trc = rc[arc]; if (trc == 0) { t1 = endn[arc]; t = x[arc]; if (t > 0) { dfct[node] -= t; dfct[t1] += t; u[arc] = t; x[arc] = 0; } else t = u[arc]; ddpos[node] -= t; ddneg[t1] -= t; } /* Increase the reduced cost on all outgoing arcs. */ trc += delprc; if ((trc < 0) && (trc > -nxtbrk)) nxtbrk = -trc; else if (trc == 0) { /* Arc goes from active to balanced. Update the rate * of dual ascent at node and at its neighbor. */ ddneg[node] += x[arc]; ddpos[endn[arc]] += x[arc]; } rc[arc] = trc; } /* Look at all arcs into node. */ for (arc = fin[node]; arc > 0; arc = nxtin[arc]) { trc = rc[arc]; if (trc == 0) { t1 = startn[arc]; t = u[arc]; if (t > 0) { dfct[node] -= t; dfct[t1] += t; x[arc] = t; u[arc] = 0; } else t = x[arc]; ddneg[t1] -= t; ddpos[node] -= t; } /* Decrease the reduced cost on all incoming arcs. */ trc -= delprc; if ((trc > 0) && (trc < nxtbrk)) nxtbrk = trc; else if (trc == 0) { /* Arc goes from inactive to balanced. Update the rate * of dual ascent at node and at its neighbor. */ ddpos[startn[arc]] += u[arc]; ddneg[node] += u[arc]; } rc[arc] = trc; } /* If price of node can be decreased further without * decreasing the dual cost (even the dual cost doesn't * increase), return to decrease the price further. */ if ((ddneg[node] <= 0) && (nxtbrk < large)) { delprc = nxtbrk; goto L63; } } } /*--------------------------------------------------------------*/ L70: /* Initialize tree data structure. */ for (i = 1; i <= n; i++) tfstou[i] = tfstin[i] = 0; for (i = 1; i <= na; i++) { tnxtin[i] = tnxtou[i] = -1; if (rc[i] == 0) { tnxtou[i] = tfstou[startn[i]]; tfstou[startn[i]] = i; tnxtin[i] = tfstin[endn[i]]; tfstin[endn[i]] = i; } } L90: /* Initialize other variables. */ feasbl = true; iter = 0; nmultinode = 0; num_augm = 0; num_ascnt = 0; num_passes = 0; numnz = n; numnz_new = 0; svitch = false; for (i = 1; i <= n; i++) mark[i] = scan[i] = false; nlabel = 0; /* RELAX4 uses an adaptive strategy to decide whether to continue * the scanning process after a multinode price change. * The threshold parameter tp and ts that control this strategy * are set in the next two lines. */ tp = 10; ts = n / 15; /* Initialize the queue of nodes with nonzero deficit. */ for (node = 1; node <= n - 1; node++) nxtqueue[node] = node + 1; nxtqueue[n] = 1; node = lastqueue = n; /*--------------------------------------------------------------*/ /* Start the relaxation algorithm. */ L100: /* Code for advancing the queue of nonzero deficit nodes. */ prevnode = node; node = nxtqueue[node]; defcit = dfct[node]; if (node == lastqueue) { numnz = numnz_new; numnz_new = 0; lastqueue = prevnode; num_passes++; } /* Code for deleting a node from the queue. */ if (defcit == 0) { nxtnode = nxtqueue[node]; if (node == nxtnode) return 0; else { nxtqueue[prevnode] = nxtnode; nxtqueue[node] = 0; node = nxtnode; goto L100; } } else posit = (defcit > 0); iter++; numnz_new++; if (posit) { /* Attempt a single node iteration from node with positive * deficit. */ pchange = false; indef = defcit; delx = 0; nb = 0; /* Check outgoing (probably) balanced arcs from node. */ for (arc = tfstou[node]; arc > 0; arc = tnxtou[arc]) { if ((rc[arc] == 0) && (x[arc] > 0)) { delx += x[arc]; nb++; save[nb] = arc; } } /* Check incoming arcs. */ for (arc = tfstin[node]; arc > 0; arc = tnxtin[arc]) { if ((rc[arc] == 0) && (u[arc] > 0)) { delx += u[arc]; nb++; save[nb] = -arc; } } /* End of initial node scan. */ L4018: /* If no price change is possible, exit. */ if (delx > defcit) { quit = (defcit < indef); goto L4016; } /* RELAX4 searches along the ascent direction for the best * price by checking the slope of the dual cost at successive * break points. First, we compute the distance to the next * break point. */ delprc = large; for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { rdcost = rc[arc]; if ((rdcost < 0) && (rdcost > -delprc)) delprc = -rdcost; } for (arc = fin[node]; arc > 0; arc = nxtin[arc]) { rdcost = rc[arc]; if ((rdcost > 0) && (rdcost < delprc)) delprc = rdcost; } /* Check if problem is infeasible. */ if ((delx < defcit) && (delprc == large)) { /* The dual cost can be decreased without bound. */ return 6; } /* Skip flow adjustment if there is no flow to modify. */ if (delx == 0) goto L4014; /* Adjust the flow on the balanced arcs incident to node to * maintain complementary slackness after the price change. */ for (j = 1; j <= nb; j++) { arc = save[j]; if (arc > 0) { node2 = endn[arc]; t1 = x[arc]; dfct[node2] += t1; if (nxtqueue[node2] == 0) { nxtqueue[prevnode] = node2; nxtqueue[node2] = node; prevnode = node2; } u[arc] += t1; x[arc] = 0; } else { narc = -arc; node2 = startn[narc]; t1 = u[narc]; dfct[node2] += t1; if (nxtqueue[node2] == 0) { nxtqueue[prevnode] = node2; nxtqueue[node2] = node; prevnode = node2; } x[narc] += t1; u[narc] = 0; } } defcit -= delx; L4014: if (delprc == large) { quit = true; goto L4019; } /* Node corresponds to a dual ascent direction. Decrease the * price of node by delprc and compute the stepsize to the next * breakpoint in the dual cost. */ nb = 0; pchange = true; dp = delprc; delprc = large; delx = 0; for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { rdcost = rc[arc] + dp; rc[arc] = rdcost; if (rdcost == 0) { nb++; save[nb] = arc; delx += x[arc]; } if ((rdcost < 0) && (rdcost > -delprc)) delprc = -rdcost; } for (arc = fin[node]; arc > 0; arc = nxtin[arc]) { rdcost = rc[arc] - dp; rc[arc] = rdcost; if (rdcost == 0) { nb++; save[nb] = -arc; delx += u[arc]; } if ((rdcost > 0) && (rdcost < delprc)) delprc = rdcost; } /* Return to check if another price change is possible. */ goto L4018; L4016: /* Perform flow augmentation at node. */ for (j = 1; j <= nb; j++) { arc = save[j]; if (arc > 0) { /* arc is an outgoing arc from node. */ node2 = endn[arc]; t1 = dfct[node2]; if (t1 < 0) { /* Decrease the total deficit by decreasing flow of * arc. */ quit = true; t2 = x[arc]; dx = defcit; if (dx > -t1) dx = -t1; if (dx > t2) dx = t2; defcit -= dx; dfct[node2] = t1 + dx; if (nxtqueue[node2] == 0) { nxtqueue[prevnode] = node2; nxtqueue[node2] = node; prevnode = node2; } x[arc] = t2 - dx; u[arc] += dx; if (defcit == 0) break; } } else { /* -arc is an incoming arc to node. */ narc = -arc; node2 = startn[narc]; t1 = dfct[node2]; if (t1 < 0) { /* Decrease the total deficit by increasing flow of * -arc. */ quit = true; t2 = u[narc]; dx = defcit; if (dx > -t1) dx = -t1; if (dx > t2) dx = t2; defcit -= dx; dfct[node2] = t1 + dx; if (nxtqueue[node2] == 0) { nxtqueue[prevnode] = node2; nxtqueue[node2] = node; prevnode = node2; } x[narc] += dx; u[narc] = t2 - dx; if (defcit == 0) break; } } } L4019: dfct[node] = defcit; /* Reconstruct the linked list of balance arcs incident to this * node. For each adjacent node, we add any newly balanced arcs * to the list, but do not bother removing formerly balanced * ones (they will be removed the next time each adjacent node * is scanned). */ if (pchange) { arc = tfstou[node]; tfstou[node] = 0; while (arc > 0) { nxtarc = tnxtou[arc]; tnxtou[arc] = -1; arc = nxtarc; } arc = tfstin[node]; tfstin[node] = 0; while (arc > 0) { nxtarc = tnxtin[arc]; tnxtin[arc] = -1; arc = nxtarc; } /* Now add the currently balanced arcs to the list for this * node (which is now empty), and the appropriate adjacent * ones. */ for (j = 1; j <= nb; j++) { arc = save[j]; if (arc < 0) arc = -arc; if (tnxtou[arc] < 0) { tnxtou[arc] = tfstou[startn[arc]]; tfstou[startn[arc]] = arc; } if (tnxtin[arc] < 0) { tnxtin[arc] = tfstin[endn[arc]]; tfstin[endn[arc]] = arc; } } } /* End of single node iteration for positive deficit node. */ } else { /* Attempt a single node iteration from node with negative * deficit. */ pchange = false; defcit = -defcit; indef = defcit; delx = 0; nb = 0; for (arc = tfstin[node]; arc > 0; arc = tnxtin[arc]) { if ((rc[arc] == 0) && (x[arc] > 0)) { delx += x[arc]; nb++; save[nb] = arc; } } for (arc = tfstou[node]; arc > 0; arc = tnxtou[arc]) { if ((rc[arc] == 0) && (u[arc] > 0)) { delx += u[arc]; nb++; save[nb] = -arc; } } L4028: if (delx >= defcit) { quit = (defcit < indef); goto L4026; } /* Compute distance to next breakpoint. */ delprc = large; for (arc = fin[node]; arc > 0; arc = nxtin[arc]) { rdcost = rc[arc]; if ((rdcost < 0) && (rdcost > -delprc)) delprc = -rdcost; } for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { rdcost = rc[arc]; if ((rdcost > 0) && (rdcost < delprc)) delprc = rdcost; } /* Check if problem is infeasible. */ if ((delx < defcit) && (delprc == large)) return 7; if (delx == 0) goto L4024; /* Flow augmentation is possible. */ for (j = 1; j <= nb; j++) { arc = save[j]; if (arc > 0) { node2 = startn[arc]; t1 = x[arc]; dfct[node2] -= t1; if (nxtqueue[node2] == 0) { nxtqueue[prevnode] = node2; nxtqueue[node2] = node; prevnode = node2; } u[arc] += t1; x[arc] = 0; } else { narc = -arc; node2 = endn[narc]; t1 = u[narc]; dfct[node2] -= t1; if (nxtqueue[node2] == 0) { nxtqueue[prevnode] = node2; nxtqueue[node2] = node; prevnode = node2; } x[narc] += t1; u[narc] = 0; } } defcit -= delx; L4024: if (delprc == large) { quit = true; goto L4029; } /* Price increase at node is possible. */ nb = 0; pchange = true; dp = delprc; delprc = large; delx = 0; for (arc = fin[node]; arc > 0; arc = nxtin[arc]) { rdcost = rc[arc] + dp; rc[arc] = rdcost; if (rdcost == 0) { nb++; save[nb] = arc; delx += x[arc]; } if ((rdcost < 0) && (rdcost > -delprc)) delprc = -rdcost; } for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { rdcost = rc[arc] - dp; rc[arc] = rdcost; if (rdcost == 0) { nb++; save[nb] = -arc; delx += u[arc]; } if ((rdcost > 0) && (rdcost < delprc)) delprc = rdcost; } goto L4028; L4026: /* Perform flow augmentation at node. */ for (j = 1; j <= nb; j++) { arc = save[j]; if (arc > 0) { /* arc is an incoming arc to node. */ node2 = startn[arc]; t1 = dfct[node2]; if (t1 > 0) { quit = true; t2 = x[arc]; dx = defcit; if (dx > t1) dx = t1; if (dx > t2) dx = t2; defcit -= dx; dfct[node2] = t1 - dx; if (nxtqueue[node2] == 0) { nxtqueue[prevnode] = node2; nxtqueue[node2] = node; prevnode = node2; } x[arc] = t2 - dx; u[arc] += dx; if (defcit == 0) break; } } else { /* -arc is an outgoing arc from node. */ narc = -arc; node2 = endn[narc]; t1 = dfct[node2]; if (t1 > 0) { quit = true; t2 = u[narc]; dx = defcit; if (dx > t1) dx = t1; if (dx > t2) dx = t2; defcit -= dx; dfct[node2] = t1 - dx; if (nxtqueue[node2] == 0) { nxtqueue[prevnode] = node2; nxtqueue[node2] = node; prevnode = node2; } x[narc] += dx; u[narc] = t2 - dx; if (defcit == 0) break; } } } L4029: dfct[node] = -defcit; /* Reconstruct the list of balanced arcs incident to node. */ if (pchange) { arc = tfstou[node]; tfstou[node] = 0; while (arc > 0) { nxtarc = tnxtou[arc]; tnxtou[arc] = -1; arc = nxtarc; } arc = tfstin[node]; tfstin[node] = 0; while (arc > 0) { nxtarc = tnxtin[arc]; tnxtin[arc] = -1; arc = nxtarc; } /* Now add the currently balanced arcs to the list for this * node (which is now empty), and the appropriate adjacent * ones. */ for (j = 1; j <= nb; j++) { arc = save[j]; if (arc <= 0) arc = -arc; if (tnxtou[arc] < 0) { tnxtou[arc] = tfstou[startn[arc]]; tfstou[startn[arc]] = arc; } if (tnxtin[arc] < 0) { tnxtin[arc] = tfstin[endn[arc]]; tfstin[endn[arc]] = arc; } } } /* End of single node iteration for a negative deficit node. */ } if (quit || (num_passes <= 3)) goto L100; /* Do a multinode iteration from node. */ nmultinode++; /* If number of nonzero deficit nodes is small, continue labeling * until a flow augmentation is done. */ svitch = (numnz < tp); /* Unmark nodes labeled earlier. */ for (j = 1; j <= nlabel; j++) { node2 = label[j]; mark[node2] = scan[node2] = false; } /* Initialize labeling. */ nlabel = 1; label[1] = node; mark[node] = true; prdcsr[node] = 0; /* Scan starting node. */ scan[node] = true; nscan = 1; dm = dfct[node]; delx = 0; for (j = 1; j <= nb; j++) { arc = save[j]; if (arc > 0) { if (posit) node2 = endn[arc]; else node2 = startn[arc]; if (!mark[node2]) { nlabel++; label[nlabel] = node2; prdcsr[node2] = arc; mark[node2] = true; delx += x[arc]; } } else { narc = -arc; if (posit) node2 = startn[narc]; else node2 = endn[narc]; if (!mark[node2]) { nlabel++; label[nlabel] = node2; prdcsr[node2] = arc; mark[node2] = true; delx += u[narc]; } } } L4120:/* Start scanning a labeled but unscanned node. */ nscan++; /* Check to see if switch needs to be set to true so to continue * scanning even after a price change. */ svitch = svitch || ((nscan > ts) && (numnz < ts)); /* Scanning will continue until either an overestimate of the * residual capacity across the cut corresponding to the scanned * set of nodes (called delx) exceeds the absolute value of the * total deficit of the scanned nodes (called dm), or else an * augmenting path is found. Arcs that are in the tree but are not * balanced are removed as part of the scanning process. */ i = label[nscan]; scan[i] = true; naugnod = 0; if (posit) { /* Scanning node i in case of positive deficit. */ prvarc = 0; arc = tfstou[i]; while (arc > 0) { /* arc is an outgoing arc from node. */ if (rc[arc] == 0) { if (x[arc] > 0) { node2 = endn[arc]; if (!mark[node2]) { /* node2 is not labeled, so add node2 to the labeled set. */ prdcsr[node2] = arc; if (dfct[node2] < 0) { naugnod++; save[naugnod] = node2; } nlabel++; label[nlabel] = node2; mark[node2] = true; delx += x[arc]; } } prvarc = arc; arc = tnxtou[arc]; } else { tmparc = arc; arc = tnxtou[arc]; tnxtou[tmparc] = -1; if (prvarc == 0) tfstou[i] = arc; else tnxtou[prvarc] = arc; } } prvarc = 0; arc = tfstin[i]; while (arc > 0) { /* arc is an incoming arc into node. */ if (rc[arc] == 0) { if (u[arc] > 0) { node2 = startn[arc]; if (!mark[node2]) { /* node2 is not labeled, so add node2 to the * labeled set. */ prdcsr[node2] = -arc; if (dfct[node2] < 0) { naugnod++; save[naugnod] = node2; } nlabel++; label[nlabel] = node2; mark[node2] = true; delx += u[arc]; } } prvarc = arc; arc = tnxtin[arc]; } else { tmparc = arc; arc = tnxtin[arc]; tnxtin[tmparc] = -1; if (prvarc == 0) tfstin[i] = arc; else tnxtin[prvarc] = arc; } } /* Correct the residual capacity of the scanned node cut. */ arc = prdcsr[i]; if (arc > 0) delx -= x[arc]; else delx -= u[-arc]; /* End of scanning of node i for positive deficit case. */ } else { /* Scanning node i for negative deficit case. */ prvarc = 0; arc = tfstin[i]; while (arc > 0) { if (rc[arc] == 0) { if (x[arc] > 0) { node2 = startn[arc]; if (!mark[node2]) { prdcsr[node2] = arc; if (dfct[node2] > 0) { naugnod++; save[naugnod] = node2; } nlabel++; label[nlabel] = node2; mark[node2] = true; delx += x[arc]; } } prvarc = arc; arc = tnxtin[arc]; } else { tmparc = arc; arc = tnxtin[arc]; tnxtin[tmparc] = -1; if (prvarc == 0) tfstin[i] = arc; else tnxtin[prvarc] = arc; } } prvarc = 0; arc = tfstou[i]; while (arc > 0) { if (rc[arc] == 0) { if (u[arc] > 0) { node2 = endn[arc]; if (!mark[node2]) { prdcsr[node2] = -arc; if (dfct[node2] > 0) { naugnod++; save[naugnod] = node2; } nlabel++; label[nlabel] = node2; mark[node2] = true; delx += u[arc]; } } prvarc = arc; arc = tnxtou[arc]; } else { tmparc = arc; arc = tnxtou[arc]; tnxtou[tmparc] = -1; if (prvarc == 0) tfstou[i] = arc; else tnxtou[prvarc] = arc; } } arc = prdcsr[i]; if (arc > 0) delx -= x[arc]; else delx -= u[-arc]; } /* Add deficit of node scanned to dm. */ dm += dfct[i]; /* Check if the set of scanned nodes correspond to a dual ascent * direction; if yes, perform a price adjustment step, otherwise * continue labeling. */ if (nscan < nlabel) { if (svitch) goto L4210; if ((delx >= dm) && (delx >= -dm)) goto L4210; } /* Try a price change. * [Note that since delx - abs(dm) is an overestimate of ascent * slope, we may occasionally try a direction that is not an * ascent direction. In this case the ascnt routines return with * quit = false, so we continue labeling nodes.] */ if (posit) { ascnt1(csa, dm, &delx, &nlabel, &feasbl, &svitch, nscan, node, &prevnode); num_ascnt++; } else { ascnt2(csa, dm, &delx, &nlabel, &feasbl, &svitch, nscan, node, &prevnode); num_ascnt++; } if (!feasbl) return 8; if (!svitch) goto L100; /* Store those newly labeled nodes to which flow augmentation is * possible. */ naugnod = 0; for (j = nscan + 1; j <= nlabel; j++) { node2 = label[j]; if (posit && (dfct[node2] < 0)) { naugnod++; save[naugnod] = node2; } else if ((!posit) && (dfct[node2] > 0)) { naugnod++; save[naugnod] = node2; } } L4210:/* Check if flow augmentation is possible. If not, return to scan * another node. */ if (naugnod == 0) goto L4120; for (j = 1; j <= naugnod; j++) { num_augm++; augnod = save[j]; if (posit) { /* Do the augmentation from node with positive deficit. */ dx = -dfct[augnod]; ib = augnod; while (ib != node) { arc = prdcsr[ib]; if (arc > 0) { if (dx > x[arc]) dx = x[arc]; ib = startn[arc]; } else { if (dx > u[-arc]) dx = u[-arc]; ib = endn[-arc]; } } if (dx > dfct[node]) dx = dfct[node]; if (dx > 0) { /* Increase (decrease) the flow of all forward (backward) * arcs in the flow augmenting path. Adjust node deficit * accordingly. */ if (nxtqueue[augnod] == 0) { nxtqueue[prevnode] = augnod; nxtqueue[augnod] = node; prevnode = augnod; } dfct[augnod] += dx; dfct[node] -= dx; ib = augnod; while (ib != node) { arc = prdcsr[ib]; if (arc > 0) { x[arc] -= dx; u[arc] += dx; ib = startn[arc]; } else { narc = -arc; x[narc] += dx; u[narc] -= dx; ib = endn[narc]; } } } } else { /* Do the augmentation from node with negative deficit. */ dx = dfct[augnod]; ib = augnod; while (ib != node) { arc = prdcsr[ib]; if (arc > 0) { if (dx > x[arc]) dx = x[arc]; ib = endn[arc]; } else { if (dx > u[-arc]) dx = u[-arc]; ib = startn[-arc]; } } if (dx > -dfct[node]) dx = -dfct[node]; if (dx > 0) { /* Update the flow and deficits. */ if (nxtqueue[augnod] == 0) { nxtqueue[prevnode] = augnod; nxtqueue[augnod] = node; prevnode = augnod; } dfct[augnod] -= dx; dfct[node] += dx; ib = augnod; while (ib != node) { arc = prdcsr[ib]; if (arc > 0) { x[arc] -= dx; u[arc] += dx; ib = endn[arc]; } else { narc = -arc; x[narc] += dx; u[narc] -= dx; ib = startn[narc]; } } } } if (dfct[node] == 0) goto L100; if (dfct[augnod] != 0) svitch = false; } /* If node still has nonzero deficit and all newly labeled nodes * have same sign for their deficit as node, we can continue * labeling. In this case, continue labeling only when flow * augmentation is done relatively infrequently. */ if (svitch && (iter > 8 * num_augm)) goto L4120; /* Return to do another relaxation iteration. */ goto L100; # undef nmultinode # undef iter # undef num_augm # undef num_ascnt # undef nsp } /*********************************************************************** * NAME * * relax4_inidat - construct linked lists for network topology * * PURPOSE * * This routine constructs two linked lists for the network topology: * one list (given by fou, nxtou) for the outgoing arcs of nodes and * one list (given by fin, nxtin) for the incoming arcs of nodes. These * two lists are required by RELAX4. * * INPUT PARAMETERS * * n = number of nodes * na = number of arcs * startn[j] = starting node for arc j, j = 1,...,na * endn[j] = ending node for arc j, j = 1,...,na * * OUTPUT PARAMETERS * * fou[i] = first arc out of node i, i = 1,...,n * nxtou[j] = next arc out of the starting node of arc j, j = 1,...,na * fin[i] = first arc into node i, i = 1,...,n * nxtin[j] = next arc into the ending node of arc j, j = 1,...,na * * WORKING PARAMETERS * * tempin[1+n], tempou[1+n] */ void relax4_inidat(struct relax4_csa *csa) { /* input parameters */ int n = csa->n; int na = csa->na; int *startn = csa->startn; int *endn = csa->endn; /* output parameters */ int *fou = csa->fou; int *nxtou = csa->nxtou; int *fin = csa->fin; int *nxtin = csa->nxtin; /* working parameters */ int *tempin = csa->label; int *tempou = csa->prdcsr; /* local variables */ int i, i1, i2; for (i = 1; i <= n; i++) { fin[i] = fou[i] = 0; tempin[i] = tempou[i] = 0; } for (i = 1; i <= na; i++) { nxtin[i] = nxtou[i] = 0; i1 = startn[i]; i2 = endn[i]; if (fou[i1] != 0) nxtou[tempou[i1]] = i; else fou[i1] = i; tempou[i1] = i; if (fin[i2] != 0) nxtin[tempin[i2]] = i; else fin[i2] = i; tempin[i2] = i; } return; } /*********************************************************************** * NAME * * ascnt1 - multi-node price adjustment for positive deficit case * * PURPOSE * * This subroutine performs the multi-node price adjustment step for * the case where the scanned nodes have positive deficit. It first * checks if decreasing the price of the scanned nodes increases the * dual cost. If yes, then it decreases the price of all scanned nodes. * There are two possibilities for price decrease: if switch = true, * then the set of scanned nodes corresponds to an elementary direction * of maximal rate of ascent, in which case the price of all scanned * nodes are decreased until the next breakpoint in the dual cost is * encountered. At this point, some arc becomes balanced and more * node(s) are added to the labeled set and the subroutine is exited. * If switch = false, then the price of all scanned nodes are decreased * until the rate of ascent becomes negative (this corresponds to the * price adjustment step in which both the line search and the * degenerate ascent iteration are implemented). * * INPUT PARAMETERS * * dm = total deficit of scanned nodes * switch = true if labeling is to continue after price change * nscan = number of scanned nodes * curnode = most recently scanned node * n = number of nodes * na = number of arcs * large = a very large integer to represent infinity (see note 3) * startn[i] = starting node for the i-th arc, i = 1,...,na * endn[i] = ending node for the i-th arc, i = 1,...,na * fou[i] = first arc leaving i-th node, i = 1,...,n * nxtou[i] = next arc leaving the starting node of j-th arc, * i = 1,...,na * fin[i] = first arc entering i-th node, i = 1,...,n * nxtin[i] = next arc entering the ending node of j-th arc, * i = 1,...,na * * UPDATED PARAMETERS * * delx = a lower estimate of the total flow on balanced arcs in * the scanned-nodes cut * nlabel = number of labeled nodes * feasbl = false if problem is found to be infeasible * prevnode = the node before curnode in queue * rc[j] = reduced cost of arc j, j = 1,...,na * u[j] = residual capacity of arc j, j = 1,...,na * x[j] = flow on arc j, j = 1,...,na * dfct[i] = deficit at node i, i = 1,...,n * label[k] = k-th node labeled, k = 1,...,nlabel * prdcsr[i] = predecessor of node i in tree of labeled nodes (0 if i * is unlabeled), i = 1,...,n * tfstou[i] = first balanced arc out of node i, i = 1,...,n * tnxtou[j] = next balanced arc out of the starting node of arc j, * j = 1,...,na * tfstin[i] = first balanced arc into node i, i = 1,...,n * tnxtin[j] = next balanced arc into the ending node of arc j, * j = 1,...,na * nxtqueue[i] = node following node i in the fifo queue (0 if node is * not in the queue), i = 1,...,n * scan[i] = true if node i is scanned, i = 1,...,n * mark[i] = true if node i is labeled, i = 1,...,n * * WORKING PARAMETERS * * save[1+na] */ static void ascnt1(struct relax4_csa *csa, int dm, int *delx, int *nlabel, int *feasbl, int *svitch, int nscan, int curnode, int *prevnode) { /* input parameters */ int n = csa->n; /* int na = csa->na; */ int large = csa->large; int *startn = csa->startn; int *endn = csa->endn; int *fou = csa->fou; int *nxtou = csa->nxtou; int *fin = csa->fin; int *nxtin = csa->nxtin; /* updated parameters */ # define delx (*delx) # define nlabel (*nlabel) # define feasbl (*feasbl) # define svitch (*svitch) # define prevnode (*prevnode) int *rc = csa->rc; int *u = csa->u; int *x = csa->x; int *dfct = csa->dfct; int *label = csa->label; int *prdcsr = csa->prdcsr; int *tfstou = csa->tfstou; int *tnxtou = csa->tnxtou; int *tfstin = csa->tfstin; int *tnxtin = csa->tnxtin; int *nxtqueue = csa->nxtqueue; char *scan = csa->scan; char *mark = csa->mark; int *save = csa->save; /* local variables */ int arc, delprc, dlx, i, j, nb, node, node2, nsave, rdcost, t1, t2, t3; /* Store the arcs between the set of scanned nodes and its * complement in save and compute delprc, the stepsize to the next * breakpoint in the dual cost in the direction of decreasing * prices of the scanned nodes. * [The arcs are stored into save by looking at the arcs incident * to either the set of scanned nodes or its complement, depending * on whether nscan > n/2 or not. This improves the efficiency of * storing.] */ delprc = large; dlx = 0; nsave = 0; if (nscan <= n / 2) { for (i = 1; i <= nscan; i++) { node = label[i]; for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { /* arc points from scanned node to an unscanned node. */ node2 = endn[arc]; if (!scan[node2]) { nsave++; save[nsave] = arc; rdcost = rc[arc]; if ((rdcost == 0) && (prdcsr[node2] != arc)) dlx += x[arc]; if ((rdcost < 0) && (rdcost > -delprc)) delprc = -rdcost; } } for (arc = fin[node]; arc > 0; arc = nxtin[arc]) { /* arc points from unscanned node to scanned node. */ node2 = startn[arc]; if (!scan[node2]) { nsave++; save[nsave] = -arc; rdcost = rc[arc]; if ((rdcost == 0) && (prdcsr[node2] != -arc)) dlx += u[arc]; if ((rdcost > 0) && (rdcost < delprc)) delprc = rdcost; } } } } else { for (node = 1; node <= n; node++) { if (scan[node]) continue; for (arc = fin[node]; arc > 0; arc = nxtin[arc]) { node2 = startn[arc]; if (scan[node2]) { nsave++; save[nsave] = arc; rdcost = rc[arc]; if ((rdcost == 0) && (prdcsr[node] != arc)) dlx += x[arc]; if ((rdcost < 0) && (rdcost > -delprc)) delprc = -rdcost; } } for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { node2 = endn[arc]; if (scan[node2]) { nsave++; save[nsave] = -arc; rdcost = rc[arc]; if ((rdcost == 0) && (prdcsr[node] != -arc)) dlx += u[arc]; if ((rdcost > 0) && (rdcost < delprc)) delprc = rdcost; } } } } /* Check if the set of scanned nodes truly corresponds to a dual * ascent direction. [Here delx + dlx is the exact sum of the flow * on arcs from the scanned set to the unscanned set plus the * (capacity - flow) on arcs from the unscanned set to the scanned * set.] If this were not the case, set switch to true and exit * subroutine. */ if ((!svitch) && (delx + dlx >= dm)) { svitch = true; return; } delx += dlx; L4: /* Check that the problem is feasible. */ if (delprc == large) { /* We can increase the dual cost without bound, so the primal * problem is infeasible. */ feasbl = false; return; } /* Decrease the prices of the scanned nodes, add more nodes to * the labeled set and check if a newly labeled node has negative * deficit. */ if (svitch) { for (i = 1; i <= nsave; i++) { arc = save[i]; if (arc > 0) { rc[arc] += delprc; if (rc[arc] == 0) { node2 = endn[arc]; if (tnxtou[arc] < 0) { tnxtou[arc] = tfstou[startn[arc]]; tfstou[startn[arc]] = arc; } if (tnxtin[arc] < 0) { tnxtin[arc] = tfstin[node2]; tfstin[node2] = arc; } if (!mark[node2]) { prdcsr[node2] = arc; nlabel++; label[nlabel] = node2; mark[node2] = true; } } } else { arc = -arc; rc[arc] -= delprc; if (rc[arc] == 0) { node2 = startn[arc]; if (tnxtou[arc] < 0) { tnxtou[arc] = tfstou[node2]; tfstou[node2] = arc; } if (tnxtin[arc] < 0) { tnxtin[arc] = tfstin[endn[arc]]; tfstin[endn[arc]] = arc; } if (!mark[node2]) { prdcsr[node2] = -arc; nlabel++; label[nlabel] = node2; mark[node2] = true; } } } } return; } else { /* Decrease the prices of the scanned nodes by delprc. Adjust * flow to maintain complementary slackness with the prices. */ nb = 0; for (i = 1; i <= nsave; i++) { arc = save[i]; if (arc > 0) { t1 = rc[arc]; if (t1 == 0) { t2 = x[arc]; t3 = startn[arc]; dfct[t3] -= t2; if (nxtqueue[t3] == 0) { nxtqueue[prevnode] = t3; nxtqueue[t3] = curnode; prevnode = t3; } t3 = endn[arc]; dfct[t3] += t2; if (nxtqueue[t3] == 0) { nxtqueue[prevnode] = t3; nxtqueue[t3] = curnode; prevnode = t3; } u[arc] += t2; x[arc] = 0; } rc[arc] = t1 + delprc; #if 0 /* by mao; 26/IV-2013 */ if (rc[arc] == 0) #else if (rc[arc] == 0 && nb < n) #endif { delx += x[arc]; nb++; prdcsr[nb] = arc; } } else { arc = -arc; t1 = rc[arc]; if (t1 == 0) { t2 = u[arc]; t3 = startn[arc]; dfct[t3] += t2; if (nxtqueue[t3] == 0) { nxtqueue[prevnode] = t3; nxtqueue[t3] = curnode; prevnode = t3; } t3 = endn[arc]; dfct[t3] -= t2; if (nxtqueue[t3] == 0) { nxtqueue[prevnode] = t3; nxtqueue[t3] = curnode; prevnode = t3; } x[arc] += t2; u[arc] = 0; } rc[arc] = t1 - delprc; #if 0 /* by mao; 26/IV-2013 */ if (rc[arc] == 0) #else if (rc[arc] == 0 && nb < n) #endif { delx += u[arc]; nb++; prdcsr[nb] = arc; } } } } if (delx <= dm) { /* The set of scanned nodes still corresponds to a dual * (possibly degenerate) ascent direction. Compute the stepsize * delprc to the next breakpoint in the dual cost. */ delprc = large; for (i = 1; i <= nsave; i++) { arc = save[i]; if (arc > 0) { rdcost = rc[arc]; if ((rdcost < 0) && (rdcost > -delprc)) delprc = -rdcost; } else { arc = -arc; rdcost = rc[arc]; if ((rdcost > 0) && (rdcost < delprc)) delprc = rdcost; } } if ((delprc != large) || (delx < dm)) goto L4; } /* Add new balanced arcs to the superset of balanced arcs. */ for (i = 1; i <= nb; i++) { arc = prdcsr[i]; if (tnxtin[arc] == -1) { j = endn[arc]; tnxtin[arc] = tfstin[j]; tfstin[j] = arc; } if (tnxtou[arc] == -1) { j = startn[arc]; tnxtou[arc] = tfstou[j]; tfstou[j] = arc; } } return; # undef delx # undef nlabel # undef feasbl # undef svitch # undef prevnode } /*********************************************************************** * NAME * * ascnt2 - multi-node price adjustment for negative deficit case * * PURPOSE * * This routine is analogous to ascnt1 but for the case where the * scanned nodes have negative deficit. */ static void ascnt2(struct relax4_csa *csa, int dm, int *delx, int *nlabel, int *feasbl, int *svitch, int nscan, int curnode, int *prevnode) { /* input parameters */ int n = csa->n; /* int na = csa->na; */ int large = csa->large; int *startn = csa->startn; int *endn = csa->endn; int *fou = csa->fou; int *nxtou = csa->nxtou; int *fin = csa->fin; int *nxtin = csa->nxtin; /* updated parameters */ # define delx (*delx) # define nlabel (*nlabel) # define feasbl (*feasbl) # define svitch (*svitch) # define prevnode (*prevnode) int *rc = csa->rc; int *u = csa->u; int *x = csa->x; int *dfct = csa->dfct; int *label = csa->label; int *prdcsr = csa->prdcsr; int *tfstou = csa->tfstou; int *tnxtou = csa->tnxtou; int *tfstin = csa->tfstin; int *tnxtin = csa->tnxtin; int *nxtqueue = csa->nxtqueue; char *scan = csa->scan; char *mark = csa->mark; int *save = csa->save; /* local variables */ int arc, delprc, dlx, i, j, nb, node, node2, nsave, rdcost, t1, t2, t3; /* Store the arcs between the set of scanned nodes and its * complement in save and compute delprc, the stepsize to the next * breakpoint in the dual cost in the direction of increasing * prices of the scanned nodes. */ delprc = large; dlx = 0; nsave = 0; if (nscan <= n / 2) { for (i = 1; i <= nscan; i++) { node = label[i]; for (arc = fin[node]; arc > 0; arc = nxtin[arc]) { node2 = startn[arc]; if (!scan[node2]) { nsave++; save[nsave] = arc; rdcost = rc[arc]; if ((rdcost == 0) && (prdcsr[node2] != arc)) dlx += x[arc]; if ((rdcost < 0) && (rdcost > -delprc)) delprc = -rdcost; } } for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { node2 = endn[arc]; if (!scan[node2]) { nsave++; save[nsave] = -arc; rdcost = rc[arc]; if ((rdcost == 0) && (prdcsr[node2] != -arc)) dlx += u[arc]; if ((rdcost > 0) && (rdcost < delprc)) delprc = rdcost; } } } } else { for (node = 1; node <= n; node++) { if (scan[node]) continue; for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { node2 = endn[arc]; if (scan[node2]) { nsave++; save[nsave] = arc; rdcost = rc[arc]; if ((rdcost == 0) && (prdcsr[node] != arc)) dlx += x[arc]; if ((rdcost < 0) && (rdcost > -delprc)) delprc = -rdcost; } } for (arc = fin[node]; arc > 0; arc = nxtin[arc]) { node2 = startn[arc]; if (scan[node2]) { nsave++; save[nsave] = -arc; rdcost = rc[arc]; if ((rdcost == 0) && (prdcsr[node] != -arc)) dlx += u[arc]; if ((rdcost > 0) && (rdcost < delprc)) delprc = rdcost; } } } } if ((!svitch) && (delx + dlx >= -dm)) { svitch = true; return; } delx += dlx; /* Check that the problem is feasible. */ L4: if (delprc == large) { feasbl = false; return; } /* Increase the prices of the scanned nodes, add more nodes to * the labeled set and check if a newly labeled node has positive * deficit. */ if (svitch) { for (i = 1; i <= nsave; i++) { arc = save[i]; if (arc > 0) { rc[arc] += delprc; if (rc[arc] == 0) { node2 = startn[arc]; if (tnxtou[arc] < 0) { tnxtou[arc] = tfstou[node2]; tfstou[node2] = arc; } if (tnxtin[arc] < 0) { tnxtin[arc] = tfstin[endn[arc]]; tfstin[endn[arc]] = arc; } if (!mark[node2]) { prdcsr[node2] = arc; nlabel++; label[nlabel] = node2; mark[node2] = true; } } } else { arc = -arc; rc[arc] -= delprc; if (rc[arc] == 0) { node2 = endn[arc]; if (tnxtou[arc] < 0) { tnxtou[arc] = tfstou[startn[arc]]; tfstou[startn[arc]] = arc; } if (tnxtin[arc] < 0) { tnxtin[arc] = tfstin[node2]; tfstin[node2] = arc; } if (!mark[node2]) { prdcsr[node2] = -arc; nlabel++; label[nlabel] = node2; mark[node2] = true; } } } } return; } else { nb = 0; for (i = 1; i <= nsave; i++) { arc = save[i]; if (arc > 0) { t1 = rc[arc]; if (t1 == 0) { t2 = x[arc]; t3 = startn[arc]; dfct[t3] -= t2; if (nxtqueue[t3] == 0) { nxtqueue[prevnode] = t3; nxtqueue[t3] = curnode; prevnode = t3; } t3 = endn[arc]; dfct[t3] += t2; if (nxtqueue[t3] == 0) { nxtqueue[prevnode] = t3; nxtqueue[t3] = curnode; prevnode = t3; } u[arc] += t2; x[arc] = 0; } rc[arc] = t1 + delprc; #if 0 /* by mao; 26/IV-2013 */ if (rc[arc] == 0) #else if (rc[arc] == 0 && nb < n) #endif { delx += x[arc]; nb++; prdcsr[nb] = arc; } } else { arc = -arc; t1 = rc[arc]; if (t1 == 0) { t2 = u[arc]; t3 = startn[arc]; dfct[t3] += t2; if (nxtqueue[t3] == 0) { nxtqueue[prevnode] = t3; nxtqueue[t3] = curnode; prevnode = t3; } t3 = endn[arc]; dfct[t3] -= t2; if (nxtqueue[t3] == 0) { nxtqueue[prevnode] = t3; nxtqueue[t3] = curnode; prevnode = t3; } x[arc] += t2; u[arc] = 0; } rc[arc] = t1 - delprc; #if 0 /* by mao; 26/IV-2013 */ if (rc[arc] == 0) #else if (rc[arc] == 0 && nb < n) #endif { delx += u[arc]; nb++; prdcsr[nb] = arc; } } } } if (delx <= -dm) { delprc = large; for (i = 1; i <= nsave; i++) { arc = save[i]; if (arc > 0) { rdcost = rc[arc]; if ((rdcost < 0) && (rdcost > -delprc)) delprc = -rdcost; } else { arc = -arc; rdcost = rc[arc]; if ((rdcost > 0) && (rdcost < delprc)) delprc = rdcost; } } if ((delprc != large) || (delx < -dm)) goto L4; } /* Add new balanced arcs to the superset of balanced arcs. */ for (i = 1; i <= nb; i++) { arc = prdcsr[i]; if (tnxtin[arc] == -1) { j = endn[arc]; tnxtin[arc] = tfstin[j]; tfstin[j] = arc; } if (tnxtou[arc] == -1) { j = startn[arc]; tnxtou[arc] = tfstou[j]; tfstou[j] = arc; } } return; # undef delx # undef nlabel # undef feasbl # undef svitch # undef prevnode } /*********************************************************************** * NAME * * auction - compute good initial flow and prices * * PURPOSE * * This subroutine uses a version of the auction algorithm for min * cost network flow to compute a good initial flow and prices for the * problem. * * INPUT PARAMETERS * * n = number of nodes * na = number of arcs * large = a very large integer to represent infinity (see note 3) * startn[i] = starting node for the i-th arc, i = 1,...,na * endn[i] = ending node for the i-th arc, i = 1,...,na * fou[i] = first arc leaving i-th node, i = 1,...,n * nxtou[i] = next arc leaving the starting node of j-th arc, * i = 1,...,na * fin[i] = first arc entering i-th node, i = 1,...,n * nxtin[i] = next arc entering the ending node of j-th arc, * i = 1,...,na * * UPDATED PARAMETERS * * rc[j] = reduced cost of arc j, j = 1,...,na * u[j] = residual capacity of arc j, j = 1,...,na * x[j] = flow on arc j, j = 1,...,na * dfct[i] = deficit at node i, i = 1,...,n * * OUTPUT PARAMETERS * * nsp = number of auction/shortest path iterations * * WORKING PARAMETERS * * p[1+n], prdcsr[1+n], save[1+na], fpushf[1+n], nxtpushf[1+na], * fpushb[1+n], nxtpushb[1+na], nxtqueue[1+n], extend_arc[1+n], * sb_level[1+n], sb_arc[1+n], path_id[1+n] * * RETURNS * * 0 = normal return * 1 = problem is found to be infeasible */ static int auction(struct relax4_csa *csa) { /* input parameters */ int n = csa->n; int na = csa->na; int large = csa->large; int *startn = csa->startn; int *endn = csa->endn; int *fou = csa->fou; int *nxtou = csa->nxtou; int *fin = csa->fin; int *nxtin = csa->nxtin; /* updated parameters */ # define crash (csa->crash) int *rc = csa->rc; int *u = csa->u; int *x = csa->x; int *dfct = csa->dfct; /* output parameters */ # define nsp (csa->nsp) /* working parameters */ int *p = csa->label; int *prdcsr = csa->prdcsr; int *save = csa->save; int *fpushf = csa->tfstou; int *nxtpushf = csa->tnxtou; int *fpushb = csa->tfstin; int *nxtpushb = csa->tnxtin; int *nxtqueue = csa->nxtqueue; int *extend_arc = csa->extend_arc; int *sb_level = csa->sb_level; int *sb_arc = csa->sb_arc; char *path_id = csa->mark; /* local variables */ int arc, bstlevel, end, eps, extarc, factor, flow, i, incr, last, lastqueue, maxcost, mincost, nas, naug, new_level, node, nolist, num_passes, nxtnode, pass, pend, pr_term, prd, prevarc, prevlevel, prevnode, pstart, pterm, rdcost, red_cost, resid, root, secarc, seclevel, start, term, thresh_dfct; /* start initialization using auction */ naug = 0; pass = 0; thresh_dfct = 0; /* factor determines by how much epsilon is reduced at each * minimization */ factor = 3; /* num_passes determines how many auction scaling phases are * performed */ num_passes = 1; /* set arc flows to satisfy cs and calculate maxcost and * mincost */ maxcost = -large; mincost = large; for (arc = 1; arc <= na; arc++) { start = startn[arc]; end = endn[arc]; rdcost = rc[arc]; if (maxcost < rdcost) maxcost = rdcost; if (mincost > rdcost) mincost = rdcost; if (rdcost < 0) { dfct[start] += u[arc]; dfct[end] -= u[arc]; x[arc] = u[arc]; u[arc] = 0; } else x[arc] = 0; } /* set initial epsilon */ if ((maxcost - mincost) >= 8) eps = (maxcost - mincost) / 8; else eps = 1; /* set initial prices to zero */ for (node = 1; node <= n; node++) p[node] = 0; /* Initialization using auction/shortest paths. */ L100: /* Start of the first scaling phase. */ pass++; if ((pass == num_passes) || (eps == 1)) crash = 0; nolist = 0; /* construct list of positive surplus nodes and queue of negative * surplus nodes */ for (node = 1; node <= n; node++) { prdcsr[node] = 0; path_id[node] = false; extend_arc[node] = 0; sb_level[node] = -large; nxtqueue[node] = node + 1; if (dfct[node] > 0) { nolist++; save[nolist] = node; } } nxtqueue[n] = 1; root = 1; prevnode = lastqueue = n; /* initialization with down iterations for negative surplus * nodes */ for (i = 1; i <= nolist; i++) { node = save[i]; nsp++; /* build the list of arcs w/ room for pushing flow and find * proper price for down iteration */ bstlevel = -large; fpushf[node] = 0; for (arc = fou[node]; arc > 0; arc = nxtou[arc]) { if (u[arc] > 0) { if (fpushf[node] == 0) { fpushf[node] = arc; nxtpushf[arc] = 0; last = arc; } else { nxtpushf[last] = arc; nxtpushf[arc] = 0; last = arc; } } if (x[arc] > 0) { new_level = p[endn[arc]] + rc[arc]; if (new_level > bstlevel) { bstlevel = new_level; extarc = arc; } } } fpushb[node] = 0; for (arc = fin[node]; arc > 0; arc = nxtin[arc]) { if (x[arc] > 0) { if (fpushb[node] == 0) { fpushb[node] = arc; nxtpushb[arc] = 0; last = arc; } else { nxtpushb[last] = arc; nxtpushb[arc] = 0; last = arc; } } if (u[arc] > 0) { new_level = p[startn[arc]] - rc[arc]; if (new_level > bstlevel) { bstlevel = new_level; extarc = -arc; } } } extend_arc[node] = extarc; p[node] = bstlevel - eps; } L200: /* Start the augmentation cycles of the new scaling phase. */ if (dfct[root] >= thresh_dfct) goto L3000; term = root; path_id[root] = true; L500: /* Main forward algorithm with root as origin. */ /* start of a new forward iteration */ pterm = p[term]; extarc = extend_arc[term]; if (extarc == 0) { /* build the list of arcs w/ room for pushing flow */ fpushf[term] = 0; for (arc = fou[term]; arc > 0; arc = nxtou[arc]) { if (u[arc] > 0) { if (fpushf[term] == 0) { fpushf[term] = arc; nxtpushf[arc] = 0; last = arc; } else { nxtpushf[last] = arc; nxtpushf[arc] = 0; last = arc; } } } fpushb[term] = 0; for (arc = fin[term]; arc > 0; arc = nxtin[arc]) { if (x[arc] > 0) { if (fpushb[term] == 0) { fpushb[term] = arc; nxtpushb[arc] = 0; last = arc; } else { nxtpushb[last] = arc; nxtpushb[arc] = 0; last = arc; } } } goto L600; } /* speculative path extension attempt */ /* note: arc > 0 means that arc is oriented from the root to the * destinations * arc < 0 means that arc is oriented from the destinations to the * root * extarc = 0 or prdarc = 0, means the extension arc or the * predecessor arc, respectively, has not been established */ if (extarc > 0) { if (u[extarc] == 0) { seclevel = sb_level[term]; goto L580; } end = endn[extarc]; bstlevel = p[end] + rc[extarc]; if (pterm >= bstlevel) { if (path_id[end]) goto L1200; term = end; prdcsr[term] = extarc; path_id[term] = true; /* if negative surplus node is found, do an augmentation */ if (dfct[term] > 0) goto L2000; /* return for another iteration */ goto L500; } } else { extarc = -extarc; if (x[extarc] == 0) { seclevel = sb_level[term]; goto L580; } start = startn[extarc]; bstlevel = p[start] - rc[extarc]; if (pterm >= bstlevel) { if (path_id[start]) goto L1200; term = start; prdcsr[term] = -extarc; path_id[term] = true; /* if negative surplus node is found, do an augmentation */ if (dfct[term] > 0) goto L2000; /* return for another iteration */ goto L500; } } L550: /* second best logic test applied to save a full node scan * if old best level continues to be best go for another * contraction */ seclevel = sb_level[term]; if (bstlevel <= seclevel) goto L800; L580: /* if second best can be used do either a contraction or start * over with a speculative extension */ if (seclevel > -large) { extarc = sb_arc[term]; if (extarc > 0) { if (u[extarc] == 0) goto L600; bstlevel = p[endn[extarc]] + rc[extarc]; } else { if (x[-extarc] == 0) goto L600; bstlevel = p[startn[-extarc]] - rc[-extarc]; } if (bstlevel == seclevel) { sb_level[term] = -large; extend_arc[term] = extarc; goto L800; } } L600: /* extension/contraction attempt was unsuccessful, so scan * terminal node */ nsp++; bstlevel = seclevel = large; for (arc = fpushf[term]; arc > 0; arc = nxtpushf[arc]) { new_level = p[endn[arc]] + rc[arc]; if (new_level < seclevel) { if (new_level < bstlevel) { seclevel = bstlevel; bstlevel = new_level; secarc = extarc; extarc = arc; } else { seclevel = new_level; secarc = arc; } } } for (arc = fpushb[term]; arc > 0; arc = nxtpushb[arc]) { new_level = p[startn[arc]] - rc[arc]; if (new_level < seclevel) { if (new_level < bstlevel) { seclevel = bstlevel; bstlevel = new_level; secarc = extarc; extarc = -arc; } else { seclevel = new_level; secarc = -arc; } } } sb_level[term] = seclevel; sb_arc[term] = secarc; extend_arc[term] = extarc; L800: /* End of node scan. */ /* if the terminal node is the root, adjust its price and change * root */ if (term == root) { p[term] = bstlevel + eps; if (p[term] >= large) { /* no path to the destination */ /* problem is found to be infeasible */ return 1; } path_id[root] = false; prevnode = root; root = nxtqueue[root]; goto L200; } /* check whether extension or contraction */ prd = prdcsr[term]; if (prd > 0) { pr_term = startn[prd]; prevlevel = p[pr_term] - rc[prd]; } else { pr_term = endn[-prd]; prevlevel = p[pr_term] + rc[-prd]; } if (prevlevel > bstlevel) { /* path extension */ if (prevlevel >= bstlevel + eps) p[term] = bstlevel + eps; else p[term] = prevlevel; if (extarc > 0) { end = endn[extarc]; if (path_id[end]) goto L1200; term = end; } else { start = startn[-extarc]; if (path_id[start]) goto L1200; term = start; } prdcsr[term] = extarc; path_id[term] = true; /* if negative surplus node is found, do an augmentation */ if (dfct[term] > 0) goto L2000; /* return for another iteration */ goto L500; } else { /* path contraction */ p[term] = bstlevel + eps; path_id[term] = false; term = pr_term; if (pr_term != root) { if (bstlevel <= pterm + eps) goto L2000; } pterm = p[term]; extarc = prd; if (prd > 0) bstlevel += eps + rc[prd]; else bstlevel += eps - rc[-prd]; /* do a second best test and if that fails, do a full node * scan */ goto L550; } L1200:/* A cycle is about to form; do a retreat sequence. */ node = term; L1600:if (node != root) { path_id[node] = false; prd = prdcsr[node]; if (prd > 0) { pr_term = startn[prd]; if (p[pr_term] == p[node] + rc[prd] + eps) { node = pr_term; goto L1600; } } else { pr_term = endn[-prd]; if (p[pr_term] == p[node] - rc[-prd] + eps) { node = pr_term; goto L1600; } } /* do a full scan and price rise at pr_term */ nsp++; bstlevel = seclevel = large; for (arc = fpushf[pr_term]; arc > 0; arc = nxtpushf[arc]) { new_level = p[endn[arc]] + rc[arc]; if (new_level < seclevel) { if (new_level < bstlevel) { seclevel = bstlevel; bstlevel = new_level; secarc = extarc; extarc = arc; } else { seclevel = new_level; secarc = arc; } } } for (arc = fpushb[pr_term]; arc > 0; arc = nxtpushb[arc]) { new_level = p[startn[arc]] - rc[arc]; if (new_level < seclevel) { if (new_level < bstlevel) { seclevel = bstlevel; bstlevel = new_level; secarc = extarc; extarc = -arc; } else { seclevel = new_level; secarc = -arc; } } } sb_level[pr_term] = seclevel; sb_arc[pr_term] = secarc; extend_arc[pr_term] = extarc; p[pr_term] = bstlevel + eps; if (pr_term == root) { prevnode = root; path_id[root] = false; root = nxtqueue[root]; goto L200; } path_id[pr_term] = false; prd = prdcsr[pr_term]; if (prd > 0) term = startn[prd]; else term = endn[-prd]; if (term == root) { prevnode = root; path_id[root] = false; root = nxtqueue[root]; goto L200; } else goto L2000; } L2000:/* End of auction/shortest path routine. */ /* do augmentation from root and correct the push lists */ incr = -dfct[root]; for (node = root;;) { extarc = extend_arc[node]; path_id[node] = false; if (extarc > 0) { node = endn[extarc]; if (incr > u[extarc]) incr = u[extarc]; } else { node = startn[-extarc]; if (incr > x[-extarc]) incr = x[-extarc]; } if (node == term) break; } path_id[term] = false; if (dfct[term] > 0) { if (incr > dfct[term]) incr = dfct[term]; } for (node = root;;) { extarc = extend_arc[node]; if (extarc > 0) { end = endn[extarc]; /* add arc to the reduced graph */ if (x[extarc] == 0) { nxtpushb[extarc] = fpushb[end]; fpushb[end] = extarc; new_level = p[node] - rc[extarc]; if (sb_level[end] > new_level) { sb_level[end] = new_level; sb_arc[end] = -extarc; } } x[extarc] += incr; u[extarc] -= incr; /* remove arc from the reduced graph */ if (u[extarc] == 0) { nas++; arc = fpushf[node]; if (arc == extarc) fpushf[node] = nxtpushf[arc]; else { prevarc = arc; arc = nxtpushf[arc]; while (arc > 0) { if (arc == extarc) { nxtpushf[prevarc] = nxtpushf[arc]; break; } prevarc = arc; arc = nxtpushf[arc]; } } } node = end; } else { extarc = -extarc; start = startn[extarc]; /* add arc to the reduced graph */ if (u[extarc] == 0) { nxtpushf[extarc] = fpushf[start]; fpushf[start] = extarc; new_level = p[node] + rc[extarc]; if (sb_level[start] > new_level) { sb_level[start] = new_level; sb_arc[start] = extarc; } } u[extarc] += incr; x[extarc] -= incr; /* remove arc from the reduced graph */ if (x[extarc] == 0) { nas++; arc = fpushb[node]; if (arc == extarc) fpushb[node] = nxtpushb[arc]; else { prevarc = arc; arc = nxtpushb[arc]; while (arc > 0) { if (arc == extarc) { nxtpushb[prevarc] = nxtpushb[arc]; break; } prevarc = arc; arc = nxtpushb[arc]; } } } node = start; } if (node == term) break; } dfct[term] -= incr; dfct[root] += incr; /* insert term in the queue if it has a large enough surplus */ if (dfct[term] < thresh_dfct) { if (nxtqueue[term] == 0) { nxtnode = nxtqueue[root]; if ((p[term] >= p[nxtnode]) && (root != nxtnode)) { nxtqueue[root] = term; nxtqueue[term] = nxtnode; } else { nxtqueue[prevnode] = term; nxtqueue[term] = root; prevnode = term; } } } /* if root has a large enough surplus, keep it in the queue and * return for another iteration */ if (dfct[root] < thresh_dfct) { prevnode = root; root = nxtqueue[root]; goto L200; } L3000:/* end of augmentation cycle */ /* Check for termination of scaling phase. If scaling phase is not * finished, advance the queue and return to take another node. */ nxtnode = nxtqueue[root]; if (root != nxtnode) { nxtqueue[root] = 0; nxtqueue[prevnode] = nxtnode; root = nxtnode; goto L200; } /* End of subproblem (scaling phase). */ /* Reduce epsilon. */ eps /= factor; if (eps < 1) eps = 1; thresh_dfct /= factor; if (eps == 1) thresh_dfct = 0; /* if another auction scaling phase remains, reset the flows & * the push lists; else reset arc flows to satisfy cs and compute * reduced costs */ if (crash == 1) { for (arc = 1; arc <= na; arc++) { start = startn[arc]; end = endn[arc]; pstart = p[start]; pend = p[end]; if (pstart > pend + eps + rc[arc]) { resid = u[arc]; if (resid > 0) { dfct[start] += resid; dfct[end] -= resid; x[arc] += resid; u[arc] = 0; } } else if (pstart < pend - eps + rc[arc]) { flow = x[arc]; if (flow > 0) { dfct[start] -= flow; dfct[end] += flow; x[arc] = 0; u[arc] += flow; } } } /* return for another phase */ goto L100; } else { crash = 1; for (arc = 1; arc <= na; arc++) { start = startn[arc]; end = endn[arc]; red_cost = rc[arc] + p[end] - p[start]; if (red_cost < 0) { resid = u[arc]; if (resid > 0) { dfct[start] += resid; dfct[end] -= resid; x[arc] += resid; u[arc] = 0; } } else if (red_cost > 0) { flow = x[arc]; if (flow > 0) { dfct[start] -= flow; dfct[end] += flow; x[arc] = 0; u[arc] += flow; } } rc[arc] = red_cost; } } return 0; # undef crash # undef nsp } /* eof */