# dsplp.f

          SUBROUTINE DSPLP (DUSRMT, MRELAS, NVARS, COSTS, PRGOPT, DATTRV,
         +   BL, BU, IND, INFO, PRIMAL, DUALS, IBASIS, WORK, LW, IWORK, LIW)
    C***BEGIN PROLOGUE  DSPLP
    C***PURPOSE  Solve linear programming problems involving at
    C            most a few thousand constraints and variables.
    C            Takes advantage of sparsity in the constraint matrix.
    C***LIBRARY   SLATEC
    C***CATEGORY  G2A2
    C***TYPE      DOUBLE PRECISION (SPLP-S, DSPLP-D)
    C***KEYWORDS  LINEAR CONSTRAINTS, LINEAR OPTIMIZATION,
    C             LINEAR PROGRAMMING, LP, SPARSE CONSTRAINTS
    C***AUTHOR  Hanson, R. J., (SNLA)
    C           Hiebert, K. L., (SNLA)
    C***DESCRIPTION
    C
    C     These are the short usage instructions; for details about
    C     other features, options and methods for defining the matrix
    C     A, see the extended usage instructions which are contained in
    C     the Long Description section below.
    C
    C   |------------|
    C   |Introduction|
    C   |------------|
    C     The subprogram DSPLP( ) solves a linear optimization problem.
    C     The problem statement is as follows
    C
    C                         minimize (transpose of costs)*x
    C                         subject to A*x=w.
    C
    C     The entries of the unknowns x and w may have simple lower or
    C     upper bounds (or both), or be free to take on any value.  By
    C     setting the bounds for x and w, the user is imposing the con-
    C     straints of the problem.  The matrix A has MRELAS rows and
    C     NVARS columns.  The vectors costs, x, and w respectively
    C     have NVARS, NVARS, and MRELAS number of entries.
    C
    C     The input for the problem includes the problem dimensions,
    C     MRELAS and NVARS, the array COSTS(*), data for the matrix
    C     A, and the bound information for the unknowns x and w, BL(*),
    C     BU(*), and IND(*).  Only the nonzero entries of the matrix A
    C     are passed to DSPLP( ).
    C
    C     The output from the problem (when output flag INFO=1) includes
    C     optimal values for x and w in PRIMAL(*), optimal values for
    C     dual variables of the equations A*x=w and the simple bounds
    C     on x in  DUALS(*), and the indices of the basic columns,
    C     IBASIS(*).
    C
    C  |------------------------------|
    C  |Fortran Declarations Required:|
    C  |------------------------------|
    C
    C     DIMENSION COSTS(NVARS),PRGOPT(*),DATTRV(*),
    C    *BL(NVARS+MRELAS),BU(NVARS+MRELAS),IND(NVARS+MRELAS),
    C    *PRIMAL(NVARS+MRELAS),DUALS(MRELAS+NVARS),IBASIS(NVARS+MRELAS),
    C    *WORK(LW),IWORK(LIW)
    C
    C     EXTERNAL DUSRMT
    C
    C     The dimensions of PRGOPT(*) and DATTRV(*) must be at least 1.
    C     The exact lengths will be determined by user-required options and
    C     data transferred to the subprogram DUSRMT( ).
    C
    C     The values of LW and LIW, the lengths of the arrays WORK(*)
    C     and IWORK(*), must satisfy the inequalities
    C
    C               LW .GE. 4*NVARS+ 8*MRELAS+LAMAT+  LBM
    C               LIW.GE.   NVARS+11*MRELAS+LAMAT+2*LBM
    C
    C     It is an error if they do not both satisfy these inequalities.
    C     (The subprogram will inform the user of the required lengths
    C     if either LW or LIW is wrong.)  The values of LAMAT and LBM
    C     nominally are
    C
    C               LAMAT=4*NVARS+7
    C     and       LBM  =8*MRELAS
    C
    C     LAMAT determines the length of the sparse matrix storage area.
    C     The value of LBM determines the amount of storage available
    C     to decompose and update the active basis matrix.
    C
    C  |------|
    C  |Input:|
    C  |------|
    C
    C     MRELAS,NVARS
    C     ------------
    C     These parameters are respectively the number of constraints (the
    C     linear relations A*x=w that the unknowns x and w are to satisfy)
    C     and the number of entries in the vector x.  Both must be .GE. 1.
    C     Other values are errors.
    C
    C     COSTS(*)
    C     --------
    C     The NVARS entries of this array are the coefficients of the
    C     linear objective function.  The value COSTS(J) is the
    C     multiplier for variable J of the unknown vector x.  Each
    C     entry of this array must be defined.
    C
    C     DUSRMT
    C     ------
    C     This is the name of a specific subprogram in the DSPLP( ) package
    C     used to define the matrix A.  In this usage mode of DSPLP( )
    C     the user places the nonzero entries of A in the
    C     array DATTRV(*) as given in the description of that parameter.
    C     The name DUSRMT must appear in a Fortran EXTERNAL statement.
    C
    C     DATTRV(*)
    C     ---------
    C     The array DATTRV(*) contains data for the matrix A as follows:
    C     Each column (numbered J) requires (floating point) data con-
    C     sisting of the value (-J) followed by pairs of values.  Each pair
    C     consists of the row index immediately followed by the value
    C     of the matrix at that entry.  A value of J=0 signals that there
    C     are no more columns.  The required length of
    C     DATTRV(*) is 2*no. of nonzeros + NVARS + 1.
    C
    C     BL(*),BU(*),IND(*)
    C     ------------------
    C     The values of IND(*) are input parameters that define
    C     the form of the bounds for the unknowns x and w.  The values for
    C     the bounds are found in the arrays BL(*) and BU(*) as follows.
    C
    C     For values of J between 1 and NVARS,
    C          if IND(J)=1, then X(J) .GE. BL(J); BU(J) is not used.
    C          if IND(J)=2, then X(J) .LE. BU(J); BL(J) is not used.
    C          if IND(J)=3, then BL(J) .LE. X(J) .LE. BU(J),(BL(J)=BU(J) ok)
    C          if IND(J)=4, then X(J) is free to have any value,
    C          and BL(J), BU(J) are not used.
    C
    C     For values of I between NVARS+1 and NVARS+MRELAS,
    C          if IND(I)=1, then W(I-NVARS) .GE. BL(I); BU(I) is not used.
    C          if IND(I)=2, then W(I-NVARS) .LE. BU(I); BL(I) is not used.
    C          if IND(I)=3, then BL(I) .LE. W(I-NVARS) .LE. BU(I),
    C          (BL(I)=BU(I) is ok).
    C          if IND(I)=4, then W(I-NVARS) is free to have any value,
    C          and BL(I), BU(I) are not used.
    C
    C     A value of IND(*) not equal to 1,2,3 or 4 is an error.  When
    C     IND(I)=3, BL(I) must be .LE. BU(I).  The condition BL(I).GT.
    C     BU(I) indicates infeasibility and is an error.
    C
    C     PRGOPT(*)
    C     ---------
    C     This array is used to redefine various parameters within DSPLP( ).
    C     Frequently, perhaps most of the time, a user will be satisfied
    C     and obtain the solutions with no changes to any of these
    C     parameters.  To try this, simply set PRGOPT(1)=1.D0.
    C
    C     For users with more sophisticated needs, DSPLP( ) provides several
    C     options that may be used to take advantage of more detailed
    C     knowledge of the problem or satisfy other utilitarian needs.
    C     The complete description of how to use this option array to
    C     utilize additional subprogram features is found under the
    C     heading  of DSPLP( ) Subprogram Options in the Extended
    C     Usage Instructions.
    C
    C     Briefly, the user should note the following value of the parameter
    C     KEY and the corresponding task or feature desired before turning
    C     to that document.
    C
    C     Value     Brief Statement of Purpose for Option
    C     of KEY
    C     ------    -------------------------------------
    C     50        Change from a minimization problem to a
    C               maximization problem.
    C     51        Change the amount of printed output.
    C               Normally, no printed output is obtained.
    C     52        Redefine the line length and precision used
    C               for the printed output.
    C     53        Redefine the values of LAMAT and LBM that
    C               were discussed above under the heading
    C               Fortran Declarations Required.
    C     54        Redefine the unit number where pages of the sparse
    C               data matrix A are stored.  Normally, the unit
    C               number is 1.
    C     55        A computation, partially completed, is
    C               being continued.  Read the up-to-date
    C               partial results from unit number 2.
    C     56        Redefine the unit number where the partial results
    C               are stored.  Normally, the unit number is 2.
    C     57        Save partial results on unit 2 either after
    C               maximum iterations or at the optimum.
    C     58        Redefine the value for the maximum number of
    C               iterations.  Normally, the maximum number of
    C               iterations is 3*(NVARS+MRELAS).
    C     59        Provide DSPLP( ) with a starting (feasible)
    C               nonsingular basis.  Normally, DSPLP( ) starts
    C               with the identity matrix columns corresponding
    C               to the vector w.
    C     60        The user has provided scale factors for the
    C               columns of A.  Normally, DSPLP( ) computes scale
    C               factors that are the reciprocals of the max. norm
    C               of each column.
    C     61        The user has provided a scale factor
    C               for the vector costs.  Normally, DSPLP( ) computes
    C               a scale factor equal to the reciprocal of the
    C               max. norm of the vector costs after the column
    C               scaling for the data matrix has been applied.
    C     62        Size parameters, namely the smallest and
    C               largest magnitudes of nonzero entries in
    C               the matrix A, are provided.  Values noted
    C               outside this range are to be considered errors.
    C     63        Redefine the tolerance required in
    C               evaluating residuals for feasibility.
    C               Normally, this value is set to RELPR,
    C               where RELPR = relative precision of the arithmetic.
    C     64        Change the criterion for bringing new variables
    C               into the basis from the steepest edge (best
    C               local move) to the minimum reduced cost.
    C     65        Redefine the value for the number of iterations
    C               between recalculating the error in the primal
    C               solution.  Normally, this value is equal to ten.
    C     66        Perform "partial pricing" on variable selection.
    C               Redefine the value for the number of negative
    C               reduced costs to compute (at most) when finding
    C               a variable to enter the basis.  Normally this
    C               value is set to NVARS.  This implies that no
    C               "partial pricing" is used.
    C     67        Adjust the tuning factor (normally one) to apply
    C               to the primal and dual error estimates.
    C     68        Pass  information to the  subprogram  DFULMT(),
    C               provided with the DSPLP() package, so that a Fortran
    C               two-dimensional array can be used as the argument
    C               DATTRV(*).
    C     69        Pass an absolute tolerance to use for the feasibility
    C               test when the usual relative error test indicates
    C               infeasibility.  The nominal value of this tolerance,
    C               TOLABS, is zero.
    C
    C
    C  |---------------|
    C  |Working Arrays:|
    C  |---------------|
    C
    C     WORK(*),LW,
    C     IWORK(*),LIW
    C     ------------
    C     The arrays WORK(*) and IWORK(*) are respectively floating point
    C     and type INTEGER working arrays for DSPLP( ) and its
    C     subprograms.  The lengths of these arrays are respectively
    C     LW and LIW.  These parameters must satisfy the inequalities
    C     noted above under the heading "Fortran Declarations Required:"
    C     It is an error if either value is too small.
    C
    C  |----------------------------|
    C  |Input/Output files required:|
    C  |----------------------------|
    C
    C     Fortran unit 1 is used by DSPLP( ) to store the sparse matrix A
    C     out of high-speed memory.  A crude
    C     upper bound for the amount of information written on unit 1
    C     is 6*nz, where nz is the number of nonzero entries in A.
    C
    C  |-------|
    C  |Output:|
    C  |-------|
    C
    C     INFO,PRIMAL(*),DUALS(*)
    C     -----------------------
    C     The integer flag INFO indicates why DSPLP( ) has returned to the
    C     user.  If INFO=1 the solution has been computed.  In this case
    C     X(J)=PRIMAL(J) and W(I)=PRIMAL(I+NVARS).  The dual variables
    C     for the equations A*x=w are in the array DUALS(I)=dual for
    C     equation number I.  The dual value for the component X(J) that
    C     has an upper or lower bound (or both) is returned in
    C     DUALS(J+MRELAS).  The only other values for INFO are .LT. 0.
    C     The meaning of these values can be found by reading
    C     the diagnostic message in the output file, or by looking for
    C     error number = (-INFO) in the Extended Usage Instructions
    C     under the heading:
    C
    C          List of DSPLP( ) Error and Diagnostic Messages.
    C
    C     BL(*),BU(*),IND(*)
    C     ------------------
    C     These arrays are output parameters only under the (unusual)
    C     circumstances where the stated problem is infeasible, has an
    C     unbounded optimum value, or both.  These respective conditions
    C     correspond to INFO=-1,-2 or -3.    See the Extended
    C     Usage Instructions for further details.
    C
    C     IBASIS(I),I=1,...,MRELAS
    C     ------------------------
    C     This array contains the indices of the variables that are
    C     in the active basis set at the solution (INFO=1).  A value
    C     of IBASIS(I) between 1 and NVARS corresponds to the variable
    C     X(IBASIS(I)).  A value of IBASIS(I) between NVARS+1 and NVARS+
    C     MRELAS corresponds to the variable W(IBASIS(I)-NVARS).
    C
    C *Long Description:
    C
    C     SUBROUTINE DSPLP(DUSRMT,MRELAS,NVARS,COSTS,PRGOPT,DATTRV,
    C    *           BL,BU,IND,INFO,PRIMAL,DUALS,IBASIS,WORK,LW,IWORK,LIW)
    C
    C     |------------|
    C     |Introduction|
    C     |------------|
    C     The subprogram DSPLP( ) solves a linear optimization problem.
    C     The problem statement is as follows
    C
    C                         minimize (transpose of costs)*x
    C                         subject to A*x=w.
    C
    C     The entries of the unknowns x and w may have simple lower or
    C     upper bounds (or both), or be free to take on any value.  By
    C     setting the bounds for x and w, the user is imposing the con-
    C     straints of the problem.
    C
    C     (The problem may also be stated as a maximization
    C     problem.  This is done by means of input in the option array
    C     PRGOPT(*).)  The matrix A has MRELAS rows and NVARS columns.  The
    C     vectors costs, x, and w respectively have NVARS, NVARS, and
    C     MRELAS number of entries.
    C
    C     The input for the problem includes the problem dimensions,
    C     MRELAS and NVARS, the array COSTS(*), data for the matrix
    C     A, and the bound information for the unknowns x and w, BL(*),
    C     BU(*), and IND(*).
    C
    C     The output from the problem (when output flag INFO=1) includes
    C     optimal values for x and w in PRIMAL(*), optimal values for
    C     dual variables of the equations A*x=w and the simple bounds
    C     on x in  DUALS(*), and the indices of the basic columns in
    C     IBASIS(*).
    C
    C  |------------------------------|
    C  |Fortran Declarations Required:|
    C  |------------------------------|
    C
    C     DIMENSION COSTS(NVARS),PRGOPT(*),DATTRV(*),
    C    *BL(NVARS+MRELAS),BU(NVARS+MRELAS),IND(NVARS+MRELAS),
    C    *PRIMAL(NVARS+MRELAS),DUALS(MRELAS+NVARS),IBASIS(NVARS+MRELAS),
    C    *WORK(LW),IWORK(LIW)
    C
    C     EXTERNAL DUSRMT (or 'NAME', if user provides the subprogram)
    C
    C     The dimensions of PRGOPT(*) and DATTRV(*) must be at least 1.
    C     The exact lengths will be determined by user-required options and
    C     data transferred to the subprogram DUSRMT( ) ( or 'NAME').
    C
    C     The values of LW and LIW, the lengths of the arrays WORK(*)
    C     and IWORK(*), must satisfy the inequalities
    C
    C               LW .GE. 4*NVARS+ 8*MRELAS+LAMAT+  LBM
    C               LIW.GE.   NVARS+11*MRELAS+LAMAT+2*LBM
    C
    C     It is an error if they do not both satisfy these inequalities.
    C     (The subprogram will inform the user of the required lengths
    C     if either LW or LIW is wrong.)  The values of LAMAT and LBM
    C     nominally are
    C
    C               LAMAT=4*NVARS+7
    C     and       LBM  =8*MRELAS
    C
    C     These values will be as shown unless the user changes them by
    C     means of input in the option array PRGOPT(*).  The value of LAMAT
    C     determines the length of the sparse matrix "staging" area.
    C     For reasons of efficiency the user may want to increase the value
    C     of LAMAT.  The value of LBM determines the amount of storage
    C     available to decompose and update the active basis matrix.
    C     Due to exhausting the working space because of fill-in,
    C     it may be necessary for the user to increase the value of LBM.
    C     (If this situation occurs an informative diagnostic is printed
    C     and a value of INFO=-28 is obtained as an output parameter.)
    C
    C  |------|
    C  |Input:|
    C  |------|
    C
    C     MRELAS,NVARS
    C     ------------
    C     These parameters are respectively the number of constraints (the
    C     linear relations A*x=w that the unknowns x and w are to satisfy)
    C     and the number of entries in the vector x.  Both must be .GE. 1.
    C     Other values are errors.
    C
    C     COSTS(*)
    C     --------
    C     The NVARS entries of this array are the coefficients of the
    C     linear objective function.  The value COSTS(J) is the
    C     multiplier for variable J of the unknown vector x.  Each
    C     entry of this array must be defined.  This array can be changed
    C     by the user between restarts.  See options with KEY=55,57 for
    C     details of checkpointing and restarting.
    C
    C     DUSRMT
    C     ------
    C     This is the name of a specific subprogram in the DSPLP( ) package
    C     that is used to define the matrix entries when this data is passed
    C     to DSPLP( ) as a linear array.  In this usage mode of DSPLP( )
    C     the user gives information about the nonzero entries of A
    C     in DATTRV(*) as given under the description of that parameter.
    C     The name DUSRMT must appear in a Fortran EXTERNAL statement.
    C     Users who are passing the matrix data with DUSRMT( ) can skip
    C     directly to the description of the input parameter DATTRV(*).
    C     Also see option 68 for passing the constraint matrix data using
    C     a standard Fortran two-dimensional array.
    C
    C     If the user chooses to provide a subprogram 'NAME'( ) to
    C     define the matrix A, then DATTRV(*) may be used to pass floating
    C     point data from the user's program unit to the subprogram
    C     'NAME'( ). The content of DATTRV(*) is not changed in any way.
    C
    C     The subprogram 'NAME'( ) can be of the user's choice
    C     but it must meet Fortran standards and it must appear in a
    C     Fortran EXTERNAL statement.  The first statement of the subprogram
    C     has the form
    C
    C          SUBROUTINE 'NAME'(I,J,AIJ, INDCAT, PRGOPT, DATTRV, IFLAG)
    C
    C     The variables I,J, INDCAT, IFLAG(10) are type INTEGER,
    C          while  AIJ, PRGOPT(*),DATTRV(*) are type REAL.
    C
    C     The user interacts with the contents of IFLAG(*) to
    C     direct the appropriate action.  The algorithmic steps are
    C     as follows.
    C
    C          Test IFLAG(1).
    C
    C             IF(IFLAG(1).EQ.1) THEN
    C
    C               Initialize the necessary pointers and data
    C               for defining the matrix A.  The contents
    C               of IFLAG(K), K=2,...,10, may be used for
    C               storage of the pointers.  This array remains intact
    C               between calls to 'NAME'( ) by DSPLP( ).
    C               RETURN
    C
    C             END IF
    C
    C             IF(IFLAG(1).EQ.2) THEN
    C
    C               Define one set of values for I,J,AIJ, and INDCAT.
    C               Each nonzero entry of A must be defined this way.
    C               These values can be defined in any convenient order.
    C               (It is most efficient to define the data by
    C               columns in the order 1,...,NVARS; within each
    C               column define the entries in the order 1,...,MRELAS.)
    C               If this is the last matrix value to be
    C               defined or updated, then set IFLAG(1)=3.
    C               (When I and J are positive and respectively no larger
    C               than MRELAS and NVARS, the value of AIJ is used to
    C               define (or update) row I and column J of A.)
    C               RETURN
    C
    C             END IF
    C
    C               END
    C
    C     Remarks:  The values of I and J are the row and column
    C               indices for the nonzero entries of the matrix A.
    C               The value of this entry is AIJ.
    C               Set INDCAT=0 if this value defines that entry.
    C               Set INDCAT=1 if this entry is to be updated,
    C                            new entry=old entry+AIJ.
    C               A value of I not between 1 and MRELAS, a value of J
    C               not between 1 and NVARS, or a value of INDCAT
    C               not equal to 0 or 1 are each errors.
    C
    C               The contents of IFLAG(K), K=2,...,10, can be used to
    C               remember the status (of the process of defining the
    C               matrix entries) between calls to 'NAME'( ) by DSPLP( ).
    C               On entry to 'NAME'( ), only the values 1 or 2 will be
    C               in IFLAG(1).  More than 2*NVARS*MRELAS definitions of
    C               the matrix elements is considered an error because
    C               it suggests an infinite loop in the user-written
    C               subprogram 'NAME'( ).  Any matrix element not
    C               provided by 'NAME'( ) is defined to be zero.
    C
    C               The REAL arrays PRGOPT(*) and DATTRV(*) are passed as
    C               arguments directly from DSPLP( ) to 'NAME'( ).
    C               The array PRGOPT(*) contains any user-defined program
    C               options.  In this usage mode the array DATTRV(*) may
    C               now contain any (type REAL) data that the user needs
    C               to define the matrix A.  Both arrays PRGOPT(*) and
    C               DATTRV(*) remain intact between calls to 'NAME'( )
    C               by DSPLP( ).
    C     Here is a subprogram that communicates the matrix values for A,
    C     as represented in DATTRV(*), to DSPLP( ).  This subprogram,
    C     called DUSRMT( ), is included as part of the DSPLP( ) package.
    C     This subprogram 'decodes' the array DATTRV(*) and defines the
    C     nonzero entries of the matrix  A for DSPLP( ) to store.  This
    C     listing is presented here as a guide and example
    C     for the users who find it necessary to write their own subroutine
    C     for this purpose.  The contents of DATTRV(*) are given below in
    C     the description of that parameter.
    C
    C     SUBROUTINE DUSRMT(I,J,AIJ, INDCAT,PRGOPT,DATTRV,IFLAG)
    C     DIMENSION PRGOPT(*),DATTRV(*),IFLAG(10)
    C
    C     IF(IFLAG(1).EQ.1) THEN
    C
    C     THIS IS THE INITIALIZATION STEP.  THE VALUES OF IFLAG(K),K=2,3,4,
    C     ARE RESPECTIVELY THE COLUMN INDEX, THE ROW INDEX (OR THE NEXT COL.
    C     INDEX), AND THE POINTER TO THE MATRIX ENTRY'S VALUE WITHIN
    C     DATTRV(*).  ALSO CHECK (DATTRV(1)=0.) SIGNIFYING NO DATA.
    C          IF(DATTRV(1).EQ.0.) THEN
    C          I = 0
    C          J = 0
    C          IFLAG(1) = 3
    C          ELSE
    C          IFLAG(2)=-DATTRV(1)
    C          IFLAG(3)= DATTRV(2)
    C          IFLAG(4)= 3
    C          END IF
    C
    C          RETURN
    C     ELSE
    C          J=IFLAG(2)
    C          I=IFLAG(3)
    C          L=IFLAG(4)
    C          IF(I.EQ.0) THEN
    C
    C     SIGNAL THAT ALL OF THE NONZERO ENTRIES HAVE BEEN DEFINED.
    C               IFLAG(1)=3
    C               RETURN
    C          ELSE IF(I.LT.0) THEN
    C
    C     SIGNAL THAT A SWITCH IS MADE TO A NEW COLUMN.
    C               J=-I
    C               I=DATTRV(L)
    C               L=L+1
    C          END IF
    C
    C          AIJ=DATTRV(L)
    C
    C     UPDATE THE INDICES AND POINTERS FOR THE NEXT ENTRY.
    C          IFLAG(2)=J
    C          IFLAG(3)=DATTRV(L+1)
    C          IFLAG(4)=L+2
    C
    C     INDCAT=0 DENOTES THAT ENTRIES OF THE MATRIX ARE ASSIGNED THE
    C     VALUES FROM DATTRV(*).  NO ACCUMULATION IS PERFORMED.
    C          INDCAT=0
    C          RETURN
    C     END IF
    C     END
    C
    C     DATTRV(*)
    C     ---------
    C     If the user chooses to use the provided subprogram DUSRMT( ) then
    C     the array DATTRV(*) contains data for the matrix A as follows:
    C     Each column (numbered J) requires (floating point) data con-
    C     sisting of the value (-J) followed by pairs of values.  Each pair
    C     consists of the row index immediately followed by the value
    C     of the matrix at that entry.  A value of J=0 signals that there
    C     are no more columns.  (See "Example of DSPLP( ) Usage," below.)
    C     The dimension of DATTRV(*) must be 2*no. of nonzeros
    C     + NVARS + 1 in this usage.  No checking of the array
    C     length is done by the subprogram package.
    C
    C     If the Save/Restore feature is in use (see options with
    C     KEY=55,57 for details of checkpointing and restarting)
    C     DUSRMT( ) can be used to redefine entries of the matrix.
    C     The matrix entries are redefined or overwritten.  No accum-
    C     ulation is performed.
    C     Any other nonzero entry of A, defined in a previous call to
    C     DSPLP( ), remain intact.
    C
    C     BL(*),BU(*),IND(*)
    C     ------------------
    C     The values of IND(*) are input parameters that define
    C     the form of the bounds for the unknowns x and w.  The values for
    C     the bounds are found in the arrays BL(*) and BU(*) as follows.
    C
    C     For values of J between 1 and NVARS,
    C          if IND(J)=1, then X(J) .GE. BL(J); BU(J) is not used.
    C          if IND(J)=2, then X(J) .LE. BU(J); BL(J) is not used.
    C          if IND(J)=3, then BL(J) .LE. X(J) .LE. BU(J),(BL(J)=BU(J) ok)
    C          if IND(J)=4, then X(J) is free to have any value,
    C          and BL(J), BU(J) are not used.
    C
    C     For values of I between NVARS+1 and NVARS+MRELAS,
    C          if IND(I)=1, then W(I-NVARS) .GE. BL(I); BU(I) is not used.
    C          if IND(I)=2, then W(I-NVARS) .LE. BU(I); BL(I) is not used.
    C          if IND(I)=3, then BL(I) .LE. W(I-NVARS) .LE. BU(I),
    C          (BL(I)=BU(I) is ok).
    C          if IND(I)=4, then W(I-NVARS) is free to have any value,
    C          and BL(I), BU(I) are not used.
    C
    C     A value of IND(*) not equal to 1,2,3 or 4 is an error.  When
    C     IND(I)=3, BL(I) must be .LE. BU(I).  The condition BL(I).GT.
    C     BU(I) indicates infeasibility and is an error.  These
    C     arrays can be changed by the user between restarts.  See
    C     options with KEY=55,57 for details of checkpointing and
    C     restarting.
    C
    C     PRGOPT(*)
    C     ---------
    C     This array is used to redefine various parameters within DSPLP( ).
    C     Frequently, perhaps most of the time, a user will be satisfied
    C     and obtain the solutions with no changes to any of these
    C     parameters.  To try this, simply set PRGOPT(1)=1.D0.
    C
    C     For users with more sophisticated needs, DSPLP( ) provides several
    C     options that may be used to take advantage of more detailed
    C     knowledge of the problem or satisfy other utilitarian needs.
    C     The complete description of how to use this option array to
    C     utilize additional subprogram features is found under the
    C     heading "Usage of DSPLP( ) Subprogram Options."
    C
    C     Briefly, the user should note the following value of the parameter
    C     KEY and the corresponding task or feature desired before turning
    C     to that section.
    C
    C     Value     Brief Statement of Purpose for Option
    C     of KEY
    C     ------    -------------------------------------
    C     50        Change from a minimization problem to a
    C               maximization problem.
    C     51        Change the amount of printed output.
    C               Normally, no printed output is obtained.
    C     52        Redefine the line length and precision used
    C               for the printed output.
    C     53        Redefine the values of LAMAT and LBM that
    C               were discussed above under the heading
    C               Fortran Declarations Required.
    C     54        Redefine the unit number where pages of the sparse
    C               data matrix A are stored.  Normally, the unit
    C               number is 1.
    C     55        A computation, partially completed, is
    C               being continued.  Read the up-to-date
    C               partial results from unit number 2.
    C     56        Redefine the unit number where the partial results
    C               are stored.  Normally, the unit number is 2.
    C     57        Save partial results on unit 2 either after
    C               maximum iterations or at the optimum.
    C     58        Redefine the value for the maximum number of
    C               iterations.  Normally, the maximum number of
    C               iterations is 3*(NVARS+MRELAS).
    C     59        Provide DSPLP( ) with a starting (feasible)
    C               nonsingular basis.  Normally, DSPLP( ) starts
    C               with the identity matrix columns corresponding
    C               to the vector w.
    C     60        The user has provided scale factors for the
    C               columns of A.  Normally, DSPLP( ) computes scale
    C               factors that are the reciprocals of the max. norm
    C               of each column.
    C     61        The user has provided a scale factor
    C               for the vector costs.  Normally, DSPLP( ) computes
    C               a scale factor equal to the reciprocal of the
    C               max. norm of the vector costs after the column
    C               scaling for the data matrix has been applied.
    C     62        Size parameters, namely the smallest and
    C               largest magnitudes of nonzero entries in
    C               the matrix A, are provided.  Values noted
    C               outside this range are to be considered errors.
    C     63        Redefine the tolerance required in
    C               evaluating residuals for feasibility.
    C               Normally, this value is set to the value RELPR,
    C               where RELPR = relative precision of the arithmetic.
    C     64        Change the criterion for bringing new variables
    C               into the basis from the steepest edge (best
    C               local move) to the minimum reduced cost.
    C     65        Redefine the value for the number of iterations
    C               between recalculating the error in the primal
    C               solution.  Normally, this value is equal to ten.
    C     66        Perform "partial pricing" on variable selection.
    C               Redefine the value for the number of negative
    C               reduced costs to compute (at most) when finding
    C               a variable to enter the basis.  Normally this
    C               value is set to NVARS.  This implies that no
    C               "partial pricing" is used.
    C     67        Adjust the tuning factor (normally one) to apply
    C               to the primal and dual error estimates.
    C     68        Pass  information to the  subprogram  DFULMT(),
    C               provided with the DSPLP() package, so that a Fortran
    C               two-dimensional array can be used as the argument
    C               DATTRV(*).
    C     69        Pass an absolute tolerance to use for the feasibility
    C               test when the usual relative error test indicates
    C               infeasibility.  The nominal value of this tolerance,
    C               TOLABS, is zero.
    C
    C
    C  |---------------|
    C  |Working Arrays:|
    C  |---------------|
    C
    C     WORK(*),LW,
    C     IWORK(*),LIW
    C     ------------
    C     The arrays WORK(*) and IWORK(*) are respectively floating point
    C     and type INTEGER working arrays for DSPLP( ) and its
    C     subprograms.  The lengths of these arrays are respectively
    C     LW and LIW.  These parameters must satisfy the inequalities
    C     noted above under the heading "Fortran Declarations Required."
    C     It is an error if either value is too small.
    C
    C  |----------------------------|
    C  |Input/Output files required:|
    C  |----------------------------|
    C
    C     Fortran unit 1 is used by DSPLP( ) to store the sparse matrix A
    C     out of high-speed memory.  This direct access file is opened
    C     within the package under the following two conditions.
    C     1. When the Save/Restore feature is used.  2. When the
    C     constraint matrix is so large that storage out of high-speed
    C     memory is required.  The user may need to close unit 1
    C     (with deletion from the job step) in the main program unit
    C     when several calls are made to DSPLP( ).  A crude
    C     upper bound for the amount of information written on unit 1
    C     is 6*nz, where nz is the number of nonzero entries in A.
    C     The unit number may be redefined to any other positive value
    C     by means of input in the option array PRGOPT(*).
    C
    C     Fortran unit 2 is used by DSPLP( ) only when the Save/Restore
    C     feature is desired.  Normally this feature is not used.  It is
    C     activated by means of input in the option array PRGOPT(*).
    C     On some computer systems the user may need to open unit
    C     2 before executing a call to DSPLP( ).  This file is type
    C     sequential and is unformatted.
    C
    C     Fortran unit=I1MACH(2) (check local setting) is used by DSPLP( )
    C     when the printed output feature (KEY=51) is used.  Normally
    C     this feature is not used.  It is activated by input in the
    C     options array PRGOPT(*).  For many computer systems I1MACH(2)=6.
    C
    C  |-------|
    C  |Output:|
    C  |-------|
    C
    C     INFO,PRIMAL(*),DUALS(*)
    C     -----------------------
    C     The integer flag INFO indicates why DSPLP( ) has returned to the
    C     user.  If INFO=1 the solution has been computed.  In this case
    C     X(J)=PRIMAL(J) and W(I)=PRIMAL(I+NVARS).  The dual variables
    C     for the equations A*x=w are in the array DUALS(I)=dual for
    C     equation number I.  The dual value for the component X(J) that
    C     has an upper or lower bound (or both) is returned in
    C     DUALS(J+MRELAS).  The only other values for INFO are .LT. 0.
    C     The meaning of these values can be found by reading
    C     the diagnostic message in the output file, or by looking for
    C     error number = (-INFO) under the heading "List of DSPLP( ) Error
    C     and Diagnostic Messages."
    C     The diagnostic messages are printed using the error processing
    C     subprogram XERMSG( ) with error category LEVEL=1.
    C     See the document "Brief Instr. for Using the Sandia Math.
    C     Subroutine Library," SAND79-2382, Nov., 1980, for further inform-
    C     ation about resetting the usual response to a diagnostic message.
    C
    C     BL(*),BU(*),IND(*)
    C     ------------------
    C     These arrays are output parameters only under the (unusual)
    C     circumstances where the stated problem is infeasible, has an
    C     unbounded optimum value, or both.  These respective conditions
    C     correspond to INFO=-1,-2 or -3.  For INFO=-1 or -3 certain comp-
    C     onents of the vectors x or w will not satisfy the input bounds.
    C     If component J of X or component I of W does not satisfy its input
    C     bound because of infeasibility, then IND(J)=-4 or IND(I+NVARS)=-4,
    C     respectively.  For INFO=-2 or -3 certain
    C     components of the vector x could not be used as basic variables
    C     because the objective function would have become unbounded.
    C     In particular if component J of x corresponds to such a variable,
    C     then IND(J)=-3.  Further, if the input value of IND(J)
    C                      =1, then BU(J)=BL(J);
    C                      =2, then BL(J)=BU(J);
    C                      =4, then BL(J)=0.,BU(J)=0.
    C
    C     (The J-th variable in x has been restricted to an appropriate
    C     feasible value.)
    C     The negative output value for IND(*) allows the user to identify
    C     those constraints that are not satisfied or those variables that
    C     would cause unbounded values of the objective function.  Note
    C     that the absolute value of IND(*), together with BL(*) and BU(*),
    C     are valid input to DSPLP( ).  In the case of infeasibility the
    C     sum of magnitudes of the infeasible values is minimized.  Thus
    C     one could reenter DSPLP( ) with these components of x or w now
    C     fixed at their present values.  This involves setting
    C     the appropriate components of IND(*) = 3, and BL(*) = BU(*).
    C
    C     IBASIS(I),I=1,...,MRELAS
    C     ------------------------
    C     This array contains the indices of the variables that are
    C     in the active basis set at the solution (INFO=1).  A value
    C     of IBASIS(I) between 1 and NVARS corresponds to the variable
    C     X(IBASIS(I)).  A value of IBASIS(I) between NVARS+1 and NVARS+
    C     MRELAS corresponds to the variable W(IBASIS(I)-NVARS).
    C
    C     Computing with the Matrix A after Calling DSPLP( )
    C     --------------------------------------------------
    C     Following the return from DSPLP( ), nonzero entries of the MRELAS
    C     by NVARS matrix A are available for usage by the user.  The method
    C     for obtaining the next nonzero in column J with a row index
    C     strictly greater than I in value, is completed by executing
    C
    C         CALL DPNNZR(I,AIJ,IPLACE,WORK,IWORK,J)
    C
    C     The value of I is also an output parameter.  If I.LE.0 on output,
    C     then there are no more nonzeroes in column J.  If I.GT.0, the
    C     output value for component number I of column J is in AIJ.  The
    C     parameters WORK(*) and IWORK(*) are the same arguments as in the
    C     call to DSPLP( ).  The parameter IPLACE is a single INTEGER
    C     working variable.
    C
    C     The data structure used for storage of the matrix A within DSPLP()
    C     corresponds to sequential storage by columns as defined in
    C     SAND78-0785.  Note that the names of the subprograms LNNZRS(),
    C     LCHNGS(),LINITM(),LLOC(),LRWPGE(), and LRWVIR() have been
    C     changed to DPNNZR(),DPCHNG(),PINITM(),IPLOC(),DPRWPG(), and
    C     DPRWVR() respectively.  The error processing subprogram LERROR()
    C     is no longer used; XERMSG() is used instead.
    C
    C  |--------------------------------|
    C  |Subprograms Required by DSPLP( )|
    C  |--------------------------------|
    C     Called by DSPLP() are DPLPMN(),DPLPUP(),DPINIT(),DPOPT(),
    C                          DPLPDM(),DPLPCE(),DPINCW(),DPLPFL(),
    C                          DPLPFE(),DPLPMU().
    C
    C     Error Processing Subprograms XERMSG(),I1MACH(),D1MACH()
    C
    C     Sparse Matrix Subprograms DPNNZR(),DPCHNG(),DPRWPG(),DPRWVR(),
    C                               PINITM(),IPLOC()
    C
    C     Mass Storage File Subprograms SOPENM(),SCLOSM(),DREADP(),DWRITP()
    C
    C     Basic Linear Algebra Subprograms DCOPY(),DASUM(),DDOT()
    C
    C     Sparse Matrix Basis Handling Subprograms LA05AD(),LA05BD(),
    C                                             LA05CD(),LA05ED(),MC20AD()
    C
    C     Vector Output Subprograms DVOUT(),IVOUT()
    C
    C     Machine-sensitive Subprograms I1MACH( ),D1MACH( ),
    C                                   SOPENM(),SCLOSM(),DREADP(),DWRITP().
    C     COMMON Block Used
    C     -----------------
    C     /LA05DD/ SMALL,LP,LENL,LENU,NCP,LROW,LCOL
    C     See the document AERE-R8269 for further details.
    C    |-------------------------|
    C    |Example of DSPLP( ) Usage|
    C    |-------------------------|
    C     PROGRAM LPEX
    C     THE OPTIMIZATION PROBLEM IS TO FIND X1, X2, X3 THAT
    C     MINIMIZE X1 + X2 + X3, X1.GE.0, X2.GE.0, X3 UNCONSTRAINED.
    C
    C     THE UNKNOWNS X1,X2,X3 ARE TO SATISFY CONSTRAINTS
    C
    C        X1 -3*X2 +4*X3 = 5
    C        X1 -2*X2     .LE.3
    C            2*X2 - X3.GE.4
    C
    C     WE FIRST DEFINE THE DEPENDENT VARIABLES
    C          W1=X1 -3*X2 +4*X3
    C          W2=X1- 2*X2
    C          W3=    2*X2 -X3
    C
    C     WE NOW SHOW HOW TO USE DSPLP( ) TO SOLVE THIS LINEAR OPTIMIZATION
    C     PROBLEM.  EACH REQUIRED STEP WILL BE SHOWN IN THIS EXAMPLE.
    C     DIMENSION COSTS(03),PRGOPT(01),DATTRV(18),BL(06),BU(06),IND(06),
    C    *PRIMAL(06),DUALS(06),IBASIS(06),WORK(079),IWORK(103)
    C
    C     EXTERNAL DUSRMT
    C     MRELAS=3
    C     NVARS=3
    C
    C     DEFINE THE ARRAY COSTS(*) FOR THE OBJECTIVE FUNCTION.
    C     COSTS(01)=1.
    C     COSTS(02)=1.
    C     COSTS(03)=1.
    C
    C     PLACE THE NONZERO INFORMATION ABOUT THE MATRIX IN DATTRV(*).
    C     DEFINE COL. 1:
    C     DATTRV(01)=-1
    C     DATTRV(02)=1
    C     DATTRV(03)=1.
    C     DATTRV(04)=2
    C     DATTRV(05)=1.
    C
    C     DEFINE COL. 2:
    C     DATTRV(06)=-2
    C     DATTRV(07)=1
    C     DATTRV(08)=-3.
    C     DATTRV(09)=2
    C     DATTRV(10)=-2.
    C     DATTRV(11)=3
    C     DATTRV(12)=2.
    C
    C     DEFINE COL. 3:
    C     DATTRV(13)=-3
    C     DATTRV(14)=1
    C     DATTRV(15)=4.
    C     DATTRV(16)=3
    C     DATTRV(17)=-1.
    C
    C     DATTRV(18)=0
    C
    C     CONSTRAIN X1,X2 TO BE NONNEGATIVE. LET X3 HAVE NO BOUNDS.
    C     BL(1)=0.
    C     IND(1)=1
    C     BL(2)=0.
    C     IND(2)=1
    C     IND(3)=4
    C
    C     CONSTRAIN W1=5,W2.LE.3, AND W3.GE.4.
    C     BL(4)=5.
    C     BU(4)=5.
    C     IND(4)=3
    C     BU(5)=3.
    C     IND(5)=2
    C     BL(6)=4.
    C     IND(6)=1
    C
    C     INDICATE THAT NO MODIFICATIONS TO OPTIONS ARE IN USE.
    C     PRGOPT(01)=1
    C
    C     DEFINE THE WORKING ARRAY LENGTHS.
    C     LW=079
    C     LIW=103
    C     CALL DSPLP(DUSRMT,MRELAS,NVARS,COSTS,PRGOPT,DATTRV,
    C    *BL,BU,IND,INFO,PRIMAL,DUALS,IBASIS,WORK,LW,IWORK,LIW)
    C
    C     CALCULATE VAL, THE MINIMAL VALUE OF THE OBJECTIVE FUNCTION.
    C     VAL=DDOT(NVARS,COSTS,1,PRIMAL,1)
    C
    C     STOP
    C     END
    C    |------------------------|
    C    |End of Example of Usage |
    C    |------------------------|
    C
    C    |-------------------------------------|
    C    |Usage of DSPLP( ) Subprogram Options.|
    C    |-------------------------------------|
    C
    C     Users frequently have a large variety of requirements for linear
    C     optimization software.  Allowing for these varied requirements
    C     is at cross purposes with the desire to keep the usage of DSPLP( )
    C     as simple as possible. One solution to this dilemma is as follows.
    C     (1) Provide a version of DSPLP( ) that solves a wide class of
    C     problems and is easy to use. (2) Identify parameters within
    C     DSPLP() that certain users may want to change.  (3) Provide a
    C     means of changing any selected number of these parameters that
    C     does not require changing all of them.
    C
    C     Changing selected parameters is done by requiring
    C     that the user provide an option array, PRGOPT(*), to DSPLP( ).
    C     The contents of PRGOPT(*) inform DSPLP( ) of just those options
    C     that are going to be modified within the total set of possible
    C     parameters that can be modified.  The array PRGOPT(*) is a linked
    C     list consisting of groups of data of the following form
    C
    C          LINK
    C          KEY
    C          SWITCH
    C          data set
    C
    C     that describe the desired options.  The parameters LINK, KEY and
    C     switch are each one word and are always required.  The data set
    C     can be comprised of several words or can be empty.  The number of
    C     words in the data set for each option depends on the value of
    C     the parameter KEY.
    C
    C     The value of LINK points to the first entry of the next group
    C     of data within PRGOPT(*).  The exception is when there are no more
    C     options to change.  In that case, LINK=1 and the values for KEY,
    C     SWITCH and data set are not referenced.  The general layout of
    C     PRGOPT(*) is as follows:
    C          ...PRGOPT(1)=LINK1 (link to first entry of next group)
    C          .  PRGOPT(2)=KEY1 (KEY to the option change)
    C          .  PRGOPT(3)=SWITCH1 (on/off switch for the option)
    C          .  PRGOPT(4)=data value
    C          .       .
    C          .       .
    C          .       .
    C          ...PRGOPT(LINK1)=LINK2 (link to first entry of next group)
    C          .  PRGOPT(LINK1+1)=KEY2 (KEY to option change)
    C          .  PRGOPT(LINK1+2)=SWITCH2 (on/off switch for the option)
    C          .  PRGOPT(LINK1+3)=data value
    C          ...     .
    C          .       .
    C          .       .
    C          ...PRGOPT(LINK)=1 (no more options to change)
    C
    C     A value of LINK that is .LE.0 or .GT. 10000 is an error.
    C     In this case DSPLP( ) returns with an error message, INFO=-14.
    C     This helps prevent using invalid but positive values of LINK that
    C     will probably extend beyond the program limits of PRGOPT(*).
    C     Unrecognized values of KEY are ignored.  If the value of SWITCH is
    C     zero then the option is turned off.  For any other value of SWITCH
    C     the option is turned on.  This is used to allow easy changing of
    C     options without rewriting PRGOPT(*).  The order of the options is
    C     arbitrary and any number of options can be changed with the
    C     following restriction.  To prevent cycling in processing of the
    C     option array PRGOPT(*), a count of the number of options changed
    C     is maintained.  Whenever this count exceeds 1000 an error message
    C     (INFO=-15) is printed and the subprogram returns.
    C
    C     In the following description of the options, the value of
    C     LATP indicates the amount of additional storage that a particular
    C     option requires.  The sum of all of these values (plus one) is
    C     the minimum dimension for the array PRGOPT(*).
    C
    C     If a user is satisfied with the nominal form of DSPLP( ),
    C     set PRGOPT(1)=1 (or PRGOPT(1)=1.D0).
    C
    C     Options:
    C
    C -----KEY = 50.  Change from a minimization problem to a maximization
    C     problem.
    C     If SWITCH=0  option is off; solve minimization problem.
    C              =1  option is on; solve maximization problem.
    C     data set =empty
    C     LATP=3
    C
    C -----KEY = 51.  Change the amount of printed output.  The nominal form
    C     of DSPLP( ) has no printed output.
    C     The first level of output (SWITCH=1) includes
    C
    C     (1) Minimum dimensions for the arrays COSTS(*),BL(*),BU(*),IND(*),
    C         PRIMAL(*),DUALS(*),IBASIS(*), and PRGOPT(*).
    C     (2) Problem dimensions MRELAS,NVARS.
    C     (3) The types of and values for the bounds on x and w,
    C         and the values of the components of the vector costs.
    C     (4) Whether optimization problem is minimization or
    C         maximization.
    C     (5) Whether steepest edge or smallest reduced cost criteria used
    C         for exchanging variables in the revised simplex method.
    C
    C     Whenever a solution has been found, (INFO=1),
    C
    C     (6) the value of the objective function,
    C     (7) the values of the vectors x and w,
    C     (8) the dual variables for the constraints A*x=w and the
    C         bounded components of x,
    C     (9) the indices of the basic variables,
    C    (10) the number of revised simplex method iterations,
    C    (11) the number of full decompositions of the basis matrix.
    C
    C     The second level of output (SWITCH=2) includes all for SWITCH=1
    C     plus
    C
    C    (12) the iteration number,
    C    (13) the column number to enter the basis,
    C    (14) the column number to leave the basis,
    C    (15) the length of the step taken.
    C
    C     The third level of output (SWITCH=3) includes all for SWITCH=2
    C     plus
    C    (16) critical quantities required in the revised simplex method.
    C          This output is rather voluminous.  It is intended to be used
    C          as a diagnostic tool in case of a failure in DSPLP( ).
    C
    C     If SWITCH=0  option is off; no printed output.
    C              =1  summary output.
    C              =2  lots of output.
    C              =3  even more output.
    C     data set =empty
    C     LATP=3
    C
    C -----KEY = 52.  Redefine the parameter, IDIGIT, which determines the
    C     format and precision used for the printed output.  In the printed
    C     output, at least ABS(IDIGIT) decimal digits per number is
    C     printed.  If IDIGIT.LT.0, 72 printing columns are used.  If
    C     IDIGIT.GT.0, 133 printing columns are used.
    C     If SWITCH=0  option is off; IDIGIT=-4.
    C              =1  option is on.
    C     data set =IDIGIT
    C     LATP=4
    C
    C -----KEY = 53.  Redefine LAMAT and LBM, the lengths of the portions of
    C     WORK(*) and IWORK(*) that are allocated to the sparse matrix
    C     storage and the sparse linear equation solver, respectively.
    C     LAMAT must be .GE. NVARS+7 and LBM must be positive.
    C     If SWITCH=0  option is off; LAMAT=4*NVARS+7
    C                                 LBM  =8*MRELAS.
    C              =1  option is on.
    C     data set =LAMAT
    C               LBM
    C     LATP=5
    C
    C -----KEY = 54. Redefine IPAGEF, the file number where the pages of the
    C     sparse data matrix are stored.  IPAGEF must be positive and
    C     different from ISAVE (see option 56).
    C     If SWITCH=0  option is off; IPAGEF=1.
    C              =1  option is on.
    C     data set =IPAGEF
    C     LATP=4
    C
    C -----KEY = 55.  Partial results have been computed and stored on unit
    C     number ISAVE (see option 56), during a previous run of
    C     DSPLP( ).  This is a continuation from these partial results.
    C     The arrays COSTS(*),BL(*),BU(*),IND(*) do not have to have
    C     the same values as they did when the checkpointing occurred.
    C     This feature makes it possible for the user to do certain
    C     types of parameter studies such as changing costs and varying
    C     the constraints of the problem.  This file is rewound both be-
    C     fore and after reading the partial results.
    C     If SWITCH=0  option is off; start a new problem.
    C              =1  option is on; continue from partial results
    C                  that are stored in file ISAVE.
    C     data set = empty
    C     LATP=3
    C
    C -----KEY = 56.  Redefine ISAVE, the file number where the partial
    C     results are stored (see option 57).  ISAVE must be positive and
    C     different from IPAGEF (see option 54).
    C     If SWITCH=0  option is off; ISAVE=2.
    C              =1  option is on.
    C     data set =ISAVE
    C     LATP=4
    C
    C -----KEY = 57.  Save the partial results after maximum number of
    C     iterations, MAXITR, or at the optimum.  When this option is on,
    C     data essential to continuing the calculation is saved on a file
    C     using a Fortran binary write operation.  The data saved includes
    C     all the information about the sparse data matrix A.  Also saved
    C     is information about the current basis.  Nominally the partial
    C     results are saved on Fortran unit 2.  This unit number can be
    C     redefined (see option 56).  If the save option is on,
    C     this file must be opened (or declared) by the user prior to the
    C     call to DSPLP( ).  A crude upper bound for the number of words
    C     written to this file is 6*nz.  Here nz= number of nonzeros in A.
    C     If SWITCH=0  option is off; do not save partial results.
    C              =1  option is on; save partial results.
    C     data set = empty
    C     LATP=3
    C
    C -----KEY = 58.  Redefine the maximum number of iterations, MAXITR, to
    C     be taken before returning to the user.
    C     If SWITCH=0  option is off; MAXITR=3*(NVARS+MRELAS).
    C              =1  option is on.
    C     data set =MAXITR
    C     LATP=4
    C
    C -----KEY = 59.  Provide DSPLP( ) with exactly MRELAS indices which
    C     comprise a feasible, nonsingular basis.  The basis must define a
    C     feasible point: values for x and w such that A*x=w and all the
    C     stated bounds on x and w are satisfied.  The basis must also be
    C     nonsingular.  The failure of either condition will cause an error
    C     message (INFO=-23 or =-24, respectively).  Normally, DSPLP( ) uses
    C     identity matrix columns which correspond to the components of w.
    C     This option would normally not be used when restarting from
    C     a previously saved run (KEY=57).
    C     In numbering the unknowns,
    C     the components of x are numbered (1-NVARS) and the components
    C     of w are numbered (NVARS+1)-(NVARS+MRELAS).  A value for an
    C     index .LE. 0 or .GT. (NVARS+MRELAS) is an error (INFO=-16).
    C     If SWITCH=0  option is off; DSPLP( ) chooses the initial basis.
    C              =1  option is on; user provides the initial basis.
    C     data set =MRELAS indices of basis; order is arbitrary.
    C     LATP=MRELAS+3
    C
    C -----KEY = 60.  Provide the scale factors for the columns of the data
    C     matrix A.  Normally, DSPLP( ) computes the scale factors as the
    C     reciprocals of the max. norm of each column.
    C     If SWITCH=0  option is off; DSPLP( ) computes the scale factors.
    C              =1  option is on; user provides the scale factors.
    C     data set =scaling for column J, J=1,NVARS; order is sequential.
    C     LATP=NVARS+3
    C
    C -----KEY = 61.  Provide a scale factor, COSTSC, for the vector of
    C     costs.  Normally, DSPLP( ) computes this scale factor to be the
    C     reciprocal of the max. norm of the vector costs after the column
    C     scaling has been applied.
    C     If SWITCH=0  option is off; DSPLP( ) computes COSTSC.
    C              =1  option is on; user provides COSTSC.
    C     data set =COSTSC
    C     LATP=4
    C
    C -----KEY = 62.  Provide size parameters, ASMALL and ABIG, the smallest
    C     and largest magnitudes of nonzero entries in the data matrix A,
    C     respectively.  When this option is on, DSPLP( ) will check the
    C     nonzero entries of A to see if they are in the range of ASMALL and
    C     ABIG.  If an entry of A is not within this range, DSPLP( ) returns
    C     an error message, INFO=-22. Both ASMALL and ABIG must be positive
    C     with ASMALL .LE. ABIG.  Otherwise,  an error message is returned,
    C     INFO=-17.
    C     If SWITCH=0  option is off; no checking of the data matrix is done
    C              =1  option is on; checking is done.
    C     data set =ASMALL
    C               ABIG
    C     LATP=5
    C
    C -----KEY = 63.  Redefine the relative tolerance, TOLLS, used in
    C     checking if the residuals are feasible.  Normally,
    C     TOLLS=RELPR, where RELPR is the machine precision.
    C     If SWITCH=0  option is off; TOLLS=RELPR.
    C              =1  option is on.
    C     data set =TOLLS
    C     LATP=4
    C
    C -----KEY = 64. Use the minimum reduced cost pricing strategy to choose
    C     columns to enter the basis.  Normally, DSPLP( ) uses the steepest
    C     edge pricing strategy which is the best local move.  The steepest
    C     edge pricing strategy generally uses fewer iterations than the
    C     minimum reduced cost pricing, but each iteration costs more in the
    C     number of calculations done.  The steepest edge pricing is
    C     considered to be more efficient.  However, this is very problem
    C     dependent.  That is why DSPLP( ) provides the option of either
    C     pricing strategy.
    C     If SWITCH=0  option is off; steepest option edge pricing is used.
    C              =1  option is on; minimum reduced cost pricing is used.
    C     data set =empty
    C     LATP=3
    C
    C -----KEY = 65.  Redefine MXITBR, the number of iterations between
    C     recalculating the error in the primal solution.  Normally, MXITBR
    C     is set to 10.  The error in the primal solution is used to monitor
    C     the error in solving the linear system.  This is an expensive
    C     calculation and every tenth iteration is generally often enough.
    C     If SWITCH=0  option is off; MXITBR=10.
    C              =1  option is on.
    C     data set =MXITBR
    C     LATP=4
    C
    C -----KEY = 66.  Redefine NPP, the number of negative reduced costs
    C     (at most) to be found at each iteration of choosing
    C     a variable to enter the basis.  Normally NPP is set
    C     to NVARS which implies that all of the reduced costs
    C     are computed at each such step.  This "partial
    C     pricing" may very well increase the total number
    C     of iterations required.  However it decreases the
    C     number of calculations at each iteration.
    C     therefore the effect on overall efficiency is quite
    C     problem-dependent.
    C
    C     if SWITCH=0 option is off; NPP=NVARS
    C              =1 option is on.
    C     data set =NPP
    C     LATP=4
    C
    C -----KEY =  67.  Redefine the tuning factor (PHI) used to scale the
    C     error estimates for the primal and dual linear algebraic systems
    C     of equations.  Normally, PHI = 1.D0, but in some environments it
    C     may be necessary to reset PHI to the range 0.001-0.01.  This is
    C     particularly important for machines with short word lengths.
    C
    C     if SWITCH = 0 option is off; PHI=1.D0.
    C               = 1 option is on.
    C     Data Set  = PHI
    C     LATP=4
    C
    C -----KEY = 68.  Used together with the subprogram DFULMT(), provided
    C     with the DSPLP() package, for passing a standard Fortran two-
    C     dimensional array containing the constraint matrix.  Thus the sub-
    C     program DFULMT must be declared in a Fortran EXTERNAL statement.
    C     The two-dimensional array is passed as the argument DATTRV.
    C     The information about the array and problem dimensions are passed
    C     in the option array PRGOPT(*).  It is an error if DFULMT() is
    C     used and this information is not passed in PRGOPT(*).
    C
    C     if SWITCH = 0 option is off; this is an error is DFULMT() is
    C                                  used.
    C               = 1 option is on.
    C     Data Set  = IA = row dimension of two-dimensional array.
    C                 MRELAS = number of constraint equations.
    C                 NVARS  = number of dependent variables.
    C     LATP = 6
    C -----KEY = 69.  Normally a relative tolerance (TOLLS, see option 63)
    C     is used to decide if the problem is feasible.  If this test fails
    C     an absolute test will be applied using the value TOLABS.
    C     Nominally TOLABS = zero.
    C     If SWITCH = 0 option is off; TOLABS = zero.
    C               = 1 option is on.
    C     Data set  = TOLABS
    C     LATP = 4
    C
    C    |-----------------------------|
    C    |Example of Option array Usage|
    C    |-----------------------------|
    C     To illustrate the usage of the option array, let us suppose that
    C     the user has the following nonstandard requirements:
    C
    C          a) Wants to change from minimization to maximization problem.
    C          b) Wants to limit the number of simplex steps to 100.
    C          c) Wants to save the partial results after 100 steps on
    C             Fortran unit 2.
    C
    C     After these 100 steps are completed the user wants to continue the
    C     problem (until completed) using the partial results saved on
    C     Fortran unit 2.  Here are the entries of the array PRGOPT(*)
    C     that accomplish these tasks.  (The definitions of the other
    C     required input parameters are not shown.)
    C
    C     CHANGE TO A MAXIMIZATION PROBLEM; KEY=50.
    C     PRGOPT(01)=4
    C     PRGOPT(02)=50
    C     PRGOPT(03)=1
    C
    C     LIMIT THE NUMBER OF SIMPLEX STEPS TO 100; KEY=58.
    C     PRGOPT(04)=8
    C     PRGOPT(05)=58
    C     PRGOPT(06)=1
    C     PRGOPT(07)=100
    C
    C     SAVE THE PARTIAL RESULTS, AFTER 100 STEPS, ON FORTRAN
    C     UNIT 2; KEY=57.
    C     PRGOPT(08)=11
    C     PRGOPT(09)=57
    C     PRGOPT(10)=1
    C
    C     NO MORE OPTIONS TO CHANGE.
    C     PRGOPT(11)=1
    C     The user makes the CALL statement for DSPLP( ) at this point.
    C     Now to restart, using the partial results after 100 steps, define
    C     new values for the array PRGOPT(*):
    C
    C     AGAIN INFORM DSPLP( ) THAT THIS IS A MAXIMIZATION PROBLEM.
    C     PRGOPT(01)=4
    C     PRGOPT(02)=50
    C     PRGOPT(03)=1
    C
    C     RESTART, USING SAVED PARTIAL RESULTS; KEY=55.
    C     PRGOPT(04)=7
    C     PRGOPT(05)=55
    C     PRGOPT(06)=1
    C
    C     NO MORE OPTIONS TO CHANGE.  THE SUBPROGRAM DSPLP( ) IS NO LONGER
    C     LIMITED TO 100 SIMPLEX STEPS BUT WILL RUN UNTIL COMPLETION OR
    C     MAX.=3*(MRELAS+NVARS) ITERATIONS.
    C     PRGOPT(07)=1
    C     The user now makes a CALL to subprogram DSPLP( ) to compute the
    C     solution.
    C    |--------------------------------------------|
    C    |End of Usage of DSPLP( ) Subprogram Options.|
    C    |--------------------------------------------|
    C
    C     |-----------------------------------------------|
    C     |List of DSPLP( ) Error and Diagnostic Messages.|
    C     |-----------------------------------------------|
    C      This section may be required to understand the meanings of the
    C     error flag =-INFO  that may be returned from DSPLP( ).
    C
    C -----1. There is no set of values for x and w that satisfy A*x=w and
    C     the stated bounds.  The problem can be made feasible by ident-
    C     ifying components of w that are now infeasible and then rede-
    C     signating them as free variables.  Subprogram DSPLP( ) only
    C     identifies an infeasible problem; it takes no other action to
    C     change this condition.  Message:
    C     DSPLP( ). THE PROBLEM APPEARS TO BE INFEASIBLE.
    C     ERROR NUMBER =         1
    C
    C     2. One of the variables in either the vector x or w was con-
    C     strained at a bound.  Otherwise the objective function value,
    C     (transpose of costs)*x, would not have a finite optimum.
    C     Message:
    C     DSPLP( ). THE PROBLEM APPEARS TO HAVE NO FINITE SOLN.
    C     ERROR NUMBER =         2
    C
    C     3.  Both of the conditions of 1. and 2. above have occurred.
    C     Message:
    C     DSPLP( ). THE PROBLEM APPEARS TO BE INFEASIBLE AND TO
    C     HAVE NO FINITE SOLN.
    C     ERROR NUMBER =         3
    C
    C -----4.  The REAL and INTEGER working arrays, WORK(*) and IWORK(*),
    C     are not long enough. The values (I1) and (I2) in the message
    C     below will give you the minimum length required.  Also redefine
    C     LW and LIW, the lengths of these arrays.  Message:
    C     DSPLP( ). WORK OR IWORK IS NOT LONG ENOUGH. LW MUST BE (I1)
    C     AND LIW MUST BE (I2).
    C               IN ABOVE MESSAGE, I1=         0
    C               IN ABOVE MESSAGE, I2=         0
    C     ERROR NUMBER =        4
    C
    C -----5. and 6.  These error messages often mean that one or more
    C     arguments were left out of the call statement to DSPLP( ) or
    C     that the values of MRELAS and NVARS have been over-written
    C     by garbage.  Messages:
    C     DSPLP( ). VALUE OF MRELAS MUST BE .GT.0. NOW=(I1).
    C               IN ABOVE MESSAGE, I1=         0
    C     ERROR NUMBER =         5
    C
    C     DSPLP( ). VALUE OF NVARS MUST BE .GT.0. NOW=(I1).
    C               IN ABOVE MESSAGE, I1=         0
    C     ERROR NUMBER =         6
    C
    C -----7.,8., and 9.  These error messages can occur as the data matrix
    C     is being defined by either DUSRMT( ) or the user-supplied sub-
    C     program, 'NAME'( ). They would indicate a mistake in the contents
    C     of DATTRV(*), the user-written subprogram or that data has been
    C     over-written.
    C     Messages:
    C     DSPLP( ). MORE THAN 2*NVARS*MRELAS ITERS. DEFINING OR UPDATING
    C     MATRIX DATA.
    C     ERROR NUMBER =        7
    C
    C     DSPLP( ). ROW INDEX (I1) OR COLUMN INDEX (I2) IS OUT OF RANGE.
    C               IN ABOVE MESSAGE, I1=         1
    C               IN ABOVE MESSAGE, I2=        12
    C     ERROR NUMBER =        8
    C
    C     DSPLP( ). INDICATION FLAG (I1) FOR MATRIX DATA MUST BE
    C     EITHER 0 OR 1.
    C               IN ABOVE MESSAGE, I1=        12
    C     ERROR NUMBER =        9
    C
    C -----10. and 11.  The type of bound (even no bound) and the bounds
    C     must be specified for each independent variable. If an independent
    C     variable has both an upper and lower bound, the bounds must be
    C     consistent.  The lower bound must be .LE. the upper bound.
    C     Messages:
    C     DSPLP( ). INDEPENDENT VARIABLE (I1) IS NOT DEFINED.
    C               IN ABOVE MESSAGE, I1=         1
    C     ERROR NUMBER =        10
    C
    C     DSPLP( ).  LOWER BOUND (R1) AND UPPER BOUND (R2) FOR INDEP.
    C     VARIABLE (I1) ARE NOT CONSISTENT.
    C               IN ABOVE MESSAGE, I1=         1
    C               IN ABOVE MESSAGE, R1=    0.
    C               IN ABOVE MESSAGE, R2=    -.1000000000E+01
    C     ERROR NUMBER =        11
    C
    C -----12. and 13.  The type of bound (even no bound) and the bounds
    C     must be specified for each dependent variable.  If a dependent
    C     variable has both an upper and lower bound, the bounds must be
    C     consistent. The lower bound must be .LE. the upper bound.
    C     Messages:
    C     DSPLP( ). DEPENDENT VARIABLE (I1) IS NOT DEFINED.
    C               IN ABOVE MESSAGE, I1=         1
    C     ERROR NUMBER =        12
    C
    C     DSPLP( ).  LOWER BOUND (R1) AND UPPER BOUND (R2) FOR DEP.
    C      VARIABLE (I1) ARE NOT CONSISTENT.
    C               IN ABOVE MESSAGE, I1=         1
    C               IN ABOVE MESSAGE, R1=    0.
    C               IN ABOVE MESSAGE, R2=    -.1000000000E+01
    C     ERROR NUMBER =        13
    C
    C -----14. - 21.  These error messages can occur when processing the
    C     option array, PRGOPT(*), supplied by the user.  They would
    C     indicate a mistake in defining PRGOPT(*) or that data has been
    C     over-written.  See heading Usage of DSPLP( )
    C     Subprogram Options, for details on how to define PRGOPT(*).
    C     Messages:
    C     DSPLP( ). THE USER OPTION ARRAY HAS UNDEFINED DATA.
    C     ERROR NUMBER =        14
    C
    C     DSPLP( ). OPTION ARRAY PROCESSING IS CYCLING.
    C     ERROR NUMBER =        15
    C
    C     DSPLP( ). AN INDEX OF USER-SUPPLIED BASIS IS OUT OF RANGE.
    C     ERROR NUMBER =        16
    C
    C     DSPLP( ). SIZE PARAMETERS FOR MATRIX MUST BE SMALLEST AND LARGEST
    C     MAGNITUDES OF NONZERO ENTRIES.
    C     ERROR NUMBER =        17
    C
    C     DSPLP( ). THE NUMBER OF REVISED SIMPLEX STEPS BETWEEN CHECK-POINTS
    C     MUST BE POSITIVE.
    C     ERROR NUMBER =        18
    C
    C     DSPLP( ). FILE NUMBERS FOR SAVED DATA AND MATRIX PAGES MUST BE
    C     POSITIVE AND NOT EQUAL.
    C     ERROR NUMBER =        19
    C
    C     DSPLP( ). USER-DEFINED VALUE OF LAMAT (I1)
    C     MUST BE .GE. NVARS+7.
    C               IN ABOVE MESSAGE, I1=         1
    C     ERROR NUMBER =         20
    C
    C     DSPLP( ). USER-DEFINED VALUE OF LBM MUST BE .GE. 0.
    C     ERROR NUMBER =         21
    C
    C -----22.  The user-option, number 62, to check the size of the matrix
    C     data has been used.  An element of the matrix does not lie within
    C     the range of ASMALL and ABIG, parameters provided by the user.
    C     (See the heading: Usage of DSPLP( ) Subprogram Options,
    C     for details about this feature.)  Message:
    C     DSPLP( ). A MATRIX ELEMENT'S SIZE IS OUT OF THE SPECIFIED RANGE.
    C     ERROR NUMBER =        22
    C
    C -----23.  The user has provided an initial basis that is singular.
    C     In this case, the user can remedy this problem by letting
    C     subprogram DSPLP( ) choose its own initial basis.  Message:
    C     DSPLP( ). A SINGULAR INITIAL BASIS WAS ENCOUNTERED.
    C     ERROR NUMBER =         23
    C
    C -----24.  The user has provided an initial basis which is infeasible.
    C     The x and w values it defines do not satisfy A*x=w and the stated
    C     bounds.  In this case, the user can let subprogram DSPLP( )
    C     choose its own initial basis.  Message:
    C     DSPLP( ). AN INFEASIBLE INITIAL BASIS WAS ENCOUNTERED.
    C     ERROR NUMBER =        24
    C
    C -----25.Subprogram DSPLP( ) has completed the maximum specified number
    C     of iterations.  (The nominal maximum number is 3*(MRELAS+NVARS).)
    C     The results, necessary to continue on from
    C     this point, can be saved on Fortran unit 2 by activating option
    C     KEY=57.  If the user anticipates continuing the calculation, then
    C     the contents of Fortran unit 2 must be retained intact.  This
    C     is not done by subprogram DSPLP( ), so the user needs to save unit
    C     2 by using the appropriate system commands.  Message:
    C     DSPLP( ). MAX. ITERS. (I1) TAKEN. UP-TO-DATE RESULTS
    C     SAVED ON FILE (I2). IF(I2)=0, NO SAVE.
    C               IN ABOVE MESSAGE, I1=       500
    C               IN ABOVE MESSAGE, I2=         2
    C     ERROR NUMBER =        25
    C
    C -----26.  This error should never happen.  Message:
    C     DSPLP( ). MOVED TO A SINGULAR POINT. THIS SHOULD NOT HAPPEN.
    C     ERROR NUMBER =        26
    C
    C -----27.  The subprogram LA05A( ), which decomposes the basis matrix,
    C     has returned with an error flag (R1).  (See the document,
    C     "Fortran subprograms for handling sparse linear programming
    C     bases", AERE-R8269, J.K. Reid, Jan., 1976, H.M. Stationery Office,
    C     for an explanation of this error.)  Message:
    C     DSPLP( ). LA05A( ) RETURNED ERROR FLAG (R1) BELOW.
    C               IN ABOVE MESSAGE, R1=    -.5000000000E+01
    C     ERROR NUMBER =        27
    C
    C -----28.  The sparse linear solver package, LA05*( ), requires more
    C     space.  The value of LBM must be increased.  See the companion
    C     document, Usage of DSPLP( ) Subprogram Options, for details on how
    C     to increase the value of LBM.  Message:
    C     DSPLP( ). SHORT ON STORAGE FOR LA05*( ) PACKAGE. USE PRGOPT(*)
    C     TO GIVE MORE.
    C     ERROR NUMBER =        28
    C
    C -----29.  The row dimension of the two-dimensional Fortran array,
    C     the number of constraint equations (MRELAS), and the number
    C     of variables (NVARS), were not passed to the subprogram
    C     DFULMT().  See KEY = 68 for details.  Message:
    C     DFULMT() OF DSPLP() PACKAGE. ROW DIM., MRELAS, NVARS ARE
    C     MISSING FROM PRGOPT(*).
    C     ERROR NUMBER =        29
    C
    C     |-------------------------------------------------------|
    C     |End of List of DSPLP( ) Error and Diagnostic Messages. |
    C     |-------------------------------------------------------|
    C***REFERENCES  R. J. Hanson and K. L. Hiebert, A sparse linear
    C                 programming subprogram, Report SAND81-0297, Sandia
    C                 National Laboratories, 1981.
    C***ROUTINES CALLED  DPLPMN, XERMSG
    C***REVISION HISTORY  (YYMMDD)
    C   811215  DATE WRITTEN
    C   890531  Changed all specific intrinsics to generic.  (WRB)
    C   890605  Corrected references to XERRWV.  (WRB)
    C   890605  Removed unreferenced labels.  (WRB)
    C   891006  Cosmetic changes to prologue.  (WRB)
    C   891006  REVISION DATE from Version 3.2
    C   891214  Prologue converted to Version 4.0 format.  (BAB)
    C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    C   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
    C   920501  Reformatted the REFERENCES section.  (WRB)
    C***END PROLOGUE  DSPLP