% % Large Scale Nonlinear Network Optimization % User's Guide % % This version: March 2nd, 1993. % \documentstyle[11pt,a4]{article} \newcommand{\R}[1]{{\bf R}^{#1}} \newcommand{\tim}[1]{\;\; \mbox{#1} \;\;} \begin{document} \begin{titlepage} \vspace*{4.5 cm} \begin{center} \begin{minipage}{8.cm} \begin{center} {\sc LSNNO, A Fortran Subroutine for \\ Large Scale Nonlinear Network \\ Optimization problems. \\ ( User's Guide ) } \\ \mbox{} \\ by D. Tuyttens \\ \mbox{} \\ Report 92/25 \hfill \today \end{center} \end{minipage} \end{center} \vspace{2.cm} \begin{center} \parbox{12.cm}{ {\bf Abstract.} This report contains the specification data and user's guide for the routine LSNNO. This routine aims at solving large scale nonlinear network optimization problems. A detailed description of the calling sequence is provided as well as a practical example of use. } \end{center} \vspace{5. cm} \begin{center} {\bf Department of Mathematics \\ Facult\'{e}s Universitaires ND de la Paix \\ B-5000 Namur, Belgium } \end{center} \end{titlepage} \section{Summary} LSNNO is a Fortran Subroutine designed to solve large scale nonlinear network optimization problems, i.e. problems of the form : \[ \begin{array}{lc} \min_{{\bf x} \in \R{n}} & f({\bf x}) \\ \mbox{subject to} & A{\bf x} = b \\ & l \leq {\bf x} \leq u, \end{array} \] where $f({\bf x})$ is a nonlinear {\em objective function}, that is twice continuously differentiable on the feasible set defined by the above constraints, where the matrix $A$ is a $m \times n$ node-arc incidence matrix, where the supply/demand vector $b$ is in $\R{m}$ and satisfies $\sum_{i=1}^m b_i = 0$, and where the vectors $l, u \in \R{n}$ are the lower and upper bounds on the flows in the oriented network associated with $A$. The routine is especially effective on problems involving a large number of variables. The objective function $f(x)$ is of the {\em partially separable} type (see Griewank and Toint - 1982), i.e. it has the form \[ f({\bf x}) = \sum_{i=1}^{nel} f_i({\bf x}_i) \hspace{2. cm} {\bf x} = ( x_1,x_2,\ldots,x_n) \] where the ${\bf x}_i$ are either small subsets of ${\bf x}$ or such that the rank of the second derivative matrix of the nonlinear element function $f_i$ is small for some other reason. The subroutine passes control back to the user to calculate function and derivative values of the element $f_i$. A considerable flexibility in the computation of the derivatives is provided, i.e, if they are not supplied by the user, the program contains an assortment of routines for computing an approximation of these quantities. If there is a linear transformation of variables such that one or more of the element functions depend on fewer transformed variables than the number of original variables ${\bf x}_i$ they involve, this may be specified. This leads to substantially improved computing time. \section{How to use the routine} \subsection{The argument list} \begin{verbatim} CALL LSNNO( NA, NN, NEL, ELVAR, ELPTR, ELST, FR, TO, + X, XST, FX, FUVAL, LFUVAL, GPTR, HPTR, + INFORM, FCALC, LFC, IFC, IT, MAXIT, EPSF, + PREC, NIPST, IWK, LIWK, WK, LWK, IPDEVC, + WHAT, FREQ, INFO, IFLAG ) \end{verbatim} \begin{description} \item[NA ] is an INTEGER variable, that must be set by the user to the number of arcs in the network; this represents also the number of variables of the problem. It is not altered by the routine. {\bf Restriction:} NA $> 0$. \item[NN ] is an INTEGER variable, that must be set by the user to the number of nodes in the network. It is not altered by the routine. {\bf Restriction:} NN $> 0$. \item[NEL] is an INTEGER variable, that must be set by the user to the number of element functions. It is not altered by the routine. {\bf Restriction:} NEL $> 0$. \item[ELVAR ] is an INTEGER array containing the indices of the variables in the first element, followed by those in the second element, etc. See \S~5 for an example. It is not altered by the routine. \item[ELPTR ] is an INTEGER array of length at least NEL+1, whose kth value is the position of the first variable of element function k, in the list ELVAR. In addition, ELPTR(NEL+1) must be equal to the total length of ELVAR plus one. See \S~5 for an example. It is not altered by the routine. \item[ELST ] is an INTEGER array of length NEL, that must be set by the user as follows. \begin{description} \item[ELST(i) = 1 : ] ith element Hessian is provided by the user, \item[ELST(i) = 2 : ] ith element Hessian is computed by differences in the gradient values, \item[ELST(i) = 3 : ] ith element Hessian is computed by differences in the function values, \item[ELST(i) = 4 : ] ith element Hessian is updated by a secant approximation ( BFGS or Rank-One formula depending on the convexity of the ith element function ) and initialized analitically, \item[ELST(i) = 5 : ] ith element Hessian is updated by a secant approximation and initialized by differences in the gradient values, \item[ELST(i) = 6 : ] ith element Hessian is updated by a secant approximation and initialized by differences in the function values, \item[ELST(i) = 7 : ] ith element Hessian is updated by a secant approximation and initialized by a correctly scaled multiple of the identity matrix. \item[ELST(i) = 8 : ] the product of the ith Hessian matrix by a vector is accomplished using Dembo's technique ( by difference in the gradient values ). In this case, the ith element Hessian matrix need not to be stored. \end{description} ELST(i) $> 0$ means that the gradient of the ith element function is supplied by the user. If ELST(i) $< 0$ ( opposite of value defined above ) the gradient of the ith element function is estimated by differences. Analytical gradients are always preferable when available. Some restrictions : \begin{description} \item[ELST(i) = -2 ] is automatically replaced by ELST(i) = -3 \item[ELST(i) = -5 ] is automatically replaced by ELST(i) = -6 \item[ELST(i) = -8 ] is not available. \end{description} Note that IT IS ALTERED by the routine and that it is meaningless on output. \item[FR ] is an INTEGER array of length NA, which must be set by the user, such that its ith component represents the origine node of arc i. It is not altered by the routine. \item[TO ] is an INTEGER array of length NA, which must be set by the user, such that its ith component represents the destination node of arc i. It is not altered by the routine. \item[X ] is a DOUBLE PRECISION array of length NA which must be set by the user to the value of the variables at the starting point. The subroutine checks if this starting point satisfies the constraints of the problem. If it is not the case, a new feasible starting point is computed using a subroutine called SIMPLX. On exit, it contains the values of the variables at the best point found (usually the solution). \item[XST ] is an INTEGER array of length NA, that must be set by the user as follows. XST(i) must be set to -1 if the ith variable is unconstrained, to 0 if it is fixed, and to 1 if it is bounded. The values of fixed variables are never changed during the minimization process. Note that IT IS ALTERED by the routine and that it is meaningless on output. \item[FX ] is DOUBLE PRECISION variable. It contains the value of the objective function at the point X. On exit, it is the best function value known. \item[FUVAL ] is a DOUBLE PRECISION array of length LFUVAL, which is used to store function and derivative information for the element functions. The user is asked to provide values for these functions and/or derivatives, evaluated at the point X, at specified locations within FUVAL, when control is returned to the calling program with INFORM $> 0$ (see INFORM). The layout of FUVAL is shown in figure~1. \setlength{\unitlength}{0.002 cm} \begin{figure}[htbp] \begin{center} \begin{picture}(6000,1600)(10,10) \put(10,1050){\framebox(2000,750){}} \put(2010,1050){\framebox(2000,750){}} \put(4010,1050){\framebox(2000,750){}} \put(610,1300){\shortstack{element\\ values}} \put(2610,1300){\shortstack{element\\ gradients}} \put(4610,1300){\shortstack{element\\ Hessians}} \put(10,700){\vector(0,1){300}} \put(2010,700){\vector(0,1){300}} \put(4010,700){\vector(0,1){300}} \put(6010,700){\vector(0,1){300}} \put(10,500){FUVAL(1)} \put(2010,500){GPTR(1)} \put(2010,250){FUVAL(NEL+1)} \put(4010,500){HPTR(1)} \put(4010,250){GPTR(NEL+1)} \put(6010,500){LFUVAL} \put(6010,250){HPTR(NEL+1)-1} \end{picture} \caption{Partitioning of the workspace array FUVAL} \end{center} \end{figure} The first segment of FUVAL contains the values of the element functions; the next two segments contain their gradients and Hessians matrices. \item[LFUVAL ] is an INTEGER variable that must be set by the user to the length of FUVAL in the calling program. It is not altered by the routine. \item[GPTR ] is an INTEGER array of length NEL+1, which need not be set on entry. The array is set in LSNNO so that its ith value ($ 1 \leq i \leq NEL $) gives the position in the array FUVAL of the first component of the gradient of the ith element function, {\em with respect to its internal variables} (see FUVAL and \S~2.2). Note that GPTR(1) = NEL+1 and that GPTR(NEL+1) = HPTR(1). \item[HPTR ] is an INTEGER array of length NEL+1, which need not be set on entry. The array is set in LSNNO so that its ith value ($ 1 \leq i \leq NEL $) gives the position in the array FUVAL of the first component of the Hessian of the ith element function, {\em with respect to its internal variables} (see FUVAL and \S~2.2). Only the upper triangular part of the Hessian matrix is stored and the storage is by columns. Note that HPTR(1) = GPTR(NEL+1) and that HPTR(NEL+1) = LFUVAL+1. \item[INFORM ] is an INTEGER variable that is set within LSNNO, and passed back to the user to prompt further action ( INFORM $> 0$ ). Possible values of INFORM are detailled on \S~2.2. On initial entry, INFORM must be set to 0. \item[FCALC ] is an INTEGER array of length NEL, which need not be set on entry. The LFC first components are set in LSNNO to the indices of the element functions for which the user must provide the value and/or derivatives at X. Precisely what element function is required at X is under control of the variable INFORM (see \S~2.2). \item[LFC ] is an INTEGER variable which need not be set on entry. It is set in LSNNO to the number of element functions for which information is needed at X. Note that LFC $\leq$ NEL. \item[IFC ] is an INTEGER array of length 3, which need not be set by the user and will contain, on exit, the actual number of calls to the function and derivatives that were made during the routine's execution. \begin{description} \item[IFC(1) : ] number of element function evaluations, \item[IFC(2) : ] number of element gradient evaluations, \item[IFC(3) : ] number of element Hessians evaluations. \end{description} \item[IT ] is an INTEGER array of length 3, which need not be set by the user and will contain, on exit, the actual number of iterations performed by the routine. \begin{description} \item[IT(1) : ] number of major iterations, \item[IT(2) : ] number of minor iterations, \item[IT(3) : ] number of C.G. iterations. \end{description} \item[MAXIT ] is an INTEGER variable that must be set by the user to the maximum number of minor iterations that the algorithm will be allowed to perform. It is not altered by the user. \item[EPSF ] is a DOUBLE PRECISION variable, that must be set by the user to a measure of the accuracy ( infinity norm of the reduced gradient ) required to stop the minimization. It is not altered by the routine. \item[PREC ] is a LOGICAL variable, that must be set to .true. if the user whishes to use a preconditioned conjugate gradient scheme to solve the linear systems of equations which arise at every iteration of the minimization algorithm ( A diagonal preconditioner is provided by the routine ). Otherwise no preconditioner is used. It is not altered by the routine. \item[NIPST ] is a LOGICAL variable, that must be set to .true. if the user whishes to use the concept of independent superbasic sets in order to exploit the parrallelism in both the constraints and the cost function structure. Otherwise no decomposition is achieved. It is not altered by the routine. \item[IWK ] is an INTEGER array of length LIWK used as workspace by LSNNO. \item[LIWK ] is an INTEGER variable which must be set by the user to the length of the working array IWK. It is not altered by the routine. \item[WK ] is an DOUBLE PRECISION array of length LWK used as workspace by LSNNO. \item[LWK ] is an INTEGER variable which must be set by the user to the length of the working array WK. It is not altered by the routine. \item[IPDEVC ] is an INTEGER variable, that must be set by the user to the output device unit number for printing of messages. It is not altered by the routine. \item[WHAT ] is an INTEGER variable, that must be set by the user to specify the amount of output generated when output occurs. \begin{description} \item[ WHAT = 0 : ] no output is generated. \item[ WHAT = 1 : ] only the problem specifications and the statistics of the execution will be performed. \item[ WHAT = 2 : ] in addition to the information output with WHAT = 1, a simple one line description of each minor iteration is given. This includes the minor iteration number, the pivoting iteration number, the current function value, the (infinity and two) norm of the reduced gradient, the total number of function, gradient and Hessian evaluations, the number of conjugate-gradient iterations that have been performed and the number of superbasics variables at current minor iteration. At each new major iteration, the following informations are given : the major iteration number, the number of de-activated variables and the number of independent superbasic sets. \item[ WHAT = 3 : ] the output is more detailed that in the case WHAT = 2 and the layout is quite different. More informations are given about the number of variables( basic, superbasic, quasi-active ), the conjugate-gradient procedure( residual norm, conditionning estimate ), the linesearch( basic and superbasic feasible steplength, steplength used, number of trials, reduction in the objective function ), the Hessian update( number of BFGS and RANK-ONE updates if a secant method is used). \item[ WHAT = 4 : ] in addition to the information output with WHAT = 3, the vector of the current iterate is printed. \end{description} This argument is not altered by the routine. \item[FREQ ] is an INTEGER variable, that must be set by the user to specify the frequency of output by the routine LSNNO. \begin{description} \item[ FREQ $\leq 0$ : ] no iteration output is generated (except warning messages), \item[ FREQ $> 0$ : ] output every FREQ iteration (if WHAT $\geq 2$). \end{description} This argument is not altered by the routine. \item[INFO ] is an INTEGER variable. It need not be set by the user and contains , on exit under an error condition, information about the error. See \S~2.5 for further details. \item[IFLAG ] is an INTEGER variable, that need not be set by the user and contains, on exit, the exit condition of the routine. If this flag is positive, an error has been detected by the routine. See \S~2.5 for further details. \end{description} \subsection{Reverse communication} LSNNO uses the reverse communication concept. In this case, no special calling sequence is to be implemented by the user to provide a routine that evaluates the element functions ( and possibly their gradients and Hessians). An example is shown in \S~5. When a return is made from LSNNO with INFORM $> 0$, LSNNO is asking the user for further information. The user should compute the required information and re-enter LSNNO with INFORM unchanged ( the value of INFORM indicates where a return is made in LSNNO). Possible values of INFORM and the information required are : \begin{description} \item[Inform = ] 1 \\ Only function values are required, \item[Inform = ] 2 \\ Only gradient values are required, \item[Inform = ] 3 \\ Only Hessian values are required, \item[Inform = ] 4 \\ Gradient and Hessian values are required. \end{description} \subsection{Internal variables and subroutine RANGE} We now turn to describe the way in which the user can pass to the algorithm information about a linear transformation that reduces the number of independent variables upon which each element function depends. For example the function \[ (x_1+x_2)\sin(x_1-x_3) \] has three {\em elemental} variables $x_1,x_2$ and $x_3$. However, this example illustrates an additional point. After the change of variables $v_1=x_1+x_2$ and $v_2=x_1-x_3$ the function can be written as the product, \[ v_1\sin(v_2) \] The variables $v_1$ and $v_2$ are known as {\em internal} variables for the function, and are obtained as {\em linear combinations} of the elemental variables. In this example, the linear transformation from the elemental variables to the internal ones is defined by \[ \left( \begin{array}{c} v_1 \\ v_2 \end{array} \right) = \left( \begin{array}{ccc} 1 & 1 & 0 \\ 1 & 0 & -1 \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) \] In general the transformation will be of the form \[ v = R x \] and this transformation is {\em useful} if the matrix R has fewer rows than columns. The Hessian matrix can be held as \[ R^TCR \] where $C$ is a square symmetric matrix of order the number of internal variables. The user should therefore provide a routine call RANGE, that defines the transformation between internal and elemental variables for the considered function. This routine has the following argument list : \begin{verbatim} SUBROUTINE RANGE( IEL, ELDIM, INTDIM, W1, W2, MODE ) \end{verbatim} \begin{description} \item[IEL ] is an INTEGER variable , that contains the index of the element function considered. It must be left unchanged by the routine. \item[ELDIM ] is an INTEGER variable, containing the number of elemental variables in the element. It must be left unchanged by the routine. \item[INTDIM ] is an INTEGER variable that must be set by RANGE to the number of internal variables in the element. \item[W1 ] is a DOUBLE PRECISION array containing the input vector $w_1$(see MODE). It must be left unchanged by the routine. \item[W2 ] is a DOUBLE PRECISION array that must be set by the routine to the result vector $w_2$, related to R and $w_1$(see MODE). \item[MODE ] is an INTEGER variable, that contains a code for the work to be accomplished by the routine : \begin{description} \item[MODE = 0 : ] the user should define INTDIM to the number of internal variables of the considered element. Note that $INTDIM \leq ELDIM$. If INTDIM = ELDIM, the routine will never be called with another value of MODE. \item[MODE = 1 : ] the user should provide the vector $w_2$ such that $Rw_1=w_2$. \item[MODE = 2 : ] the user should provide the vector $w_2$ such that $R^Tw_1=w_2$. \item[MODE = 3 : ] the user should provide the vector $w_2$ with the smallest norm such that $Rw_2=w_1$. Equivalently, $w_2$ is the result of the application of the Moore-Penrose generalized inverse of R to $w_1$, i.e., \[ w_2 =R^T(RR^T)^{-1}w_1 \] It is worth noticing that the vector $w_1$ always belongs to the range of R. \end{description} \end{description} Finally, we note that the user is not forced to use this internal representation of the element functions. LSNNO will still work but at expense of extra storage and computational effort. This situation is nevertheless not recommended, especially when the gradients are not available analytically, or when the Hessian matrix is updated by a secant formula. It may, in these cases, result in much slower convergence, or possibly in no convergence at all. An example of the use of RANGE is shown in \S~5. \subsection{Functions XLOWER, XUPPER and RHS} The fact that the ith variable of the problem is bounded is signalled to the routine LSNNO by XST(i) being equal to 1. LSNNO will then need to know the actual values of the bounds on this variables. This information is provided by two user-supplied double precision functions XLOWER(i) and XUPPER(i). XLOWER(i) returns the value of the lower bound on the ith variable, while XUPPER(i) returns the value of the upper bound on this variable. Their specification is as follows : \begin{verbatim} DOUBLE PRECISION FUNCTION XLOWER(I) DOUBLE PRECISION FUNCTION XUPPER(I) \end{verbatim} The same technique is used to obtain the values of the right hand side vector b of the network constraints, i.e. the supply/demand values on the nodes. The double precision function RHS(i) returns this information for the ith node and is specified as follows : \begin{verbatim} DOUBLE PRECISION FUNCTION RHS(I) \end{verbatim} An example of the use of these functions is shown in \S~5. \subsection{Error messages} If IFLAG is positive on return from LSNNO, an error has been detected. The variable IFLAG is set within LSNNO to an appropriate error index and INFO to a (sometimes) meaningful value. A complete list of these errors, together with the associated value of IFLAG and INFO is given below. \begin{description} \item[IFLAG = 1 :] The total available work space provided in the vector WK, of length LWK, is too small. INFO = minimum necessary length of the vector WK. \item[IFLAG = 2 :] The total available work space provided in the vector IWK, of length LIWK, is too small. INFO = minimum necessary length of the vector IWK. \item[IFLAG = 3 :] The total available work space provided in the vector FUVAL, of length LFUVAL, is too small. INFO = minimum necessary length of the vector FUVAL. \item[IFLAG = 4 :] The bounds on the kth variable are inconsistent. INFO = k. \item[IFLAG = 5 :] The initial value of the fixed variable k is inconsistent with the bounds. INFO = k. \item[IFLAG = 6 :] The status of element function IEL corresponds to a facility that is not available at present (ELST(IEL)=-8). INFO = IEL. \item[IFLAG = 7 :] Maximum number of iterations executed before obtaining an initial feasible solution. Too much accuracy is required or the parameter MXITER in routine SIMPLX is to small. INFO = MXITER. \item[IFLAG = 8 :] The problem is infeasible. The set of feasible points is empty. \item[IFLAG = 9 :] The linear subproblem solved for searching an initial feasible solution is unbounded. \item[IFLAG = 10 :] The graph is not connected. INFO is meaningless. \item[IFLAG = 11 :] The number of independent superbasic sets is larger than the maximun predicted (INFO). Increase the parameter MAXINS in subroutine ANALYS. \item[IFLAG = 12 :] The algorithm is stopped because the difference between two successive function values along the current linesearch is insignificant compared to the noise on these function values. INFO is meaningless. \item[IFLAG = 13 :] More than MAXIT minor iterations have been performed. Recheck the datas of the problem and the derivatives, or increase MAXIT. INFO = MAXIT. \item[The following cases should never occur.] \item[IFLAG = 14 :] The routine SIMPLX is stopped because the linear function is increasing. \item[IFLAG = 15 :] Before a maximal basis spanning tree is searched, the total number of variables partitioned in the fixed, superbasic and nonbasic sets (INFO) differs from ARCS, i.e. the total number of variables of the problem. \item[IFLAG = 16 :] The number of basic variables (INFO) is incorrect. \item[IFLAG = 17 :] The total number of superbasics variables partitioned in the differents insets (INFO) differs from the total number of superbasics variables of the problem. \item[IFLAG = 18 :] The directional derivative at the beginning of a linesearch is non-negative. INFO is meaningless. \end{description} \subsection{Common} The following common blocks are used in LSNNO, they are not to be set by the user. The variables in common block PRBDIM contain the fundamental specifications of the problem. \begin{verbatim} COMMON / PRBDIM / ARCS, NODES, ELEM \end{verbatim} \begin{description} \item[ARCS ] is an INTEGER variable that is set to the number of oriented arcs in the network, it is also the number of variables of the problem. \item[NODES ] is an INTEGER variable that is set to the number of nodes in the network. \item[ELEM ] is an INTEGER variable that is set to the number of element functions. \end{description} The variables in common block PRBMCH contain the machine dependent constants. \begin{verbatim} COMMON / PRBMCH / EPSMCH, HUGE, TINY, TOL \end{verbatim} \begin{description} \item[EPSMCH ] is a DOUBLE PRECISION variable that is set to the flaoting point precision. \item[HUGE ] is a DOUBLE PRECISION variable that is set to the biggest floating number. \item[TINY ] is a DOUBLE PRECISION variable that is set to the smallest floating number. \item[TOL ] is a DOUBLE PRECISION variable that is set to the smallest difference between two numbers for considering them to be different. \end{description} The variables in common block PRBCPU provide CPU times information during the minimization. \begin{verbatim} COMMON / PRBCPU / XOCPU, TOTCPU, CGCPU, PRDCPU, LNSCPU, UPDCPU \end{verbatim} \begin{description} \item[X0CPU ] is a REAL variable which gives the CPU time spent in searching an initial feasible solution. \item[TOTCPU ] is a REAL variable which gives the total CPU time spent by the minimization process (Excluding X0CPU). \item[CGCPU ] is a REAL variable which gives the CPU time spent in the Conjugate Gradient procedure. \item[PRDCPU ] is a REAL variable which gives the CPU time spent by the matrix-vector products in C.G. \item[LNSCPU ] is a REAL variable which gives the CPU time spent in the linesearch. \item[UPDCPU ] is a REAL variable which gives the CPU time spent in updating the second derivatives. \end{description} The variables in common block PRBCST contain the constants used in the software, namely $0, 1, 2, 3, 0.5, 10^{-1}, 10^{-2}, 10^{-4}$. \begin{verbatim} COMMON / PRBCST / ZERO, ONE, TWO, THREE, HALF, TENM1, TENM2, TENM4 \end{verbatim} \section{General information} \begin{description} \item[Use of common : ] see \S 2.6. \item[Workspace : ] Provided by the user (see arguments IWK and WK). \item[Machine dependent routine : ] call SECOND(time). This subroutine call is machine dependent and uses only one REAL argument. It returns the CPU time used from the begin of the execution. \item[Input/Output : ] No input; output on device number IPDEVC. \item[Restrictions : ] NA $> 0$, NN $> 0$, NEL $> 0$. \item[Portability : ] Fortran 77. \end{description} \section{Method} The basic method that is implemented by the routine LSNNO is described in detail by Toint and Tuyttens (1989). The software itself is described in Toint and Tuyttens (1992). The concept of partial separability was first suggested by Griewank and Toint (1982). The extensions of this concept to network constraints is new. The method used is iterative. At each iteration, a quadratic model of the objective function is constructed. The step is determined by using a truncated conjugate gradient algorithm, followed by a linesearch. The strategy for treating bound constraints is based on the usual projection device. This technique is similar to that used by Bertsekas (1982). The method combines specialized data structure (see Bradley, Brown and Graves 1977), MINOS-like partitioning of the variables (see Murthag and Saunders 1978) and the concept of independent superbasic sets (Escudero 1986) for handling the network constraints. \begin{description} \item[References] \item[] D.P. Bertsekas, "Projected Newton methods for optimization problems with simple constraints", SIAM Journal of Control and Optimization, vol.20, pp.221--246, 1980. \item[] G.H. Bradley, G.G. Brown and G.W. Graves, "Design and implementation of large scale transshipment algorithms", Management Science, vol.24, pp.1--34, 1977. \item[] L.F. Escudero, "Performance evaluation of independent superbasic sets on nonlinear replicated networks", European Journal of Operational Research, vol.23, pp.343--355,1986. \item[] A. Griewank and Ph.L. Toint, "Partitioned variable metric updates for large structured optimization problems", Numerische Mathematik, vol.39, pp.119--137, 1982. \item[] B. Murtagh and M. Saunders, "Large scale linearly constrained optimization", Mathematical Programming, vol.14,pp.41--72,1978. \item[] Ph.L. Toint and D. Tuyttens, "On large scale nonlinear network optimization", Mathematical Programming, Series B, 48(1), pp. 125--159, 1990. \item[] Ph.L. Toint and D. Tuyttens, "LSNNO: a Fortran subroutine for solving large-scale nonlinear network optimization problems'', ACM Transactions on Mathematical Software, 18(3), pp. 308--328, 1992. \end{description} \section{Example of use} We now consider the small example problem \[ \begin{array}{lc} {\rm minimize} & x_2e^{x_1+x_3}+(x_3x_4)^2+(x_3-x_5)^2 \\ \mbox{subject to} & \left\{ \begin{array}{cccccccccccr} & x_1 & + & x_2 & & & & & & & = & 10 \\ - & x_1 & & & - & x_3 & + & x_4 & & & = & 0 \\ & & - & x_2 & + & x_3 & & & + & x_5 & = & 0 \\ & & & & & & - & x_4 & - & x_5 & = & -10 \end{array} \right. \\ & 2 \leq x_1 \leq 4 \\ & 6 \leq x_2 \leq 8 \\ & 0 \leq x_3 \leq 5 \end{array} \] The directed network corresponding to the four flow conservation equations is described in figure 2. \setlength{\unitlength}{0.002 cm} \begin{figure}[htbp] \begin{center} \begin{picture}(6000,1600)(10,10) % % make the circles for the nodes % \put(1000,800){\circle{250}} \put(3000,390){\circle{250}} \put(3000,1300){\circle{250}} \put(5000,800){\circle{250}} % % number of the nodes % \put(1000,800){\makebox(0,0){1}} \put(3000,390){\makebox(0,0){3}} \put(3000,1300){\makebox(0,0){2}} \put(5000,800){\makebox(0,0){4}} % % supply/demand % \put(500,800){\vector(1,0){300}} \put(400,800){\makebox(0,0)[r]{+10}} \put(5200,800){\vector(1,0){300}} \put(5600,800){\makebox(0,0)[l]{-10}} % % oriented arcs % \put(3000,525){\vector(0,1){600}} \put(1125,825){\vector(4,1){1750}} \put(1125,775){\vector(4,-1){1750}} \put(3125,1275){\vector(4,-1){1750}} \put(3125,315){\vector(4,1){1750}} % % number of arcs % \put(3100,800){\makebox(0,0){3}} \put(2000,1200){\makebox(0,0){1}} \put(2000,400){\makebox(0,0){2}} \put(4000,1200){\makebox(0,0){4}} \put(4000,400){\makebox(0,0){5}} \end{picture} \caption{Network diagram for example test} \end{center} \end{figure} The network has 4 nodes and 5 directed arcs. We can then set NA = 5 and NN = 4. The network is represented by the two vectors FR (the origins of the arcs), TO (the endpoints of the arcs) and the function RHS (the supply/demand on the nodes), as follows : \[ \begin{array}{cccccc} I & 1 & 2 & 3 & 4 & 5 \\ FR(I) & 1 & 1 & 3 & 2 & 3 \\ TO(I) & 2 & 3 & 2 & 4 & 4 \end{array} \] and \def\baselinestrech{0.8} {\small {\tt \begin{verbatim} DOUBLE PRECISION FUNCTION RHS(I) INTEGER I IF( I.EQ.1 ) THEN RHS = 10.0D0 ELSE IF( I.EQ.4 ) THEN RHS = -10.D0 ELSE RHS = 0.0D0 ENDIF RETURN END \end{verbatim} } } \def\baselinestrech{1.2} We note that $\sum_{i=1}^{NN} b(i) = 0$. The variables 1, 2 and 3 are bounded, but the variables 4 and 5 are unbounded. We build the vector XST as follows : \[ \begin{array}{cccccc} I & 1 & 2 & 3 & 4 & 5 \\ XST(I) & 1 & 1 & 1 & -1 & -1 \end{array} \] The considered function depends on the 5 variables and we consider the following partition into three elements. We set NEL = 3. \[ \begin{array}{ccccc} x_2e^{x_1+x_3} & + & (x_3x_4)^2 & + & (x_3-x_5)^2 \\ \uparrow & & \uparrow & & \uparrow \\ \tim{element 1 } & & \tim{element 2} & & \tim{element 3} \end{array} \] The first element involves variables 1, 2 and 3, the second element involves variables 3 and 4, and the third element involves the variables 3 and 5. We now build the vector ELVAR corresponding to the problem : \[ \begin{array}{cccccccc} I & 1 & 2 & 3 & 4 & 5 & 6 & 7 \\ ELVAR(I) & 1 & 2 & 3 & 3 & 4 & 3 & 5 \\ & \uparrow & & & \uparrow & & \uparrow & \\ & el.1 & & & el.2 & & el.3 & \end{array} \] In this vector, we now locate the position of the first variable of each element, and build the vector ELPTR as follows : \[ \begin{array}{ccccc} I & 1 & 2 & 3 & 4 \\ ELPTR(I) & 1 & 4 & 6 & 8 \end{array} \] We observe that ELPTR(NEL+1) is set to the total length of ELVAR plus 1, as required. The analytical gradients and Hessians are available for the three element functions, and we build the vector ELST as follows : \[ \begin{array}{cccc} I & 1 & 2 & 3 \\ ELST(I) & 1 & 1 & 1 \end{array} \] We observe that element 1 may be reduced to dependence on two variables with the matrix \[ R_1 = \left( \begin{array}{ccc} 1 & 0 & 1 \\ 0 & 1 & 0 \end{array} \right) \] and that element 3 may be reduced to dependence on one variable with the matrix \[ R_3 = \left( \begin{array}{cc} 1 & -1 \end{array} \right) \] We now give an example of a subroutine that would describe the function of our simple problem~: \def\baselinestrech{0.8} {\small {\tt \begin{verbatim} SUBROUTINE ELFNCT( FUVAL, X, NA, NEL, FCALC, LFC, + DER, ELVAR, ELPTR, GPTR, HPTR ) INTEGER NA, NEL, LFC, FCALC(1), DER, + ELVAR(1), ELPTR(1), GPTR(1), HPTR(1) DOUBLE PRECISION FUVAL(1), X(1) INTEGER I, IEL, IGSTRT, IHSTRT DOUBLE PRECISION TEMP1, TEMP2, TEMP3 DO 10 I = 1 , LFC IEL = FCALC(I) IGSTRT = GPTR(IEL) - 1 IHSTRT = HPTR(IEL) - 1 TEMP1 = DEXP(X(1) + X(3)) TEMP2 = X(3) * X(4) TEMP3 = X(3) - X(5) IF( DER.EQ.0 ) THEN IF( IEL.EQ.1 ) THEN FUVAL( IEL ) = X(2) * TEMP1 ELSE IF( IEL.EQ.2 ) THEN FUVAL( IEL ) = TEMP2**2.0D0 ELSE FUVAL( IEL ) = TEMP3**2.0D0 ENDIF ELSE IF( DER.EQ.1 ) THEN IF( IEL.EQ.1 ) THEN FUVAL( IGSTRT + 1 ) = TEMP1 * X(2) FUVAL( IGSTRT + 2 ) = TEMP1 ELSE IF( IEL.EQ.2 ) THEN FUVAL( IGSTRT + 1 ) = 2.0D0 * X(4) * TEMP2 FUVAL( IGSTRT + 2 ) = 2.0D0 * X(3) * TEMP2 ELSE FUVAL( IGSTRT + 1 ) = 2.0D0 * TEMP3 ENDIF ELSE IF( IEL.EQ.1 ) THEN FUVAL( IHSTRT + 1 ) = TEMP1 * X(2) FUVAL( IHSTRT + 2 ) = TEMP1 FUVAL( IHSTRT + 3 ) = 0.0D0 ELSE IF( IEL.EQ.2 ) THEN FUVAL( IHSTRT + 1 ) = 2.0D0 * X(4) * X(4) FUVAL( IHSTRT + 2 ) = 4.0D0 * TEMP2 FUVAL( IHSTRT + 3 ) = 2.0D0 * X(3) * X(3) ELSE FUVAL( IHSTRT + 1 ) = 2.0D0 ENDIF ENDIF 10 CONTINUE RETURN END \end{verbatim} } } \def\baselinestrech{1.2} Note that the derivatives are taken with respect to the internal variables, this means that we must transform from elemental to internal variables before evaluating the derivatives. We now write down the routine RANGE associated with our example. \def\baselinestrech{0.8} {\small {\tt \begin{verbatim} SUBROUTINE RANGE( IEL, ELDIM, INTDIM, W1, W2, MODE ) INTEGER IEL, ELDIM, INTDIM, MODE DOUBLE PRECISION W1(1), W2(1) IF( MODE.EQ.0 ) THEN IF( IEL.EQ.1 ) THEN INTDIM = 2 ELSE IF( IEL.EQ.2 ) THEN INTDIM = ELDIM ELSE INTDIM = 1 ENDIF RETURN ENDIF GOTO( 1, 2, 3 ) MODE C C R*w1=w2 C 1 CONTINUE IF( IEL.EQ.1 ) THEN W2(1) = W1(1) + W1(3) W2(2) = W1(2) ELSE W2(1) = W1(1) - W1(2) ENDIF RETURN C C R'*w1 = w2 C 2 CONTINUE IF( IEL.EQ.1 ) THEN W2(1) = W1(1) W2(2) = W1(2) W2(3) = W1(1) ELSE W2(1) = W1(1) W2(1) = -W1(1) ENDIF RETURN C C R*w2=w1 C 3 CONTINUE IF( IEL.EQ.1 ) THEN W2(1) = 0.5D0 * W1(1) W2(2) = W1(2) W2(3) = 0.5D0 * W1(1) ELSE W2(1) = 0.5D0 * W1(1) W2(2) = - W2(1) ENDIF RETURN C END \end{verbatim} } } \def\baselinestrech{1.2} As one can see, to write the routine RANGE is not difficult. The last thing we have to set up for our example is the bound specification. In the case where variables are unbounded, we have to build an artificial lower bound and an artificial upper bound. This gives the following functions : \def\baselinestrech{0.8} {\small {\tt \begin{verbatim} DOUBLE PRECISION FUNCTION XLOWER(I) INTEGER I IF( I.EQ.1 ) THEN XLOWER = 2.0D0 ELSE IF( I.EQ.2 ) THEN XLOWER = 6.0D0 ELSE IF( I.EQ.3 ) THEN XLOWER = 0.0D0 ELSE XLOWER = -1.0D20 ENDIF RETURN END DOUBLE PRECISION FUNCTION XUPPER(I) INTEGER I IF( I.EQ.1 ) THEN XUPPER = 4.0D0 ELSE IF( I.EQ.2 ) THEN XUPPER = 8.0D0 ELSE IF( I.EQ.3 ) THEN XUPPER = 5.0D0 ELSE XUPPER = 1.0D20 ENDIF RETURN END \end{verbatim} } } \def\baselinestrech{1.2} This completes the supply information to LSNNO. The problem may now be solved using the following program. We start with the initial point equal to the origin and terminate when the infinity norm of the reduced gradient is smaller than $10^{-5}$, using an unpreconditioned conjugate gradient scheme. \def\baselinestrech{0.8} {\small {\tt \begin{verbatim} PROGRAM EXMAIN INTEGER NAMAX, NELMAX, NA, NN, NEL, INFORM, + LELVAR, LELPTR, LFUVAL, LGPTR, + LHPTR, LFCALC, LFC, LIWK, LWK, MAXIT, + IPDEVC, WHAT, FREQ, INFO, IFLAG PARAMETER ( NAMAX = 5 , NELMAX = 3 ) PARAMETER ( LIWK = 10000 , LWK = 10000 ) PARAMETER ( LELVAR = 7 , LELPTR = NELMAX + 1 ) PARAMETER ( LGPTR = NELMAX + 1 , LHPTR = NELMAX + 1 ) PARAMETER ( LFUVAL = 10000 , LFCALC = NELMAX ) PARAMETER ( MAXIT = 100 , IPDEVC = 6 ) PARAMETER ( FREQ = 1 , WHAT = 2 ) INTEGER ELVAR(LELVAR), ELPTR(LELPTR), ELST(NELMAX), + XST(NAMAX), FR(NAMAX), TO(NAMAX), IT(3), + GPTR(LGPTR), HPTR(LHPTR), IWK(LIWK), + FCALC(LFCALC), IFC(3), I LOGICAL PREC, NIPST DOUBLE PRECISION X(NAMAX), FX, FUVAL(LFUVAL), + EPSF, WK(LWK) DATA NA / 5 /, NN / 4 /, NEL / 3 / DATA FR / 1, 1, 3, 2, 3 / DATA TO / 2, 3, 2, 4, 4 / DATA ELVAR / 1, 2, 3, 3, 4, 3, 5 / DATA ELPTR / 1, 4, 6, 8 / DATA ELST / 1, 1, 1 / DATA X / 0.0D0, 0.0D0, 0.0D0, 0.0D0, 0.0D0 / DATA XST / 1, 1, 1, -1, -1 / DATA EPSF / 1.D-5 / DATA PREC / .FALSE. /, NIPST / .FALSE. / INFORM = 0 * * SOLVING THE PROBLEN WITH LSNNO. * 100 CONTINUE CALL LSNNO( NA, NN, NEL, ELVAR, ELPTR, ELST, FR, TO, + X, XST, FX, FUVAL, LFUVAL, GPTR, HPTR, + INFORM, FCALC, LFC, IFC, IT, MAXIT, EPSF, + PREC, NIPST, IWK, LIWK, WK, LWK, IPDEVC, + WHAT, FREQ, INFO, IFLAG ) * IF( INFORM.GT.0 ) THEN * IF( INFORM.EQ.1 ) THEN CALL ELFNCT( FUVAL, X, NA, NEL, FCALC, LFC, 0, + ELVAR, ELPTR, GPTR, HPTR ) ENDIF * IF( INFORM.EQ.2 .OR. INFORM.EQ.4 ) THEN CALL ELFNCT( FUVAL, X, NA, NEL, FCALC, LFC, 1, + ELVAR, ELPTR, GPTR, HPTR ) ENDIF * IF( INFORM.EQ.3 .OR. INFORM.EQ.4 ) THEN CALL ELFNCT( FUVAL, X, NA, NEL, FCALC, LFC, 2, + ELVAR, ELPTR, GPTR, HPTR ) ENDIF * GO TO 100 * ELSE WRITE(IPDEVC,1000) DO 200 I = 1 , NA WRITE(IPDEVC,1010) I,X(I) 200 CONTINUE ENDIF * STOP 1000 FORMAT(/, ' VARIABLE NAME VALUE ', + /, ' ------------- ----- ') 1010 FORMAT(7X,'X',I1,4X,D12.5) END \end{verbatim} } } \def\baselinestrech{1.2} This produces the following output. \def\baselinestrech{0.8} {\small {\tt \begin{verbatim} ********** PROBLEM SPECIFICATIONS ********** NUMBER OF ARCS 5 NUMBER OF NODES 4 NUMBER OF ELEMENTS 3 INTEGER SPACE USED 180 DOUBLE SPACE USED 56 MACHINE PRECISION 0.22204D-15 ********** FLOW FEASIBILITY CHECKING ********** * FLOW IS INFEASIBLE * NORM OF INFEASIBILITY = 0.10000D+02 SUM OF THE RHS VALUES = 0.00000D+00 * FLOW IS FEASIBLE * NORM OF INFEASIBILITY = 0.00000D+00 SUM OF THE RHS VALUES = 0.00000D+00 >>>>>>>>>>>>>>>>>>>>>>>> LSNNO - OPTIMIZATION PHASE <<<<<<<<<<<<<<<<<<<<<<<< BASIS NOT FULLY FREE - 2 NONBASIC VARIABLES. SUPERBASIC SET IS EMPTY. >> MAJOR ITERATION 1 - DEACTIVATION OF 2 VARIABLE(S) - 1 IN-SET(S) << min it piv fx rgmax rgnor #f #g #h cgit #su ------------------------------------------------------------------------------ - - - 0.50660D+05 0.418D+05 0.531D+05 1.0 1.0 1.0 - 2 1 1 - 0.22222D+05 0.185D+05 0.237D+05 2.0 2.0 2.0 2 2 2 2 - 0.95944D+04 0.813D+04 0.104D+05 3.0 3.0 3.0 4 2 3 3 - 0.40516D+04 0.351D+04 0.446D+04 4.0 4.0 4.0 6 2 >> MAJOR ITERATION 2 - DEACTIVATION OF 2 VARIABLE(S) - 1 IN-SET(S) << min it piv fx rgmax rgnor #f #g #h cgit #su ------------------------------------------------------------------------------ - - - 0.40516D+04 0.351D+04 0.446D+04 4.0 4.0 4.0 - 2 4 1 - 0.16562D+04 0.149D+04 0.187D+04 5.0 5.0 5.0 8 2 5 2 - 0.36359D+03 0.000D+00 0.000D+00 6.0 6.0 6.0 15 2 >> MAJOR ITERATION 3 - DEACTIVATION OF 1 VARIABLE(S) - 1 IN-SET(S) << min it piv fx rgmax rgnor #f #g #h cgit #su ------------------------------------------------------------------------------ - - - 0.36359D+03 0.218D+03 0.218D+03 6.0 6.0 6.0 - 1 6 1 - 0.18960D+03 0.100D+03 0.100D+03 7.0 7.0 7.0 16 1 7 2 - 0.12311D+03 0.000D+00 0.000D+00 8.0 8.0 8.0 17 1 >>>>>>>>>>>>>>>>>>>>>>>> LSNNO - END OF OPTIMIZATION <<<<<<<<<<<<<<<<<<<<<<<< * OPTIMAL SOLUTION FOUND * * FLOW IS FEASIBLE * NORM OF INFEASIBILITY = 0.88818D-15 SUM OF THE RHS VALUES = 0.00000D+00 * STATISTICS : -------------- Strict complementarity is not verified. Problem specifications number of arcs : 5 number of nodes : 4 number of elements : 3 Function value at starting point : 0.506595D+05 at optimum point : 0.123112D+03 Iterations counts Major iterations : 3 minor iterations : 7 pivoting iterations : 0 C.G. iterations : 17 average dimension : 1.7 function calls : 8.0 gradient calls : 8.0 Hessian calls : 8.0 Variables status at optimum superbasics at optimum : 0 actives at optimum : 3 de-activated variables : 3 activated variables : 5 final tolerance required : 0.100000D-04 reduced gradient inf. norm : 0.000000D+00 reduced gradient l2 norm : 0.000000D+00 C.P.U. time for initial solution : 0.0 optimization : 0.0 congugate gradient : 0.0 matrix vector products : 0.0 linesearch : 0.0 Hessian updates : 0.0 ********** END OF STATISTICS ********** VARIABLE NAME VALUE ------------- ----- X1 0.20000D+01 X2 0.80000D+01 X3 0.00000D+00 X4 0.20000D+01 X5 0.80000D+01 \end{verbatim} } } \def\baselinestrech{1.2} \end{document} .