 
C           RICHARD J. HANSON, DIV. 1642
C           SANDIA NATIONAL LABS.
C           ALBUQUERQUE, NEW MEXICO 87185
C     SUBROUTINE SPLP(USRMAT,MRELAS,NVARS,COSTS,PRGOPT,DATTRV,
C    *           BL,BU,IND,INFO,PRIMAL,DUALS,IBASIS,WORK,LW,IWORK,LIW)
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.)
C
C   \------------\
C   \INTRODUCTION\
C   \------------\
C     THE SUBPROGRAM SPLP( ) 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 SPLP( ).
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 USRMAT
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 USRMAT( ).
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     USRMAT
C     ------
C     THIS IS THE NAME OF A SPECIFIC SUBPROGRAM IN THE SPLP( ) PACKAGE
C     USED TO DEFINE THE MATRIX A.  IN THIS USAGE MODE OF SPLP( )
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 USRMAT 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 SPLP( ).
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.E0.
C
C     FOR USERS WITH MORE SOPHISTICATED NEEDS, SPLP( ) 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 SPLP( ) 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 DICUSSED 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 25.
C     55        A COMPUTATION, PARTIALLY COMPLETED, IS
C               BEING CONTINUED.  READ THE UP-TO-DATE
C               PARTIAL RESULTS FROM UNIT NUMBER 26.
C     56        REDEFINE THE UNIT NUMBER WHERE THE PARTIAL RESULTS
C               ARE STORED.  NORMALLY, THE UNIT NUMBER IS 26.
C     57        SAVE PARTIAL RESULTS ON UNIT 26 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 SPLP( ) WITH A STARTING (FEASIBLE)
C               NONSINGULAR BASIS.  NORMALLY, SPLP( ) 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, SPLP( ) 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, SPLP( ) 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
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 SPLP( ) AND ITS
C     SUBPROGRAMS.  THE LENGTHS OF THESE ARRAYS ARE RESPECTIVELY
C     LW AND LIW.  THESE PARAMETERS MUST SATISFY THE INEQUALITES
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 25 IS USED BY SPLP( ) 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 25
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 SPLP( ) 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 SPLP( ) 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     \\
C     \SAMPLE PROBLEM:\
C     \\
C
      PROGRAM LPEX
C     PROGRAM LPEX(OUTPUT,TAPE6=OUTPUT)
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 SPLP( ) TO SOLVE THIS LINEAR OPTIMIZATION
C     PROBLEM.  EACH REQUIRED STEP WILL BE SHOWN IN THIS EXAMPLE.
      DIMENSION COSTS(03),PRGOPT(04),DATTRV(18),BL(06),BU(06),IND(06),
     *PRIMAL(06),DUALS(06),IBASIS(06),WORK(079),IWORK(103)
C
      EXTERNAL USRMAT
      MRELAS=3
      NVARS=3
C
C     DEFINE THE ARRAY COSTS(*) FOR THE OBJECTIVE FUNCTION.
      COSTS(01)=1.
      COSTS(02)=1.
      COSTS(03)=1.
C
C     PLACE THE NONZERO INFORMATION ABOUT THE MATRIX IN DATTRV(*).
      DATTRV(01)=-1
      DATTRV(02)=1
      DATTRV(03)=1.
      DATTRV(04)=2
      DATTRV(05)=1.
C
      DATTRV(06)=-2
      DATTRV(07)=1
      DATTRV(08)=-3.
      DATTRV(09)=2
      DATTRV(10)=-2.
      DATTRV(11)=3
      DATTRV(12)=2.
C
      DATTRV(13)=-3
      DATTRV(14)=1
      DATTRV(15)=4.
      DATTRV(16)=3
      DATTRV(17)=-1.
C
      DATTRV(18)=0
C
C     CONSTRAIN X1,X2 TO BE NONNEGATIVE. LET X3 HAVE NO BOUNDS.
      BL(1)=0.
      IND(1)=1
      BL(2)=0.
      IND(2)=1
      IND(3)=4
C
C     CONSTRAIN W1=5,W2.LE.3, AND W3.GE.4.
      BL(4)=5.
      BU(4)=5.
      IND(4)=3
      BU(5)=3.
      IND(5)=2
      BL(6)=4.
      IND(6)=1
C
C     INDICATE THAT NO ADDITIONAL OPTIONS ARE IN USE.
      PRGOPT(01)=4
      PRGOPT(02)=51
      PRGOPT(03)=1
      PRGOPT(04)=1
C
C     DEFINE THE WORKING ARRAY LENGTHS.
      LW=079
      LIW=103
      CALL SPLP(USRMAT,MRELAS,NVARS,COSTS,PRGOPT,DATTRV,
     *BL,BU,IND,INFO,PRIMAL,DUALS,IBASIS,WORK,LW,IWORK,LIW)
      STOP
      END
 
 
