Fortran subroutines for network (cid:13)ow optimization using an interior point algorithm

We describe Fortran subroutines for network (cid:13)ow optimization using an interior point network (cid:13)ow algorithm, that, together with a Fortran language driver, make up PDNET. The algorithm is described in detail and its implementation is outlined. Usage of the package is described and some computational experiments are reported. Source code for the software can be downloaded at http://www.research.att.com/~mgcr/ pdnet .


Introduction
Given a directed graph ( , ) G = V E , where V is a set of m vertices and E a set of n edges, let ( , ) i j denote a directed edge from vertex i to vertex j . The minimum cost network flow problem can be formulated as ( , )  Most often, the problem data are assumed to be integer. In matrix notation, the above network flow problem can be formulated as a primal linear program of the form where c is a m n × vector whose elements are ij c , A is the m n × vertex-edge incidence matrix of the graph ( , ) G = V E , i.e. for each edge ( , ) i j in E there is an associated column in matrix A with exactly two nonzero entries: an entry 1 in row i and an entry 1 − in row j ; b , x , and u are defined as above, and s is an n-dimensional vector of upper bound slacks. Furthermore, an appropriate variable change allows us to assume that the lower bounds l are zero, without loss of generality. The dual of (3) is where y is the m -dimensional vector of dual variables and w and z are n -dimensional vectors of dual slacks.
If graph G is disconnected and has p connected components, there are exactly p redundant flow conservation constraints, which are sometimes removed from the problem formulation. Without loss of generality, we rule out trivially infeasible problems by assuming 0, 1, , where k V is the set of vertices for the k -th component of G .
If it is required that the flow ij x be integer, (2) is replaced with , integer, ( , ) In the description of the algorithm, we assume without loss of generality, that 0 ij l = for all ( , ) i j ∈E and that 0 c ≠ . A simple change of variables is done in the subroutines to transform the original problem into an equivalent one with 0 ij l = for all ( , ) i j ∈E . The flow is transformed back to the original problem upon termination. The case where 0 c = is a simple feasibility problem, and is handled by solving a maximum flow problem [1].
Before concluding this introduction, we present some notation and outline the remainder of the paper. We denote the i -th column of A by i A , the i -th row of A by i A i and a submatrix of A formed by columns with indices in set S by S A . Let n x ∈ R . We denote by X the n n × diagonal matrix having the elements of x in the diagonal. The Euclidean or 2-norm is denoted by · ‖‖ .
This paper describes Fortran subroutines used in an implementation of PDNET, an interior point network flow method introduced in Portugal, Resende, Veiga & Júdice [12]. The paper is organized as follows. In Section 2 we review the truncated infeasible-primal feasible-dual interior point method for linear programming. The implementation of this algorithm to handle network flow problems is described in Section 3. Section 4 describes the subroutines and their usage. Computational results, comparing PDNET with the network optimizer in CPLEX 10, are reported in Section 5. Concluding remarks are made in Section 6.

Truncated primal-infeasible dual-feasible interior point algorithm
In this section, we recall the interior point algorithm implemented in PDNET. Let x y s w z ∆ ∆ ∆ ∆ ∆ is obtained as the solution of the linear system of equations where e is a vector of ones of appropriate order, k r is such that Primal and dual steps are taken in the direction ( , , , where p α and d α are step-sizes in the primal and dual spaces, respectively, and are given by The solution of the linear system (5) is obtained in two steps. First, we compute the k y ∆ component of the direction as the solution of the system of normal equations where k Θ is given by The remaining components of the direction are then recovered by

Implementation
We discuss how the truncated primal-infeasible dual-feasible algorithm can be used for solving network flow problems. For ease of discussion, we assume, without loss of generality, that the graph G is connected. However, disconnected graphs are handled by PDNET.

Computing the Newton direction
Since the exact solution of (10) can be computationally expensive, a preconditioned conjugate gradient (PCG) algorithm is used to compute approximately an interior point search direction at each iteration. The PCG solves the linear system where M is a positive definite matrix and k Θ is given by (11), and The aim is to make the preconditioned matrix , and improve the efficiency of the conjugate gradient algorithm by reducing the number of iterations it takes to find a feasible direction.
Pseudo-code for the preconditioned conjugate gradient algorithm implemented in PDNET is presented in Figure 1. The matrix-vector multiplications in line 7 are of the form and can be carried out without forming k A A Θ explicitly. PDNET uses as its initial direction 0 y ∆ the direction y ∆ produced in the previous call to the conjugate gradient algorithm, i.e.
during the previous interior point iteration. The first time pcg is called, we assume The preconditioned residual is computed in lines 3 and 11 when the system of linear equations is solved. PDNET uses primal-dual variants of the diagonal and spanning tree preconditioners described in [15,16].    O m divisions. This preconditioner has been shown to be effective during the initial interior point iterations [11,14,15,16,18].
In the spanning tree preconditioner [16], one identifies a maximal spanning tree of the graph G , using as weights the diagonal elements of the current scaling matrix, where e is a unit n -vector. An exact maximal spanning tree is computed with the Fibonacci heap variant of Prim's algorithm [13], as described in [1]. At the k -th interior point iteration, let be the indices of the edges of the maximal spanning tree. The spanning tree preconditioner is For simplicity of notation, we include in A heuristic is used to select the preconditioner. The initial selection is the diagonal preconditioner, since it tends to outperform the other preconditioners during the initial interior point iterations. The number of conjugate gradients taken at each interior point iteration is monitored. If the number of conjugate gradient iterations exceeds / 4 m , the current computation of the direction is discarded, and a new conjugate gradient computation is done with the spanning tree preconditioner. The diagonal preconditioner is not used again. The diagonal preconditioner is limited to at most 30 interior point iterations. If at iteration 30 the diagonal preconditioner is still in effect, at iteration 31 the spanning tree preconditioner is triggered. Also, as a safeguard, a hard limit of 1000 conjugate gradient iterations is imposed.
To determine when the approximate direction and tightens the tolerance after each interior point iteration k as follows: has the complexity of one conjugate gradient iteration and should not be carried out every conjugate gradient iteration. Since ( ) Since, on network linear programs, the preconditioned conjugate gradient method finds good directions in few iterations, this estimate is quite accurate in practice. Since it is inexpensive, it is computed at each conjugate gradient iteration.

Stopping criteria for interior point method
In [15], two stopping criteria for the interior point method were used. The first, called the primal-basic (PB) stopping rule, uses the spanning tree computed for the tree preconditioner.
If the network flow problem has a unique solution, the edges of the tree converge to the optimal basic sequence of the problem. Let T be the index set of the edges of the tree, and define to the index set of edges that are fixed to their upper bounds. If the solution * x T of the linear system integer, then * x T has only integer components. Optimality of * x T can be verified by computing a lower bound on the optimal objective function value. This can be done with a strategy introduced independently in [15] and [10,17]. Denote by * A tentative optimal dual solution * y (having a possibly better objective value than the current dual interior point solution k y ) can be found by orthogonally projecting k y onto the supporting affine space of the dual face complementary to * x T . In an attempt to preserve dual feasibility, we compute * y as the solution of the least squares problem Resende & Veiga [15] describe a ( ) O m operation procedure to compute this projection. A feasible dual solution * * * ( , , ) y z w is built by adjusting the dual slacks. Let x s and * * * ( , , ) y w z are optimal primal and dual solutions, respectively. If the data is integer and x s is a primal optimal (integer) solution.
To apply the second stopping procedure of [15], called the maximum flow (MF) stopping criterion, an indicator function to partition the edge set into active and inactive (fixed at upper or lower bounds) is needed. In PDNET, the indicator used is the so-called primal-dual indicator, studied by Gay [5] and El-Bakry, Tapia & Zhang [4]. Let ξ be a small tolerance. Edge i is classified as inactive at its lower bound if 1 and Edge i is classified as inactive at its upper bound if 1 and The remaining edges are set active. In PDNET, ξ is initially set to 3 10 − and this tolerance is We select a tentative optimal dual face F as a maximum weighted spanning forest limited to the active edges as determined by the indicator. The edge weights used in PDNET are those of the scaling matrix k Θ .
As in the PB indicator, we project the current dual interior solution k y orthogonally onto F . Once the projected dual solution * y is at hand, we attempt to find a feasible flow * x complementary to * y . A refined tentative optimal face is selected by redefining the set of active edges as where r ε is a small tolerance ( By considering only the active edges, a restricted network is built. Flow on this network must satisfy Clearly, from the flow balance constraints (16) In addition, for each edge ( , ) i j ∈ F there is an associated capacity ij u . Let The if and only if x F is a feasible flow for the restricted network [15]. Therefore, finding a feasible flow for the restricted network involves the solution of a maximum flow problem. Furthermore, if the data is integer, this feasible flow is integer, as we can select a maximum flow algorithm that provides an integer solution.
Since this stopping criterion involves the solution of a maximum flow problem, it should not be triggered until the interior point algorithm is near the optimal solution. The criterion is triggered at iteration k , when k µ µ < ε occurs for first time. The choice 1 µ = ε used in PDNET is appropriate for the set of test problems considered here. In a more general purpose implementation, a scale invariant criterion is desirable. All subsequent iterations test this stopping rule. In PDNET, the implementation of Goldfarb & Grigoriadis [6] of Dinic's algorithm is used to solve the maximum flow problems.

Other implementation issues
To conclude this section, we make some remarks on other important implementation issues of the primal-infeasible, dual-feasible algorithm, namely the starting solution, the adjustment of parameter k µ , and the primal and dual stepsizes.
The stepsize parameters p and d are both set to 0.995 throughout the iterations, slightly more conservative than as suggested by [9].

Fortran subroutines
The current implementation PDNET consists of a collection of core subroutines with additional subsidiary modules written in Fortran [7]. The software distribution associated with this paper provides fully functional utilities by including Fortran reference implementations for the main program, providing default parameter settings and supplying additional routines for data input and output functions. In this section, we describe the usage of the PDNET framework with comments on user extensions. Specific instructions for building and installing PDNET from the source code are provided with the software distribution. Table 1 lists all modules provided in the PDNET software distribution. Fortran functions for reading network data We adopted several design and programming style guidelines in implementing the PDNET routines in Fortran. We required that all arrays and variables be passed as subroutine arguments, avoiding the use of COMMON blocks. The resulting code is expected to compile and run without modification under most software development and computing environments.

PDNET core subroutines
The PDNET core subroutines are invoked via a single interface provided by subroutine pdnet(). Following the specifications listed in Table 2, the calling program must provide data via the input arguments and allocate the internal and output arrays appropriately as described in Subsection 4.3 in file fdriver.f90.
We provide reference implementations of the PDNET main program which also serve as guides for developing custom applications invoking the PDNET core subroutines. In Subsection 4.2, we discuss in detail the input and output routines used in the reference implementations. Subsection 4.4 discusses the setting of parameters in PDNET. In addition, the core subroutines call an external function for maximum flow computation, which is provided in a subsidiary module, with its interface discussed in Subsection 4.5.

Data input
Programs invoking the PDNET subroutines must supply an instance for the minimum-cost flow problem. As presented in Table 3, the data structure includes the description of the underlying graph, node netflow values, arc costs, capacities and lower bounds. All node and arc attributes are integer valued.   The reference implementations of the main program PDNET reads network data from an input file by invoking functions available in the Fortran module pdnet_read.f90. These functions build PDNET input data structures from data files in the DIMACS format [2] with instances of minimum cost flows problems. As illustrated in the PDNET module fdriver.f90, we provide specialized interfaces used in Fortran programs using dynamic memory allocation.

Memory allocation
In PDNET, total memory utilization is carefully managed by allocating each individual array to a temporary vector passed as argument to internal PDNET functions. Furthermore, input and output arrays, presented in Table 4, are passed as arguments to subroutine pdnet() and must be allocated by the calling procedure.

Parameter setting
PDNET execution is controlled by a variety of parameters, selecting features of the underlying algorithm and the interface with the calling program. A subset of these parameters is exposed to the calling procedure via a PDNET core subroutine. The remaining parameters are set at compile time with default values assigned in inside module pdnet_default.f90.
The run time parameters listed in Tables 5 and 6 are set with Fortran function pdnet_setintparm(). The integer parameters are assigned to components of vector info and double precision parameters are assigned to vector dinfo.

Maximum flow computation
PDNET includes a maximum flow computation module called pdnet_maxflow, featuring the implementation of Goldfarb and Grigoriadis [6] of Dinic's algorithm, that is used to check the maximum flow stopping criterion. Furthermore, a modification of this module, called pdnet_feasible, is called by the module pdnet_checkfeas() after reading the network file to compute a maximum flow on the network, therefore checking for infeasibility.

Running PDNET
Module fdriver inputs the network structure and passes control to module pdnet(). This subroutine starts by reading the control parameters from the vector info, with pdnet_getinfo(). Correctness of the input is checked with pdnet_check(), and data structure created with pdnet_datastruct(). Internal parameters for methods are set with pdnet_default(), and additional structures are created with pdnet_buildstruct(). Subroutines pdnet_transform() and pdnet_perturb() are called in order to shift the lower bounds to zero and to transform the data into double precision, respectively. Subroutines pdnet_probdata() and pdnet_checkfeas() check the problem characteristics and inquire if the network has enough capacity to transport the proposed amount of commodity. The primal-dual main loop is then started and the right-hand-side of the search direction system is computed by pdnet_comprhs(). The maximum spanning tree is computed by pdnet_heap(), and its optimality is tested by pdnet_optcheck(). Under certain conditions ( 1 µ < ), the maxflow stopping criterion is invoked by subroutine pdnet_checkmaxflow(). If required, preconditioner switch takes place, followed by a call to pdnet_precconjgrd() to solve the Newton direction linear system using the chosen preconditioner. A summary of the iteration is printed by pdnet_printout(). Primal and dual updates are made by pdnet_updatesol() and stopping criteria check takes place, before returning to the start of the iteration loop.
We now present a usage example. Let us consider the problem of finding the minimum cost flow on the network represented in Figure 2. In this figure, the integers close to the nodes represent the node's produced (positive value) or consumed (negative value) net flows, and the integers close to the arcs stand for their unit cost. Furthermore the capacity of each arc is bounded between 0 and 10. In Figure 3, we show the DIMACS format representation of this problem, stored in file test.min. Furthermore, Figures 4 and 5 show the beginning and the end of the printout given by PDNET, respectively.

Computational results
A preliminary version of PDNET was tested extensively with results reported in [12]. In this section, we report on a limited experiment with the version of PDNET in the software distribution.
The experiments were done on a PC with an Intel Pentium IV processor running at 2 Ghz and 2 Gb of main memory. The operating system is Ubuntu Linux 7.10 (kernel 2.6.22). The code was compiled on the Intel Fortran compiler version 10.0 using the flags -O3. CPU times in seconds were computed by calling the Fortran 90 function cpu_time(). The test problems are instances of the classes mesh, grid and netgen_lo of minimum-cost network flow problems, taken from the First DIMACS Implementation Challenge [3]. The specifications of the mesh instances generated are presented in Table 7. The specifications used in the GRIDGEN generator to build the grid problems are displayed in Table 8. Finally, the instances of the test set netgen_lo were generated according to the guidelines stated in the computational study of Resende and Veiga [15], using the NETGEN generator following the specifications presented in Table 9. All these generators can be downloaded from the FTP site dimacs.rutgers.edu.  Instances of increasing dimension were considered. The three netgen_lo instances presented were generated considering 9 x = , 13 x = and 15 x = respectively. In Table 10, the number of iterations (IT) and the CPU time (CPU) of PDNET and CPLEX 10 are compared. The reported results show that PDNET tends to outperform CPLEX in the larger instances of the mesh and netgen_lo set problems, but fails to do so in grid set problems. For testset netgen_lo the difference is quite noticeable for the largest instance. We observed that in this problem about 90% of the CPU time and iteration count in CPLEX 10 was spent computing a feasible solution. The results in this experience and those reported in [12] show that this code is quite competitive with CPLEX and other efficient network flow codes for large-scale problems.

Concluding remarks
In this paper, we describe a Fortran implementation of PDNET, a primal-infeasible dualfeasible interior point method for solving large--scale linear network flow problems. The subroutines are described, and directions for usage are given. A number of technical features of the implementation, which enable the user to control several aspects of the program execution, are also presented. Some computational experience with a number of test problems from the DIMACS collection is reported. These results illustrate the efficiency of PDNET for the solution of linear network flow problems. Source code for the software is available for download at http://www.research.att.com/~mgcr/pdnet.