Logo Search packages:      
Sourcecode: fgarch version File versions  Download package

GarchModellingB.f

C SQP

C PART I:    QSQP.F
C PART II:   MQSUBS.F
C PART III:  PQSUBS.F  



C ##############################################################################
C PART I: QSQP.F


************************************************************************
* SUBROUTINE PSQPN              ALL SYSTEMS                   97/01/22
* PURPOSE :
* EASY TO USE SUBROUTINE FOR GENERAL NONLINEAR PROGRAMMING PROBLEMS.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NB  CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED.
*         NB>0-SIMPLE BOUNDS ACCEPTED.
*  II  NC  NUMBER OF LINEAR CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE
*         X(I) IS UNBOUNDED. IX(I)=1-LOVER BOUND XL(I).LE.X(I).
*         IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND
*         XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  CF(NC+1)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*         IC(KC)=0-CONSTRAINT CF(KC) IS NOT USED. IC(KC)=1-LOVER
*         CONSTRAINT CL(KC).LE.CF(KC). IC(KC)=2-UPPER CONSTRAINT
*         CF(KC).LE.CU(KC). IC(KC)=3-TWO SIDE CONSTRAINT
*         CL(KC).LE.CF(KC).LE.CU(KC). IC(KC)=5-EQUALITY CONSTRAINT
*         CF(KC).EQ.CL(KC).
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  IA  IA(NF)  AUXILIARY ARRAY.
*  RA  RA((NF+NC+7)*NF+3*NC+1)  AUXILIARY ARRAY.
*  II  IPAR(7)  INTEGER PAREMETERS:
*      IPAR(1)  PRINT SPECIFICATION. IPAR(1)=0-NO PRINT.
*         ABS(IPAR(1))=1-PRINT OF FINAL RESULTS.
*         ABS(IPAR(1))=2-PRINT OF FINAL RESULTS AND ITERATIONS.
*         IPAR(1)>0-BASIC FINAL RESULTS. IPAR(1)<0-EXTENDED FINAL
*         RESULTS.
*      IPAR(2)  MAXIMUM NUMBER OF ITERATIONS.
*      IPAR(3)  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
*      IPAR(4)  SCALING OF THE BFGS UPDATE. IPAR(4)=1-NO SCALING.
*         IPAR(4)=2-SCALING IN THE FIRST ITERATION.
*         IPAR(4)=3-CONTROLLED SCALING.
*      IPAR(5)  CORRECTION OF THE BFGS UPDATE IF A NEGATIVE CURVATURE
*         OCCURS. IPAR(5)=1-NO CORRECTION. IPAR(5)=2-POWELL'S
*         CORRECTION.
*      IPAR(6)  RESTART AFTER UNSUCCESSFUL UPDATE. IPAR(6)=0-RESTART
*         SUPPRESSED. IPAR(6)=1-RESTART PERFORMED.
*      IPAR(7)  INTERPOLATION IN LINE SEARCH. IPAR(7)=1-BISECTION.
*         IPAR(7)=2-TWO POINT QUADRATIC INTERPOLATION. IPAR(7)=3-THREE
*         POINT QUADRATIC INTERPOLATION. IPAR(7)=4-THREE POINT CUBIC
*         INTERPOLATION.
*  RI  RPAR(7)  REAL PARAMETERS:
*      RPAR(1)  MAXIMUM STEPSIZE.
*      RPAR(2)  TOLERANCE FOR CHANGE OF VARIABLES.
*      RPAR(3)  TOLERANCE FOR CONSTRAINT VIOLATIONS.
*      RPAR(4)  TOLERANCE FOR THE GRADIENT OF THE LAGRANGIAN FUNCTION.
*      RPAR(5)  TOLERANCE FOR A DESCENT DIRECTION.
*      RPAR(6)  TOLERANCE FOR A FUNCTION DECREASE IN THE LINE SEARCH.
*      RPAR(7)  PENALTY COEFFICIENT.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE OF THE LAGRANGIAN FUNCTION.
*  RO  CMAX  MAXIMUM CONSTRAINT VIOLATION.
*  IO  ITERM  CAUSE OF TERMINATION.
*
* VARIABLES IN COMMON /STATSQP/ (STATISTICS) :
*  IO  NRES  NUMBER OF RESTARTS.
*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITION.
*  IO  NREM  NUMBER OF CONSTRAINT DELETIONS.
*  IO  NADD  NUMBER OF CONSTRAINT ADDITIONS.
*  IO  NIT  NUMBER OF ITERATIONS.
*  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
*  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
*  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
*
* SUBPROGRAMS USED :
*  S   PSQP  RECURSIVE QUADRATIC PROGRAMMING METHOD WITH THE BFGS
*         VARIABLE METRIC UPDATE.
*
* EXTERNAL SUBROUTINES :
*  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS A VECTOR OF VARIABLES AND FF IS THE
*         VALUE OF THE OBJECTIVE FUNCTION.
*  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS A VECTOR OF VARIABLES AND GC(NF) IS
*         THE GRADIENT OF THE OBJECTIVE FUNCTION.
*  SE  CON  COMPUTATION OF THE VALUE OF THE CONSTRAINT FUNCTION.
*         CALLING SEQUENCE: CALL CON(NF,KC,X,FC) WHERE NF IS THE
*         NUMBER OF VARIABLES, KC IS THE INDEX OF THE CONSTRAINT
*         FUNCTION, X(NF) IS A VECTOR OF VARIABLES AND FC IS THE
*         VALUE OF THE CONSTRAINT FUNCTION.
*  SE  DCON  COMPUTATION OF THE GRADIENT OF THE CONSTRAINT FUNCTION.
*         CALLING SEQUENCE: CALL DCON(NF,KC,X,GC) WHERE NF IS THE
*         NUMBER OF VARIABLES, KC IS THE INDEX OF THE CONSTRAINT
*         FUNCTION, X(NF) IS A VECTOR OF VARIABLES AND GC(NF) IS THE
*         GRADIENT OF THE CONSTRAINT FUNCTION.
*
      SUBROUTINE PSQPN(NF,NB,NC,X,IX,XL,XU,CF,IC,CL,CU,IA,RA,IPAR,RPAR,
     +                 F,GMAX,CMAX,ITERM)
*
*     POINTERS FOR AUXILIARY ARRAYS
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION F,CMAX,GMAX
      INTEGER          ITERM,NB,NC,NF
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION CF(*),CL(*),CU(*),RA(*),RPAR(7),X(*),XL(*),XU(*)
      INTEGER          IA(*),IC(*),IPAR(7),IX(*)
*     ..
*     .. Scalars in Common ..
      INTEGER          NADD,NDEC,NFG,NFH,NFV,NIT,NREM,NRES
*     ..
*     .. Local Scalars ..
      INTEGER          LCFD,LCFO,LCG,LCP,LCR,LCZ,LG,LGC,LGF,LGO,LH,LIA,
     +                 LS,LXO
*     ..
*     .. External Subroutines ..
      EXTERNAL         PSQP
*     ..
*     .. Common blocks ..
      COMMON           /STATSQP/NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
*     ..
      LCG = 1
      LCFO = LCG + NF*NC
      LCFD = LCFO + NC + 1
      LGC = LCFD + NC
      LCR = LGC + NF
      LCZ = LCR + NF*(NF+1)/2
      LCP = LCZ + NF
      LGF = LCP + NC
      LG = LGF + NF
      LH = LG + NF
      LS = LH + NF* (NF+1)/2
      LXO = LS + NF
      LGO = LXO + NF
      LIA = 1
      CALL PSQP(NF,NB,NC,X,IX,XL,XU,CF,IC,CL,CU,RA,RA(LCFO),RA(LCFD),
     +          RA(LGC),IA,RA(LCR),RA(LCZ),RA(LCP),RA(LGF),RA(LG),
     +          RA(LH),RA(LS),RA(LXO),RA(LGO),RPAR(1),RPAR(2),RPAR(3),
     +          RPAR(4),RPAR(5),RPAR(6),RPAR(7),CMAX,GMAX,F,IPAR(1),
     +          IPAR(2),IPAR(3),IPAR(4),IPAR(5),IPAR(6),IPAR(7),ITERM)
      RETURN

      END
************************************************************************
* SUBROUTINE PSQP               ALL SYSTEMS                   97/01/22
* PURPOSE :
* RECURSIVE QUADRATIC PROGRAMMING METHOD WITH THE BFGS VARIABLE METRIC
* UPDATE FOR GENERAL NONLINEAR PROGRAMMING PROBLEMS.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NB  CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED.
*         NB>0-SIMPLE BOUNDS ACCEPTED.
*  II  NC  NUMBER OF LINEAR CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE
*         X(I) IS UNBOUNDED. IX(I)=1-LOVER BOUND XL(I).LE.X(I).
*         IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND
*         XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RO  CF(NC+1)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*         IC(KC)=0-CONSTRAINT CF(KC) IS NOT USED. IC(KC)=1-LOVER
*         CONSTRAINT CL(KC).LE.CF(KC). IC(KC)=2-UPPER CONSTRAINT
*         CF(KC).LE.CU(KC). IC(KC)=3-TWO SIDE CONSTRAINT
*         CL(KC).LE.CF(KC).LE.CU(KC). IC(KC)=5-EQUALITY CONSTRAINT
*         CF(KC).EQ.CL(KC).
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RA  CFO(NA)  VECTOR CONTAINING SAVED VALUES OF THE CONSTRAINT
*         FUNCTIONS.
*  RA  CFD(NA)  VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT
*         FUNCTIONS.
*  RA  GC(NF)  GRADIENT OF THE SELECTED CONSTRAINT FUNCTION.
*  IO  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RO  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RO  CZ(NF)  VECTOR OF LAGRANGE MULTIPLIERS.
*  RO  GF(NF)  GRADIENT OF THE MODEL FUNCTION.
*  RO  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RU  H(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OR INVERSION OF THE
*         HESSIAN MATRIX APPROXIMATION.
*  RO  S(NF)  DIRECTION VECTOR.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RI  GO(NF)  GRADIENTS DIFFERENCE.
*  RI  XMAX  MAXIMUM STEPSIZE.
*  RI  TOLX  TOLERANCE FOR CHANGE OF VARIABLES.
*  RI  TOLC  TOLERANCE FOR CONSTRAINT VIOLATIONS.
*  RI  TOLG  TOLERANCE FOR THE GRADIENT OF THE LAGRANGIAN FUNCTION.
*  RI  TOLD  TOLERANCE FOR A DESCENT DIRECTION.
*  RI  TOLS  TOLERANCE FOR A FUNCTION DECREASE IN THE LINE SEARCH.
*  RI  RPF  PENALTY COEFFICIENT.
*  RO  CMAX  MAXIMUM CONSTRAINT VIOLATION.
*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE OF THE LAGRANGIAN FUNCTION.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*         FUNCTIONS.
*  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
*         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
*         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
*         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
*         RESULTS.
*  II  MIT  MAXIMUN NUMBER OF ITERATIONS.
*  II  MFV  MAXIMUN NUMBER OF FUNCTION EVALUATIONS.
*  II  IEXT  TYPE OF OBJECTIVE FUNCTION. IEXT<0-MAXIMUM OF POSITIVE
*         VALUES. IEXT=0-MAXIMUM OF ABSOLUTE VALUES. IEXT>0-MAXIMUM
*         OF NEGATIVE VALUES.
*  II  MET  SELECTION OF SELF SCALING. MET=1-SELF SCALING SUPPRESSED.
*         MET=2 INITIAL SELF SCALING. MET=3 CONTROLLED SELF SCALING.
*  II  MEC  CORRECTION IF THE NEGATIVE CURVATURE OCCURS.
*         MEC=1-CORRECTION SUPPRESSED. MEC=2-POWELL'S CORRECTION.
*  II  MER  RESTART AFTER UNSUCCESSFUL UPDATE. MER=0-RESTART
*         SUPPRESSED. MER=1-RESTART PERFORMED.
*  II  MES  INTERPOLATION METHOD SELECTION. MES=1-BISECTION. MES=2-TWO
*         POINT QUADRATIC INTERPOLATION. MES=3-THREE POINT QUADRATIC
*         INTERPOLATION. MES=4-THREE POINT CUBIC INTERPOLATION.
*  IO  ITERM  CAUSE OF TERMINATION.
*
* VARIABLES IN COMMON /STATSQP/ (STATISTICS) :
*  IO  NRES  NUMBER OF RESTARTS.
*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITION.
*  IO  NREM  NUMBER OF CONSTRAINT DELETIONS.
*  IO  NADD  NUMBER OF CONSTRAINT ADDITIONS.
*  IO  NIT  NUMBER OF ITERATIONS.
*  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
*  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
*  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
*
* SUBPROGRAMS USED :
*  S   PC1F01  COMPUTATION OF THE VALUE AND THE GRADIENT OF THE
*         CONSTRAINT FUNCTION.
*  S   PF1F01  COMPUTATION OF THE VALUE AND THE GRADIENT OF THE
*         OBJECTIVE FUNCTION.
*  S   PLQDB1  GENERAL QUADRATIC PROGRAMMING SUBROUTINE BASED ON THE
*         GOLDFARB-IDNANI DUAL METHOD.
*  S   PLNEWS  IDENTIFICATION OF ACTIVE SIMPLE BOUNDS.
*  S   PLREDL  TRANSFORMATION OF THE INCOMPATIBLE QUADRATIC PROGRAMMING
*         SUBPROBLEM.
*  S   PP0AF8  COMPUTATION OF VALUE OF THE AUGMENTED LAGRANGIAN
*         FUNCTION.
*  S   PPSET2  COMPUTATION OF THE NEW PENALTY PARAMETERS.
*  S   PS0L02  LINE SEARCH USING ONLY FUNCTION VALUES.
*  S   PYTRND  DETERMINATION OF DIFFERENCES FOR VARIABLE METRIC
*         UPDATES.
*  S   PUDBG1  VARIABLE METRIC UPDATE AFTER GILL-MURRAY DECOMPOSITION.
*  S   MXDSMI  SYMMETRIC MATRIX IS REPLACED BY THE UNIT MATRIX.
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVINA  ABSOLUTE VALUES OF ELEMENTS OF AN INTEGER VECTOR.
*  RF  MXVMAX  L-INFINITY NORM OF A VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
* EXTERNAL SUBROUTINES :
*  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS A VECTOR OF VARIABLES AND FF IS THE
*         VALUE OF THE OBJECTIVE FUNCTION.
*  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION.
*         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER
*         OF VARIABLES, X(NF) IS A VECTOR OF VARIABLES AND GC(NF) IS
*         THE GRADIENT OF THE OBJECTIVE FUNCTION.
*  SE  CON  COMPUTATION OF THE VALUE OF THE CONSTRAINT FUNCTION.
*         CALLING SEQUENCE: CALL CON(NF,KC,X,FC) WHERE NF IS THE
*         NUMBER OF VARIABLES, KC IS THE INDEX OF THE CONSTRAINT
*         FUNCTION, X(NF) IS A VECTOR OF VARIABLES AND FC IS THE
*         VALUE OF THE CONSTRAINT FUNCTION.
*  SE  DCON  COMPUTATION OF THE GRADIENT OF THE CONSTRAINT FUNCTION.
*         CALLING SEQUENCE: CALL DCON(NF,KC,X,GC) WHERE NF IS THE
*         NUMBER OF VARIABLES, KC IS THE INDEX OF THE CONSTRAINT
*         FUNCTION, X(NF) IS A VECTOR OF VARIABLES AND GC(NF) IS THE
*         GRADIENT OF THE CONSTRAINT FUNCTION.
*
* METHOD :
* RECURSIVE QUADRATIC PROGRAMMING METHOD WITH THE BFGS VARIABLE METRIC
* UPDATE.
*
      SUBROUTINE PSQP(NF,NB,NC,X,IX,XL,XU,CF,IC,CL,CU,CG,CFO,CFD,GC,
     +                ICA,CR,CZ,CP,GF,G,H,S,XO,GO,XMAX,TOLX,TOLC,TOLG,
     +                TOLD,TOLS,RPF,CMAX,GMAX,F,IPRNT,MIT,MFV,MET,MEC,
     +                MER,MES,ITERM)
*     .. Scalar Arguments ..
      DOUBLE PRECISION F,CMAX,GMAX,RPF,TOLC,TOLD,TOLG,TOLS,TOLX,XMAX
      INTEGER          IPRNT,ITERM,MET,MEC,MER,MES,MFV,MIT,NB,NC,NF
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION CF(*),CFD(*),CFO(*),CG(*),CL(*),CP(*),CR(*),
     +                 CZ(*),CU(*),G(*),GC(*),GF(*),GO(*),H(*),S(*),
     +                 X(*),XL(*),XO(*),XU(*)
      INTEGER          IC(*),ICA(*),IX(*)
*     ..
*     .. Scalars in Common ..
      INTEGER          NADD,NDEC,NFG,NFH,NFV,NIT,NREM,NRES
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION ALF1,ALF2,CMAXO,DMAX,EPS7,EPS9,ETA0,ETA2,ETA9,
     +                 FMAX,FMIN,FO,GNORM,P,PO,R,RMAX,RMIN,RO,SNORM,
     +                 TOLB,TOLF,UMAX,RP,FP,PP,FF,FC
      INTEGER          I,IDECF,IEXT,IREST,ITERD,ITERL,ITERH,ITERQ,ITERS,
     +                 KBC,KBF,KC,KD,KIT,LD,MRED,MTESF,MTESX,N,NTESF,K,
     +                 NTESX,IEST,INITS,KTERS,MAXST,ISYS,MFP,NRED,IPOM,
     +                 LDS
*     ..
*     .. External Functions ..
      DOUBLE PRECISION MXVDOT, MXVMAX
      EXTERNAL         MXVDOT, MXVMAX
*     ..
*     .. External Subroutines ..
      EXTERNAL         PC1F01,PF1F01,PLNEWS,PLQDB1,PLREDL,PP0AF8,PPSET2,
     +                 PS0L02,PUDBG1,PYTRND,MXDSMI,MXVCOP,MXVINA,MXVSET
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC        ABS,MAX,MIN,SQRT
*     ..
*     .. Common blocks ..
      COMMON           /STATSQP/NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
*     ..
CR      IF (ABS(IPRNT).GT.1) WRITE (6,FMT='(1X,''ENTRY TO PSQP :'')')
*
*     INITIATION
*
      KBF = 0
      KBC = 0
      IF (NB.GT.0) KBF = 2
      IF (NC.GT.0) KBC = 2
      NRES = 0
      NDEC = 0
      NREM = 0
      NADD = 0
      NIT = 0
      NFV = 0
      NFG = 0
      NFH = 0
      ISYS = 0
      IEST = 0
      IEXT = 0
      MTESX = 2
      MTESF = 2
      INITS = 1
      ITERM = 0
      ITERS = 0
      ITERD = 0
      ITERQ = 0
      MRED = 20
      IREST = 1
      ITERS = 2
      KTERS = 5
      IDECF = 1
      ETA0 = 1.0D-15
      ETA2 = 1.0D-15
      ETA9 = 1.0D60
      EPS7 = 1.0D-15
      EPS9 = 1.0D-8
      ALF1 = 1.0D-10
      ALF2 = 1.0D10
      FMAX = 1.0D60
      FMIN = -FMAX
      TOLB = -FMAX
      DMAX = ETA9
      TOLF = 1.0D-16
      IF (XMAX.LE.0.0D0) XMAX = 1.0D+3
      IF (TOLX.LE.0.0D0) TOLX = 1.0D-16
      IF (TOLG.LE.0.0D0) TOLG = 1.0D-6
      IF (TOLC.LE.0.0D0) TOLC = 1.0D-6
      IF (TOLD.LE.0.0D0) TOLD = 1.0D-8
      IF (TOLS.LE.0.0D0) TOLS = 1.0D-4
      IF (RPF.LE.0.0D0) RPF = 1.0D-4
      IF (MET.LE.0) MET = 2
      IF (MEC.LE.0) MEC = 2
      IF (MES.LE.0) MES = 1
      IF (MIT.LE.0) MIT = 1000
      IF (MFV.LE.0) MFV = 2000
      KD = 1
      LD = -1
      KIT = 0
      CALL MXVSET(NC,0.0D0,CP)
*
*     INITIAL OPERATIONS WITH SIMPLE BOUNDS
*
      IF (KBF.GT.0) THEN
          DO 20 I = 1,NF
              IF ((IX(I).EQ.3.OR.IX(I).EQ.4) .AND. XU(I).LE.XL(I)) THEN
                  XU(I) = XL(I)
                  IX(I) = 5

              ELSE IF (IX(I).EQ.5 .OR. IX(I).EQ.6) THEN
                  XL(I) = X(I)
                  XU(I) = X(I)
                  IX(I) = 5
              END IF

              IF (IX(I).EQ.1 .OR. IX(I).EQ.3) X(I) = MAX(X(I),XL(I))
              IF (IX(I).EQ.2 .OR. IX(I).EQ.3) X(I) = MIN(X(I),XU(I))
   20     CONTINUE
      END IF
*
*     INITIAL OPERATIONS WITH GENERAL CONSTRAINTS
*
      IF (KBC.GT.0) THEN
          K = 0
          DO 30 KC = 1,NC
              IF ((IC(KC).EQ.3.OR.IC(KC).EQ.4) .AND.
     +            CU(KC).LE.CL(KC)) THEN
                  CU(KC) = CL(KC)
                  IC(KC) = 5

              ELSE IF (IC(KC).EQ.5 .OR. IC(KC).EQ.6) THEN
                  CU(KC) = CL(KC)
                  IC(KC) = 5
              END IF

              K = K + NF
   30     CONTINUE
      END IF

      IF (KBF.GT.0) THEN
          DO 31 I = 1,NF
              IF (IX(I).GE.5) IX(I) = -IX(I)
              IF (IX(I).LE.0) THEN

              ELSE IF ((IX(I).EQ.1.OR.IX(I).EQ.3).AND.X(I).LE.XL(I))
     +                THEN
                  X(I) = XL(I)

              ELSE IF ((IX(I).EQ.2.OR.IX(I).EQ.3).AND.X(I).GE.XU(I))
     +                THEN
                  X(I) = XU(I)
              END IF

              CALL PLNEWS(X,IX,XL,XU,EPS9,I,ITERL)
              IF (IX(I).GT.10) IX(I) = 10 - IX(I)
   31     CONTINUE
      END IF

      FO = FMIN
      GMAX = ETA9
      DMAX = ETA9
   40 CONTINUE
      LDS=LD
      CALL PF1F01(NF,X,GF,GF,FF,F,KD,LD,IEXT)
      LD=LDS
      CALL PC1F01(NF,NC,X,FC,CF,CL,CU,IC,GC,CG,CMAX,KD,LD)
      CF(NC+1)=F
CR      IF (ABS(IPRNT).GT.1) WRITE (6,FMT=
CR     +'(1X,''NIT='',I4,2X,''NFV='',I4,2X,''NFG='',I4,2X,''F='',
CR     +  D12.6,2X,''C='',D8.2,2X,''G='',D8.2)')
CR     +  NIT,NFV,NFG,F,CMAX,GMAX

      if (iprnt.eq.1) then
            call dblepr("LLH improved to:", -1, f, 1)
            call dblepr("With X:", -1, x, nf)
            
      endif
*
*     START OF THE ITERATION WITH TESTS FOR TERMINATION.
*
      IF (ITERM.LT.0) GO TO 80
      IF (ITERS.EQ.0) GO TO 50
      IF (F.LE.TOLB) THEN
          ITERM = 3
          GO TO 80

      END IF

      IF (DMAX.LE.TOLX) THEN
          ITERM = 1
          NTESX = NTESX + 1
          IF (NTESX.GE.MTESX) GO TO 80

      ELSE
          NTESX = 0
      END IF

   50 IF (NIT.GE.MIT) THEN
          ITERM = 12
          GO TO 80

      END IF

      IF (NFV.GE.MFV) THEN
          ITERM = 11
          GO TO 80

      END IF

      ITERM = 0
      NIT = NIT + 1
   60 CONTINUE
*
*     RESTART
*
      N = NF
      IF (IREST.GT.0) THEN
          CALL MXDSMI(N,H)
          LD = MIN(LD,1)
          IDECF = 1
          IF (KIT.LT.NIT) THEN
              NRES = NRES + 1
              KIT = NIT

          ELSE
              ITERM = -10
              IF (ITERS.LT.0) ITERM = ITERS - 5
              GO TO 80

          END IF

      END IF
*
*     DIRECTION DETERMINATION USING A QUADRATIC PROGRAMMING PROCEDURE
*
      CALL MXVCOP(NC+1,CF,CFO)
      MFP=2
      IPOM=0
   61 CONTINUE
      CALL PLQDB1(NF,NC,X,IX,XL,XU,CF,CFD,IC,ICA,CL,CU,CG,CR,CZ,G,GF,
     &            H,S,MFP,KBF,KBC,IDECF,ETA2,ETA9,EPS7,EPS9,UMAX,GMAX,
     &            N,ITERQ)
      IF (ITERQ.LT.0) THEN
          IF (IPOM.LT.10) THEN
              IPOM=IPOM+1
              CALL PLREDL(NC,CF,IC,CL,CU,KBC)
              GO TO 61
          END IF
          ITERD=ITERQ-10
          GO TO 62
      END IF
      IPOM=0
      ITERD=1
      GMAX=MXVMAX(NF,G)
      GNORM=SQRT(MXVDOT(NF,G,G))
      SNORM=SQRT(MXVDOT(NF,S,S))
   62 CONTINUE
      IF (ITERD.LT.0) ITERM = ITERD
      IF (ITERM.NE.0) GO TO 80
      CALL MXVCOP(NC+1,CFO,CF)
*
*     TEST FOR SUFFICIENT DESCENT
*
      P = MXVDOT(NF,G,S)
      IREST = 1
      IF (SNORM.LE.0.0D0) THEN

      ELSE IF (P+TOLD*GNORM*SNORM.LE.0.0D0) THEN
          IREST = 0
      END IF

      IF (IREST.EQ.0) THEN
          NRED = 0
          RMIN = ALF1*GNORM/SNORM
          RMAX = MIN(ALF2*GNORM/SNORM,XMAX/SNORM)

      ELSE
          GO TO 60

      END IF

      IF (GMAX.LE.TOLG.AND.CMAX.LE.TOLC) THEN
          ITERM=4
          GO TO 80
      ENDIF
      CALL PPSET2(NF,N,NC,ICA,CZ,CP)
      CALL MXVINA(NC,IC)
      CALL PP0AF8(NF,N,NC,CF,IC,ICA,CL,CU,CZ,RPF,FC,F)
*
*     PREPARATION OF LINE SEARCH
*
      RO = 0.0D0
      FO = F
      PO = P
      CMAXO = CMAX
      CALL MXVCOP(NF,X,XO)
      CALL MXVCOP(NF,G,GO)
      CALL MXVCOP(NF,GF,CR)
      CALL MXVCOP(NC+1,CF,CFO)
*
*     LINE SEARCH WITHOUT DIRECTIONAL DERIVATIVES
*
11170 CONTINUE
      CALL PS0L02(R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX,
     & TOLS,KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS,
     & MES,ISYS)
      GOTO (11174,11172) ISYS+1
11172 CONTINUE
      CALL MXVDIR(NF,R,S,XO,X)
      LDS=LD
      CALL PF1F01(NF,X,GF,G,FF,F,KD,LD,IEXT)
      LD=LDS
      CALL PC1F01(NF,NC,X,FC,CF,CL,CU,IC,GC,CG,CMAX,KD,LD)
      CF(NC+1)=F
      CALL PP0AF8(NF,N,NC,CF,IC,ICA,CL,CU,CZ,RPF,FC,F)
      GOTO 11170
11174 CONTINUE
      KD=1
*
*     DECISION AFTER UNSUCCESSFUL LINE SEARCH
*
      IF (ITERS.LE.0) THEN
          R = 0.0D0
          F = FO
          P = PO
          CALL MXVCOP(NF,XO,X)
          CALL MXVCOP(NF,CR,GF)
          CALL MXVCOP(NC+1,CFO,CF)
          IREST = 1
          LD = KD
          GO TO 60
      END IF
*
*     COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE
*     FUNCTION TOGETHER WITH THE VALUES AND THE GRADIENTS OF THE
*     APPROXIMATED FUNCTIONS
*
      IF (KD.GT.LD) THEN
          LDS=LD
          CALL PF1F01(NF,X,GF,GF,FF,F,KD,LD,IEXT)
          LD=LDS
          CALL PC1F01(NF,NC,X,FC,CF,CL,CU,IC,GC,CG,CMAX,KD,LD)
      END IF
*
*     PREPARATION OF VARIABLE METRIC UPDATE
*
      CALL MXVCOP(NF,GF,G)
      CALL PYTRND(NF,N,X,XO,ICA,CG,CZ,G,GO,R,F,FO,P,PO,CMAX,CMAXO,
     +            DMAX,KD,LD,ITERS)
*
*     VARIABLE METRIC UPDATE
*
      CALL PUDBG1(N,H,G,S,XO,GO,R,PO,NIT,KIT,ITERH,MET,MEC)
      IF (MER.GT.0.AND.ITERH.GT.0) IREST=1
*
*     END OF THE ITERATION
*
      GO TO 40

   80 CONTINUE
CR      IF (IPRNT.GT.1 .OR. IPRNT.LT.0) WRITE (6,
CR     +    FMT='(1X,''EXIT FROM PSQP :'')')
CR      IF (IPRNT.NE.0) WRITE (6,FMT=
CR     +'(1X,''NIT='',I4,2X,''NFV='',I4,2X,''NFG='',I4,2X,''F='',
CR     +  D12.6,2X,''C='',D8.2,2X,''G='',D8.2,2X,''ITERM='',I3)')
CR     +  NIT,NFV,NFG,F,CMAX,GMAX,ITERM
CR      IF (IPRNT.LT.0) WRITE (6,FMT='(1X,''X ='',5D15.7:/(4X,5D15.7))')
CR     +    (X(I),I=1,NF)
      RETURN

      END

      
C ##############################################################################
C PART II: MQSUBS.F


* SUBROUTINE MXDCMM               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX A
* BY A VECTOR X.
*
* PARAMETERS :
*  II  N  NUMBER OF ROWS OF THE MATRIX A.
*  II  M  NUMBER OF COLUMNS OF THE MATRIX A.
*  RI  A(N*M)  RECTANGULAR MATRIX STORED COLUMNWISE IN THE
*         ONE-DIMENSIONAL ARRAY.
*  RI  X(M)  INPUT VECTOR.
*  RO  Y(N)  OUTPUT VECTOR EQUAL TO A*X.
*
* SUBPROGRAMS USED :
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE MXDCMM(N,M,A,X,Y)
      INTEGER N,M
      DOUBLE PRECISION A(*),X(*),Y(*)
      INTEGER J,K
      CALL MXVSET(N,0.0D 0,Y)
      K=0
      DO 1 J=1,M
      CALL MXVDIR(N,X(J),A(K+1),Y,Y)
      K=K+N
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXDPGB                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A DENSE SYMMETRIC
* POSITIVE DEFINITE MATRIX A+E USING THE FACTORIZATION A+E=L*D*TRANS(L)
* OBTAINED BY THE SUBROUTINE MXDPGF.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXDPGF.
*  RU  X(N)  ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR
*         EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR
*         EQUATIONS.
*  II  JOB  OPTION. IF JOB=0 THEN X:=(A+E)**(-1)*X. IF JOB>0 THEN
*         X:=L**(-1)*X. IF JOB<0 THEN X:=TRANS(L)**(-1)*X.
*
* METHOD :
* BACK SUBSTITUTION
*
      SUBROUTINE MXDPGB(N,A,X,JOB)
C     .. Scalar Arguments ..
      INTEGER JOB,N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(*),X(*)
C     ..
C     .. Local Scalars ..
      INTEGER I,II,IJ,J
C     ..
      IF (JOB.GE.0) THEN
*
*     PHASE 1 : X:=L**(-1)*X
*
          IJ = 0
          DO 20 I = 1,N
              DO 10 J = 1,I - 1
                  IJ = IJ + 1
                  X(I) = X(I) - A(IJ)*X(J)
   10         CONTINUE
              IJ = IJ + 1
   20     CONTINUE
      END IF

      IF (JOB.EQ.0) THEN
*
*     PHASE 2 : X:=D**(-1)*X
*
          II = 0
          DO 30 I = 1,N
              II = II + I
              X(I) = X(I)/A(II)
   30     CONTINUE
      END IF

      IF (JOB.LE.0) THEN
*
*     PHASE 3 : X:=TRANS(L)**(-1)*X
*
          II = N* (N-1)/2
          DO 50 I = N - 1,1,-1
              IJ = II
              DO 40 J = I + 1,N
                  IJ = IJ + J - 1
                  X(I) = X(I) - A(IJ)*X(J)
   40         CONTINUE
              II = II - I
   50     CONTINUE
      END IF

      RETURN

      END
* SUBROUTINE MXDPGF                ALL SYSTEMS                89/12/01
* PORTABILITY : ALL SYSTEMS
* 89/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* FACTORIZATION A+E=L*D*TRANS(L) OF A DENSE SYMMETRIC POSITIVE DEFINITE
* MATRIX A+E WHERE D AND E ARE DIAGONAL POSITIVE DEFINITE MATRICES AND
* L IS A LOWER TRIANGULAR MATRIX. IF A IS SUFFICIENTLY POSITIVE
* DEFINITE THEN E=0.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2)  ON INPUT A GIVEN DENSE SYMMETRIC (USUALLY POSITIVE
*         DEFINITE) MATRIX A STORED IN THE PACKED FORM. ON OUTPUT THE
*         COMPUTED FACTORIZATION A+E=L*D*TRANS(L).
*  IO  INF  AN INFORMATION OBTAINED IN THE FACTORIZATION PROCESS. IF
*         INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE AND E=0. IF
*         INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE AND E>0. IF
*         INF>0 THEN A IS INDEFINITE AND INF IS AN INDEX OF THE
*         MOST NEGATIVE DIAGONAL ELEMENT USED IN THE FACTORIZATION
*         PROCESS.
*  RU  ALF  ON INPUT A DESIRED TOLERANCE FOR POSITIVE DEFINITENESS. ON
*         OUTPUT THE MOST NEGATIVE DIAGONAL ELEMENT USED IN THE
*         FACTORIZATION PROCESS (IF INF>0).
*  RO  TAU  MAXIMUM DIAGONAL ELEMENT OF THE MATRIX E.
*
* METHOD :
* P.E.GILL, W.MURRAY : NEWTON TYPE METHODS FOR UNCONSTRAINED AND
* LINEARLY CONSTRAINED OPTIMIZATION, MATH. PROGRAMMING 28 (1974)
* PP. 311-350.
*
      SUBROUTINE MXDPGF(N,A,INF,ALF,TAU)
C     .. Scalar Arguments ..
      DOUBLE PRECISION ALF,TAU
      INTEGER INF,N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(*)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION BET,DEL,GAM,RHO,SIG,TOL
      INTEGER I,IJ,IK,J,K,KJ,KK,L
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,MAX
C     ..
      L = 0
      INF = 0
      TOL = ALF
*
*     ESTIMATION OF THE MATRIX NORM
*
      ALF = 0.0D0
      BET = 0.0D0
      GAM = 0.0D0
      TAU = 0.0D0
      KK = 0
      DO 20 K = 1,N
          KK = KK + K
          BET = MAX(BET,ABS(A(KK)))
          KJ = KK
          DO 10 J = K + 1,N
              KJ = KJ + J - 1
              GAM = MAX(GAM,ABS(A(KJ)))
   10     CONTINUE
   20 CONTINUE
      BET = MAX(TOL,BET,GAM/N)
*      DEL = TOL*BET
      DEL = TOL*MAX(BET,1.0D0)
      KK = 0
      DO 60 K = 1,N
          KK = KK + K
*
*     DETERMINATION OF A DIAGONAL CORRECTION
*
          SIG = A(KK)
          IF (ALF.GT.SIG) THEN
              ALF = SIG
              L = K
          END IF

          GAM = 0.0D0
          KJ = KK
          DO 30 J = K + 1,N
              KJ = KJ + J - 1
              GAM = MAX(GAM,ABS(A(KJ)))
   30     CONTINUE
          GAM = GAM*GAM
          RHO = MAX(ABS(SIG),GAM/BET,DEL)
          IF (TAU.LT.RHO-SIG) THEN
              TAU = RHO - SIG
              INF = -1
          END IF
*
*     GAUSSIAN ELIMINATION
*
          A(KK) = RHO
          KJ = KK
          DO 50 J = K + 1,N
              KJ = KJ + J - 1
              GAM = A(KJ)
              A(KJ) = GAM/RHO
              IK = KK
              IJ = KJ
              DO 40 I = K + 1,J
                  IK = IK + I - 1
                  IJ = IJ + 1
                  A(IJ) = A(IJ) - A(IK)*GAM
   40         CONTINUE
   50     CONTINUE
   60 CONTINUE
      IF (L.GT.0 .AND. ABS(ALF).GT.DEL) INF = L
      RETURN

      END
* FUNCTION MXDPGP                  ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE NUMBER MXDPGP=TRANS(X)*D**(-1)*Y WHERE D IS A
* DIAGONAL MATRIX IN THE FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
* SUBROUTINE MXDPGF.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXDPGF.
*  RI  X(N)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RR  MXDPGP  COMPUTED NUMBER MXDPGP=TRANS(X)*D**(-1)*Y.
*
      DOUBLE PRECISION FUNCTION MXDPGP(N,A,X,Y)
C     .. Scalar Arguments ..
      INTEGER N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(N* (N+1)/2),X(N),Y(N)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION TEMP
      INTEGER I,J
C     ..
      J = 0
      TEMP = 0.0D0
      DO 10 I = 1,N
          J = J + I
          TEMP = TEMP + X(I)*Y(I)/A(J)
   10 CONTINUE
      MXDPGP = TEMP
      RETURN

      END
* SUBROUTINE MXDPGS                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SCALING OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX A+E USING THE
* FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE SUBROUTINE MXDPGF.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXDPGF.
*  RI  ALF  SCALING FACTOR.
*
      SUBROUTINE MXDPGS(N,A,ALF)
C     .. Scalar Arguments ..
      DOUBLE PRECISION ALF
      INTEGER N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(*)
C     ..
C     .. Local Scalars ..
      INTEGER I,J
C     ..
      J = 0
      DO 10 I = 1,N
          J = J + I
          A(J) = A(J)*ALF
   10 CONTINUE
      RETURN

      END
* SUBROUTINE MXDPGU                ALL SYSTEMS                89/12/01
* PORTABILITY : ALL SYSTEMS
* 89/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* CORRECTION OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX A+E IN THE
* FACTORED FORM A+E=L*D*TRANS(L) OBTAINED BY THE SUBROUTINE MXDPGF.
* THE CORRECTION IS DEFINED AS A+E:=A+E+ALF*X*TRANS(X) WHERE ALF IS A
* GIVEN SCALING FACTOR AND X IS A GIVEN VECTOR.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
*         SUBROUTINE MXDPGF.
*  RI  ALF  SCALING FACTOR IN THE CORRECTION TERM.
*  RI  X(N)  VECTOR IN THE CORRECTION TERM.
*  RA  Y(N) AUXILIARY VECTOR.
*
* METHOD :
* P.E.GILL, W.MURRAY, M.SAUNDERS: METHODS FOR COMPUTING AND MODIFYING
* THE LDV FACTORS OF A MATRIX, MATH. OF COMP. 29 (1974) PP. 1051-1077.
*
      SUBROUTINE MXDPGU(N,A,ALF,X,Y)
C     .. Parameters ..
      DOUBLE PRECISION ZERO,ONE,FOUR,CON
      PARAMETER (ZERO=0.0D0,ONE=1.0D0,FOUR=4.0D0,CON=1.0D-8)
C     ..
C     .. Scalar Arguments ..
      DOUBLE PRECISION ALF
      INTEGER N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION A(*),X(*),Y(*)
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION B,D,P,R,T,TO
      INTEGER I,II,IJ,J
C     ..
C     .. External Subroutines ..
      EXTERNAL MXVSCL
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC SQRT
C     ..
      IF (ALF.GE.ZERO) THEN
*
*     FORWARD CORRECTION IN CASE WHEN THE SCALING FACTOR IS NONNEGATIVE
*
          ALF = SQRT(ALF)
          CALL MXVSCL(N,ALF,X,Y)
          TO = ONE
          II = 0
          DO 30 I = 1,N
              II = II + I
              D = A(II)
              P = Y(I)
              T = TO + P*P/D
              R = TO/T
              A(II) = D/R
              B = P/ (D*T)
              IF (A(II).LE.FOUR*D) THEN
*
*     AN EASY FORMULA FOR LIMITED DIAGONAL ELEMENT
*
                  IJ = II
                  DO 10 J = I + 1,N
                      IJ = IJ + J - 1
                      D = A(IJ)
                      Y(J) = Y(J) - P*D
                      A(IJ) = D + B*Y(J)
   10             CONTINUE

              ELSE
*
*     A MORE COMPLICATE BUT NUMERICALLY STABLE FORMULA FOR UNLIMITED
*     DIAGONAL ELEMENT
*
                  IJ = II
                  DO 20 J = I + 1,N
                      IJ = IJ + J - 1
                      D = A(IJ)
                      A(IJ) = R*D + B*Y(J)
                      Y(J) = Y(J) - P*D
   20             CONTINUE
              END IF

              TO = T
   30     CONTINUE

      ELSE
*
*     BACKWARD CORRECTION IN CASE WHEN THE SCALING FACTOR IS NEGATIVE
*
          ALF = SQRT(-ALF)
          CALL MXVSCL(N,ALF,X,Y)
          TO = ONE
          IJ = 0
          DO 50 I = 1,N
              D = Y(I)
              DO 40 J = 1,I - 1
                  IJ = IJ + 1
                  D = D - A(IJ)*Y(J)
   40         CONTINUE
              Y(I) = D
              IJ = IJ + 1
              TO = TO - D*D/A(IJ)
   50     CONTINUE
          IF (TO.LE.ZERO) TO = CON
          II = N* (N+1)/2
          DO 70 I = N,1,-1
              D = A(II)
              P = Y(I)
              T = TO + P*P/D
              A(II) = D*TO/T
              B = -P/ (D*TO)
              TO = T
              IJ = II
              DO 60 J = I + 1,N
                  IJ = IJ + J - 1
                  D = A(IJ)
                  A(IJ) = D + B*Y(J)
                  Y(J) = Y(J) + P*D
   60         CONTINUE
              II = II - I
   70     CONTINUE
      END IF

      RETURN

      END
* SUBROUTINE MXDPRB                ALL SYSTEMS                89/12/01
* PORTABILITY : ALL SYSTEMS
* 89/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A DENSE SYMMETRIC
* POSITIVE DEFINITE MATRIX A USING THE FACTORIZATION A=TRANS(R)*R.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2) FACTORIZATION A=TRANS(R)*R.
*  RU  X(N)  ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR
*         EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR
*         EQUATIONS.
*  II  JOB  OPTION. IF JOB=0 THEN X:=A**(-1)*X. IF JOB>0 THEN
*         X:=TRANS(R)**(-1)*X. IF JOB<0 THEN X:=R**(-1)*X.
*
* METHOD :
* BACK SUBSTITUTION
*
      SUBROUTINE MXDPRB(N,A,X,JOB)
*     .. Scalar Arguments ..
      INTEGER          JOB,N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(*),X(*)
*     ..
*     .. Local Scalars ..
      INTEGER          I,II,IJ,J
*     ..
      IF (JOB.GE.0) THEN
*
*     PHASE 1 : X:=TRANS(R)**(-1)*X
*
          IJ = 0
          DO 20 I = 1,N
              DO 10 J = 1,I - 1
                  IJ = IJ + 1
                  X(I) = X(I) - A(IJ)*X(J)
   10         CONTINUE
              IJ = IJ + 1
              X(I) = X(I)/A(IJ)
   20     CONTINUE
      END IF

      IF (JOB.LE.0) THEN
*
*     PHASE 2 : X:=R**(-1)*X
*
          II = N* (N+1)/2
          DO 40 I = N,1,-1
              IJ = II
              DO 30 J = I + 1,N
                  IJ = IJ + J - 1
                  X(I) = X(I) - A(IJ)*X(J)
   30         CONTINUE
              X(I) = X(I)/A(II)
              II = II - I
   40     CONTINUE
      END IF

      RETURN

      END
* SUBROUTINE MXDPRC                ALL SYSTEMS                92/12/01
* PORTABILITY : ALL SYSTEMS
* 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* CORRECTION OF A SINGULAR DENSE SYMMETRIC POSITIVE SEMIDEFINITE MATRIX
* A DECOMPOSED AS A=TRANS(R)*R.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  IO  INF  AN INFORMATION OBTAINED IN THE CORRECTION PROCESS. IF
*         INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE. IF
*         INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE.
*         PROCESS.
*  RI  TOL  DESIRED TOLERANCE FOR POSITIVE DEFINITENESS.
*
      SUBROUTINE MXDPRC(N,A,INF,TOL)
      INTEGER N
      REAL*8  A(N*(N+1)/2),TOL
      REAL*8 TOL1,TEMP
      INTEGER L,I
      REAL*8 ZERO
      PARAMETER (ZERO=0.0D 0)
      INF=0
      TOL1=SQRT(TOL)
      TEMP=TOL1
      DO 3 I=1,N*(N+1)/2
      TEMP=MAX(TEMP,ABS(A(I)))
    3 CONTINUE
      TEMP=TEMP*TOL1
      L=0
      DO 1 I=1,N
      L=L+I
      IF (ABS(A(L)).LE.TEMP) THEN
      A(L)=SIGN(TEMP,A(L))
      INF=-1
      ENDIF
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXDPRM                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A GIVEN VECTOR X BY A DENSE SYMMETRIC POSITIVE
* DEFINITE MATRIX A USING THE FACTORIZATION A=TRANS(R)*R.
*
* PARAMETERS :
*  II  N ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2) FACTORIZATION A=TRANS(R)*R.
*  RU  X(N)  ON INPUT THE GIVEN VECTOR. ON OUTPUT THE RESULT OF
*         MULTIPLICATION.
*  II  JOB  OPTION. IF JOB=0 THEN X:=A*X. IF JOB>0 THEN X:=R*X.
*         IF JOB<0 THEN X:=TRANS(R)*X.
*
      SUBROUTINE MXDPRM( N, A, X, JOB)
      INTEGER N, JOB
      REAL*8  A(N*(N+1)/2), X(N)
      INTEGER  I, J, II, IJ
      IF (JOB .GE. 0)  THEN
C
C     PHASE 1 : X:=R*X
C
      II = 0
      DO 3  I = 1, N
      II = II + I
      X(I) = A(II) * X(I)
      IJ = II
      DO 2  J = I+1, N
      IJ = IJ + J - 1
      X(I) = X(I) + A(IJ) * X(J)
    2 CONTINUE
    3 CONTINUE
      ENDIF
      IF (JOB .LE. 0)  THEN
C
C     PHASE 2 : X:=TRANS(R)*X
C
      IJ = N*(N+1)/2
      DO 6  I = N, 1, -1
      X(I) = A(IJ) * X(I)
      DO 5  J = I-1, 1, -1
      IJ = IJ - 1
      X(I) = X(I) + A(IJ) * X(J)
    5 CONTINUE
      IJ = IJ - 1
    6 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE MXDRGR               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* PLANE ROTATION IS APPLIED TO A ROWWISE STORED DENSE RECTANGULAR
* MATRIX A.
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX A.
*  RU  A(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
*         ONE-DIMENSIONAL ARRAY.
*  II  K  FIRST INDEX OF THE PLANE ROTATION.
*  II  L  SECOND INDEX OF THE PLANE ROTATION.
*  RI  CK  DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  RI  CL  OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  II  IER  TYPE OF THE PLANE ROTATION. IER=0-GENERAL PLANE ROTATION.
*         IER=1-PERMUTATION. IER=2-TRANSFORMATION SUPPRESSED.
*
* SUBPROGRAMS USED :
*  S   MXVROT  PLANE ROTATION APPLIED TO TWO ELEMENTS.
*
      SUBROUTINE MXDRGR(N,A,K,L,CK,CL,IER)
      INTEGER N,K,L,IER
      DOUBLE PRECISION A(*),CK,CL
      INTEGER I,IK,IL
      IF (IER.NE.0.AND.IER.NE.1) RETURN
      IK=(K-1)*N
      IL=(L-1)*N
      DO 1 I=1,N
      IK=IK+1
      IL=IL+1
      CALL MXVROT(A(IK),A(IL),CK,CL,IER)
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXDRMD               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A BY
* A VECTOR X AND ADDITION OF A SCALED VECTOR ALF*Y.
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX A.
*  II  M  NUMBER OF ROWS OF THE MATRIX A.
*  RI  A(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
*         ONE-DIMENSIONAL ARRAY.
*  RI  X(N)  INPUT VECTOR.
*  RI  ALF  SCALING FACTOR.
*  RI  Y(M)  INPUT VECTOR.
*  RO  Z(M)  OUTPUT VECTOR EQUAL TO A*X+ALF*Y.
*
      SUBROUTINE MXDRMD(N,M,A,X,ALF,Y,Z)
      INTEGER N,M
      REAL*8 A(M*N),X(N),ALF,Y(M),Z(M)
      REAL*8 TEMP
      INTEGER I,J,K
      K=0
      DO 2 J=1,M
      TEMP=ALF*Y(J)
      DO 1 I=1,N
      TEMP=TEMP+A(K+I)*X(I)
    1 CONTINUE
      Z(J)=TEMP
      K=K+N
    2 CONTINUE
      RETURN
      END
* SUBROUTINE MXDRMI               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* ROWWISE STORED DENSE RECTANGULAR MATRIX A IS SET TO BE A PART OF THE
* UNIT MATRIX.
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX A.
*  II  M  NUMBER OF ROWS OF THE MATRIX A.
*  RO  A(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE ONE-DIMENSIONAL
*          ARRAY. THIS MATRIX IS SET TO TRANS([I,0]).
*
      SUBROUTINE MXDRMI(N,M,A)
      INTEGER N,M
      REAL*8 A(M*N)
      INTEGER I,J,K
      REAL*8 ZERO,ONE
      PARAMETER (ZERO=0.0D 0,ONE=1.0D 0)
      K=0
      DO 2 J=1,M
      DO 1 I=1,N
      A(I+K)=ZERO
      IF (I.EQ.J) A(I+K)=ONE
    1 CONTINUE
      K=K+N
    2 CONTINUE
      RETURN
      END
* SUBROUTINE MXDRMM               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A BY
* A VECTOR X.
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX A.
*  II  M  NUMBER OF ROWS OF THE MATRIX A.
*  RI  A(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
*         ONE-DIMENSIONAL ARRAY.
*  RI  X(N)  INPUT VECTOR.
*  RO  Y(M)  OUTPUT VECTOR EQUAL TO A*X.
*
      SUBROUTINE MXDRMM(N,M,A,X,Y)
      INTEGER N,M
      DOUBLE PRECISION A(*),X(*),Y(*)
      DOUBLE PRECISION TEMP
      INTEGER I,J,K
      K=0
      DO 2 J=1,M
      TEMP=0.0D 0
      DO 1 I=1,N
      TEMP=TEMP+A(K+I)*X(I)
    1 CONTINUE
      Y(J)=TEMP
      K=K+N
    2 CONTINUE
      RETURN
      END
* FUNCTION  MXDRMN               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* EUCLIDEAN NORM OF A PART OF THE I-TH COLUMN OF A ROWWISE STORED DENSE
* RECTANGULAR MATRIX A IS COMPUTED.
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX A.
*  II  M  NUMBER OF ROWS OF THE MATRIX A.
*  RI  A(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
*         ONE-DIMENSIONAL ARRAY.
*  II  I  INDEX OF THE COLUMN WHOSE NORM IS COMPUTED.
*  II  J  INDEX OF THE FIRST ELEMENT FROM WHICH THE NORM IS COMPUTED.
*
      FUNCTION MXDRMN(N,M,A,I,J)
      INTEGER N,M,I,J
      REAL*8 A(M*N),MXDRMN
      REAL*8 POM,DEN
      INTEGER K,L
      REAL*8 ZERO
      PARAMETER (ZERO=0.0D 0)
      DEN=ZERO
      L=(J-1)*N
      DO 1 K=J,M
      DEN=MAX(DEN,ABS(A(L+I)))
      L=L+N
    1 CONTINUE
      POM=ZERO
      IF (DEN.GT.ZERO) THEN
      L=(J-1)*N
      DO 2 K=J,M
      POM=POM+(A(L+I)/DEN)**2
      L=L+N
    2 CONTINUE
      ENDIF
      MXDRMN=DEN*SQRT(POM)
      RETURN
      END
* SUBROUTINE MXDRMV               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* K-TH COLUMN OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A IS COPIED
* TO THE VECTOR X.
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX A.
*  II  M  NUMBER OF ROWS OF THE MATRIX A.
*  RI  A(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
*         ONE-DIMENSIONAL ARRAY.
*  RO  X(M)  OUTPUT VECTOR SUCH THAT X(J)=A(J,K) FOR ALL J.
*  II  K  INDEX OF THE ROW BEING COPIED TO THE OUTPUT VECTOR.
*
      SUBROUTINE MXDRMV(N,M,A,X,K)
      INTEGER N,M,K
      DOUBLE PRECISION A(*),X(*)
      INTEGER I,J
      IF (K.LT.1.OR.K.GT.N) RETURN
      I=K
      DO 1 J=1,M
      X(J)=A(I)
      I=I+N
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXDRQF               ALL SYSTEMS                92/12/01
* PORTABILITY : ALL SYSTEMS
* 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* QR DECOMPOSITION OF ROWWISE STORED DENSE RECTANGULAR MATRIX Q USING
* HOUSEHOLDER TRANSFORMATIONS WITHOUT PIVOTING.
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX Q.
*  II  M  NUMBER OF ROWS OF THE MATRIX Q.
*  RU  Q(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
*         ONE-DIMENSIONAL ARRAY.
*  RO  R(N*(N+1)/2)  UPPER TRIANGULAR MATRIX STORED IN THE PACKED FORM.
*
* SUBPROGRAMS USED :
*  S   MXDRMN  EUCLIDEAN NORM OF A PART OF THE ROWWISE STORED
*         RECTANGULAR MATRIX COLUMN.
*
* METHOD :
* P.A.BUSSINGER, G.H.GOLUB : LINEAR LEAST SQUARES SOLUTION BY
* HOUSEHOLDER TRANSFORMATION. NUMER. MATH. 7 (1965) 269-276.
*
      SUBROUTINE MXDRQF(N,M,Q,R)
      INTEGER N,M
      REAL*8 Q(M*N),R(N*(N+1)/2)
      REAL*8 ALF,POM,MXDRMN
      INTEGER I,J,K,L,JP,KP,NM
      REAL*8 ZERO,ONE
      PARAMETER (ZERO=0.0D 0,ONE=1.0D 0)
      NM=MIN(N,M)
C
C     QR DECOMPOSITION
C
      L=0
      KP=0
      DO 6 K=1,NM
      POM=MXDRMN(N,M,Q,K,K)
      IF (Q(KP+K).NE.ZERO) POM=SIGN(POM,Q(KP+K))
      R(L+K)=-POM
      JP=0
      DO 1 J=1,K-1
      R(L+J)=Q(JP+K)
      Q(JP+K)=ZERO
      JP=JP+N
    1 CONTINUE
      IF (POM.NE.ZERO) THEN
C
C     HOUSEHOLDER TRANSFORMATION
C
      DO 2 J=K,M
      Q(JP+K)=Q(JP+K)/POM
      JP=JP+N
    2 CONTINUE
      Q(KP+K)=Q(KP+K)+ONE
      DO 5 I=K+1,N
      ALF=ZERO
      JP=KP
      DO 3 J=K,M
      ALF=ALF+Q(JP+K)*Q(JP+I)
      JP=JP+N
    3 CONTINUE
      ALF=ALF/Q(KP+K)
      JP=KP
      DO 4 J=K,M
      Q(JP+I)=Q(JP+I)-ALF*Q(JP+K)
      JP=JP+N
    4 CONTINUE
    5 CONTINUE
      ENDIF
      L=L+K
      KP=KP+N
    6 CONTINUE
C
C     EXPLICIT FORMULATION OF THE ORTHOGONAL MATRIX
C
      KP=N*N
      DO 11 K=N,1,-1
      KP=KP-N
      IF (Q(KP+K).NE.ZERO) THEN
      DO 9 I=K+1,N
      ALF=ZERO
      JP=KP
      DO 7 J=K,M
      ALF=ALF+Q(JP+K)*Q(JP+I)
      JP=JP+N
    7 CONTINUE
      ALF=ALF/Q(KP+K)
      JP=KP
      DO 8 J=K,M
      Q(JP+I)=Q(JP+I)-ALF*Q(JP+K)
      JP=JP+N
    8 CONTINUE
    9 CONTINUE
      JP=KP
      DO 10 J=K,M
      Q(JP+K)=-Q(JP+K)
      JP=JP+N
   10 CONTINUE
      ENDIF
      Q(KP+K)=Q(KP+K)+ONE
   11 CONTINUE
      RETURN
      END
* SUBROUTINE MXDRQU               ALL SYSTEMS                92/12/01
* PORTABILITY : ALL SYSTEMS
* 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* UPDATE OF A QR DECOMPOSITION. THIS QR DECOMPOSITION IS UPDATED
* BY THE RULE Q*R:=Q*R+ALF*X*TRANS(Y).
*
* PARAMETERS :
*  II  N  NUMBER OF COLUMNS OF THE MATRIX Q.
*  II  M  NUMBER OF ROWS OF THE MATRIX Q.
*  RU  Q(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
*         ONE-DIMENSIONAL ARRAY (PART OF THE ORTHOGONAL MATRIX).
*  RU  R(N*(N+1)/2)  UPPER TRIANGULAR MATRIX STORED IN A PACKED FORM.
*  RI  ALF  SCALAR PARAMETER.
*  RI  X(M)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RA  Z(N)  AUXILIARY VECTOR.
*  IO  INF  INFORMATION. IF INF=0 THEN X LIES IN THE COLUMN SPACE OF Q.
*         IF INF=1 THEN X DOES NOT LIE IN THE COLUMN SPACE OF Q.
*
* SUBPROGRAMS USED :
*  RF  MXVNOR  EUCLIDEAN NORM OF A VECTOR.
*  S   MXVORT  COMPUTATION OF THE PLANE ROTATION MATRIX.
*  S   MXVROT  PLANE ROTATION IS APPLIED TO TWO NUMBERS.
*
* METHOD :
* J.W.DANIEL, W.B.GRAGG, L.KAUFMAN, G.W.STEWARD : REORTHOGONALIZATION
* AND STABLE ALGORITHMS FOR UPDATING THE GRAM-SCHMIDT QR FACTORIZATION.
* MATHEMATICS OF COMPUTATION 30 (1976) 772-795.
      SUBROUTINE MXDRQU(N,M,Q,R,ALF,X,Y,Z,INF)
      INTEGER N,M,INF
      REAL*8 Q(M*N),R(N*(N+1)/2),ALF,X(M),Y(N),Z(N)
      REAL*8 CK,CL,ZK,ZL,MXVNOR
      INTEGER J,K,L,KJ,KK,IER
      REAL*8 ONE,CON
      PARAMETER (ONE=1.0D 0,CON=1.0D-10)
      INF=0
      KK=N*(N+1)/2
C
C     COMPUTATION OF THE VECTOR TRANS(Q)*X
C
      CALL MXDCMM(N,M,Q,X,Z)
      IF (M.GT.N) THEN
C
C     IF X DOES NOT LIE IN THE COLUMN SPACE OF Q WE HAVE TO USE
C     A SUBPROBLEM WHOSE DIMENSION IS BY ONE GREATER (INF=1).
C
      ZK=MXVNOR(M,X)
      CALL MXDRMD(N,M,Q,Z,-ONE,X,X)
      ZL=MXVNOR(M,X)
      IF (ZL.GT.CON*ZK) THEN
      INF=1
      CALL MXVSCL(M,-ONE/ZL,X,X)
      CALL MXVORT(Z(N),ZL,CK,CL,IER)
      IF (IER.EQ.0.OR.IER.EQ.1) THEN
      CALL MXVROT(R(KK),ZL,CK,CL,IER)
      KJ=N
      DO 1 J=1,M
      CALL MXVROT(Q(KJ),X(J),CK,CL,IER)
      KJ=KJ+N
    1 CONTINUE
      ENDIF
      ENDIF
      ENDIF
C
C     APPLICATION OF PLANE ROTATIONS TO THE VECTOR Z SO THAT
C     TRANS(Q1)*Z=E1 WHERE Q1 IS AN ORTHOGONAL MATRIX (ACCUMULATION OF
C     THE PLANE ROTATIONS) AND E1 IS THE FIRST COLUMN OF THE UNIT
C     MATRIX. AT THE SAME TIME BOTH THE UPPER HESSENBERG MATRIX
C     TRANS(Q1)*R AND THE ORTHOGONAL MATRIX Q*Q1 ARE CONSTRUCTED SO THAT
C     Q*Q1*R1=Q*Q1*(TRANS(Q1)*R+ALF*E1*TRANS(Y)) WHERE R1 IS AN UPPER
C     HESSENBERG MATRIX.
C
      DO 4 L=N,2,-1
      K=L-1
      KK=KK-L
      CALL MXVORT(Z(K),Z(L),CK,CL,IER)
      IF (IER.EQ.0.OR.IER.EQ.1) THEN
      CALL MXVROT(R(KK),Z(L),CK,CL,IER)
      KJ=KK
      DO 2 J=L,N
      KJ=KJ+J-1
      CALL MXVROT(R(KJ),R(KJ+1),CK,CL,IER)
    2 CONTINUE
      KJ=K
      DO 3 J=1,M
      CALL MXVROT(Q(KJ),Q(KJ+1),CK,CL,IER)
      KJ=KJ+N
    3 CONTINUE
      ENDIF
    4 CONTINUE
      Z(1)=ALF*Z(1)
      KJ=1
      DO 5 J=1,N
      R(KJ)=R(KJ)+Z(1)*Y(J)
      KJ=KJ+J
    5 CONTINUE
C
C     APPLICATION OF PLANE ROTATIONS TO THE UPPER HESSENBERG MATRIX R1
C     GIVEN ABOVE SO THAT R2=TRANS(Q2)*R1 WHERE Q2 IS AN ORTHOGONAL
C     MATRIX (ACCUMULATION OF THE PLANE ROTATIONS) AND R2 IS AN UPPER
C     TRIANGULAR MATRIX. WE OBTAIN THE NEW QR DECOMPOSITION Q*Q1*Q2*R2.
C
      KK=1
      DO 8 L=2,N
      K=L-1
      CALL MXVORT(R(KK),Z(L),CK,CL,IER)
      IF (IER.EQ.0.OR.IER.EQ.1) THEN
      KJ=KK
      DO 6 J=L,N
      KJ=KJ+J-1
      CALL MXVROT(R(KJ),R(KJ+1),CK,CL,IER)
    6 CONTINUE
      KJ=K
      DO 7 J=1,M
      CALL MXVROT(Q(KJ),Q(KJ+1),CK,CL,IER)
      KJ=KJ+N
    7 CONTINUE
      ENDIF
      KK=KK+L
    8 CONTINUE
C
C     BACK TRANSFORMATION OF THE GREATER SUBPROBLEM IF INF=1.
C
      IF (INF.EQ.1) THEN
      CALL MXVORT(R(KK),ZL,CK,CL,IER)
      IF (IER.EQ.0.OR.IER.EQ.1) THEN
      KJ=N
      DO 9 J=1,M
      CALL MXVROT(Q(KJ),X(J),CK,CL,IER)
      KJ=KJ+N
    9 CONTINUE
      ENDIF
      ENDIF
      RETURN
      END
* SUBROUTINE MXDSMI                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DENSE SYMMETRIC MATRIX A IS SET TO THE UNIT MATRIX WITH THE SAME
* ORDER.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RO  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM
*         WHICH IS SET TO THE UNIT MATRIX (I.E. A:=I).
*
      SUBROUTINE MXDSMI(N,A)
*     .. Scalar Arguments ..
      INTEGER          N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(*)
*     ..
*     .. Local Scalars ..
      INTEGER          I,M
*     ..
      M = N* (N+1)/2
      DO 10 I = 1,M
          A(I) = 0.0D0
   10 CONTINUE
      M = 0
      DO 20 I = 1,N
          M = M + I
          A(M) = 1.0D0
   20 CONTINUE
      RETURN

      END
* SUBROUTINE MXDSMM                ALL SYSTEMS                89/12/01
* PORTABILITY : ALL SYSTEMS
* 89/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MULTIPLICATION OF A DENSE SYMMETRIC MATRIX A BY A VECTOR X.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  RI  X(N)  INPUT VECTOR.
*  RO  Y(N)  OUTPUT VECTOR EQUAL TO  A*X.
*
      SUBROUTINE MXDSMM(N,A,X,Y)
*     .. Scalar Arguments ..
      INTEGER          N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(*),X(*),Y(*)
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION TEMP
      INTEGER          I,J,K,L
*     ..
      K = 0
      DO 30 I = 1,N
          TEMP = 0.0D0
          L = K
          DO 10 J = 1,I
              L = L + 1
              TEMP = TEMP + A(L)*X(J)
   10     CONTINUE
          DO 20 J = I + 1,N
              L = L + J - 1
              TEMP = TEMP + A(L)*X(J)
   20     CONTINUE
          Y(I) = TEMP
          K = K + I
   30 CONTINUE
      RETURN

      END
* SUBROUTINE MXDSMR               ALL SYSTEMS                92/12/01
* PORTABILITY : ALL SYSTEMS
* 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* PLANE ROTATION IS APPLIED TO A DENSE SYMMETRIC MATRIX A. THE CASE
* K=L+1 IS REQUIRED.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2) DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  II  K  FIRST INDEX OF PLANE ROTATION.
*  II  L  SECOND INDEX OF PLANE ROTATION.
*  RO  CK  DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  RO  CL  OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  IO  IER  INFORMATION ON THE TRANSFORMATION. IER<0-K OR L OUT OF
*         RANGE. IER=0-PLANE ROTATION. IER=1-PERMUTATION.
*         IER=2-TRANSFORMATION SUPPRESSED.
*
* SUBPROGRAMS USED :
*  S   MXVROT  PLANE ROTATION IS APPLIED TO TWO NUMBERS.
*
      SUBROUTINE MXDSMR(N,A,K,L,CK,CL,IER)
      INTEGER N,K,L,IER
      DOUBLE PRECISION A(*),CK,CL
      DOUBLE PRECISION AKK,AKL,ALL,CKK,CKL,CLL
      INTEGER J,KJ,LJ,KK,KL,LL
      IF (IER.NE.0.AND.IER.NE.1) RETURN
      IF (K.NE.L+1) THEN
      IER=-1
      RETURN
      ENDIF
      LJ=L*(L-1)/2
      DO 1 J=1,N
      IF (J.LE.L) THEN
      LJ=LJ+1
      KJ=LJ+L
      ELSE
      LJ=LJ+J-1
      KJ=LJ+1
      ENDIF
      IF (J.NE.K.AND.J.NE.L) THEN
      CALL MXVROT(A(KJ),A(LJ),CK,CL,IER)
      ENDIF
    1 CONTINUE
      IF (IER.EQ.0) THEN
      CKK=CK**2
      CKL=CK*CL
      CLL=CL**2
      LL=L*(L+1)/2
      KL=LL+L
      KK=LL+K
      AKL=(CKL+CKL)*A(KL)
      AKK=CKK*A(KK)+CLL*A(LL)+AKL
      ALL=CLL*A(KK)+CKK*A(LL)-AKL
      A(KL)=(CLL-CKK)*A(KL)+CKL*(A(KK)-A(LL))
      A(KK)=AKK
      A(LL)=ALL
      ELSE
      LL=L*(L+1)/2
      KK=LL+K
      AKK=A(KK)
      A(KK)=A(LL)
      A(LL)=AKK
      ENDIF
      RETURN
      END
* SUBROUTINE MXDSMS                ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SCALING OF A DENSE SYMMETRIC MATRIX.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM
*         WHICH IS SCALED BY THE VALUE ALF (I.E. A:=ALF*A).
*  RI  ALF  SCALING FACTOR.
*
      SUBROUTINE MXDSMS( N, A, ALF)
      INTEGER N
      REAL*8  A(N*(N+1)/2), ALF
      INTEGER I,M
      M = N * (N+1) / 2
      DO 1  I = 1, M
      A(I) = A(I) * ALF
    1 CONTINUE
      RETURN
      END
* SUBROUTINE MXDSMU                ALL SYSTEMS                89/12/01
C PORTABILITY : ALL SYSTEMS
C 89/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* UPDATE OF A DENSE SYMMETRIC MATRIX A. THIS UPDATE IS DEFINED AS
* A:=A+ALF*X*TRANS(X) WHERE ALF IS A GIVEN SCALING FACTOR AND X IS
* A GIVEN VECTOR.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RU  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  RI  ALF  SCALING FACTOR IN THE CORRECTION TERM.
*  RI  X(N)  VECTOR IN THE CORRECTION TERM.
*
      SUBROUTINE MXDSMU( N, A, ALF, X)
      INTEGER N
      REAL*8  A(N*(N+1)/2), X(N), ALF
      REAL*8 TEMP
      INTEGER  I, J, K
      K = 0
      DO 2  I = 1, N
      TEMP = ALF * X(I)
      DO 1  J = 1, I
      K = K + 1
      A(K) = A(K) + TEMP * X(J)
    1 CONTINUE
    2 CONTINUE
      RETURN
      END
* SUBROUTINE MXDSMV                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* K-TH ROW OF A DENSE SYMMETRIC MATRIX A IS COPIED TO THE VECTOR X.
*
* PARAMETERS :
*  II  N  ORDER OF THE MATRIX A.
*  RI  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM.
*  RO  X(N)  OUTPUT VECTOR.
*  II  K  INDEX OF COPIED ROW.
*
      SUBROUTINE MXDSMV(N,A,X,K)
*     .. Scalar Arguments ..
      INTEGER          K,N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(*),X(*)
*     ..
*     .. Local Scalars ..
      INTEGER          I,L
*     ..
      L = K* (K-1)/2
      DO 10 I = 1,N
          IF (I.LE.K) THEN
              L = L + 1

          ELSE
              L = L + I - 1
          END IF

          X(I) = A(L)
   10 CONTINUE
      RETURN

      END
* SUBROUTINE MXVCOP                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COPYING OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RO  Y(N)  OUTPUT VECTOR WHERE Y:= X.
*
      SUBROUTINE MXVCOP(N,X,Y)
C     .. Scalar Arguments ..
      INTEGER N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION X(*),Y(*)
C     ..
C     .. Local Scalars ..
      INTEGER I
C     ..
      DO 10 I = 1,N
          Y(I) = X(I)
   10 CONTINUE
      RETURN

      END
* SUBROUTINE MXVDIF                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VECTOR DIFFERENCE.
*
* PARAMETERS :
*  RI  X(N)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RO  Z(N)  OUTPUT VECTOR WHERE Z:= X - Y.
*
      SUBROUTINE MXVDIF(N,X,Y,Z)
C     .. Scalar Arguments ..
      INTEGER N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION X(*),Y(*),Z(*)
C     ..
C     .. Local Scalars ..
      INTEGER I
C     ..
      DO 10 I = 1,N
          Z(I) = X(I) - Y(I)
   10 CONTINUE
      RETURN

      END
* SUBROUTINE MXVDIR                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VECTOR AUGMENTED BY THE SCALED VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  A  SCALING FACTOR.
*  RI  X(N)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RO  Z(N)  OUTPUT VECTOR WHERE Z:= Y + A*X.
*
      SUBROUTINE MXVDIR(N,A,X,Y,Z)
C     .. Scalar Arguments ..
      DOUBLE PRECISION A
      INTEGER N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION X(*),Y(*),Z(*)
C     ..
C     .. Local Scalars ..
      INTEGER I
C     ..
      DO 10 I = 1,N
          Z(I) = Y(I) + A*X(I)
   10 CONTINUE
      RETURN

      END
* FUNCTION MXVDOT                  ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DOT PRODUCT OF TWO VECTORS.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RR  MXVDOT  VALUE OF DOT PRODUCT MXVDOT=TRANS(X)*Y.
*
      DOUBLE PRECISION
     +  FUNCTION MXVDOT(N,X,Y)
*     .. Scalar Arguments ..
      INTEGER          N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION X(*),Y(*)
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION TEMP
      INTEGER          I
*     ..
      TEMP = 0.0D0
      DO 10 I = 1,N
          TEMP = TEMP + X(I)*Y(I)
   10 CONTINUE
      MXVDOT = TEMP
      RETURN

      END
* SUBROUTINE MXVINA             ALL SYSTEMS                   90/12/01
* PORTABILITY : ALL SYSTEMS
* 90/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* ELEMENTS OF THE INTEGER VECTOR ARE REPLACED BY THEIR ABSOLUTE VALUES.
*
* PARAMETERS :
*  II  N DIMENSION OF THE INTEGER VECTOR.
*  IU  IX(N)  INTEGER VECTOR WHICH IS UPDATED SO THAT IX(I):=ABS(IX(I))
*         FOR ALL I.
*
      SUBROUTINE MXVINA(N,IX)
C     .. Scalar Arguments ..
      INTEGER N
C     ..
C     .. Array Arguments ..
      INTEGER IX(*)
C     ..
C     .. Local Scalars ..
      INTEGER I
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS
C     ..
      DO 10 I = 1,N
          IX(I) = ABS(IX(I))
          IF (IX(I).GT.10) IX(I) = IX(I) - 10
   10 CONTINUE
      RETURN

      END
* SUBROUTINE MXVIND               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* CHANGE OF THE INTEGER VECTOR ELEMENT FOR THE CONSTRAINT ADDITION.
*
* PARAMETERS :
*  IU  IX(N)  INTEGER VECTOR.
*  II  I  INDEX OF THE CHANGED ELEMENT.
*  II JOB  CHANGE SPECIFICATION. IS JOB.EQ.0 THEN IX(I)=10-IX(I).
*
      SUBROUTINE MXVIND(IX,I,JOB)
      INTEGER IX(*),I,JOB
      IF (JOB.EQ.0) IX(I)=10-IX(I)
      RETURN
      END
* SUBROUTINE MXVINV               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* CHANGE OF THE INTEGER VECTOR ELEMENT FOR THE CONSTRAINT ADDITION.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  IU  IX(N)  INTEGER VECTOR.
*  II  I  INDEX OF THE CHANGED ELEMENT.
*  II  JOB  CHANGE SPECIFICATION
*
      SUBROUTINE MXVINV(IX,I,JOB)
*     .. Scalar Arguments ..
      INTEGER          I,JOB
*     ..
*     .. Array Arguments ..
      INTEGER          IX(*)
*     ..
      IF ((IX(I).EQ.3.OR.IX(I).EQ.5) .AND. JOB.LT.0) IX(I) = IX(I) + 1
      IF ((IX(I).EQ.4.OR.IX(I).EQ.6) .AND. JOB.GT.0) IX(I) = IX(I) - 1
      IX(I) = -IX(I)
      RETURN

      END
* FUNCTION MXVMAX               ALL SYSTEMS                   91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* L-INFINITY NORM OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RR  MXVMAX  L-INFINITY NORM OF THE VECTOR X.
*
      DOUBLE PRECISION FUNCTION MXVMAX(N,X)
C     .. Scalar Arguments ..
      INTEGER N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION X(*)
C     ..
C     .. Local Scalars ..
      INTEGER I
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,MAX
C     ..
      MXVMAX = 0.0D0
      DO 10 I = 1,N
          MXVMAX = MAX(MXVMAX,ABS(X(I)))
   10 CONTINUE
      RETURN

      END
* SUBROUTINE MXVNEG                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* CHANGE THE SIGNS OF VECTOR ELEMENTS.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RO  Y(N)  OUTPUT VECTOR WHERE Y:= - X.
*
      SUBROUTINE MXVNEG(N,X,Y)
C     .. Scalar Arguments ..
      INTEGER N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION X(*),Y(*)
C     ..
C     .. Local Scalars ..
      INTEGER I
C     ..
      DO 10 I = 1,N
          Y(I) = -X(I)
   10 CONTINUE
      RETURN

      END
* FUNCTION  MXVNOR               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* EUCLIDEAN NORM OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RR  MXVNOR  EUCLIDEAN NORM OF X.
*
      FUNCTION MXVNOR(N,X)
      INTEGER N
      REAL*8 X(N),MXVNOR
      REAL*8 POM,DEN
      INTEGER I
      REAL*8 ZERO
      PARAMETER (ZERO=0.0D 0)
      DEN=ZERO
      DO 1 I=1,N
      DEN=MAX(DEN,ABS(X(I)))
    1 CONTINUE
      POM=ZERO
      IF (DEN.GT.ZERO) THEN
      DO 2 I=1,N
      POM=POM+(X(I)/DEN)**2
    2 CONTINUE
      ENDIF
      MXVNOR=DEN*SQRT(POM)
      RETURN
      END
* SUBROUTINE MXVORT               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF AN ELEMENTARY ORTHOGONAL MATRIX FOR PLANE ROTATION.
*
* PARAMETERS :
*  RU  XK  FIRST VALUE FOR PLANE ROTATION (XK IS TRANSFORMED TO
*         SQRT(XK**2+XL**2))
*  RU  XL  SECOND VALUE FOR PLANE ROTATION (XL IS TRANSFORMED TO
*         ZERO)
*  RO  CK  DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  RO  CL  OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  IO  IER  INFORMATION ON THE TRANSFORMATION. IER=0-GENERAL PLANE
*         ROTATION. IER=1-PERMUTATION. IER=2-TRANSFORMATION SUPPRESSED.
*
      SUBROUTINE MXVORT(XK,XL,CK,CL,IER)
*     .. Scalar Arguments ..
      DOUBLE PRECISION CK,CL,XK,XL
      INTEGER          IER
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION DEN,POM
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC        ABS,SQRT
*     ..
      IF (XL.EQ.0.0D0) THEN
          IER = 2

      ELSE IF (XK.EQ.0.0D0) THEN
          XK = XL
          XL = 0.0D0
          IER = 1

      ELSE
          IF (ABS(XK).GE.ABS(XL)) THEN
              POM = XL/XK
              DEN = SQRT(1.0D0+POM*POM)
              CK = 1.0D0/DEN
              CL = POM/DEN
              XK = XK*DEN

          ELSE
              POM = XK/XL
              DEN = SQRT(1.0D0+POM*POM)
              CL = 1.0D0/DEN
              CK = POM/DEN
              XK = XL*DEN
          END IF

          XL = 0.0D0
          IER = 0
      END IF

      RETURN

      END
* SUBROUTINE MXVROT               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* PLANE ROTATION IS APPLIED TO TWO VALUES.
*
* PARAMETERS :
*  RU  XK  FIRST VALUE FOR PLANE ROTATION.
*  RU  XL  SECOND VALUE FOR PLANE ROTATION.
*  RI  CK  DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  RI  CL  OFF-DIAGONAL ELEMENT OF THE ELEMENTARY ORTHOGONAL MATRIX.
*  II  IER  INFORMATION ON THE TRANSFORMATION. IER=0-GENERAL PLANE
*         ROTATION. IER=1-PERMUTATION. IER=2-TRANSFORMATION SUPPRESSED.
*
      SUBROUTINE MXVROT(XK,XL,CK,CL,IER)
*     .. Scalar Arguments ..
      DOUBLE PRECISION CK,CL,XK,XL
      INTEGER          IER
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION YK,YL
*     ..
      IF (IER.EQ.0) THEN
          YK = XK
          YL = XL
          XK = CK*YK + CL*YL
          XL = CL*YK - CK*YL

      ELSE IF (IER.EQ.1) THEN
          YK = XK
          XK = XL
          XL = YK
      END IF

      RETURN

      END
* SUBROUTINE MXVSAV                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DIFFERENCE OF TWO VECTORS RETURNED IN THE SUBSTRACTED ONE.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RU  Y(N)  UPDATE VECTOR WHERE Y:= X - Y.
*
      SUBROUTINE MXVSAV(N,X,Y)
*     .. Scalar Arguments ..
      INTEGER          N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION X(*),Y(*)
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION TEMP
      INTEGER          I
*     ..
      DO 10 I = 1,N
          TEMP = Y(I)
          Y(I) = X(I) - Y(I)
          X(I) = TEMP
   10 CONTINUE
      RETURN

      END
* SUBROUTINE MXVSCL                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SCALING OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RI  A  SCALING FACTOR.
*  RO  Y(N)  OUTPUT VECTOR WHERE Y:= A*X.
*
      SUBROUTINE MXVSCL(N,A,X,Y)
C     .. Scalar Arguments ..
      DOUBLE PRECISION A
      INTEGER N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION X(*),Y(*)
C     ..
C     .. Local Scalars ..
      INTEGER I
C     ..
      DO 10 I = 1,N
          Y(I) = A*X(I)
   10 CONTINUE
      RETURN

      END
* SUBROUTINE MXVSET                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* A SCALAR IS SET TO ALL THE ELEMENTS OF A VECTOR.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  A  INITIAL VALUE.
*  RO  X(N)  OUTPUT VECTOR SUCH THAT X(I)=A FOR ALL I.
*
      SUBROUTINE MXVSET(N,A,X)
C     .. Scalar Arguments ..
      DOUBLE PRECISION A
      INTEGER N
C     ..
C     .. Array Arguments ..
      DOUBLE PRECISION X(*)
C     ..
C     .. Local Scalars ..
      INTEGER I
C     ..
      DO 10 I = 1,N
          X(I) = A
   10 CONTINUE
      RETURN

      END
* SUBROUTINE MXVSUM                ALL SYSTEMS                88/12/01
* PORTABILITY : ALL SYSTEMS
* 88/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SUM OF TWO VECTORS.
*
* PARAMETERS :
*  II  N  VECTOR DIMENSION.
*  RI  X(N)  INPUT VECTOR.
*  RI  Y(N)  INPUT VECTOR.
*  RO  Z(N)  OUTPUT VECTOR WHERE Z:= X + Y.
*
      SUBROUTINE MXVSUM(N,X,Y,Z)
*     .. Scalar Arguments ..
      INTEGER          N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION X(*),Y(*),Z(*)
*     ..
*     .. Local Scalars ..
      INTEGER          I
*     ..
      DO 10 I = 1,N
          Z(I) = X(I) + Y(I)
   10 CONTINUE
      RETURN

      END


C ##############################################################################
C PART III: PQSUBS


* SUBROUTINE PA0GS1             ALL SYSTEMS                 97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE:
* NUMERICAL COMPUTATION OF THE GRADIENT OF THE MODEL FUNCTION.
*
* PARAMETERS :
*  II  N  NUMBER OF VARIABLES.
*  II  KA  INDEF OF THE APPROXIMATED FUNCTION.
*  RI  X(N)  VECTOR OF VARIABLES.
*  RO  GA(N)  GRADIENT OF THE APPROXIMATED FUNCTION.
*  RI  FA  VALUE OF THE APPROXIMATED FUNCTION.
*  RI  ETA1  PRECISION OF THE COMPUTED VALUES.
*  IU  NAV  NUMBER OF APPROXIMATED FUNCTION EVALUATIONS.
*
      SUBROUTINE PA0GS1(N,KA,X,GA,FA,ETA1,NAV)
      INTEGER N,KA,NAV
      REAL*8 X(*),GA(*),FA,ETA1
      REAL*8 XSTEP,XTEMP,FTEMP,ETA
      INTEGER IVAR
      ETA=SQRT(ETA1)
      FTEMP=FA
      DO 4 IVAR=1,N
C
C     STEP SELECTION
C
      XSTEP=1.0D 0
      XSTEP=ETA*MAX(ABS(X(IVAR)),XSTEP)*SIGN(1.0D 0,X(IVAR))
      XTEMP=X(IVAR)
      X(IVAR)=X(IVAR)+XSTEP
      XSTEP=X(IVAR)-XTEMP
      NAV=NAV+1
      CALL FUN(N,KA,X,FA)
C
C     NUMERICAL DIFFERENTIATION
C
      GA(IVAR)=(FA-FTEMP)/XSTEP
      X(IVAR)=XTEMP
    4 CONTINUE
      FA=FTEMP
      RETURN
      END
* SUBROUTINE PA1SQ1             ALL SYSTEMS                 97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE FUNCTION
* WHICH IS DEFINED AS A SUM OF SQUARES.
*
* PARAMETERS:
*  II  N  NUMBER OF VARIABLES.
*  RI  X(N)  VECTOR OF VARIABLES.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  AF(N)  VALUES OF THE APPROXIMATED FUNCTIONS.
*  RI  GA(NF)  GRADIENT OF THE APPROXIMATED FUNCTION.
*  RI  AG(N*N)  RECTANGULAR MATRIX WHICH IS USED FOR THE DIRECTION
*         VECTOR DETERMINATION.
*  RO  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RI  ETA1  PRECISION OF THE COMPUTES FUNCTION VALUES.
*  II  KDA  DEGREE OF COMPUTED DERIVATIVES.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IO  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  IU  NFV  NUMBER OF OBJECTIVE FUNCTION VALUES COMPUTED.
*  IU  NFG  NUMBER OF OBJECTIVE FUNCTION GRADIENTS COMPUTED.
*
* SUBPROGRAMS USED :
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PA1SQ1(N,X,F,AF,GA,AG,G,ETA1,KDA,KD,LD,NFV,NFG)
      INTEGER N,KDA,KD,LD,NFV,NFG
      REAL*8 X(*),F,AF(*),GA(*),AG(*),G(*),ETA1
      INTEGER KA,NAV
      REAL*8 FA
      IF (KD.LE.LD) RETURN
      IF (KD.GE.0.AND.LD.LT.0) THEN
      F=0.0D0
      NFV=NFV+1
      ENDIF
      IF (KD.GE.1.AND.LD.LT.1) THEN
      CALL MXVSET(N,0.0D0,G)
      IF (KDA.GT.0) NFG=NFG+1
      ENDIF
      NAV=0
      DO 30 KA=1,N
      IF (KD.LT.0) GO TO 30
      IF (LD.GE.0) THEN
      FA = AF(KA)
      GO TO 10
      ELSE
      CALL FUN(N,KA,X,FA)
      AF(KA) = FA
      ENDIF
      IF (LD.GE.0) GO TO 10
      F=F+FA*FA
   10 IF (KD.LT.1) GO TO 30
      IF (KDA.GT.0) THEN
      CALL DFUN(N,KA,X,GA)
      ELSE
      CALL PA0GS1(N,KA,X,GA,FA,ETA1,NAV)
      ENDIF
      CALL MXVDIR(N,FA,GA,G,G)
      CALL MXVCOP(N,GA,AG((KA-1)*N+1))
   30 CONTINUE
      NFV=NFV+NAV/N
      IF (KD.GE.0.AND.LD.LT.0) F=0.5D0*F
      LD=KD
      RETURN
      END
* SUBROUTINE PC1F01             ALL SYSTEMS                 97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE VALUE AND THE GRADIENT OF THE CONSTRAINT FUNCTION.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RI  FC  VALUE OF THE SELECTED CONSTRAINT FUNCTION.
*  RI  CF(NC)  VECTOR CONTAINING VALUES OF CONSTRAINT FUNCTIONS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  GC(NF)  GRADIENT OF THE SELECTED CONSTRAINT FUNCTION.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE GRADIENTS OF CONSTRAINT
*         FUNCTIONS.
*  RO  CMAX  MAXIMUM CONSTRAINT VIOLATION.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  II  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*
* SUBPROGRAMS USED :
*  S   MXVCOP  COPYING OF A VECTOR.
*
      SUBROUTINE PC1F01(NF,NC,X,FC,CF,CL,CU,IC,GC,CG,CMAX,KD,LD)
*     .. Scalar Arguments ..
      DOUBLE PRECISION FC,CMAX
      INTEGER          KD,LD,NC,NF
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION CF(*),CG(*),CL(*),CU(*),GC(*),X(*)
      INTEGER          IC(*)
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION POM,TEMP
      INTEGER          KC
*     ..
*     .. External Subroutines ..
      EXTERNAL         CON,DCON,MXVCOP
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC        MAX,MIN
*     ..

      IF (KD.LE.LD) RETURN
      IF (LD.LT.0) CMAX = 0.0D0
      DO 20 KC = 1,NC
          IF (KD.LT.0) GO TO 20
          IF (LD.GE.0) THEN
              FC = CF(KC)
              GO TO 10

          ELSE
              CALL CON(NF,KC,X,FC)
              CF(KC) = FC
          END IF

          IF (IC(KC).GT.0) THEN
              POM = 0.0D0
              TEMP = CF(KC)
              IF (IC(KC).EQ.1.OR.IC(KC).GE.3) POM = MIN(POM,
     +                                              TEMP - CL(KC))
              IF (IC(KC).EQ.2.OR.IC(KC).GE.3) POM = MIN(POM,
     +                                              CU(KC) - TEMP)
              IF (POM.LT.0.0D0) THEN
                  CMAX = MAX(CMAX,-POM)
              END IF

          END IF

   10     IF (KD.LT.1) GO TO 20
          IF (LD.GE.1) THEN
              CALL MXVCOP(NF,CG((KC - 1) * NF + 1),GC)
              GO TO 20

          ELSE
              CALL DCON(NF,KC,X,GC)
              CALL MXVCOP(NF,GC,CG((KC - 1) * NF + 1))
          END IF

   20 CONTINUE
      LD = KD
      RETURN

      END
* SUBROUTINE PF1F01                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE FUNCTION.
*
* PARAMETERS:
*  II  NF  NUMBER OF VARIABLES.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RI  GF(NF)  GRADIENT OF THE MODEL FUNCTION.
*  RO  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RI  FF  VALUE OF THE MODEL FUNCTION.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  II  KD  DEGREE OF REQUIRED DERIVATIVES.
*  II  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  II  IEXT  TYPE OF EXTREMUM. IEXT=0-MINIMUM. IEXT=1-MAXIMUM.
*
* SUBPROGRAMS USED :
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN.
*
      SUBROUTINE PF1F01(NF,X,GF,G,FF,F,KD,LD,IEXT)
*     .. Scalar Arguments ..
      DOUBLE PRECISION F,FF
      INTEGER          IEXT,KD,LD,NF
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION GF(*),G(*),X(*)
*     ..
*     .. Scalars in Common ..
      INTEGER          NADD,NDEC,NFG,NFH,NFV,NIT,NREM,NRES
*     ..
*     .. External Subroutines ..
      EXTERNAL         FUN,DFUN
*     ..
*     .. Common blocks ..
      COMMON           /STATSQP/NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
*     ..

      IF (KD.LE.LD) RETURN
      IF (LD.GE.0) GO TO 10
      NFV = NFV + 1
      CALL OBJ(NF,X,FF)
      IF (IEXT.LE.0) THEN
          F =  FF

      ELSE
          F = -FF

      END IF
   10 IF (KD.LT.1) GO TO 20
      IF (LD.GE.1) GO TO 20
      NFG = NFG + 1
      CALL DOBJ(NF,X,GF)
      IF (IEXT.GT.0) THEN
          CALL MXVNEG(NF,GF,G)
      END IF

   20 LD = KD
      RETURN

      END
* SUBROUTINE PLADB0               ALL SYSTEMS                97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* NEW LINEAR CONSTRAINT OR A NEW SIMPLE BOUND IS ADDED TO THE
* ACTIVE SET.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  IU  N  ACTUAL NUMBER OF VARIABLES.
*  IU  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RU  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RO  GMAX  MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER.
*  II  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IU  NADD  NUMBER OF CONSTRAINT ADDITIONS.
*  IO  IER  ERROR INDICATOR.
*
* SUBPROGRAMS USED :
*  S   PLADR0  CORRECTION OF KERNEL OF THE ORTHOGONAL PROJECTION
*         AFTER CONSTRAINT ADDITION.
*  S   MXDRMM  PREMULTIPLICATION OF A VECTOR BY A ROWWISE STORED DENSE
*         RECTANGULAR MATRIX.
*  S   MXDRMV  COPY OF THE SELECTED COLUMN OF A ROWWISE STORED DENSE
*         RECTANGULAR MATRIX.
*  S   MXDRGR  PLANE ROTATION OF A TRANSPOSED DENSE RECTANGULAR MATRIX.
*  S   MXVORT  DETERMINATION OF AN ELEMENTARY ORTHOGONAL MATRIX FOR
*         PLANE ROTATION.
*
      SUBROUTINE PLADB0(NF,N,ICA,CG,CR,CZ,S,EPS7,GMAX,UMAX,INEW,NADD,
     & IER)
      INTEGER NF,N,ICA(*),INEW,NADD,IER
      DOUBLE PRECISION CG(*),CR(*),CZ(*),S(*),EPS7,GMAX,UMAX
      DOUBLE PRECISION CK,CL
      INTEGER K,L,N1
      CALL PLADR0(NF,N,ICA,CG,CR,S,EPS7,GMAX,UMAX,INEW,NADD,IER)
      IF (IER.NE.0) RETURN
      IF (N.GT.0) THEN
      N1=N+1
      IF (INEW.GT.0) THEN
      CALL MXDRMM(NF,N1,CZ,CG((INEW-1)*NF+1),S)
      ELSE
      CALL MXDRMV(NF,N1,CZ,S,-INEW)
      ENDIF
      DO 1 L=1,N
      K=L+1
      CALL MXVORT(S(K),S(L),CK,CL,IER)
      CALL MXDRGR(NF,CZ,K,L,CK,CL,IER)
      IF (IER.LT.0) RETURN
    1 CONTINUE
      ENDIF
      IER=0
      RETURN
      END
* SUBROUTINE PLADB4               ALL SYSTEMS                97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* NEW LINEAR CONSTRAINT OR A NEW SIMPLE BOUND IS ADDED TO THE ACTIVE
* SET. TRANSFORMED HESSIAN MATRIX APPROXIMATION OR ITS INVERSION
* IS UPDATED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  IU  N  ACTUAL NUMBER OF VARIABLES.
*  IU  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RU  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RU  H(NF*(NF+1)/2)  TRANSFORMED HESSIAN MATRIX APPROXIMATION OR
*         ITS INVERSION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RO  GMAX  MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER.
*  IU  IDECF  DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION.
*         IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION.
*         IDECF=10-DIAGONAL MATRIX.
*  II  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IU  NADD  NUMBER OF CONSTRAINT ADDITIONS.
*  IO  IER  ERROR INDICATOR.
*
* SUBPROGRAMS USED :
*  S   PLADR0  CORRECTION OF KERNEL OF THE ORTHOGONAL PROJECTION
*         AFTER CONSTRAINT ADDITION.
*  S   MXDRMM  PREMULTIPLICATION OF A VECTOR BY A ROWWISE STORED DENSE
*         RECTANGULAR MATRIX.
*  S   MXDRMV  COPY OF THE SELECTED COLUMN OF A ROWWISE STORED DENSE
*         RECTANGULAR MATRIX.
*  S   MXDRGR  PLANE ROTATION OF A TRANSPOSED DENSE RECTANGULAR MATRIX.
*         RECTANGULAR MATRIX.
*  S   MXDSMR  PLANE ROTATION OF A DENSE SYMMETRIC MATRIX.
*  S   MXVORT  DETERMINATION OF AN ELEMENTARY ORTHOGONAL MATRIX FOR
*         PLANE ROTATION.
*
      SUBROUTINE PLADB4(NF,N,ICA,CG,CR,CZ,H,S,EPS7,GMAX,UMAX,IDECF,
     & INEW,NADD,IER)
      INTEGER NF,N,ICA(*),IDECF,INEW,NADD,IER
      DOUBLE PRECISION CG(*),CR(*),CZ(*),H(*),S(*),EPS7,GMAX,UMAX
      DOUBLE PRECISION CK,CL
      INTEGER I,J,K,L,N1
      IF (IDECF.NE.0.AND.IDECF.NE.9) THEN
      IER=-2
      RETURN
      ENDIF
      CALL PLADR0(NF,N,ICA,CG,CR,S,EPS7,GMAX,UMAX,INEW,NADD,IER)
      IF (IER.NE.0) RETURN
      IF (N.GT.0) THEN
      N1=N+1
      IF (INEW.GT.0) THEN
      CALL MXDRMM(NF,N1,CZ,CG((INEW-1)*NF+1),S)
      ELSE
      CALL MXDRMV(NF,N1,CZ,S,-INEW)
      ENDIF
      DO 1 L=1,N
      K=L+1
      CALL MXVORT(S(K),S(L),CK,CL,IER)
      CALL MXDRGR(NF,CZ,K,L,CK,CL,IER)
      CALL MXDSMR(N1,H,K,L,CK,CL,IER)
      IF (IER.LT.0) RETURN
    1 CONTINUE
      IF (IDECF.EQ.9) THEN
      L=N*(N+1)/2
      IF (H(L+N1).NE.0.0D 0) THEN
      CL=1.0D 0/H(L+N1)
      K=0
      DO 3 I=1,N
      CK=CL*H(L+I)
      DO 2 J=1,I
      K=K+1
      H(K)=H(K)-CK*H(L+J)
    2 CONTINUE
    3 CONTINUE
      ENDIF
      ENDIF
      ENDIF
      IER=0
      RETURN
      END
* SUBROUTINE PLADR0               ALL SYSTEMS                97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* TRIANGULAR DECOMPOSITION OF KERNEL OF THE ORTHOGONAL PROJECTION
* IS UPDATED AFTER CONSTRAINT ADDITION.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  IU  N  ACTUAL NUMBER OF VARIABLES.
*  IU  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RU  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RO  GMAX  MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER.
*  II  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IU  NADD  NUMBER OF CONSTRAINT ADDITIONS.
*  IO  IER  ERROR INDICATOR.
*
* SUBPROGRAMS USED :
*  S   MXSPRB  SPARSE BACK SUBSTITUTION.
*  S   MXVCOP  COPYING OF A VECTOR.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*
      SUBROUTINE PLADR0(NF,N,ICA,CG,CR,S,EPS7,GMAX,UMAX,INEW,NADD,IER)
      INTEGER NF,N,ICA(*),INEW,NADD,IER
      DOUBLE PRECISION CG(*),CR(*),S(*),EPS7,GMAX,UMAX
      DOUBLE PRECISION MXVDOT
      INTEGER NCA,NCR,I,J,K,L
      IER=0
      IF (N.LE.0) IER=2
      IF (INEW.EQ.0) IER=3
      IF (IER.NE.0) RETURN
      NCA=NF-N
      NCR=NCA*(NCA+1)/2
      IF (INEW.GT.0) THEN
      CALL MXVCOP(NF,CG((INEW-1)*NF+1),S)
      GMAX=MXVDOT(NF,CG((INEW-1)*NF+1),S)
      DO 1 J=1,NCA
      L=ICA(J)
      IF (L.GT.0) THEN
      CR(NCR+J)=MXVDOT(NF,CG((L-1)*NF+1),S)
      ELSE
      I=-L
      CR(NCR+J)=S(I)
      ENDIF
    1 CONTINUE
      ELSE
      K=-INEW
      GMAX=1.0D 0
      DO 2 J=1,NCA
      L=ICA(J)
      IF (L.GT.0) THEN
      CR(NCR+J)=CG((L-1)*NF+K)*GMAX
      ELSE
      CR(NCR+J)=0.0D 0
      ENDIF
    2 CONTINUE
      ENDIF
      IF (NCA.EQ.0) THEN
      UMAX=GMAX
      ELSE
      CALL MXDPRB(NCA,CR,CR(NCR+1),1)
      UMAX=GMAX-MXVDOT(NCA,CR(NCR+1),CR(NCR+1))
      ENDIF
      IF (UMAX.LE.EPS7*GMAX) THEN
      IER=1
      RETURN
      ELSE
      N=N-1
      NCA=NCA+1
      NCR=NCR+NCA
      ICA(NCA)=INEW
      CR(NCR)=SQRT(UMAX)
      NADD=NADD+1
      ENDIF
      RETURN
      END
* SUBROUTINE PLLPB1             ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE INITIAL FEASIBLE POINT AND THE LINEAR PROGRAMMING
* SUBROUTINE.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NC  NUMBER OF LINEAR CONSTRAINTS.
*  RU  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RO  XO(NF)  SAVED VECTOR OF VARIABLES.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RU  CF(NF)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  RA  CFD(NF)  VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT
*         FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  IO  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RO  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RO  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RO  GO(NF)  SAVED GRADIENT OF THE OBJECTIVE FUNCTION.
*  RA  S(NF)  DIRECTION VECTOR.
*  II  MFP  TYPE OF FEASIBLE POINT. MFP=1-ARBITRARY FEASIBLE POINT.
*         MFP=2-OPTIMUM FEASIBLE POINT. MFP=3-REPEATED SOLUTION.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  RI  ETA9  MAXIMUM FOR REAL NUMBERS.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RI  EPS9  TOLERANCE FOR ACTIVITY OF CONSTRAINTS.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER.
*  RO  GMAX  MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE.
*  IO  N  DIMENSION OF THE MANIFOLD DEFINED BY ACTIVE CONSTRAINTS.
*  IO  ITERL  TYPE OF FEASIBLE POINT. ITERL=1-ARBITRARY FEASIBLE POINT.
*         ITERL=2-OPTIMUM FEASIBLE POINT. ITERL=-1 FEASIBLE POINT DOES
*         NOT EXISTS. ITERL=-2 OPTIMUM FEASIBLE POINT DOES NOT EXISTS.
*
* SUBPROGRAMS USED :
*  S   PLINIT  DETERMINATION OF INITIAL POINT SATISFYING SIMPLE BOUNDS.
*  S   PLMAXL  MAXIMUM STEPSIZE USING LINEAR CONSTRAINTS.
*  S   PLMAXS  MAXIMUM STEPSIZE USING SIMPLE BOUNDS.
*  S   PLMAXT  MAXIMUM STEPSIZE USING TRUST REGION BOUNDS.
*  S   PLNEWL  IDENTIFICATION OF ACTIVE LINEAR CONSTRAINTS.
*  S   PLNEWS  IDENTIFICATION OF ACTIVE SIMPLE BOUNDS.
*  S   PLNEWT  IDENTIFICATION OF ACTIVE TRUST REGION BOUNDS.
*  S   PLDIRL  NEW VALUES OF CONSTRAINT FUNCTIONS.
*  S   PLDIRS  NEW VALUES OF VARIABLES.
*  S   PLSETC  INITIAL VALUES OF CONSTRAINT FUNCTIONS.
*  S   PLSETG  DETERMINATION OF THE FIRST PHASE GRADIENT VECTOR.
*  S   PLTRBG  DETERMINATION OF LAGRANGE MULTIPLIERS AND COMPUTATION
*  S   PLADB0  CONSTRAINT ADDITION.
*  S   PLRMB0  CONSTRAINT DELETION.
*  S   MXDCMM  PREMULTIPLICATION OF A VECTOR BY A DENSE RECTANGULAR
*         MATRIX STORED BY COLUMNS.
*  S   MXDRMM  PREMULTIPLICATION OF A VECTOR BY A DENSE RECTANGULAR
*         MATRIX STORED BY ROWS.
*  S   MXDSMI  DETERMINATION OF THE INITIAL UNIT DENSE SYMMETRIC
*         MATRIX.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVINA  ABSOLUTE VALUES OF ELEMENTS OF AN INTEGER VECTOR.
*  S   MXVINC  UPDATE OF AN INTEGER VECTOR.
*  S   MXVIND  CHANGE OF THE INTEGER VECTOR FOR CONSTRAINT ADDITION.
*  S   MXVINT  CHANGE OF THE INTEGER VECTOR FOR TRUST REGION BOUND
*         ADDITION.
*  S   MXVMUL  DIAGONAL PREMULTIPLICATION OF A VECTOR.
*  S   MXVNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PLLPB1(NF,NC,X,IX,XO,XL,XU,CF,CFD,IC,ICA,CL,CU,CG,
     & CR,CZ,G,GO,S,MFP,KBF,KBC,ETA9,EPS7,EPS9,UMAX,GMAX,N,ITERL)
      INTEGER NF,NC,IX(*),IC(*),ICA(*),MFP,KBF,KBC,N,ITERL
      DOUBLE PRECISION X(*),XO(*),XL(*),XU(*),CF(*),CFD(*),CL(*),
     & CU(*),CG(*),CR(*),CZ(*),G(*),GO(*),S(*),ETA9,EPS7,EPS9,
     & UMAX,GMAX
      DOUBLE PRECISION POM,CON,DMAX
      INTEGER NCA,NCR,NCZ,IPOM,I,K,IOLD,INEW,IER,KREM,KC,NRED
      INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      COMMON /STATSQP/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      CON=ETA9
*
*     INITIATION
*
      CALL MXVCOP(NF,X,XO)
      IPOM=0
      NRED=0
      KREM=0
      ITERL=1
      DMAX=0.0D 0
      IF (MFP.EQ.3) GO TO 5
      IF (KBF.GT.0) CALL MXVINA(NF,IX)
*
*     SHIFT OF VARIABLES FOR SATISFYING SIMPLE BOUNDS
*
      CALL PLINIT(NF,X,IX,XL,XU,EPS9,KBF,INEW,ITERL)
      IF (ITERL.LT.0) THEN
      GO TO 6
      ENDIF
      N=0
      NCA=0
      NCZ=0
      DO 1 I=1,NF
      IF (KBF.GT.0.AND.IX(I).LT.0) THEN
      NCA=NCA+1
      ICA(NCA)=-I
      ELSE
      N=N+1
      CALL MXVSET(NF,0.0D 0,CZ(NCZ+1))
      CZ(NCZ+I)=1.0D 0
      NCZ=NCZ+NF
      ENDIF
    1 CONTINUE
      CALL MXDSMI(NCA,CR)
      IF (NC.GT.0) THEN
      CALL MXDRMM(NF,NC,CG,X,CF)
*
*     ADDITION OF ACTIVE CONSTRAINTS AND INITIAL CHECK OF FEASIBILITY
*
      CALL MXVINA(NC,IC)
C      IF (NF.GT.N) CALL PLSETC(NF,NC,X,XO,CF,IC,CG,S)
      DO 2 KC=1,NC
      IF (IC(KC).NE.0) THEN
      INEW=0
      CALL PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW)
      CALL PLADB0(NF,N,ICA,CG,CR,CZ,S,EPS7,GMAX,UMAX,INEW,NADD,IER)
      CALL MXVIND(IC,KC,IER)
      IF (IC(KC).LT.-10) IPOM=1
      ENDIF
    2 CONTINUE
      ENDIF
    3 IF (IPOM.EQ.1) THEN
*
*     CHECK OF FEASIBILITY AND UPDATE OF THE FIRST PHASE OBJECTIVE
*     FUNCTION
*
      CALL PLSETG(NF,NC,IC,CG,GO,INEW)
      IF (INEW.EQ.0) IPOM=0
      ENDIF
      IF (IPOM.EQ.0.AND.ITERL.EQ.0) THEN
*
*     FEASIBILITY ACHIEVED
*
      ITERL=1
      CALL MXVCOP(NF,G,GO)
      IF (MFP.EQ.1) GO TO 6
      ENDIF
*
*     LAGRANGE MULTIPLIERS AND REDUCED GRADIENT DETERMINATION
*
    5 CALL PLTRBG(NF,N,NC,IX,IC,ICA,CG,CR,CZ,GO,S,EPS7,GMAX,UMAX,IOLD)
      INEW=0
      IF (GMAX.EQ.0.0D 0) THEN
*
*     OPTIMUM ON A LINEAR MANIFOLD OBTAINED
*
      IF (IOLD.EQ.0) THEN
      IF (IPOM.EQ.0) THEN
*
*     OPTIMAL SOLUTION ACHIEVED
*
      ITERL= 2
      GO TO 6
      ELSE
      IPOM=0
      DO 22 KC=1,NC
      IF (IC(KC).LT.-10) THEN
      INEW=0
      CALL PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW)
      IF (IC(KC).LT.-10) IPOM=1
      ENDIF
   22 CONTINUE
      IF (IPOM.EQ.0) THEN
*
*     OPTIMAL SOLUTION ACHIEVED
*
      CALL MXVCOP(NF,GO,G)
      ITERL= 2
      GO TO 6
      ELSE
*
*     FEASIBLE SOLUTION DOES NOT EXIST
*
      CALL MXVCOP(NF,GO,G)
      ITERL=-1
      GO TO 6
      ENDIF
      ENDIF
      ELSE
*
*     CONSTRAINT DELETION
*
      CALL PLRMB0(NF,N,ICA,CG,CR,CZ,GO,S,IOLD,KREM,NREM,IER)
      KC=ICA(NF-N+1)
      IF (KC.GT.0) THEN
      IC(KC)=-IC(KC)
      ELSE
      K=-KC
      IX(K)=-IX(K)
      ENDIF
      DMAX=0.0D 0
      GO TO 5
      ENDIF
      ELSE
*
*     DIRECTION DETERMINATION
*
      NCA=NF-N
      NCR=NCA*(NCA+1)/2
      CALL MXDCMM(NF,N,CZ,S,CR(NCR+1))
      CALL MXVNEG(NF,CR(NCR+1),S)
*
*     STEPSIZE SELECTION
*
      POM=CON
      CALL PLMAXL(NF,NC,CF,CFD,IC,CL,CU,CG,S,POM,KBC,KREM,INEW)
      CALL PLMAXS(NF,X,IX,XL,XU,S,POM,KBF,KREM,INEW)
      IF (INEW.EQ.0) THEN
      IF (IPOM.EQ.0) THEN
*
*     BOUNDED SOLUTION DOES NOT EXIST
*
      ITERL=-2
      ELSE
*
*     FEASIBLE SOLUTION DOES NOT EXIST
*
      ITERL=-3
      ENDIF
      GO TO 6
      ELSE
*
*     STEP REALIZATION
*
      CALL PLDIRS(NF,X,IX,S,POM,KBF)
      CALL PLDIRL(NC,CF,CFD,IC,POM,KBC)
*
*     CONSTRAINT ADDITION
*
      IF (INEW.GT.0) THEN
      KC=INEW
      INEW=0
      CALL PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW)
      CALL PLADB0(NF,N,ICA,CG,CR,CZ,S,EPS7,GMAX,UMAX,INEW,NADD,IER)
      CALL MXVIND(IC,KC,IER)
      ELSE IF (INEW+NF.GE.0) THEN
      I=-INEW
      INEW=0
      CALL PLNEWS(X,IX,XL,XU,EPS9,I,INEW)
      CALL PLADB0(NF,N,ICA,CG,CR,CZ,S,EPS7,GMAX,UMAX,INEW,NADD,IER)
      CALL MXVIND(IX,I,IER)
      ENDIF
      DMAX=POM
      IF (DMAX.GT.0.0D 0) NRED=NRED+1
      GO TO 3
      ENDIF
      ENDIF
    6 CONTINUE
      RETURN
      END
* SUBROUTINE PLRMB0               ALL SYSTEMS                92/12/01
* PORTABILITY : ALL SYSTEMS
* 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* OLD LINEAR CONSTRAINT OR AN OLD SIMPLE BOUND IS REMOVED FROM THE
* ACTIVE SET.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RU  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RU  GN(NF)  TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION.
*  II  IOLD  INDEX OF THE OLD ACTIVE CONSTRAINT.
*  IO  KREM  AUXILIARY VARIABLE.
*  IU  NREM NUMBER OF CONSTRAINT DELETION.
*  IO  IER  ERROR INDICATOR.
*
* SUBPROGRAMS USED :
*  S   PLRMR0  CORRECTION OF KERNEL OF THE ORTHOGONAL PROJECTION
*         AFTER CONSTRAINT DELETION.
*  S   MXDPRB  BACK SUBSTITUTION.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*  S   MXVMUL  DIAGONAL PREMULTIPLICATION OF A VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PLRMB0(NF,N,ICA,CG,CR,CZ,G,GN,IOLD,KREM,NREM,IER)
      INTEGER NF,N,ICA(*),IOLD,KREM,NREM,IER
      DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*)
      DOUBLE PRECISION MXVDOT
      INTEGER NCA,NCR,NCZ,I,J,KC
      IER=0
      IF (N.EQ.NF) IER=2
      IF (IOLD.EQ.0) IER=3
      IF (IER.NE.0) RETURN
      NCA=NF-N
      NCR=NCA*(NCA-1)/2
      NCZ=N*NF
      CALL PLRMR0(NF,ICA,CR,CZ(NCZ+1),N,IOLD,KREM,IER)
      CALL MXVSET(NCA,0.0D 0,CZ(NCZ+1))
      CZ(NCZ+NCA)=1.0D 0
      CALL MXDPRB(NCA,CR,CZ(NCZ+1),-1)
      CALL MXVCOP(NCA,CZ(NCZ+1),CR(NCR+1))
      N=N+1
      CALL MXVSET(NF,0.0D 0,CZ(NCZ+1))
      DO 4 J=1,NCA
      KC=ICA(J)
      IF (KC.GT.0) THEN
      CALL MXVDIR(NF,CR(NCR+J),CG((KC-1)*NF+1),CZ(NCZ+1),CZ(NCZ+1))
      ELSE
      I=-KC
      CZ(NCZ+I)=CZ(NCZ+I)+CR(NCR+J)
      ENDIF
    4 CONTINUE
      GN(N)=MXVDOT(NF,CZ(NCZ+1),G)
      NREM=NREM+1
      IER=0
      RETURN
      END
* SUBROUTINE PLQDB1             ALL SYSTEMS                   97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  IO  N  DIMENSION OF THE MANIFOLD DEFINED BY ACTIVE CONSTRAINTS.
*  II  NC  NUMBER OF LINEAR CONSTRAINTS.
*  RU  X(NF)   VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  II  IXA(NF)  VECTOR CONTAINING INFORMATION ON TRUST REGION ACTIVITY.
*  RI  XN(NF)  VECTOR OF SCALING FACTORS.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  CF(NF)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  RO  CFD(NC)  VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT
*            FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RO  CZ(NF)  VECTOR OF LAGRANGE MULTIPLIERS.
*  RO  G(NF)  GRADIENT OF THE LAGRANGIAN FUNCTION.
*  RO  GO(NF)  SAVED GRADIENT OF THE OBJECTIVE FUNCTION.
*  RU  H(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OR INVERSION OF THE
*         HESSIAN MATRIX APPROXIMATION.
*  RI  S(NF)  DIRECTION VECTOR.
*  RI  ETA2  TOLERANCE FOR POSITIVE DEFINITENESS OF THE HESSIAN MATRIX.
*  RI  ETA9  MAXIMUM FOR REAL NUMBERS.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RI  EPS9  TOLERANCE FOR ACTIVITY OF CONSTRAINTS.
*  RU  XDEL  TRUST REGION BOUND.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER.
*  RO  GMAX  MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE.
*  II  MFP  TYPE OF FEASIBLE POINT. MFP=1-ARBITRARY FEASIBLE POINT.
*         MFP=2-OPTIMUM FEASIBLE POINT. MFP=3-REPEATED SOLUTION.
*
* COMMON DATA :
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  II  NORMF  SCALING SPECIFICATION. NORMF=0-NO SCALING PERFORMED.
*         NORMF=1-SCALING FACTORS ARE DETERMINED AUTOMATICALLY.
*         NORMF=2-SCALING FACTORS ARE SUPPLIED BY USER.
*  IU  IDECF  DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION.
*         IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION.
*         IDECF=10-DIAGONAL MATRIX.
*  IU  NDECF  NUMBER OF DECOMPOSITIONS.
*  IO  ITERQ  TYPE OF FEASIBLE POINT. ITERQ=1-ARBITRARY FEASIBLE POINT.
*         ITERQ=2-OPTIMUM FEASIBLE POINT. ITERQ=-1 FEASIBLE POINT DOES
*         NOT EXISTS. ITERQ=-2 OPTIMUM FEASIBLE POINT DOES NOT EXISTS.
*
* SUBPROGRAMS USED :
*  S   PLMINS  DETERMINATION OF THE NEW ACTIVE SIMPLE BOUND.
*  S   PLMINL  DETERMINATION OF THE NEW ACTIVE LINEAR CONSTRAINT.
*  S   PLMINT  DETERMINATION OF THE NEW ACTIVE TRUST REGION BOUND.
*  S   PLADR1  ADDITION OF A NEW ACTIVE CONSTRAINT.
*  S   PLRMR0  CONSTRAIN DELETION.
*  S   PLSOB1  TRANSFORMATION OF THE LOCAL SOLUTION TO THE SOLUTION
*         OF THE ORIGINAL QP PROBLEM.
*  S   MXDPGF  GILL-MURRAY DECOMPOSITION OF A DENSE SYMMETRIC MATRIX.
*  S   MXDPGB  BACK SUBSTITUTION AFTER GILL-MURRAY DECOMPOSITION.
*  S   MXDPRB  BACK SUBSTITUTION.
*  S   MXDSMM  MATRIX VECTOR PRODUCT.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S   MXVINA  ABSOLUTE VALUES OF ELEMENTS OF AN INTEGER VECTOR.
*  S   MXVINV  CHANGE OF AN INTEGER VECTOR AFTER CONSTRAINT ADDITION.
*  S   MXVNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN.
*
      SUBROUTINE PLQDB1(NF,NC,X,IX,XL,XU,CF,CFD,IC,ICA,CL,CU,CG,CR,CZ,
     & G,GO,H,S,MFP,KBF,KBC,IDECF,ETA2,ETA9,EPS7,EPS9,UMAX,GMAX,N,ITERQ)
      INTEGER NF,NC,IX(*),IC(*),ICA(*),MFP,KBF,KBC,IDECF,N,ITERQ
      REAL*8 X(*),XL(*),XU(*),CF(*),CFD(*),CL(*),CU(*),CG(*),CR(*),
     & CZ(*),G(*),GO(*),H(*),S(*),ETA2,ETA9,EPS7,EPS9,UMAX,GMAX
      REAL*8 CON,TEMP,STEP,STEP1,STEP2,DMAX,PAR,SNORM
      INTEGER NCA,NCR,I,J,K,IOLD,JOLD,INEW,JNEW,KNEW,INF,IER,KREM,KC,
     & NRED
      INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      COMMON /STATSQP/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      CON=ETA9
      IF (IDECF.LT.0) IDECF=1
      IF (IDECF.EQ.0) THEN
C
C     GILL-MURRAY DECOMPOSITION
C
      TEMP=ETA2
      CALL MXDPGF(NF,H,INF,TEMP,STEP)
      NDEC=NDEC+1
      IDECF=1
      ENDIF
      IF (IDECF.GE.2.AND.IDECF.LE.8) THEN
      ITERQ=-10
      RETURN
      ENDIF
C
C     INITIATION
C
      NRED=0
      JOLD=0
      JNEW=0
      ITERQ=0
      DMAX=0.0D0
      IF (MFP.EQ.3) GO TO 1
      N=NF
      NCA=0
      NCR=0
      IF (KBF.GT.0) CALL MXVINA(NF,IX)
      IF (KBC.GT.0) CALL MXVINA(NC,IC)
C
C     DIRECTION DETERMINATION
C
    1 CALL MXVNEG(NF,GO,S)
      DO 2 J=1,NCA
      KC=ICA(J)
      IF (KC.GT.0) THEN
      CALL MXVDIR(NF,CZ(J),CG((KC-1)*NF+1),S,S)
      ELSE
      K=-KC
      S(K)=S(K)+CZ(J)
      ENDIF
    2 CONTINUE
      CALL MXVCOP(NF,S,G)
      IF (IDECF.EQ.1) THEN
      CALL MXDPGB(NF,H,S,0)
      ELSE
      CALL MXDSMM(NF,H,G,S)
      ENDIF
      IF (ITERQ.EQ.3) GO TO 7
C
C     CHECK OF FEASIBILITY
C
      INEW=0
      PAR=0.0D0
      CALL PLMINN(NF,NC,CF,CFD,IC,CL,CU,CG,S,EPS9,PAR,KBC,INEW,KNEW)
      CALL PLMINS(NF,IX,X,XL,XU,S,KBF,INEW,KNEW,EPS9,PAR)
      IF (INEW.EQ.0) THEN
C
C     SOLUTION ACHIEVED
C
      CALL MXVNEG(NF,G,G)
      ITERQ=2
      GO TO 7
      ELSE
      SNORM=0.0D0
      ENDIF
    4 IER=0
C
C     STEPSIZE DETERMINATION
C
      CALL PLADR1(NF,N,NC,ICA,CG,CR,H,S,G,EPS7,GMAX,UMAX,IDECF,INEW,
     & NADD,IER,1)
      CALL MXDPRB(NCA,CR,G,-1)
      IF (KNEW.LT.0) CALL MXVNEG(NCA,G,G)
C
C     PRIMAL STEPSIZE
C
      IF (IER.NE.0) THEN
      STEP1=CON
      ELSE
      STEP1=-PAR/UMAX
      ENDIF
C
C     DUAL STEPSIZE
C
      IOLD=0
      STEP2=CON
      DO 5 J=1,NCA
      KC=ICA(J)
      IF (KC.GE.0) THEN
      K=IC(KC)
      ELSE
      I=-KC
      K=IX(I)
      ENDIF
      IF (K.LE.-5) THEN
      ELSE IF ((K.EQ.-1.OR.K.EQ.-3.).AND.G(J).LE.0.0D0) THEN
      ELSE IF ((K.EQ.-2.OR.K.EQ.-4.).AND.G(J).GE.0.0D0) THEN
      ELSE
      TEMP=CZ(J)/G(J)
      IF (STEP2.GT.TEMP) THEN
      IOLD=J
      STEP2=TEMP
      ENDIF
      ENDIF
    5 CONTINUE
C
C     FINAL STEPSIZE
C
      STEP=MIN(STEP1,STEP2)
      IF (STEP.GE.CON) THEN
C
C     FEASIBLE SOLUTION DOES NOT EXIST
C
      ITERQ=-1
      GO TO 7
      ENDIF
C
C     NEW LAGRANGE MULTIPLIERS
C
      DMAX=STEP
      CALL MXVDIR(NCA,-STEP,G,CZ,CZ)
      SNORM=SNORM+SIGN(1,KNEW)*STEP
      PAR=PAR-(STEP/STEP1)*PAR
      IF (STEP.EQ.STEP1) THEN
      IF (N.LE.0) THEN
C
C     IMPOSSIBLE SITUATION
C
      ITERQ=-5
      GO TO 7
      ENDIF
C
C     CONSTRAINT ADDITION
C
      IF (IER.EQ.0) THEN
      N=N-1
      NCA=NCA+1
      NCR=NCR+NCA
      CZ(NCA)=SNORM
      ENDIF
      IF (INEW.GT.0) THEN
      KC=INEW
      CALL MXVINV(IC,KC,KNEW)
      ELSE IF (ABS(KNEW).EQ.1) THEN
      I=-INEW
      CALL MXVINV(IX,I,KNEW)
      ELSE
      I=-INEW
      IF (KNEW.GT.0) IX(I)=-3
      IF (KNEW.LT.0) IX(I)=-4
      ENDIF
      NRED=NRED+1
      NADD=NADD+1
      JNEW=INEW
      JOLD=0
      GO TO 1
      ELSE
C
C     CONSTRAINT DELETION
C
      DO 6 J=IOLD,NCA-1
      CZ(J)=CZ(J+1)
    6 CONTINUE
      CALL PLRMF0(NF,NC,IX,IC,ICA,CR,IC,G,N,IOLD,KREM,IER)
      NCR=NCR-NCA
      NCA=NCA-1
      JOLD=IOLD
      JNEW=0
      IF (KBC.GT.0) CALL MXVINA(NC,IC)
      IF (KBF.GT.0) CALL MXVINA(NF,IX)
      DO 8 J=1,NCA
      KC=ICA(J)
      IF (KC.GT.0) THEN
      IC(KC)=-IC(KC)
      ELSE
      KC=-KC
      IX(KC)=-IX(KC)
      ENDIF
    8 CONTINUE
      GO TO 4
      ENDIF
    7 CONTINUE
      RETURN
      END
* SUBROUTINE PLADR1               ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* TRIANGULAR DECOMPOSITION OF KERNEL OF THE GENERAL PROJECTION
* IS UPDATED AFTER CONSTRAINT ADDITION.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  IU  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC  NUMBER OF LINEARIZED CONSTRAINTS.
*  IU  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  H(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OR INVERSION OF THE
*         HESSIAN MATRIX APPROXIMATION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RO  G(NF)  VECTOR USED IN THE DUAL RANGE SPACE QUADRATIC PROGRAMMING
*         METHOD.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RO  GMAX  MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER.
*  RO  E  AUXILIARY VARIABLE.
*  RI  T  AUXILIARY VARIABLE.
*  IU  IDECF  DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION.
*         IDECF=1-GILL-MURRAY DECOMPOSITION. IDECF=9-INVERSION.
*         IDECF=10-DIAGONAL MATRIX.
*  II  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IU  NADD  NUMBER OF CONSTRAINT ADDITIONS.
*  IO  IER  ERROR INDICATOR.
*  II  JOB  SPECIFICATION OF COMPUTATION. OUTPUT VECTOR G IS NOT OR IS
*         COMPUTED IN CASE WHEN N.LE.0 IF JOB=0 OR JOB=1 RESPECTIVELY.
*
* SUBPROGRAMS USED :
*  S   MXDPGB  BACK SUBSTITUTION.
*  S   MXDPRB  BACK SUBSTITUTION.
*  S   MXDSMM  MATRIX-VECTOR PRODUCT.
*  S   MXDSMV  COPYING OF A ROW OF DENSE SYMMETRIC MATRIX.
*  S   MXVCOP  COPYING OF A VECTOR.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*
      SUBROUTINE PLADR1(NF,N,NC,ICA,CG,CR,H,S,G,EPS7,GMAX,UMAX,IDECF,
     & INEW,NADD,IER,JOB)
      INTEGER NF,N,NC,ICA(*),IDECF,INEW,NADD,IER,JOB
      REAL*8 CG(*),CR(*),H(*),S(*),G(*),EPS7,GMAX,UMAX
      REAL*8 MXVDOT
      INTEGER NCA,NCR,JCG,J,K,L
      IER=0
      IF (JOB.EQ.0.AND.N.LE.0) IER=2
      IF (INEW.EQ.0) IER=3
      IF (IDECF.NE.1.AND.IDECF.NE.9) IER=-2
      IF (IER.NE.0) RETURN
      NCA=NF-N
      NCR=NCA*(NCA+1)/2
      IF (INEW.GT.0) THEN
      JCG=(INEW-1)*NF+1
      IF (IDECF.EQ.1) THEN
      CALL MXVCOP(NF,CG(JCG),S)
      CALL MXDPGB(NF,H,S,0)
      ELSE
      CALL MXDSMM(NF,H,CG(JCG),S)
      ENDIF
      GMAX=MXVDOT(NF,CG(JCG),S)
      ELSE
      K=-INEW
      IF (IDECF.EQ.1) THEN
      CALL MXVSET(NF,0.0D0,S)
      S(K)=1.0D 0
      CALL MXDPGB(NF,H,S,0)
      ELSE
      CALL MXDSMV(NF,H,S,K)
      ENDIF
      GMAX=S(K)
      ENDIF
      DO 1 J=1,NCA
      L=ICA(J)
      IF (L.GT.0) THEN
      G(J)=MXVDOT(NF,CG((L-1)*NF+1),S)
      ELSE
      L=-L
      G(J)=S(L)
      ENDIF
    1 CONTINUE
      IF (N.EQ.0) THEN
      CALL MXDPRB(NCA,CR,G,1)
      UMAX=0.0D0
      IER=2
      RETURN
      ELSE IF (NCA.EQ.0) THEN
      UMAX=GMAX
      ELSE
      CALL MXDPRB(NCA,CR,G,1)
      UMAX=GMAX-MXVDOT(NCA,G,G)
      CALL MXVCOP(NCA,G,CR(NCR+1))
      ENDIF
      IF (UMAX.LE.EPS7*GMAX) THEN
      IER=1
      RETURN
      ELSE
      NCA=NCA+1
      NCR=NCR+NCA
      ICA(NCA)=INEW
      CR(NCR)=SQRT(UMAX)
      IF (JOB.EQ.0) THEN
      N=N-1
      NADD=NADD+1
      ENDIF
      ENDIF
      RETURN
      END
* SUBROUTINE PLDIRL               ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE NEW VALUES OF THE CONSTRAINT FUNCTIONS.
*
* PARAMETERS :
*  II  NC  NUMBER OF CONSTRAINTS.
*  RU  CF(NF)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  RI  CFD(NF)  VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT
*         FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  STEP  CURRENT STEPSIZE.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*
      SUBROUTINE PLDIRL(NC,CF,CFD,IC,STEP,KBC)
      INTEGER NC,IC(*),KBC
      DOUBLE PRECISION CF(*),CFD(*),STEP
      INTEGER KC
      IF (KBC.GT.0) THEN
      DO 1 KC=1,NC
      IF (IC(KC).GE.0.AND.IC(KC).LE.10) THEN
      CF(KC)=CF(KC)+STEP*CFD(KC)
      ELSE IF (IC(KC).LT.-10) THEN
      CF(KC)=CF(KC)+STEP*CFD(KC)
      ENDIF
    1 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PLDIRS               ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE NEW VECTOR OF VARIABLES.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  RU  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RI  S(NF)  DIRECTION VECTOR.
*  RI  STEP  CURRENT STEPSIZE.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*
      SUBROUTINE PLDIRS(NF,X,IX,S,STEP,KBF)
      INTEGER NF,IX(*),KBF
      DOUBLE PRECISION X(*),S(*),STEP
      INTEGER I
      DO 1 I=1,NF
      IF (KBF.LE.0) THEN
      X(I)=X(I)+STEP*S(I)
      ELSE IF (IX(I).GE.0.AND.IX(I).LE.10) THEN
      X(I)=X(I)+STEP*S(I)
      ELSE IF (IX(I).LT.-10) THEN
      X(I)=X(I)+STEP*S(I)
      ENDIF
    1 CONTINUE
      RETURN
      END
* SUBROUTINE PLINIT             ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE INITIAL POINT WHICH SATISFIES SIMPLE BOUNDS.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  RU  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  EPS9  TOLERANCE FOR ACTIVE CONSTRAINTS.
*  IO  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IO  IND  INDICATOR. IF IND.NE.0 THEN TRUST REGION BOUNDS CANNOT
*         BE SATISFIED.
*
* SUBPROGRAMS USED :
*  S   PLNEWS  TEST ON ACTIVITY OF A GIVEN SIMPLE BOUND.
*
      SUBROUTINE PLINIT(NF,X,IX,XL,XU,EPS9,KBF,INEW,IND)
      INTEGER NF,IX(*),KBF,INEW,IND
      DOUBLE PRECISION X(*),XL(*),XU(*),EPS9
      INTEGER I
      IND=0
      IF (KBF.GT.0) THEN
      DO 1 I=1,NF
      CALL PLNEWS(X,IX,XL,XU,EPS9,I,INEW)
      IF (IX(I).LT.5) THEN
      ELSE IF (IX(I).EQ.5) THEN
      IX(I)=-5
      ELSE IF (IX(I).EQ.11.OR.IX(I).EQ.13) THEN
      X(I)=XL(I)
      IX(I)=10-IX(I)
      ELSE IF (IX(I).EQ.12.OR.IX(I).EQ.14) THEN
      X(I)=XU(I)
      IX(I)=10-IX(I)
      ENDIF
    1 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PLMAXL               ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE MAXIMUM STEPSIZE USING LINEAR CONSTRAINTS.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CURRENT LINEAR CONSTRAINTS.
*  RI  CF(NF)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCYIONS.
*  RO  CFD(NF)  VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT
*         FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  S(NF)  DIRECTION VECTOR.
*  RO  STEP  MAXIMUM STEPSIZE.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  II  KREM  INDICATION OF LINEARLY DEPENDENT GRADIENTS.
*  IO  INEW  INDEX OF THE NEW ACTIVE FUNCTION.
*
* SUBPROGRAMS USED :
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*
      SUBROUTINE PLMAXL(NF,NC,CF,CFD,IC,CL,CU,CG,S,STEP,KBC,KREM,INEW)
      INTEGER NF,NC,IC(*),KBC,KREM,INEW
      DOUBLE PRECISION CF(*),CFD(*),CL(*),CU(*),CG(*),S(*),STEP
      DOUBLE PRECISION TEMP,MXVDOT
      INTEGER JCG,KC
      IF (KBC.GT.0) THEN
      JCG=1
      DO 1 KC=1,NC
      IF (KREM.GT.0.AND.IC(KC).GT.10) IC(KC)=IC(KC)-10
      IF (IC(KC).GT.0.AND.IC(KC).LE.10) THEN
      TEMP=MXVDOT(NF,CG(JCG),S)
      CFD(KC)=TEMP
      IF (TEMP.LT.0.0D 0) THEN
      IF (IC(KC).EQ.1.OR.IC(KC).GE.3) THEN
      TEMP=(CL(KC)-CF(KC))/TEMP
      IF (TEMP.LE.STEP) THEN
      INEW=KC
      STEP=TEMP
      ENDIF
      ENDIF
      ELSE IF (TEMP.GT.0.0D 0) THEN
      IF (IC(KC).EQ.2.OR.IC(KC).GE.3) THEN
      TEMP=(CU(KC)-CF(KC))/TEMP
      IF (TEMP.LE.STEP) THEN
      INEW=KC
      STEP=TEMP
      ENDIF
      ENDIF
      ENDIF
      ELSE IF (IC(KC).LT.-10) THEN
      TEMP=MXVDOT(NF,CG(JCG),S)
      CFD(KC)=TEMP
      IF (TEMP.GT.0.0D 0) THEN
      IF (IC(KC).EQ.-11.OR.IC(KC).EQ.-13.OR.IC(KC).EQ.-15) THEN
      TEMP=(CL(KC)-CF(KC))/TEMP
      IF (TEMP.LE.STEP) THEN
      INEW=KC
      STEP=TEMP
      ENDIF
      ENDIF
      ELSE IF (TEMP.LT.0.0D 0) THEN
      IF (IC(KC).EQ.-12.OR.IC(KC).EQ.-14.OR.IC(KC).EQ.-16) THEN
      TEMP=(CU(KC)-CF(KC))/TEMP
      IF (TEMP.LE.STEP) THEN
      INEW=KC
      STEP=TEMP
      ENDIF
      ENDIF
      ENDIF
      ENDIF
      JCG=JCG+NF
    1 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PLMAXS               ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE MAXIMUM STEPSIZE USING THE SIMPLE BOUNDS
* FOR VARIABLES.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  S(NF)  DIRECTION VECTOR.
*  RO  STEP  MAXIMUM STEPSIZE.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  IO  KREM  INDICATION OF LINEARLY DEPENDENT GRADIENTS.
*  IO  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*
      SUBROUTINE PLMAXS(NF,X,IX,XL,XU,S,STEP,KBF,KREM,INEW)
      INTEGER NF,IX(*),KBF,KREM,INEW
      DOUBLE PRECISION X(*),XL(*),XU(*),S(*),STEP
      DOUBLE PRECISION TEMP
      INTEGER I
      IF (KBF.GT.0) THEN
      DO 1 I=1,NF
      IF (KREM.GT.0.AND.IX(I).GT.10) IX(I)=IX(I)-10
      IF (IX(I).GT.0.AND.IX(I).LE.10) THEN
      IF (S(I).LT.0.0D 0) THEN
      IF (IX(I).EQ.1.OR.IX(I).GE.3) THEN
      TEMP=(XL(I)-X(I))/S(I)
      IF (TEMP.LE.STEP) THEN
      INEW=-I
      STEP=TEMP
      ENDIF
      ENDIF
      ELSE IF (S(I).GT.0.0D 0) THEN
      IF (IX(I).EQ.2.OR.IX(I).GE.3) THEN
      TEMP=(XU(I)-X(I))/S(I)
      IF (TEMP.LE.STEP) THEN
      INEW=-I
      STEP=TEMP
      ENDIF
      ENDIF
      ENDIF
      ENDIF
    1 CONTINUE
      ENDIF
      KREM=0
      RETURN
      END
* SUBROUTINE PLNEWL             ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* TEST ON ACTIVITY OF A GIVEN LINEAR CONSTRAINT.
*
* PARAMETERS :
*  II  KC  INDEX OF A GIVEN CONSTRAINT.
*  RI  CF(NC)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  IU  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  EPS9  TOLERANCE FOR ACTIVE CONSTRAINTS.
*  IO  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*
      SUBROUTINE PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW)
      INTEGER KC,IC(*),INEW
      DOUBLE PRECISION CF(*),CL(*),CU(*),EPS9
      DOUBLE PRECISION TEMP
      IF (IC(KC).LT.-10) IC(KC)=-IC(KC)-10
      IF (IC(KC).LE.0) THEN
      ELSE IF (IC(KC).EQ.1) THEN
      TEMP=EPS9*MAX(ABS(CL(KC)),1.0D 0)
      IF (CF(KC).GT.CL(KC)+TEMP) THEN
      ELSE IF (CF(KC).GE.CL(KC)-TEMP) THEN
      IC(KC)=11
      INEW=KC
      ELSE
      IC(KC)=-11
      ENDIF
      ELSE IF (IC(KC).EQ.2) THEN
      TEMP=EPS9*MAX(ABS(CU(KC)),1.0D 0)
      IF (CF(KC).LT.CU(KC)-TEMP) THEN
      ELSE IF (CF(KC).LE.CU(KC)+TEMP) THEN
      IC(KC)=12
      INEW=KC
      ELSE
      IC(KC)=-12
      ENDIF
      ELSE IF (IC(KC).EQ.3.OR.IC(KC).EQ.4) THEN
      TEMP=EPS9*MAX(ABS(CL(KC)),1.0D 0)
      IF (CF(KC).GT.CL(KC)+TEMP) THEN
      TEMP=EPS9*MAX(ABS(CU(KC)),1.0D 0)
      IF (CF(KC).LT.CU(KC)-TEMP) THEN
      ELSE IF (CF(KC).LE.CU(KC)+TEMP) THEN
      IC(KC)=14
      INEW=KC
      ELSE
      IC(KC)=-14
      ENDIF
      ELSE IF (CF(KC).GE.CL(KC)-TEMP) THEN
      IC(KC)=13
      INEW=KC
      ELSE
      IC(KC)=-13
      ENDIF
      ELSE IF (IC(KC).EQ.5.OR.IC(KC).EQ.6) THEN
      TEMP=EPS9*MAX(ABS(CL(KC)),1.0D 0)
      IF (CF(KC).GT.CL(KC)+TEMP) THEN
      TEMP=EPS9*MAX(ABS(CU(KC)),1.0D 0)
      IF (CF(KC).LT.CU(KC)-TEMP) THEN
      ELSE IF (CF(KC).LE.CU(KC)+TEMP) THEN
      IC(KC)=16
      INEW=KC
      ELSE
      IC(KC)=-16
      ENDIF
      ELSE IF (CF(KC).GE.CL(KC)-TEMP) THEN
      IC(KC)=15
      INEW=KC
      ELSE
      IC(KC)=-15
      ENDIF
      ENDIF
      RETURN
      END
* SUBROUTINE PLMINN             ALL SYSTEMS                   97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE NEW ACTIVE LINEAR CONSTRAINT.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  RI  CF(NC)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  RO  CFD(NC)  VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT
*            FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  S(NF)  DIRECTION VECTOR.
*  RI  EPS9  TOLERANCE FOR ACTIVE CONSTRAINTS.
*  RA  PAR  AUXILIARY VARIABLE.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  IO  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IO  KNEW  SIGNUM OF THE NEW ACTIVE NORMAL.
*
* SUBPROGRAMS USED :
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*
      SUBROUTINE PLMINN(NF,NC,CF,CFD,IC,CL,CU,CG,S,EPS9,PAR,KBC,INEW,
     & KNEW)
      INTEGER NF,NC,IC(*),KBC,INEW,KNEW
      REAL*8 CF(*),CFD(*),CL(*),CU(*),CG(*),S(*),EPS9,PAR
      REAL*8 TEMP,POM,MXVDOT
      INTEGER JCG,KC
      IF (KBC.GT.0) THEN
      JCG=1
      DO 1 KC=1,NC
      IF (IC(KC).GT.0) THEN
      TEMP=MXVDOT(NF,CG(JCG),S)
      CFD(KC)=TEMP
      TEMP=CF(KC)+TEMP
      IF (IC(KC).EQ.1.OR.IC(KC).GE.3) THEN
      POM=TEMP-CL(KC)
      IF (POM.LT.MIN(PAR,-EPS9*MAX(ABS(CL(KC)),1.0D 0))) THEN
      INEW=KC
      KNEW= 1
      PAR=POM
      ENDIF
      ENDIF
      IF (IC(KC).EQ.2.OR.IC(KC).GE.3) THEN
      POM=CU(KC)-TEMP
      IF (POM.LT.MIN(PAR,-EPS9*MAX(ABS(CU(KC)),1.0D 0))) THEN
      INEW=KC
      KNEW=-1
      PAR=POM
      ENDIF
      ENDIF
      ENDIF
      JCG=JCG+NF
    1 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PLMINS             ALL SYSTEMS                   91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF THE NEW ACTIVE SIMPLE BOUND.
*
* PARAMETERS :
*  II  NF DECLARED NUMBER OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RI  XO(NF)  SAVED VECTOR OF VARIABLES.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  S(NF)  DIRECTION VECTOR.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  IO  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IO  KNEW  SIGNUM OF THE NEW NORMAL.
*  RI  EPS9  TOLERANCE FOR ACTIVE CONSTRAINTS.
*  RA  PAR  AUXILIARY VARIABLE.
*
      SUBROUTINE PLMINS(NF,IX,XO,XL,XU,S,KBF,INEW,KNEW,EPS9,PAR)
*     .. Scalar Arguments ..
      DOUBLE PRECISION EPS9,PAR
      INTEGER          INEW,KBF,KNEW,NF
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION S(*),XL(*),XO(*),XU(*)
      INTEGER          IX(*)
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION POM,TEMP
      INTEGER          I
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC        ABS,MAX,MIN
*     ..
      IF (KBF.GT.0) THEN
          DO 10 I = 1,NF
              IF (IX(I).GT.0) THEN
                  TEMP = 1.0D0
                  IF (IX(I).EQ.1 .OR. IX(I).GE.3) THEN
                      POM = XO(I) + S(I)*TEMP - XL(I)
                      IF (POM.LT.MIN(PAR,-EPS9*MAX(ABS(XL(I)),
     +                    TEMP))) THEN
                          INEW = -I
                          KNEW = 1
                          PAR = POM
                      END IF

                  END IF

                  IF (IX(I).EQ.2 .OR. IX(I).GE.3) THEN
                      POM = XU(I) - S(I)*TEMP - XO(I)
                      IF (POM.LT.MIN(PAR,-EPS9*MAX(ABS(XU(I)),
     +                    TEMP))) THEN
                          INEW = -I
                          KNEW = -1
                          PAR = POM
                      END IF

                  END IF

              END IF

   10     CONTINUE
      END IF

      RETURN

      END
* SUBROUTINE PLNEWS             ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* TEST ON ACTIVITY OF A GIVEN SIMPLE BOUND.
*
* PARAMETERS :
*  RI  X(NF)  VECTOR OF VARIABLES.
*  IU  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  EPS9  TOLERANCE FOR ACTIVE CONSTRAINTS.
*  II  I  INDEX OF TESTED SIMPLE BOUND.
*  IO  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*
      SUBROUTINE PLNEWS(X,IX,XL,XU,EPS9,I,INEW)
      INTEGER IX(*),I,INEW
      DOUBLE PRECISION X(*),XL(*),XU(*),EPS9
      DOUBLE PRECISION TEMP
      TEMP=1.0D 0
      IF (IX(I).LE.0) THEN
      ELSE IF (IX(I).EQ.1) THEN
      IF (X(I).LE.XL(I)+EPS9*MAX(ABS(XL(I)),TEMP)) THEN
      IX(I)=11
      INEW=-I
      ENDIF
      ELSE IF (IX(I).EQ.2) THEN
      IF (X(I).GE.XU(I)-EPS9*MAX(ABS(XU(I)),TEMP)) THEN
      IX(I)=12
      INEW=-I
      ENDIF
      ELSE IF (IX(I).EQ.3.OR.IX(I).EQ.4) THEN
      IF (X(I).LE.XL(I)+EPS9*MAX(ABS(XL(I)),TEMP)) THEN
      IX(I)=13
      INEW=-I
      ENDIF
      IF (X(I).GE.XU(I)-EPS9*MAX(ABS(XU(I)),TEMP)) THEN
      IX(I)=14
      INEW=-I
      ENDIF
      ENDIF
      RETURN
      END
* SUBROUTINE PLREDL               ALL SYSTEMS                   98/12/01
C PORTABILITY : ALL SYSTEMS
C 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* TRANSFORMATION OF THE INCOMPATIBLE QUADRATIC PROGRAMMING SUBPROBLEM.
*
* PARAMETERS :
*  II  NC  NUMBER OF CURRENT LINEAR CONSTRAINTS.
*  RI  CF(NF)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*
      SUBROUTINE PLREDL(NC,CF,IC,CL,CU,KBC)
      INTEGER NC,IC(NC),KBC,K
      REAL*8 CF(*),CL(*),CU(*)
      REAL*8 TEMP
      INTEGER KC
      IF (KBC.GT.0) THEN
      DO 1 KC=1,NC
      K=IC(KC)
      IF (ABS(K).EQ.1.OR.ABS(K).EQ.3.OR.ABS(K).EQ.4) THEN
      TEMP=(CF(KC)-CL(KC))
      IF (TEMP.LT.0) CF(KC)=CL(KC)+0.1D 0*TEMP
      ENDIF
      IF (ABS(K).EQ.2.OR.ABS(K).EQ.3.OR.ABS(K).EQ.4) THEN
      TEMP=(CF(KC)-CU(KC))
      IF (TEMP.GT.0) CF(KC)=CU(KC)+0.1D 0*TEMP
      ENDIF
      IF (ABS(K).EQ.5.OR.ABS(K).EQ.6) THEN
      TEMP=(CF(KC)-CL(KC))
      CF(KC)=CL(KC)+0.1D 0*TEMP
      ENDIF
    1 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PLRMF0             ALL SYSTEMS                   91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* OPERATIONS AFTER CONSTRAINT DELETION.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  II  IA(NA)  VECTOR CONTAINING TYPES OF DEVIATIONS.
*  IU  IAA(NF+1)  VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS.
*  RU  AR((NF+1)*(NF+2)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RA  S(NF+1)  AUXILIARY VECTOR.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  IOLD  INDEX OF THE OLD ACTIVE CONSTRAINT.
*  IO  KREM  AUXILIARY VARIABLE.
*  IO  IER  ERROR INDICATOR.
*
* SUBPROGRAMS USED :
*  S   PLRMR0  CORRECTION OF KERNEL OF THE ORTHOGONAL PROJECTION
*         AFTER CONSTRAINT DELETION.
*
      SUBROUTINE PLRMF0(NF,NC,IX,IA,IAA,AR,IC,S,N,IOLD,KREM,IER)
*     .. Scalar Arguments ..
      INTEGER          IER,IOLD,KREM,N,NC,NF
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION AR(*),S(*)
      INTEGER          IA(*),IAA(*),IC(*),IX(*)
*     ..
*     .. Scalars in Common ..
      INTEGER          NADD,NDEC,NFG,NFH,NFV,NIT,NREM,NRES
*     ..
*     .. Local Scalars ..
      INTEGER          L
*     ..
*     .. External Subroutines ..
      EXTERNAL         PLRMR0
*     ..
*     .. Common blocks ..
      COMMON           /STATSQP/NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
*     ..
      CALL PLRMR0(NF,IAA,AR,S,N,IOLD,KREM,IER)
      N = N + 1
      NREM = NREM + 1
      L = IAA(NF-N+1)
      IF (L.GT.NC) THEN
          L = L - NC
          IA(L) = -IA(L)

      ELSE IF (L.GT.0) THEN
          IC(L) = -IC(L)

      ELSE
          L = -L
          IX(L) = -IX(L)
      END IF

      RETURN

      END
* SUBROUTINE PLRMR0               ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* TRIANGULAR DECOMPOSITION OF KERNEL OF THE ORTHOGONAL PROJECTION IS
* UPDATED AFTER CONSTRAINT DELETION.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  IU  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RU  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RA  G(NF)  AUXILIARY VECTOR.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  IOLD  INDEX OF THE OLD ACTIVE CONSTRAINT.
*  IO  KREM  AUXILIARY VARIABLE.
*  IO  IER  ERROR INDICATOR.
*
* SUBPROGRAMS USED :
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVORT  DETERMINATION OF AN ELEMENTARY ORTHOGONAL MATRIX FOR
*         PLANE ROTATION.
*  S   MXVROT  PLANE ROTATION OF A VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PLRMR0(NF,ICA,CR,G,N,IOLD,KREM,IER)
*     .. Scalar Arguments ..
      INTEGER          IER,IOLD,KREM,N,NF
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION CR(*),G(*)
      INTEGER          ICA(*)
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION CK,CL
      INTEGER          I,J,K,KC,L,NCA
*     ..
*     .. External Subroutines ..
      EXTERNAL         MXVCOP,MXVORT,MXVROT,MXVSET
*     ..
      NCA = NF - N
      IF (IOLD.LT.NCA) THEN
          K = IOLD* (IOLD-1)/2
          KC = ICA(IOLD)
          CALL MXVCOP(IOLD,CR(K+1),G)
          CALL MXVSET(NCA-IOLD,0.0D0,G(IOLD+1))
          K = K + IOLD
          DO 20 I = IOLD + 1,NCA
              K = K + I
              CALL MXVORT(CR(K-1),CR(K),CK,CL,IER)
              CALL MXVROT(G(I-1),G(I),CK,CL,IER)
              L = K
              DO 10 J = I,NCA - 1
                  L = L + J
                  CALL MXVROT(CR(L-1),CR(L),CK,CL,IER)
   10         CONTINUE
   20     CONTINUE
          K = IOLD* (IOLD-1)/2
          DO 30 I = IOLD,NCA - 1
              L = K + I
              ICA(I) = ICA(I+1)
              CALL MXVCOP(I,CR(L+1),CR(K+1))
              K = L
   30     CONTINUE
          ICA(NCA) = KC
          CALL MXVCOP(NCA,G,CR(K+1))
      END IF

      KREM = 1
      RETURN

      END
* SUBROUTINE PLSETC             ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF INITIAL VALUES OF THE CONSTRAINT FUNCTIONS.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CURRENT LINEAR CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RI  XO(NF)  SAVED VECTOR OF VARIABLES.
*  RU  CF(NF)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCYIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  CG(NF*MCL)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RA  S(NF)  AUXILIARY VECTOR.
*
* SUBPROGRAMS USED :
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*  S   MXVMUL  DIAGONAL PREMULTIPLICATION OF A VECTOR.
*
      SUBROUTINE PLSETC(NF,NC,X,XO,CF,IC,CG,S)
      INTEGER NF,NC,IC(*)
      DOUBLE PRECISION X(*),XO(*),CF(*),CG(*),S(*)
      DOUBLE PRECISION MXVDOT
      INTEGER JCG,KC
      CALL MXVDIF(NF,X,XO,S)
      JCG=0
      DO 1 KC=1,NC
      IF (IC(KC).NE.0) CF(KC)=CF(KC)+MXVDOT(NF,S,CG(JCG+1))
      JCG=JCG+NF
    1 CONTINUE
      RETURN
      END
* SUBROUTINE PLSETG             ALL SYSTEMS                   97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* GRADIENT DETERMINATION IN THE FIRST PHASE OF LP SUBROUTINE.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RO  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  IO  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*
* SUBPROGRAMS USED :
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PLSETG(NF,NC,IC,CG,G,INEW)
      INTEGER NF,NC,IC(*),INEW
      DOUBLE PRECISION CG(*),G(*)
      INTEGER KC
      CALL MXVSET(NF,0.0D 0,G)
      INEW=0
      DO 4 KC=1,NC
      IF (IC(KC).GE.-10) THEN
      ELSE IF (IC(KC).EQ.-11.OR.IC(KC).EQ.-13.OR.IC(KC).EQ.-15) THEN
      CALL MXVDIR(NF,-1.0D 0,CG((KC-1)*NF+1),G,G)
      INEW=1
      ELSE IF (IC(KC).EQ.-12.OR.IC(KC).EQ.-14.OR.IC(KC).EQ.-16) THEN
      CALL MXVDIR(NF, 1.0D 0,CG((KC-1)*NF+1),G,G)
      INEW=1
      ENDIF
    4 CONTINUE
      RETURN
      END
* SUBROUTINE PLTLAG               ALL SYSTEMS                97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER IS
* COMPUTED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC  NUMBER OF LINEARIZED CONSTRAINTS.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  II  IA(NA)  VECTOR CONTAINING TYPES OF DEVIATIONS.
*  II  IAA(NF+1)  VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS.
*  RI  AZ(NF+1)  VECTOR OF LAGRANGE MULTIPLIERS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  EPS7  TOLERANCE FOR LINEAR AND QUADRATIC PROGRAMMING.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
*  IO  IOLD  INDEX OF THE REMOVED CONSTRAINT.
*
      SUBROUTINE PLTLAG(NF,N,NC,IX,IA,IAA,AZ,IC,EPS7,UMAX,IOLD)
      INTEGER NF,N,NC,IX(*),IA(*),IAA(*),IC(*),IOLD
      DOUBLE PRECISION AZ(*),EPS7,UMAX
      DOUBLE PRECISION TEMP
      INTEGER NAA,J,K,L
      IOLD=0
      UMAX=0.0D 0
      NAA=NF-N
      DO 2 J=1,NAA
      TEMP=AZ(J)
      L=IAA(J)
      IF (L.GT.NC) THEN
      L=L-NC
      K=IA(L)
      ELSE IF (L.GT.0) THEN
      K=IC(L)
      ELSE
      L=-L
      K=IX(L)
      ENDIF
      IF (K.LE.-5) THEN
      ELSE IF ((K.EQ.-1.OR.K.EQ.-3).AND.UMAX+TEMP.GE.0.0D 0) THEN
      ELSE IF ((K.EQ.-2.OR.K.EQ.-4).AND.UMAX-TEMP.GE.0.0D 0) THEN
      ELSE
      IOLD=J
      UMAX=ABS(TEMP)
      ENDIF
    2 CONTINUE
      IF (UMAX.LE.EPS7) IOLD=0
      RETURN
      END
* SUBROUTINE PLTRBG               ALL SYSTEMS                97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* GRADIENT OF THE OBJECTIVE FUNCTION IS SCALED AND REDUCED. LAGRANGE
* MULTIPLIERS ARE DETERMINED. TEST VALUES GMAX AND UMAX ARE COMPUTED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CURRENT LINEAR CONSTRAINTS.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE. VECTOR CZ(1,NF) CONTAINS LAGRANGE
*         MULTIPLIERS BEING DETERMINED.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RO  GN(NF)  TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION IF IT IS
*         NONZERO.
*  RI  EPS7  TOLERANCE FOR LINEAR AND QUADRATIC PROGRAMMING.
*  RO  GMAX  NORM OF THE TRANSFORMED GRADIENT.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
*  IO  IOLD  INDEX OF THE REMOVED CONSTRAINT.
*
* SUBPROGRAMS USED :
*  S   PLVLAG  GRADIENT IS PREMULTIPLIED BY THE MATRIX WHOSE COLUMNS
*         ARE NORMALS OF THE ACTIVE CONSTRAINTS.
*  S   PLTLAG  COMPUTATION OF THE MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE
*         LAGRANGE MULTIPLIER.
*  S   MXDRMM  PREMULTIPLICATION OF A VECTOR BY A ROWWISE STORED DENSE
*         RECTANGULAR MATRIX.
*  S   MXDPRB  BACK SUBSTITUTION AFTER A CHOLESKI DECOMPOSITION.
*  RF  MXVMAX  L-INFINITY NORM OF A VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PLTRBG(NF,N,NC,IX,IC,ICA,CG,CR,CZ,G,GN,EPS7,GMAX,UMAX,
     & IOLD)
      INTEGER NF,N,NC,IX(*),IC(*),ICA(*),IOLD
      DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*),EPS7,GMAX,UMAX
      DOUBLE PRECISION MXVMAX
      INTEGER NCA,NCZ
      GMAX=0.0D 0
      IF (N.GT.0) THEN
      CALL MXDRMM(NF,N,CZ,G,GN)
      GMAX=MXVMAX(N,GN)
      ENDIF
      IF (NF.GT.N.AND.GMAX.LE.EPS7) THEN
      NCA=NF-N
      NCZ=N*NF
      CALL PLVLAG(NF,N,NC,ICA,CG,CG,G,CZ(NCZ+1))
      CALL MXDPRB(NCA,CR,CZ(NCZ+1),0)
      CALL PLTLAG(NF,N,NC,IX,IC,ICA,CZ(NCZ+1),IC,EPS7,UMAX,IOLD)
      IF (UMAX.LE.EPS7) IOLD=0
      CALL MXVSET(N,0.0D 0,GN)
      GMAX=0.0D 0
      ELSE
      IOLD=0
      UMAX=0.0D 0
      ENDIF
      RETURN
      END
* SUBROUTINE PLVLAG               ALL SYSTEMS                97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* GRADIENT OF THE OBJECTIVE FUNCTION IS PREMULTIPLIED BY TRANSPOSE
* OF THE MATRIX WHOSE COLUMNS ARE NORMALS OF CURRENT ACTIVE CONSTRAINTS
* AND GRADIENTS OF CURRENT ACTIVE FUNCTIONS.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC  NUMBER OF LINEARIZED CONSTRAINTS.
*  II  IAA(NF+1)  VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS.
*  RI  AG(NF*NA)  VECTOR CONTAINING SCALING PARAMETERS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RO  GN(NF+1)  OUTPUT VECTOR.
*
* SUBPROGRAMS USED :
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*
      SUBROUTINE PLVLAG(NF,N,NC,IAA,AG,CG,G,GN)
      INTEGER NF,N,NC,IAA(*)
      DOUBLE PRECISION AG(*),CG(*),G(*),GN(*)
      DOUBLE PRECISION MXVDOT
      INTEGER NAA,J,L
      NAA=NF-N
      DO 1 J=1,NAA
      L=IAA(J)
      IF (L.GT.NC) THEN
      L=L-NC
      GN(J)=MXVDOT(NF,AG((L-1)*NF+1),G)
      ELSE IF (L.GT.0) THEN
      GN(J)=MXVDOT(NF,CG((L-1)*NF+1),G)
      ELSE
      L=-L
      GN(J)=G(L)
      ENDIF
    1 CONTINUE
      RETURN
      END
* SUBROUTINE PNINT1                ALL SYSTEMS                91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH WITH DIRECTIONAL
* DERIVATIVES.
*
* PARAMETERS :
*  RI  RL  LOWER VALUE OF THE STEPSIZE PARAMETER.
*  RI  RU  UPPER VALUE OF THE STEPSIZE PARAMETER.
*  RI  FL  VALUE OF THE OBJECTIVE FUNCTION FOR R=RL.
*  RI  FU  VALUE OF THE OBJECTIVE FUNCTION FOR R=RU.
*  RI  PL  DIRECTIONAL DERIVATIVE FOR R=RL.
*  RI  PU  DIRECTIONAL DERIVATIVE FOR R=RU.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER OBTAINED.
*  II  MODE  MODE OF LINE SEARCH.
*  II  MTYP  METHOD SELECTION. MTYP=1-BISECTION. MTYP=2-QUADRATIC
*         INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE).
*         MTYP=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL
*         DERIVATIVES). MTYP=4-CUBIC INTERPOLATION. MTYP=5-CONIC
*         INTERPOLATION.
*  IO  MERR  ERROR INDICATOR. MERR=0 FOR NORMAL RETURN.
*
* METHOD :
* EXTRAPOLATION OR INTERPOLATION WITH STANDARD MODEL FUNCTIONS.
*
      SUBROUTINE PNINT1(RL,RU,FL,FU,PL,PU,R,MODE,MTYP,MERR)
      REAL*8 RL, RU, FL, FU, PL, PU, R
      INTEGER MODE,MTYP,MERR,NTYP
      REAL*8 A,B,C,D,DIS,DEN
      REAL*8 C1L,C1U,C2L,C2U,C3L
      PARAMETER (C1L=1.1D 0,C1U=1.0D 3,C2L=1.0D-2,C2U=0.9D 0,
     & C3L=0.1D 0)
      MERR=0
      IF (MODE.LE.0) RETURN
      IF (PL.GE.0.0D 0) THEN
      MERR=2
      RETURN
      ELSE IF (RU.LE.RL) THEN
      MERR=3
      RETURN
      ENDIF
      DO 1 NTYP=MTYP,1,-1
      IF (NTYP.EQ.1) THEN
C
C     BISECTION
C
      IF (MODE.EQ.1) THEN
      R=4.0D 0*RU
      RETURN
      ELSE
      R=0.5D 0*(RL+RU)
      RETURN
      ENDIF
      ELSE IF (NTYP.EQ.MTYP) THEN
      A = (FU-FL)/(PL*(RU-RL))
      B = PU/PL
      ENDIF
      IF (NTYP.EQ.2) THEN
C
C     QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH ONE DIRECTIONAL
C     DERIVATIVE
C
      DEN = 2.0D 0*(1.0D 0-A)
      ELSE IF (NTYP.EQ.3) THEN
C
C     QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH TWO DIRECTIONAL
C     DERIVATIVES
C
      DEN = 1.0D 0 - B
      ELSE IF (NTYP.EQ.4) THEN
C
C     CUBIC EXTRAPOLATION OR INTERPOLATION
C
      C = B - 2.0D 0*A + 1.0D 0
      D = B - 3.0D 0*A + 2.0D 0
      DIS = D*D - 3.0D 0*C
      IF (DIS.LT.0.0D 0) GO TO 1
      DEN = D + SQRT(DIS)
      ELSE IF (NTYP.EQ.5) THEN
C
C     CONIC EXTRAPOLATION OR INTERPOLATION
C
      DIS = A*A - B
      IF (DIS.LT.0.0D 0) GO TO 1
      DEN = A + SQRT(DIS)
      IF (DEN.LE.0.0D 0) GO TO 1
      DEN = 1.0D 0 - B*(1.0D 0/DEN)**3
      ENDIF
      IF (MODE.EQ.1.AND.DEN.GT.0.0D 0.AND.DEN.LT.1.0D 0) THEN
C
C     EXTRAPOLATION ACCEPTED
C
      R = RL + (RU-RL)/DEN
      R = MAX(R,C1L*RU)
      R = MIN(R,C1U*RU)
      RETURN
      ELSE IF (MODE.EQ.2.AND.DEN.GT.1.0D 0) THEN
C
C     INTERPOLATION ACCEPTED
C
      R = RL + (RU-RL)/DEN
      IF (RL.EQ.0.0D 0) THEN
      R = MAX(R,RL+C2L*(RU-RL))
      ELSE
      R = MAX(R,RL+C3L*(RU-RL))
      ENDIF
      R = MIN(R,RL+C2U*(RU-RL))
      RETURN
      ENDIF
    1 CONTINUE
      END
* SUBROUTINE PNINT3                ALL SYSTEMS                91/12/01
* PORTABILITY : ALL SYSTEMS
* 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH WITHOUT DIRECTIONAL
* DERIVATIVES.
*
* PARAMETERS :
*  RI  RO  INITIAL VALUE OF THE STEPSIZE PARAMETER.
*  RI  RL  LOWER VALUE OF THE STEPSIZE PARAMETER.
*  RI  RU  UPPER VALUE OF THE STEPSIZE PARAMETER.
*  RI  RI  INNER VALUE OF THE STEPSIZE PARAMETER.
*  RI  FO  VALUE OF THE OBJECTIVE FUNCTION FOR R=RO.
*  RI  FL  VALUE OF THE OBJECTIVE FUNCTION FOR R=RL.
*  RI  FU  VALUE OF THE OBJECTIVE FUNCTION FOR R=RU.
*  RI  FI  VALUE OF THE OBJECTIVE FUNCTION FOR R=RI.
*  RO  PO  INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER OBTAINED.
*  II  MODE  MODE OF LINE SEARCH.
*  II  MTYP  METHOD SELECTION. MTYP=1-BISECTION. MTYP=2-TWO POINT
*         QUADRATIC INTERPOLATION. MTYP=2-THREE POINT QUADRATIC
*         INTERPOLATION.
*  IO  MERR  ERROR INDICATOR. MERR=0 FOR NORMAL RETURN.
*
* METHOD :
* EXTRAPOLATION OR INTERPOLATION WITH STANDARD MODEL FUNCTIONS.
*
      SUBROUTINE PNINT3(RO,RL,RU,RI,FO,FL,FU,FI,PO,R,MODE,MTYP,MERR)
*     .. Parameters ..
      DOUBLE PRECISION ZERO,HALF,ONE,TWO,THREE,C1L,C1U,C2L,C2U,C3L
      PARAMETER        (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0,
     +                 THREE=3.0D0,C1L=1.1D0,C1U=1.0D3,C2L=1.0D-2,
     +                 C2U=0.9D0,C3L=1.0D-1)
*     ..
*     .. Scalar Arguments ..
      DOUBLE PRECISION FI,FL,FO,FU,PO,R,RI,RL,RO,RU
      INTEGER          MERR,MODE,MTYP
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION AI,AL,AU,DEN,DIS
      INTEGER          NTYP
      LOGICAL          L1,L2
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC        MAX,MIN,SQRT
*     ..
      MERR = 0
      IF (MODE.LE.0) RETURN
      IF (PO.GE.ZERO) THEN
          MERR = 2
          RETURN

      ELSE IF (RU.LE.RL) THEN
          MERR = 3
          RETURN

      END IF

      L1 = RL .LE. RO
      L2 = RI .LE. RL
      DO 10 NTYP = MTYP,1,-1
          IF (NTYP.EQ.1) THEN
*
*     BISECTION
*
              IF (MODE.EQ.1) THEN
                  R = TWO*RU
                  RETURN

              ELSE IF (RI-RL.LE.RU-RI) THEN
                  R = HALF* (RI+RU)
                  RETURN

              ELSE
                  R = HALF* (RL+RI)
                  RETURN

              END IF

          ELSE IF (NTYP.EQ.MTYP .AND. L1) THEN
              IF (.NOT.L2) AI = (FI-FO)/ (RI*PO)
              AU = (FU-FO)/ (RU*PO)
          END IF

          IF (L1 .AND. (NTYP.EQ.2.OR.L2)) THEN
*
*     TWO POINT QUADRATIC EXTRAPOLATION OR INTERPOLATION
*
              IF (AU.GE.ONE) GO TO 10
              R = HALF*RU/ (ONE-AU)

          ELSE IF (.NOT.L1 .OR. .NOT.L2 .AND. NTYP.EQ.3) THEN
*
*     THREE POINT QUADRATIC EXTRAPOLATION OR INTERPOLATION
*
              AL = (FI-FL)/ (RI-RL)
              AU = (FU-FI)/ (RU-RI)
              DEN = AU - AL
              IF (DEN.LE.ZERO) GO TO 10
              R = RI - HALF* (AU* (RI-RL)+AL* (RU-RI))/DEN

          ELSE IF (L1 .AND. .NOT.L2 .AND. NTYP.EQ.4) THEN
*
*     THREE POINT CUBIC EXTRAPOLATION OR INTERPOLATION
*
              DIS = (AI-ONE)* (RU/RI)
              DEN = (AU-ONE)* (RI/RU) - DIS
              DIS = AU + AI - DEN - TWO* (ONE+DIS)
              DIS = DEN*DEN - THREE*DIS
              IF (DIS.LT.ZERO) GO TO 10
              DEN = DEN + SQRT(DIS)
              IF (DEN.EQ.ZERO) GO TO 10
              R = (RU-RI)/DEN

          ELSE
              GO TO 10

          END IF

          IF (MODE.EQ.1 .AND. R.GT.RU) THEN
*
*     EXTRAPOLATION ACCEPTED
*
              R = MAX(R,C1L*RU)
              R = MIN(R,C1U*RU)
              RETURN

          ELSE IF (MODE.EQ.2 .AND. R.GT.RL .AND. R.LT.RU) THEN
*
*     INTERPOLATION ACCEPTED
*
              IF (RI.EQ.ZERO .AND. NTYP.NE.4) THEN
                  R = MAX(R,RL+C2L* (RU-RL))

              ELSE
                  R = MAX(R,RL+C3L* (RU-RL))
              END IF

              R = MIN(R,RL+C2U* (RU-RL))
              IF (R.EQ.RI) GO TO 10
              RETURN

          END IF

   10 CONTINUE
      END
* SUBROUTINE PNSTEP                ALL SYSTEMS                89/12/01
* PORTABILITY : ALL SYSTEMS
* 89/01/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DETERMINATION OF A SCALING FACTOR FOR THE BOUNDARY STEP.
*
* PARAMETERS :
*  RI  DEL  MAXIMUM STEPSIZE.
*  RI  A  INPUT PARAMETER.
*  RI  B  INPUT PARAMETER.
*  RI  C  INPUT PARAMETER.
*  RO  ALF  SCALING FACTOR FOR THE BOUNDARY STEP SUCH THAT
*         A**2+2*B*ALF+C*ALF**2=DEL**2.
*
      SUBROUTINE PNSTEP(DEL,A,B,C,ALF)
      REAL*8  DEL, A, B, C, ALF
      REAL*8  DEN, DIS
      REAL*8  ZERO
      PARAMETER (ZERO = 0.0D 0)
      ALF = ZERO
      DEN = (DEL+A) * (DEL-A)
      IF (DEN .LE. ZERO) RETURN
      DIS = B*B + C*DEN
      IF (B .GE. ZERO) THEN
      ALF = DEN / (SQRT(DIS) + B)
      ELSE
      ALF = (SQRT(DIS) - B) / C
      ENDIF
      RETURN
      END
* SUBROUTINE PP0AF8             ALL SYSTEMS                 97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF VALUE OF THE AUGMENTED LAGRANGIAN FUNCTION.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  N  DIMENSION OF THE CONSTRAINT NULL SPACE.
*  II  NC  NUMBER OF CONSTRAINTS.
*  RI  CF(NC+1)  VECTOR CONTAINING VALUES OF THE CONSTRAINTS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CZ(NC)  VECTOR OF LAGRANGE MULTIPLIERS.
*  RI  RPF  PENALTY COEFFICIENT.
*  RO  FC  VALUE OF THE PENALTY TERM.
*  RO  F  VALUE OF THE PENALTY FUNCTION.
*
      SUBROUTINE PP0AF8(NF,N,NC,CF,IC,ICA,CL,CU,CZ,RPF,FC,F)
      INTEGER NF,N,NC,IC(*),ICA(*)
      REAL*8 CF(*),CL(*),CU(*),CZ(*),RPF,FC,F
      REAL*8 POM,TEMP
      INTEGER J,KC
      FC=0.0D0
      DO 1 KC=1,NC
      IF (IC(KC).GT.0) THEN
      POM=0.0D0
      TEMP=CF(KC)
      IF (IC(KC).EQ.1.OR.IC(KC).GE.3) POM=MIN(POM,TEMP-CL(KC))
      IF (IC(KC).EQ.2.OR.IC(KC).GE.3) POM=MIN(POM,CU(KC)-TEMP)
      FC=FC+RPF*ABS(POM)
      ENDIF
    1 CONTINUE
      DO 2 J=1,NF-N
      KC=ICA(J)
      IF (KC.GT.0) THEN
      POM=0.0D0
      TEMP=CF(KC)
      IF (IC(KC).EQ.1.OR.IC(KC).EQ.3.OR.IC(KC).EQ.5)
     & POM=MIN(POM,TEMP-CL(KC))
      IF (IC(KC).EQ.2.OR.IC(KC).EQ.4.OR.IC(KC).EQ.6)
     & POM=MAX(POM,TEMP-CU(KC))
      FC=FC-CZ(J)*POM
      ENDIF
    2 CONTINUE
      F=CF(NC+1)+FC
      RETURN
      END
* SUBROUTINE PPSET2             ALL SYSTEMS                   97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* COMPUTATION OF THE NEW PENALTY PARAMETERS.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CZ(NF)  VECTOR OF LAGRANGE MULTIPLIERS.
*  RI  CP(NC)  VECTOR CONTAINING PENALTY PARAMETERS.
*
      SUBROUTINE PPSET2(NF,N,NC,ICA,CZ,CP)
      INTEGER NF,N,NC,ICA(*)
      REAL*8 CZ(*),CP(*)
      REAL*8 TEMP
      INTEGER J,L,KC
      DO 1 KC=1,NC
      CP(KC)=0.5D 0*CP(KC)
    1 CONTINUE
      DO 2 J=1,NF-N
      L=ICA(J)
      IF (L.GT.0) THEN
      TEMP=ABS(CZ(J))
      CP(L)=MAX(TEMP,CP(L)+0.5D 0*TEMP)
      ENDIF
    2 CONTINUE
      RETURN
      END
* SUBROUTINE PS0G01                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SIMPLE SEARCH WITH TRUST REGION UPDATE.
*
* PARAMETERS :
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  INITIAL VALUE OF THE OBJECTIVE FUNCTION.
*  RI  PO  INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  PP  QUADRATIC PART OF THE PREDICTED FUNCTION VALUE.
*  RU  XDEL  TRUST REGION BOUND.
*  RO  XDELO  PREVIOUS TRUST REGION BOUND.
*  RI  XMAX MAXIMUM STEPSIZE.
*  RI  RMAX  MAXIMUM VALUE OF THE STEPSIZE PARAMETER.
*  RI  SNORM  EUCLIDEAN NORM OF THE DIRECTION VECTOR.
*  RI  BET1  LOWER BOUND FOR STEPSIZE REDUCTION.
*  RI  BET2  UPPER BOUND FOR STEPSIZE REDUCTION.
*  RI  GAM1  LOWER BOUND FOR STEPSIZE EXPANSION.
*  RI  GAM2  UPPER BOUND FOR STEPSIZE EXPANSION.
*  RI  EPS4  FIRST TOLERANCE FOR RATIO DF/DFPRED. STEP BOUND IS
*         DECREASED IF DF/DFPRED<EPS4.
*  RI  EPS5  SECOND TOLERANCE FOR RATIO DF/DFPRED. STEP BOUND IS
*         INCREASED IF IT IS ACTIVE AND DF/DFPRED>EPS5.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IO  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  IU  IDIR INDICATOR FOR DIRECTION DETERMINATION.
*         IDIR=0-BASIC DETERMINATION. IDIR=1-DETERMINATION
*         AFTER STEPSIZE REDUCTION. IDIR=2-DETERMINATION AFTER
*         STEPSIZE EXPANSION.
*  IO  ITERS  TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-STEP
*         BOUND WAS DECREASED. ITERS=2-STEP BOUND WAS UNCHANGED.
*         ITERS=3-STEP BOUND WAS INCREASED. ITERS=6-FIRST STEPSIZE.
*  II  ITERD TERMINATION INDICATOR. ITERD<0-BAD DECOMPOSITION.
*         ITERD=0-DESCENT DIRECTION. ITERD=1-NEWTON LIKE STEP.
*         ITERD=2-INEXACT NEWTON LIKE STEP. ITERD=3-BOUNDARY STEP.
*         ITERD=4-DIRECTION WITH THE NEGATIVE CURVATURE.
*         ITERD=5-MARQUARDT STEP.
*  IO  MAXST  MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM
*         STEPSIZE WAS NOT OR WAS REACHED.
*  IO  NRED  ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  II  MRED  MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  II  KTERS  TERMINATION SELECTION. KTERS=1-NORMAL TERMINATION.
*         KTERS=6-FIRST STEPSIZE.
*  II  MES1  SWITCH FOR EXTRAPOLATION. MES1=1-CONSTANT INCREASING OF
*         THE INTERVAL. MES1=2-EXTRAPOLATION SPECIFIED BY THE PARAMETER
*         MES. MES1=3 SUPPRESSED EXTRAPOLATION.
*  II  MES2  SWITCH FOR TERMINATION. MES2=1-NORMAL TERMINATION.
*         MES2=2-TERMINATION AFTER AT LEAST TWO STEPS (ASYMPTOTICALLY
*         PERFECT LINE SEARCH).
*  II  MES3  SAFEGUARD AGAINST ROUNDING ERRORS. MES3=0-SAFEGUARD
*         SUPPRESSED. MES3=1-FIRST LEVEL OF SAFEGUARD. MES3=2-SECOND
*         LEVEL OF SAFEGUARD.
*  IU  ISYS  CONTROL PARAMETER.
*
* COMMON DATA :
*
* METHOD :
* G.A.SCHULTZ, R.B.SCHNABEL, R.H.BYRD: A FAMILY OF TRUST-REGION-BASED
* ALGORITHMS FOR UNCONSTRAINED MINIMIZATION WITH STRONG GLOBAL
* CONVERGENCE PROPERTIES, SIAM J. NUMER.ANAL. 22 (1985) PP. 47-67.
*
      SUBROUTINE PS0G01(R,F,FO,PO,PP,XDEL,XDELO,XMAX,RMAX,SNORM,
     & BET1,BET2,GAM1,GAM2,EPS4,EPS5,KD,LD,IDIR,ITERS,ITERD,
     & MAXST,NRED,MRED,KTERS,MES1,MES2,MES3,ISYS)
      INTEGER KD,LD,LDA,LDC,IDIR,ITERS,ITERD,MAXST,NRED,MRED,KTERS,
     & MES1,MES2,MES3,ISYS
      REAL*8 R,F,FO,PO,PP,XDEL,XDELO,XMAX,RMAX,SNORM,BET1,BET2,GAM1,
     & GAM2,EPS4,EPS5
      REAL*8 DF,DFPRED
      INTEGER NRED1,NRED2
      SAVE NRED1,NRED2
      GO TO (1,2) ISYS+1
    1 CONTINUE
      IF (IDIR.EQ.0) THEN
      NRED1=0
      NRED2=0
      ENDIF
      IDIR=0
      XDELO=XDEL
C
C     COMPUTATION OF THE NEW FUNCTION VALUE
C
      R=MIN(1.0D 0,RMAX)
      KD= 0
      LD=-1
      ISYS=1
      RETURN
    2 CONTINUE
      IF(KTERS.LT.0.OR.KTERS.GT.5) THEN
      ITERS=6
      ELSE
      DF=FO-F
      DFPRED=-R*(PO+R*PP)
      IF (DF.LT.EPS4*DFPRED) THEN
C
C     STEP IS TOO LARGE, IT HAS TO BE REDUCED
C
      IF (MES1.EQ.1) THEN
      XDEL=BET2*SNORM
      ELSE IF (MES1.EQ.2) THEN
      XDEL=BET2*MIN(0.5D 0*XDEL,SNORM)
      ELSE
      XDEL=0.5D 0*PO*SNORM/(PO+DF)
      XDEL=MAX(XDEL,BET1*SNORM)
      XDEL=MIN(XDEL,BET2*SNORM)
      ENDIF
      ITERS=1
      IF (MES3.LE.1) THEN
      NRED2=NRED2+1
      ELSE
      IF (ITERD.GT.2) NRED2=NRED2+1
      ENDIF
      ELSE IF (DF.LE.EPS5*DFPRED) THEN
C
C     STEP IS SUITABLE
C
      ITERS=2
      ELSE
C
C     STEP IS TOO SMALL, IT HAS TO BE ENLARGED
C
      IF (MES2.EQ.2) THEN
      XDEL=MAX(XDEL,GAM1*SNORM)
      ELSE IF (ITERD.GT.2) THEN
      XDEL=GAM1*XDEL
      ENDIF
      ITERS=3
      ENDIF
      XDEL=MIN(XDEL,XMAX,GAM2*SNORM)
      IF (FO.LE.F) THEN
      IF (NRED1.GE.MRED) THEN
      ITERS=-1
      ELSE
      IDIR=1
      ITERS=0
      NRED1=NRED1+1
      ENDIF
      ENDIF
      ENDIF
      MAXST=0
      IF (XDEL.GE.XMAX) MAXST=1
      IF (MES3.EQ.0) THEN
      NRED=NRED1
      ELSE
      NRED=NRED2
      ENDIF
      ISYS=0
      RETURN
      END
* SUBROUTINE PS0L02                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
*  EXTENDED LINE SEARCH WITHOUT DIRECTIONAL DERIVATIVES.
*
* PARAMETERS :
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RO  RO  INITIAL VALUE OF THE STEPSIZE PARAMETER.
*  RO  RP  PREVIOUS VALUE OF THE STEPSIZE PARAMETER.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  INITIAL VALUE OF THE OBJECTIVE FUNCTION.
*  RO  FP  PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
*  RI  PO  INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  PP  PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  FMIN  LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FMAX  UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
*  RI  RMIN  MINIMUM VALUE OF THE STEPSIZE PARAMETER
*  RI  RMAX  MAXIMUM VALUE OF THE STEPSIZE PARAMETER
*  RI  TOLS  TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE
*         CHANGE OF THE FUNCTION VALUE).
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IO  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  II  NIT  ACTUAL NUMBER OF ITERATIONS.
*  II  KIT  NUMBER OF THE ITERATION AFTER LAST RESTART.
*  IO  NRED  ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  II  MRED  MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  IO  MAXST  MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM
*         STEPSIZE WAS NOT OR WAS REACHED.
*  II  IEST  LOWER BOUND SPECIFICATION. IEST=0 OR IEST=1 IF LOWER BOUND
*         IS NOT OR IS GIVEN.
*  II  INITS  CHOICE OF THE INITIAL STEPSIZE. INITS=0-INITIAL STEPSIZE
*         IS SPECIFIED IN THE CALLING PROGRAM. INITS=1-UNIT INITIAL
*         STEPSIZE. INITS=2-COMBINED UNIT AND QUADRATICALLY ESTIMATED
*         INITIAL STEPSIZE. INITS=3-QUADRATICALLY ESTIMATED INITIAL
*         STEPSIZE.
*  IO  ITERS  TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT
*         LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY
*         STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE.
*         ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE.
*         ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION.
*         ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL
*         DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION.
*  II  KTERS  TERMINATION SELECTION. KTERS=1-PERFECT LINE SEARCH.
*         KTERS=2-GOLDSTEIN STEPSIZE. KTERS=3-CURRY STEPSIZE.
*         KTERS=4-EXTENDED CURRY STEPSIZE. KTERS=5-ARMIJO STEPSIZE.
*         KTERS=6-FIRST STEPSIZE.
*  II  MES  METHOD SELECTION. MES=1-BISECTION. MES=2-QUADRATIC
*         INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE).
*         MES=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL
*         DERIVATIVES). MES=4-CUBIC INTERPOLATION. MES=5-CONIC
*         INTERPOLATION.
*  IU  ISYS  CONTROL PARAMETER.
*
* SUBPROGRAM USED :
*  S   PNINT3  EXTRAPOLATION OR INTERPOLATION WITHOUT DIRECTIONAL
*         DERIVATIVES.
*
* METHOD :
* SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH EXTENDED TERMINATION
* CRITERIA.
*
      SUBROUTINE PS0L02(R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX,
     & TOLS,KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS,
     & MES,ISYS)
      INTEGER KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS,
     & MES,ISYS
      REAL*8 R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX,TOLS
      REAL*8 RL,FL,RU,FU,RI,FI,RTEMP,TOL
      INTEGER MTYP,MERR,MODE,INIT1,MES1,MES2
      LOGICAL L1,L2,L3,L4,L6,L7
      PARAMETER(TOL=1.0D-4)
      SAVE MTYP
      SAVE RL,FL,RU,FU,RI,FI
      GO TO (1,3) ISYS+1
    1 CONTINUE
      MES1=2
      MES2=2
      ITERS=0
      IF (PO.GE.0.0D 0) THEN
      R=0.0D 0
      ITERS=-2
      GO TO 4
      ENDIF
      IF(RMAX.LE.0.0D 0) THEN
      ITERS= 0
      GO TO 4
      ENDIF
C
C     INITIAL STEPSIZE SELECTION
C
      IF (INITS.GT.0) THEN
      RTEMP=FMIN-F
      ELSE IF (IEST.EQ.0) THEN
      RTEMP=F-FP
      ELSE
      RTEMP=MAX(F-FP,1.0D 1*(FMIN-F))
      ENDIF
      INIT1=ABS(INITS)
      RP=0.0D 0
      FP=FO
      PP=PO
      IF (INIT1.EQ.0) THEN
      ELSE IF (INIT1.EQ.1.OR.INITS.GE.1.AND.IEST.EQ.0) THEN
      R=1.0D 0
      ELSE IF (INIT1.EQ.2) THEN
      R=MIN(1.0D 0,4.0D 0*RTEMP/PO)
      ELSE IF (INIT1.EQ.3) THEN
      R=MIN(1.0D 0, 2.0D 0*RTEMP/PO)
      ELSE IF (INIT1.EQ.4) THEN
      R=2.0D 0*RTEMP/PO
      ENDIF
      RTEMP=R
      R=MAX(R,RMIN)
      R=MIN(R,RMAX)
      MODE=0
      RU=0.0D 0
      FU=FO
      RI=0.0D 0
      FI=FO
C
C     NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION)
C
    2 CALL PNINT3(RO,RL,RU,RI,FO,FL,FU,FI,PO,R,MODE,MTYP,MERR)
      IF (MERR.GT.0) THEN
      ITERS=-MERR
      GO TO 4
      ELSE IF (MODE.EQ.1) THEN
      NRED=NRED-1
      R=MIN(R,RMAX)
      ELSE IF (MODE.EQ.2) THEN
      NRED=NRED+1
      ENDIF
C
C     COMPUTATION OF THE NEW FUNCTION VALUE
C
      KD= 0
      LD=-1
      ISYS=1
      RETURN
    3 CONTINUE
      IF (ITERS.NE.0) GO TO 4
      IF (F.LE.FMIN) THEN
      ITERS=7
      GO TO 4
      ELSE
      L1=R.LE.RMIN.AND.NIT.NE.KIT
      L2=R.GE.RMAX
      L3=F-FO.LE.TOLS*R*PO.OR.F-FMIN.LE.(FO-FMIN)/1.0D 1
      L4=F-FO.GE.(1.0D 0-TOLS)*R*PO.OR.MES2.EQ.2.AND.MODE.EQ.2
      L6=RU-RL.LE.TOL*RU.AND.MODE.EQ.2
      L7=MES2.LE.2.OR.MODE.NE.0
      MAXST=0
      IF (L2) MAXST=1
      ENDIF
C
C     TEST ON TERMINATION
C
      IF (L1.AND..NOT.L3) THEN
      ITERS=0
      GO TO 4
      ELSE IF (L2.AND..NOT.F.GE.FU) THEN
      ITERS=7
      GO TO 4
      ELSE IF (L6) THEN
      ITERS=1
      GO TO 4
      ELSE IF (L3.AND.L7.AND.KTERS.EQ.5) THEN
      ITERS=5
      GO TO 4
      ELSE IF (L3.AND.L4.AND.L7.AND.(KTERS.EQ.2.OR.KTERS.EQ.3.OR.
     *         KTERS.EQ.4)) THEN
      ITERS=2
      GO TO 4
      ELSE IF (KTERS.LT.0.OR.KTERS.EQ.6.AND.L7) THEN
      ITERS=6
      GOTO 4
      ELSE IF(ABS(NRED).GE.MRED) THEN
      ITERS=-1
      GO TO 4
      ELSE
      RP=R
      FP=F
      MODE=MAX(MODE,1)
      MTYP=ABS(MES)
      IF(F.GE.FMAX) MTYP=1
      ENDIF
      IF (MODE.EQ.1) THEN
C
C     INTERVAL CHANGE AFTER EXTRAPOLATION
C
      RL=RI
      FL=FI
      RI=RU
      FI=FU
      RU=R
      FU=F
      IF (F.GE.FI) THEN
      NRED=0
      MODE=2
      ELSE IF ( MES1 .EQ. 1) THEN
      MTYP=1
      ENDIF
C
C     INTERVAL CHANGE AFTER INTERPOLATION
C
      ELSE IF (R.LE.RI) THEN
      IF (F.LE.FI) THEN
      RU=RI
      FU=FI
      RI=R
      FI=F
      ELSE
      RL=R
      FL=F
      ENDIF
      ELSE
      IF (F.LE.FI) THEN
      RL=RI
      FL=FI
      RI=R
      FI=F
      ELSE
      RU=R
      FU=F
      ENDIF
      ENDIF
      GO TO 2
    4 ISYS=0
      RETURN
      END
* SUBROUTINE PS1L01                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
*  STANDARD LINE SEARCH WITH DIRECTIONAL DERIVATIVES.
*
* PARAMETERS :
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RO  RO  INITIAL VALUE OF THE STEPSIZE PARAMETER.
*  RO  RP  PREVIOUS VALUE OF THE STEPSIZE PARAMETER.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  INITIAL VALUE OF THE OBJECTIVE FUNCTION.
*  RO  FP  PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
*  RO  P  VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  PO  INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  PP  PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  FMIN  LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FMAX  UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
*  RI  RMIN  MINIMUM VALUE OF THE STEPSIZE PARAMETER
*  RI  RMAX  MAXIMUM VALUE OF THE STEPSIZE PARAMETER
*  RI  TOLS  TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE
*         CHANGE OF THE FUNCTION VALUE).
*  RI  TOLP  TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE
*         CHANGE OF THE DIRECTIONAL DERIVATIVE).
*  RO  PAR1  PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC
*         UPDATES.
*  RO  PAR2  PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC
*         UPDATES.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IO  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  II  NIT  ACTUAL NUMBER OF ITERATIONS.
*  II  KIT  NUMBER OF THE ITERATION AFTER LAST RESTART.
*  IO  NRED  ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  II  MRED  MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  IO  MAXST  MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM
*         STEPSIZE WAS NOT OR WAS REACHED.
*  II  IEST  LOWER BOUND SPECIFICATION. IEST=0 OR IEST=1 IF LOWER BOUND
*         IS NOT OR IS GIVEN.
*  II  INITS  CHOICE OF THE INITIAL STEPSIZE. INITS=0-INITIAL STEPSIZE
*         IS SPECIFIED IN THE CALLING PROGRAM. INITS=1-UNIT INITIAL
*         STEPSIZE. INITS=2-COMBINED UNIT AND QUADRATICALLY ESTIMATED
*         INITIAL STEPSIZE. INITS=3-QUADRATICALLY ESTIMATED INITIAL
*         STEPSIZE.
*  IO  ITERS  TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT
*         LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY
*         STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE.
*         ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE.
*         ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION.
*         ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL
*         DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION.
*  II  KTERS  TERMINATION SELECTION. KTERS=1-PERFECT LINE SEARCH.
*         KTERS=2-GOLDSTEIN STEPSIZE. KTERS=3-CURRY STEPSIZE.
*         KTERS=4-EXTENDED CURRY STEPSIZE. KTERS=5-ARMIJO STEPSIZE.
*         KTERS=6-FIRST STEPSIZE.
*  II  MES  METHOD SELECTION. MES=1-BISECTION. MES=2-QUADRATIC
*         INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE).
*         MES=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL
*         DERIVATIVES). MES=4-CUBIC INTERPOLATION. MES=5-CONIC
*         INTERPOLATION.
*  IU  ISYS  CONTROL PARAMETER.
*
* SUBPROGRAM USED :
*  S   PNINT1  EXTRAPOLATION OR INTERPOLATION WITH DIRECTIONAL
*         DERIVATIVES.
*
* METHOD :
* SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH STANDARD TERMINATION
* CRITERIA.
*
      SUBROUTINE PS1L01(R,RO,RP,F,FO,FP,P,PO,PP,FMIN,FMAX,RMIN,RMAX,
     & TOLS,TOLP,PAR1,PAR2,KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,
     & ITERS,KTERS,MES,ISYS)
      INTEGER KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS,
     & MES,ISYS
      REAL*8 R,RO,RP,F,FO,FP,P,PO,PP,FMIN,FMAX,RMIN,RMAX,TOLS,TOLP,
     & PAR1,PAR2
      REAL*8 RL,FL,PL,RU,FU,PU,RTEMP
      INTEGER MTYP,MERR,MODE,INIT1,MES1,MES2,MES3
      LOGICAL L1,L2,L3,L5,L7,M1,M2,M3
      REAL*8 CON,CON1
      PARAMETER (CON=1.0D-2,CON1=1.0D-13)
      SAVE MTYP
      SAVE RL,FL,PL,RU,FU,PU
      GO TO (1,3) ISYS+1
    1 CONTINUE
      MES1=2
      MES2=2
      MES3=2
      ITERS=0
      IF (PO.GE.0.0D 0) THEN
      R=0.0D 0
      ITERS=-2
      GO TO 4
      ENDIF
      IF(RMAX.LE.0.0D 0) THEN
      ITERS=0
      GO TO 4
      ENDIF
C
C     INITIAL STEPSIZE SELECTION
C
      IF (INITS.GT.0) THEN
      RTEMP=FMIN-F
      ELSE IF (IEST.EQ.0) THEN
      RTEMP=F-FP
      ELSE
      RTEMP=MAX(F-FP,1.0D 1*(FMIN-F))
      ENDIF
      INIT1=ABS(INITS)
      RP=0.0D 0
      FP=FO
      PP=PO
      IF (INIT1.EQ.0) THEN
      ELSE IF (INIT1.EQ.1.OR.INITS.GE.1.AND.IEST.EQ.0) THEN
      R=1.0D 0
      ELSE IF (INIT1.EQ.2) THEN
      R=MIN(1.0D 0,4.0D 0*RTEMP/PO)
      ELSE IF (INIT1.EQ.3) THEN
      R=MIN(1.0D 0, 2.0D 0*RTEMP/PO)
      ELSE IF (INIT1.EQ.4) THEN
      R=2.0D 0*RTEMP/PO
      ENDIF
      R=MAX(R,RMIN)
      R=MIN(R,RMAX)
      MODE=0
      RU=0.0D 0
      FU=FO
      PU=PO
C
C     NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION)
C
    2 CALL PNINT1(RL,RU,FL,FU,PL,PU,R,MODE,MTYP,MERR)
      IF (MERR.GT.0) THEN
      ITERS=-MERR
      GO TO 4
      ELSE IF (MODE.EQ.1) THEN
      NRED=NRED-1
      R=MIN(R,RMAX)
      ELSE IF (MODE.EQ.2) THEN
      NRED=NRED+1
      ENDIF
C
C     COMPUTATION OF THE NEW FUNCTION VALUE AND THE NEW DIRECTIONAL
C     DERIVATIVE
C
      KD= 1
      LD=-1
      ISYS=1
      RETURN
    3 CONTINUE
      IF (MODE.EQ.0) THEN
      PAR1=P/PO
      PAR2=F-FO
      ENDIF
      IF (ITERS.NE.0) GO TO 4
      IF (F.LE.FMIN) THEN
      ITERS=7
      GO TO 4
      ELSE
      L1=R.LE.RMIN.AND.NIT.NE.KIT
      L2=R.GE.RMAX
      L3=F-FO.LE.TOLS*R*PO
      L5=P.GE.TOLP*PO.OR.MES2.EQ.2.AND.MODE.EQ.2
      L7=MES2.LE.2.OR.MODE.NE.0
      M1=.FALSE.
      M2=.FALSE.
      M3=L3
      IF(MES3.GE.1) THEN
      M1=ABS(P).LE.CON*ABS(PO).AND.FO-F.GE.(CON1/CON)*ABS(FO)
      L3=L3.OR.M1
      ENDIF
      IF(MES3.GE.2) THEN
      M2=ABS(P).LE.0.5D 0*ABS(PO).AND.ABS(FO-F).LE.2.0D 0*CON1*ABS(FO)
      L3=L3.OR.M2
      ENDIF
      MAXST=0
      IF (L2) MAXST=1
      ENDIF
C
C     TEST ON TERMINATION
C
      IF (L1.AND..NOT.L3) THEN
      ITERS=0
      GO TO 4
      ELSE IF (L2.AND.L3.AND..NOT.L5)  THEN
      ITERS=7
      GO TO 4
      ELSE IF (M3.AND.MES1.EQ.3) THEN
      ITERS=5
      GO TO 4
      ELSE IF (L3.AND.L5.AND.L7) THEN
      ITERS=4
      GO TO 4
      ELSE IF (KTERS.LT.0.OR.KTERS.EQ.6.AND.L7) THEN
      ITERS=6
      GOTO 4
      ELSE IF (ABS(NRED).GE.MRED) THEN
      ITERS=-1
      GO TO 4
      ELSE
      RP=R
      FP=F
      PP=P
      MODE=MAX(MODE,1)
      MTYP=ABS(MES)
      IF(F.GE.FMAX) MTYP=1
      ENDIF
      IF (MODE.EQ.1) THEN
C
C     INTERVAL CHANGE AFTER EXTRAPOLATION
C
      RL=RU
      FL=FU
      PL=PU
      RU=R
      FU=F
      PU=P
      IF (.NOT.L3) THEN
      NRED=0
      MODE=2
      ELSE IF ( MES1 .EQ. 1) THEN
      MTYP=1
      ENDIF
      ELSE
C
C     INTERVAL CHANGE AFTER INTERPOLATION
C
      IF (.NOT.L3) THEN
      RU=R
      FU=F
      PU=P
      ELSE
      RL=R
      FL=F
      PL=P
      ENDIF
      ENDIF
      GO TO 2
    4 ISYS=0
      RETURN
      END
* SUBROUTINE PUDBG1                ALL SYSTEMS                92/12/01
* PORTABILITY : ALL SYSTEMS
* 92/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VARIABLE METRIC UPDATE OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX
* USING THE FACTORIZATION B=L*D*TRANS(L).
*
* PARAMETERS :
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  RU  H(M)  FACTORIZATION B=L*D*TRANS(L) OF A POSITIVE
*         DEFINITE APPROXIMATION OF THE HESSIAN MATRIX.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RI  GO(NF)  GRADIENTS DIFFERENCE.
*  RI  R  VALUE OF THE STEPSIZE PARAMETER.
*  RI  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
*  II  NIT  ACTUAL NUMBER OF ITERATIONS.
*  II  KIT  NUMBER OF THE ITERATION AFTER LAST RESTART.
*  IO  ITERH  TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION.
*         ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS.
*  II  MET  SELECTION OF SELF SCALING. MET=1-SELF SCALING SUPPRESSED.
*         MET=2 INITIAL SELF SCALING.
*  II  MOD  CORRECTION IF THE NEGATIVE CURVATURE OCCURS.
*         MOD=1-CORRECTION SUPPRESSED. MOD=2-POWELL'S CORRECTION.
*
* SUBPROGRAMS USED :
*  S   MXDPGU  CORRECTION OF A DENSE SYMMETRIC POSITIVE DEFINITE
*         MATRIX IN THE FACTORED FORM B=L*D*TRANS(L).
*  S   MXDPGS  SCALING OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX
*         IN THE FACTORED FORM B=L*D*TRANS(L).
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  RF  MXVDOT  DOT PRODUCT OF VECTORS.
*  S   MXVSCL  SCALING OF A VECTOR.
*
* METHOD :
* BFGS VARIABLE METRIC METHOD.
*
      SUBROUTINE PUDBG1(N,H,G,S,XO,GO,R,PO,NIT,KIT,ITERH,MET,MOD)
*     .. Scalar Arguments ..
      DOUBLE PRECISION PO,R
      INTEGER          ITERH,KIT,MET,MOD,N,NIT
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION G(*),GO(*),H(*),S(*),XO(*)
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION A,B,C,GAM,PAR,DIS
      LOGICAL          L1,L3
*     ..
*     .. External Functions ..
      DOUBLE PRECISION MXVDOT,MXDPGP
      EXTERNAL         MXVDOT,MXDPGP
*     ..
*     .. External Subroutines ..
      EXTERNAL         MXDPGS,MXDPGU,MXVDIF,MXVSCL
*     ..
      L1 = MET .GE. 3 .OR. MET .EQ. 2 .AND. NIT .EQ. KIT
      L3 = .NOT. L1
*
*     DETERMINATION OF THE PARAMETERS B, C
*
      B = MXVDOT(N,XO,GO)
      A = 0.0D0
      IF (L1) THEN
          CALL MXVCOP(N,GO,S)
          CALL MXDPGB(N,H,S,1)
          A=MXDPGP(N,H,S,S)
          IF (A.LE.0.0D0) THEN
              ITERH=1
              RETURN

          END IF

      END IF

      CALL MXVDIF(N,GO,G,S)
      CALL MXVSCL(N,R,S,S)
      C = -R*PO
      IF (C.LE.0.0D0) THEN
          ITERH = 3
          RETURN

      END IF

      IF (MOD.GT.1) THEN
          IF (B.LE.1.0D-4*C) THEN
*
*     POWELL'S CORRECTION
*
              DIS=(1.0D0-0.1D0)*C/(C-B)
              CALL MXVDIF(N,GO,S,GO)
              CALL MXVDIR(N,DIS,GO,S,GO)
              B=C+DIS*(B-C)
              IF (L1) A=C+2.0D0*(1.0D0-DIS)*(B-C)+DIS*DIS*(A-C)
          ENDIF
      ELSE
          IF (B.LE.1.0D-4*C) THEN
              ITERH = 2
              RETURN

          ENDIF

      ENDIF

      IF (L1) THEN
*
*     DETERMINATION OF THE PARAMETER GAM (SELF SCALING)
*
          PAR = C/B
          GAM = PAR
          IF (MET.GT.1) THEN
              IF (NIT.NE.KIT) THEN
                  L3=GAM.LT.0.5D0.OR.GAM.GT.4.0D0
              ENDIF

          ENDIF

      ENDIF
      IF (L3) THEN
          GAM = 1.0D0
          PAR = GAM

      END IF
*
*     BFGS UPDATE
*
      CALL MXDPGU(N,H,PAR/B,GO,XO)
      CALL MXDPGU(N,H,-1.0D0/C,S,XO)
      ITERH = 0
      IF (GAM.EQ.1.0D0) RETURN
      CALL MXDPGS(N,H,1.0D0/GAM)
      RETURN

      END
* SUBROUTINE PUDBI1                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VARIABLE METRIC UPDATES.
*
* PARAMETERS :
*  II  N  NUMBER OF VARIABLES.
*  RU  H(N*(N+1)/2)  UPDATED APPROXIMATION OF THE INVERSE HESSIAN
*         MATRIX.
*  RA  S(N)  AUXILIARY VECTOR.
*  RI  XO(N)  VECTOR OF VARIABLES DIFFERENCE.
*  RI  GO(N)  GRADIENT DIFFERENCE.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RI  PO  INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  PAR1  PARAMETER FOR CONTROL SCALING.
*  RO  PAR2  PARAMETER FOR CONTROL SCALING.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  INITIAL VALUE OF THE OBJECTIVE FUNCTION.
*  RI  P  CURRENT VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  ETA9  THE MAXIMUM REAL NUMBER.
*  II  NIT  NUMBER OF ITERATIONS.
*  II  KIT  INDEX OF THE ITERATION WITH THE LAST RESTART.
*  II  MET  VARIABLE METRIC UPDATE.
*  II  MET1  SCALING STRATEGY.
*  II  MET2  CORRECTION RULE.
*  IU  IDECF  DECOMPOSITION INDICATOR.
*  II  ITERD  TERMINATION INDICATOR. ITERD<0-BAD DECOMPOSITION.
*         ITERD=0-DESCENT DIRECTION. ITERD=1-NEWTON LIKE STEP.
*         ITERD=2-INEXACT NEWTON LIKE STEP. ITERD=3-BOUNDARY STEP.
*         ITERD=4-DIRECTION WITH THE NEGATIVE CURVATURE.
*         ITERD=5-MARQUARDT STEP.
*  IO  ITERH  UPDATE INDICATOR. ITERH=0-SUCCESSFUL UPDATE.
*         ITERH>0-UNSUCCESSFUL UPDATE.
*
* METHOD :
* VARIOUS VARIABLE METRIC UPDATES INCLUDING BFGS UPDATE.
*
      SUBROUTINE PUDBI1(N,H,S,XO,GO,R,PO,PAR1,PAR2,F,FO,P,ETA9,NIT,KIT,
     & MET,MET1,MET2,IDECF,ITERD,ITERH)
      INTEGER N,NIT,KIT,MET,MET1,MET2,IDECF,ITERD,ITERH
      REAL*8 H(*),S(*),XO(*),GO(*),R,PO,ETA9
      REAL*8 PAR1,PAR2
      REAL*8 F,FO,P
      REAL*8 AA,CC
      REAL*8 MXVDOT
      REAL*8 DIS,POM,POM3,POM4,A,B,C,GAM,RHO,PAR
      REAL*8 DEN
      LOGICAL L1,L2,L3
      IF (MET.LE.0) GO TO 22
      IF (IDECF.NE.9) THEN
      ITERH=-1
      GO TO 22
      ENDIF
      L1=ABS(MET1).GE.3.OR.ABS(MET1).EQ.2.AND.NIT.EQ.KIT
      L3=.NOT.L1
*
*     DETERMINATION OF THE PARAMETERS A, B, C
*
      B=MXVDOT(N,XO,GO)
      IF (B.LE.0.0D 0) THEN
      ITERH=2
      GO TO 22
      ENDIF
      CALL MXDSMM(N,H,GO,S)
      A=MXVDOT(N,GO,S)
      IF (A.LE.0.0D 0) THEN
      ITERH=1
      GO TO 22
      ENDIF
      IF(MET.NE.1.OR.L1) THEN
      IF (ITERD.NE.1) THEN
      C=0.0D 0
      ELSE
      C=-R*PO
      IF (C.LE.0.0D 0) THEN
      ITERH=3
      GO TO 22
      ENDIF
      ENDIF
      ELSE
      C=0.0D 0
      ENDIF
*
*     DETERMINATION OF THE PARAMETER RHO (NONQUADRATIC PROPERTIES)
*
      IF (MET2.EQ.1) THEN
      RHO=1.0D 0
      ELSE IF (FO-F+P.EQ.0) THEN
      RHO=1.0D 0
      ELSE
      RHO=0.5D 0*B/(FO-F+P)
      ENDIF
      IF(RHO.LE.1.0D-2) RHO=1.0D 0
      IF(RHO.GE.1.0D 2) RHO=1.0D 0
      AA=A/B
      CC=C/B
      IF (L1) THEN
*
*     DETERMINATION OF THE PARAMETER GAM (SELF SCALING)
*
      PAR=A/B
      POM3=0.7D 0
      POM4=6.0D 0
      GAM=RHO/PAR
      IF (NIT.NE.KIT) THEN
      IF (MET1.EQ.3) THEN
      L2=PAR2.LE.0.0D 0
      L3=L2.AND.ABS(PAR1).LE.0.2D 0
      L3=L3.OR.(.NOT.L2.AND.GAM.GT.1.0D 0)
      L3=L3.OR.(L2.AND.PAR1.LT.0.0D 0.AND.GAM.GT.1.0D 0)
      L3=L3.OR.(L2.AND.PAR1.GT.0.0D 0.AND.GAM.LT.1.0D 0)
      L3=L3.OR.GAM.LT.POM3
      L3=L3.OR.GAM.GT.POM4
      ELSE IF (MET1.EQ.4) THEN
      L3=GAM.LT.POM3.OR.GAM.GT.POM4
      ENDIF
      ENDIF
      ENDIF
      IF (L3) THEN
      GAM=1.0D 0
      PAR=RHO/GAM
      ENDIF
      IF (MET.EQ.1) GO TO 17
*
*     NEW UPDATE
*
      POM=1.0D 0/(AA*CC)
      IF (POM.LT.1.0D 0) THEN
      POM=MAX(1.0D-15,(SQRT(C/A)-POM)/(1.0D 0-POM))
      GO TO 20
      ENDIF
   17 CONTINUE
*
*     BFGS UPDATE
*
      POM=1.0D 0
      DIS=PAR+AA
      CALL MXVDIR(N,-DIS,XO,S,XO)
      DIS=1.0D 0/(B*DIS)
      CALL MXDSMU(N,H,DIS,XO)
      CALL MXDSMU(N,H,-DIS,S)
      GO TO 21
   20 CONTINUE
*
*     GENERAL UPDATE
*
      DEN=PAR+POM*AA
      DIS=POM/DEN
      CALL MXDSMU(N,H,(PAR*DIS-1.0D 0)/A,S)
      CALL MXVDIR(N,-DIS,S,XO,S)
      CALL MXDSMU(N,H,DEN/B,S)
   21 CONTINUE
      IF (GAM.EQ.1.0D 0) GO TO 22
*
*     SCALING
*
      CALL MXDSMS(N,H,GAM)
   22 CONTINUE
      ITERH=0
      RETURN
      END
* SUBROUTINE PUDBQ1                ALL SYSTEMS                97/12/01
C PORTABILITY : ALL SYSTEMS
C 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* BROYDEN GOOD UPDATE OF A RECTANGULAR MATRIX AFTER THE QR
* DECOMPOSITION.
*
* PARAMETERS :
*  II  N  NUMBER OF VARIABLES.
*  II  NA  NUMBER OF EQUATIONS.
*  RU  H(N*(N+1)/2)  UPDATED UPPER TRIANGULAR MATRIX.
*  RI  ETA2  PARAMETER WHICH CONTROLS A NONSINGULARITY
*  RU  AG(N*NA)  UPDATED RECTANGULAR MATRIX.
*  RI  AF(NA)  RIGHT HAND SIDES OF NONLINEAR EQUATIONS.
*  RA  S(N)  AUXILIARY VECTOR.
*  RI  XO(N)  VECTOR OF VARIABLES DIFFERENCE.
*  RI  AFO(NA)  RIGHT HAND SIDES DIFFERENCE.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  II  MET  VARIABLE METRIC UPDATE.
*  IO  ITERH  UPDATE INDICATOR. ITERH=0-SUCCESSFUL UPDATE.
*         ITERH>0-UNSUCCESSFUL UPDATE.
*  IU  IDECA  DECOMPOSITION INDICATOR.
*  II  NDECA  NUMBER OF DECOMPOSITIONS.
*
* METHOD :
* VARIOUS VARIABLE METRIC UPDATES INCLUDING BFGS UPDATE.
*
      SUBROUTINE PUDBQ1(N,NA,H,ETA2,AG,AF,S,XO,AFO,R,MET,ITERH,IDECA,
     & NDECA)
      INTEGER N,NA,MET,INF,ITERH,IDECA,NDECA
      REAL*8 H(*),ETA2,AG(*),AF(*),S(*),XO(*),AFO(*),R
      REAL*8 DEN,MXVDOT
      IF (MET.LE.0) RETURN
      IF (IDECA.EQ.0) THEN
*
*     QR DECOMPOSITION
*
      DEN=ETA2
      CALL MXDRQF(N,NA,AG,H)
      CALL MXDPRC(N,H,INF,DEN)
      NDECA=NDECA+1
      IDECA=1
      ELSE IF (IDECA.NE.1) THEN
      ITERH=-1
      GOTO 22
      ENDIF
*
*     THE GOOD BROYDEN UPDATE
*
      DEN=MXVDOT(N,XO,XO)
      IF (DEN.LE.0.0D 0)  THEN
      ITERH=2
      GOTO 22
      ENDIF
      CALL MXVCOP(N,XO,S)
   21 CONTINUE
      CALL MXVNEG(N,XO,XO)
      CALL MXDPRM(N,H,XO,1)
      CALL MXDRMD(N,NA,AG,XO,1.0D 0,AFO,AFO)
      CALL MXDRQU(N,NA,AG,H,1.0D 0/DEN,AFO,S,XO,INF)
      ITERH=0
   22 CONTINUE
      RETURN
      END
* SUBROUTINE PUDRV1                ALL SYSTEMS                97/12/01
* PORTABILITY : ALL SYSTEMS
* 97/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DRIVER FOR HYBRID QUASI-NEWTON UPDATES.
*
* PARAMETERS:
*  RI  R  VALUE OF THE STEPSIZE PARAMETER.
*  RI  FO  PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
*  RI  F  CURRENT VALUE OF THE OBJECTIVE FUNCTION.
*  RI  PO  PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE.
*  II  IPOM1  UPDATE SELECTION.
*  II  IPOM2  METHOD SELECTION.
*  IO  NRED  ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
*  II  IREST  RESTART SPECIFICATION. IF IREST=0 DOES NOT HOLD THEN A
*         RESTART IS PERFORMED.
*
      SUBROUTINE PUDRV1(R,FO,F,PO,IPOM1,IPOM2,NRED,IREST)
      INTEGER IPOM1,IPOM2,NRED,IREST
      REAL*8 R,FO,F,PO
      REAL*8 POM
      REAL*8 CON1,CON2
      PARAMETER(CON1=1.0D-1,CON2=1.0D-2)
      POM=(FO-F)/FO
      GO TO (1,2,3,4) IPOM2
    1 IREST=1
      IF (NRED.LE.0) THEN
      IPOM1=2
      IREST=0
      ELSE
      IPOM1=0
      ENDIF
      GO TO 5
    2 IREST=1
      IF(POM.GE.CON2) THEN
      IPOM1=0
      ELSE IF (F-FO.LE.R*PO) THEN
      IPOM1=0
      ELSE
      IPOM1=1
      IREST=0
      ENDIF
      GO TO 5
    3 IREST=1
      IF (NRED.LE.0) THEN
      IF (IPOM1.NE.1) THEN
      IPOM1=2
      IREST=0
      ELSE
      IPOM1=0
      ENDIF
      ELSE IF (POM.GE.CON2) THEN
      IPOM1=0
      ELSE IF (IPOM1.NE.2) THEN
      IPOM1=1
      IREST=0
      ELSE
      IPOM1=0
      ENDIF
      GO TO 5
    4 IREST=1
      IPOM1=0
    5 CONTINUE
      RETURN
      END
* SUBROUTINE PYADB4             ALL SYSTEMS                   98/12/01
* PORTABILITY : ALL SYSTEMS
* 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* NEW LINEAR CONSTRAINTS OR NEW SIMPLE BOUNDS ARE ADDED TO THE ACTIVE
* SET. GILL-MURRAY FACTORIZATION OF THE TRANSFORMED HESSIAN MATRIX
* APPROXIMATION IS UPDATED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  IU  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC  NUMBER OF LINEARIZED CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  IU  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  CF(NC)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCTIONS.
*  RI  CFD(NC) VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT FUNCTIONS.
*  IU  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  IU  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RU  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RU  H(NF*(NF+1)/2)  GILL-MURRAY FACTORIZATION OF THE TRANSFORMED
*         HESSIAN MATRIX APPROXIMATION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  RI  R  VALUE OF THE STEPSIZE PARAMETER.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RI  EPS9  TOLERANCE FOR ACTIVE CONSTRAINTS.
*  RO  GMAX  MAXIMUM ABSOLUTE VALUE OF A PARTIAL DERIVATIVE.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF A NEGATIVE LAGRANGE MULTIPLIER.
*  II  KBF  TYPE OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. KBF=1-ONE
*         SIDED SIMPLE BOUNDS. KBF=2-TWO SIDED SIMPLE BOUNDS.
*  II  KBC  TYPE OF CONSTRAINTS. KBC=0-NO CONSTRAINTS. KBC=1-CONSTRAINTS
*         WITH ONE SIDED BOUNDS. KBC=2-CONSTRAINTS WITH TWO SIDED
*         BOUNDS.
*  IU  INEW  INDEX OF THE NEW ACTIVE CONSTRAINT.
*  IO  IER  ERROR INDICATOR.
*  IO  ITERM  TERMINATION INDICATOR.
*
* COMMON DATA :
*  IU  NADD  NUMBER OF CONSTRAINT ADDITIONS.
*
* SUBPROGRAMS USED :
*  S   PLADB4  ADDITION OF A NEW ACTIVE CONSTRAINT.
*  S   PLNEWS  IDENTIFICATION OF ACTIVE UPPER BOUNDS.
*  S   PLNEWL  IDENTIFICATION OF ACTIVE LINEAR CONSTRAINRS.
*  S   PLDIRL  NEW VALUES OF CONSTRAINT FUNCTIONS.
*  S   MXVIND  CHANGE OF THE INTEGER VECTOR FOR CONSTRAINT ADDITION.
*
      SUBROUTINE PYADB4(NF,N,NC,X,IX,XL,XU,CF,CFD,IC,ICA,CL,CU,CG,
     & CR,CZ,H,S,R,EPS7,EPS9,GMAX,UMAX,KBF,KBC,INEW,IER,ITERM)
      INTEGER NF,N,NC,IX(*),IC(*),ICA(*),KBF,KBC,INEW,IER,ITERM
      DOUBLE PRECISION X(*),XL(*),XU(*),CF(*),CFD(*),CL(*),CU(*),
     & CG(*),CR(*),CZ(*),H(*),S(*),R,EPS7,EPS9,GMAX,UMAX
      INTEGER I,J,K,L,IJ,IK,KC,KJ,KK,LL
      DOUBLE PRECISION DEN,TEMP
      INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      COMMON /STATSQP/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      IF (KBC.GT.0) THEN
      IF (R.NE.0.0D 0) CALL PLDIRL(NC,CF,CFD,IC,R,KBC)
      IF (INEW.NE.0) THEN
      IF (KBF.GT.0) THEN
      DO 1 I=1,NF
      INEW=0
      CALL PLNEWS(X,IX,XL,XU,EPS9,I,INEW)
      CALL PLADB4(NF,N,ICA,CG,CR,CZ,H,S,EPS7,GMAX,UMAX,9,INEW,
     & NADD,IER)
      CALL MXVIND(IX,I,IER)
      IF (IER.LT.0) THEN
      ITERM=-15
      RETURN
      ENDIF
    1 CONTINUE
      ENDIF
      DO 2 KC=1,NC
      INEW=0
      CALL PLNEWL(KC,CF,IC,CL,CU,EPS9,INEW)
      CALL PLADB4(NF,N,ICA,CG,CR,CZ,H,S,EPS7,GMAX,UMAX,9,INEW,
     & NADD,IER)
      CALL MXVIND(IC,KC,IER)
      IF (IER.LT.0) THEN
      ITERM=-15
      RETURN
      ENDIF
    2 CONTINUE
      ENDIF
      ELSE IF (KBF.GT.0) THEN
      K=0
      DO 7 L=1,NF
      IF (IX(L).GE.0) K=K+1
      INEW=0
      CALL PLNEWS(X,IX,XL,XU,EPS9,L,INEW)
      IF (INEW.NE.0) THEN
      IX(L)=10-IX(L)
      KK=K*(K-1)/2
      DEN=H(KK+K)
      IF (DEN.NE.0.0D 0) THEN
      IJ=0
      KJ=KK
      DO 4 J=1,N
      IF (J.LE.K) THEN
      KJ=KJ+1
      ELSE
      KJ=KJ+J-1
      ENDIF
      IF (J.NE.K) TEMP=H(KJ)/DEN
      IK=KK
      DO 3 I=1,J
      IF (I.LE.K) THEN
      IK=IK+1
      ELSE
      IK=IK+I-1
      ENDIF
      IJ=IJ+1
      IF (I.NE.K.AND.J.NE.K) H(IJ)=H(IJ)+TEMP*H(IK)
    3 CONTINUE
    4 CONTINUE
      ENDIF
      LL=KK+K
      DO 6 I=K+1,N
      DO 5 J=1,I
      LL=LL+1
      IF (J.NE.K) THEN
      KK=KK+1
      H(KK)=H(LL)
      ENDIF
    5 CONTINUE
    6 CONTINUE
      N=N-1
      ENDIF
    7 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PYFUT1                ALL SYSTEMS                98/12/01
C PORTABILITY : ALL SYSTEMS
C 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* TERMINATION CRITERIA AND TEST ON RESTART.
*
* PARAMETERS :
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  RI  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
*  RI  UMAX  MAXIMUN ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
*  RO  GMAX  NORM OF THE TRANSFORMED GRADIENT.
*  RI  DMAX  MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
*  RI  TOLX  LOWER BOUND FOR STEPLENGTH.
*  RI  TOLF  LOWER BOUND FOR FUNCTION DECREASE.
*  RI  TOLB  LOWER BOUND FOR FUNCTION VALUE.
*  RI  TOLG  LOWER BOUND FOR GRADIENT.
*  II  KD  DEGREE OF REQUIRED DERIVATIVES.
*  IU  NIT  ACTUAL NUMBER OF ITERATIONS.
*  II  KIT  NUMBER OF THE ITERATION AFTER RESTART.
*  II  MIT  MAXIMUM NUMBER OF ITERATIONS.
*  IU  NFV  ACTUAL NUMBER OF COMPUTED FUNCTION VALUES.
*  II  MFV  MAXIMUM NUMBER OF COMPUTED FUNCTION VALUES.
*  IU  NFG  ACTUAL NUMBER OF COMPUTED GRADIENT VALUES.
*  II  MFG  MAXIMUM NUMBER OF COMPUTED GRADIENT VALUES.
*  IU  NTESX  ACTUAL NUMBER OF TESTS ON STEPLENGTH.
*  II  MTESX  MAXIMUM NUMBER OF TESTS ON STEPLENGTH.
*  IU  NTESF  ACTUAL NUMBER OF TESTS ON FUNCTION DECREASE.
*  II  MTESF  MAXIMUM NUMBER OF TESTS ON FUNCTION DECREASE.
*  II  ITES  SYSTEM VARIBLE WHICH SPECIFIES TERMINATION. IF ITES=0
*         THEN TERMINATION IS SUPPRESSED.
*  II  IRES1  RESTART SPECIFICATION. RESTART IS PERFORMED AFTER
*         IRES1*N+IRES2 ITERATIONS.
*  II  IRES2  RESTART SPECIFICATION. RESTART IS PERFORMED AFTER
*         IRES1*N+IRES2 ITERATIONS.
*  IU  IREST  RESTART INDICATOR. RESTART IS PERFORMED IF IREST>0.
*  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
*         ITERS=0 FOR ZERO STEP.
*  IO  ITERM  TERMINATION INDICATOR. ITERM=1-TERMINATION AFTER MTESX
*         UNSUFFICIENT STEPLENGTHS. ITERM=2-TERMINATION AFTER MTESF
*         UNSUFFICIENT FUNCTION DECREASES. ITERM=3-TERMINATION ON LOWER
*         BOUND FOR FUNCTION VALUE. ITERM=4-TERMINATION ON LOWER BOUND
*         FOR GRADIENT. ITERM=11-TERMINATION AFTER MAXIMUM NUMBER OF
*         ITERATIONS. ITERM=12-TERMINATION AFTER MAXIMUM NUMBER OF
*         COMPUTED FUNCTION VALUES.
*
      SUBROUTINE PYFUT1(N,F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLB,TOLG,KD,
     & NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF,ITES,IRES1,
     & IRES2,IREST,ITERS,ITERM)
      INTEGER N,KD,NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF,
     & ITES,IRES1,IRES2,IREST,ITERS,ITERM
      REAL*8  F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLG,TOLB
      REAL*8  TEMP
      IF (ITERM.LT.0) RETURN
      IF (ITES .LE.0) GOTO 2
      IF (ITERS.EQ.0) GOTO 1
      IF (NIT.LE.0) FO=F+MIN(SQRT(ABS(F)),ABS(F)/1.0D 1)
      IF (F.LE.TOLB) THEN
      ITERM = 3
      RETURN
      ENDIF
      IF (KD.GT.0) THEN
      IF (GMAX.LE.TOLG.AND.UMAX.LE.TOLG) THEN
      ITERM = 4
      RETURN
      ENDIF
      ENDIF
      IF (NIT.LE.0) THEN
      NTESX = 0
      NTESF = 0
      ENDIF
      IF (DMAX.LE.TOLX) THEN
      ITERM = 1
      NTESX = NTESX + 1
      IF (NTESX .GE. MTESX) RETURN
      ELSE
      NTESX = 0
      ENDIF
      TEMP=ABS(FO-F)/MAX(ABS(F),1.0D 0)
      IF (TEMP.LE.TOLF) THEN
      ITERM = 2
      NTESF = NTESF+1
      IF (NTESF.GE.MTESF) RETURN
      ELSE
      NTESF = 0
      ENDIF
    1 IF (NIT.GE.MIT) THEN
      ITERM = 11
      RETURN
      ENDIF
      IF (NFV.GE.MFV) THEN
      ITERM = 12
      RETURN
      ENDIF
      IF (NFG.GE.MFG) THEN
      ITERM = 13
      RETURN
      ENDIF
    2 ITERM = 0
      IF (N.GT.0.AND.NIT-KIT.GE.IRES1*N+IRES2) THEN
      IREST=MAX(IREST,1)
      ENDIF
      NIT = NIT + 1
      RETURN
      END
* SUBROUTINE PYRMB1               ALL SYSTEMS                98/12/01
* PORTABILITY : ALL SYSTEMS
* 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* OLD LINEAR CONSTRAINT OR AN OLD SIMPLE BOUND IS REMOVED FROM THE
* ACTIVE SET. TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION AND
* TRANSFORMED HESSIAN MATRIX APPROXIMATION ARE UPDATED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  IU  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  IU  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  IU  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RU  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RU  GN(NF)  TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION.
*  RU  H(NF*(NF+1)/2)  TRANSFORMED HESSIAN MATRIX APPROXIMATION.
*  RI  EPS8  TOLERANCE FOR CONSTRAINT TO BE REMOVED.
*  RI  UMAX  MAXIMUN ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
*  RI  GMAX  NORM OF THE TRANSFORMED GRADIENT.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  II  IOLD  INDEX OF THE REMOVED CONSTRAINT.
*  IA  KOLD  AUXILIARY VARIABLE.
*  IA  KREM  AUXILIARY VARIABLE.
*  IO  IER  ERROR INDICATOR.
*  IO  ITERM  TERMINATION INDICATOR.
*
* COMMON DATA :
*  IU  NREM  NUMBER OF CONSTRAINT DELETIONS.
*
* SUBPROGRAMS USED :
*  S   PLRMB0  CONSTRAINT DELETION.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PYRMB1(NF,N,IX,IC,ICA,CG,CR,CZ,G,GN,H,EPS8,UMAX,
     & GMAX,KBF,KBC,IOLD,KOLD,KREM,IER,ITERM)
      INTEGER NF,N,IX(*),IC(*),ICA(*),KBF,KBC,IOLD,KOLD,KREM,IER,
     $ ITERM
      DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*),H(*),EPS8,UMAX,
     & GMAX
      INTEGER I,J,K,KC,L
      INTEGER NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      COMMON /STATSQP/ NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      IF (KBC.GT.0) THEN
      IF (UMAX.GT.EPS8*GMAX) THEN
      CALL PLRMB0(NF,N,ICA,CG,CR,CZ,G,GN,IOLD,KREM,NREM,IER)
      IF (IER.LT.0) THEN
      ITERM=-16
      ELSE IF (IER.GT.0) THEN
      IOLD=0
      ELSE
      K=N*(N-1)/2
      CALL MXVSET(N,0.0D 0,H(K+1))
      H(K+N)=1.0D 0
      KC=ICA(NF-N+1)
      IF (KC.GT.0) THEN
      IC(KC)=-IC(KC)
      ELSE
      K=-KC
      IX(K)=-IX(K)
      ENDIF
      ENDIF
      ELSE
      IOLD=0
      ENDIF
      ELSE IF (KBF.GT.0) THEN
      IF (UMAX.GT.EPS8*GMAX)  THEN
      IX(IOLD)=MIN(ABS(IX(IOLD)),3)
      DO 1 I=N,KOLD,-1
      GN(I+1)=GN(I)
    1 CONTINUE
      GN(KOLD)=G(IOLD)
      N=N+1
      K=N*(N-1)/2
      L=K+N
      DO 3 I=N,KOLD,-1
      DO 2 J=I,1,-1
      IF (I.NE.KOLD.AND.J.NE.KOLD) THEN
      H(L)=H(K)
      K=K-1
      L=L-1
      ELSE IF (I.EQ.KOLD.AND.J.EQ.KOLD) THEN
      H(L)=1.0D 0
      L=L-1
      ELSE
      H(L)=0.0D 0
      L=L-1
      ENDIF
    2 CONTINUE
    3 CONTINUE
      ELSE
      IOLD=0
      KOLD=0
      ENDIF
      ENDIF
      RETURN
      END
* SUBROUTINE PYTRBD             ALL SYSTEMS                   98/12/01
* PORTABILITY : ALL SYSTEMS
* 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED
* AND TRANSFORMED. TEST VALUE DMAX IS DETERMINED.
*
* PARAMETERS :
*  II  NF DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RU  GO(NF)  GRADIENTS DIFFERENCE.
*  RI  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM CURRENT
*         REDUCED SUBSPACE.
*  RU  SN(NF)  TRANSFORMED DIRECTION VECTOR.
*  RI  R  VALUE OF THE STEPSIZE PARAMETER.
*  RU  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
*  RU  P  NEW VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RU  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  DMAX  MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
*  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
*         ITERS=0 FOR ZERO STEP.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*
* SUBPROGRAMS USED :
*  S   MXDRMM  PREMULTIPLICATION OF A VECTOR BY TRANSPOSE OF A DENSE
*         RECTANGULAR MATRIX.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  S   MXVMUL  DIAGONAL PREMULTIPLICATION OF A VECTOR.
*  S   MXVSAV  DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
*         SUBSTRACTED ONE.
*  S   MXVSCL  SCALING OF A VECTOR.
*
      SUBROUTINE PYTRBD(NF,N,X,IX,XO,G,GO,CZ,SN,R,F,FO,P,PO,DMAX,ITERS,
     & KBF,KBC)
      INTEGER NF,N,IX(*),ITERS,KBF,KBC
      DOUBLE PRECISION X(*),XO(*),G(*),GO(*),CZ(*),SN(*),R,F,FO,P,PO,
     & DMAX
      INTEGER I,K
      IF (ITERS.GT.0) THEN
      CALL MXVDIF(NF,X,XO,XO)
      CALL MXVDIF(NF,G,GO,GO)
      PO=R*PO
      P=R*P
      ELSE
      F = FO
      P = PO
      CALL MXVSAV(NF,X,XO)
      CALL MXVSAV(NF,G,GO)
      ENDIF
      DMAX = 0.0D 0
      IF (KBC.GT.0) THEN
      DO 1 I=1,NF
      DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0))
    1 CONTINUE
      IF (N.GT.0) THEN
      CALL MXVSCL(N,R,SN,XO)
      CALL MXVCOP(NF,GO,SN)
      CALL MXDRMM(NF,N,CZ,SN,GO)
      ENDIF
      ELSE IF (KBF.GT.0) THEN
      K=0
      DO 2 I=1,NF
      IF (IX(I).LT.0) GO TO 2
      K=K+1
      DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0))
      XO(K)=XO(I)
      GO(K)=GO(I)
    2 CONTINUE
      ELSE
      DO 3 I=1,NF
      DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0))
    3 CONTINUE
      ENDIF
      RETURN
      END
* SUBROUTINE PYTRBG               ALL SYSTEMS                98/12/01
* PORTABILITY : ALL SYSTEMS
* 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* GRADIENT OF THE OBJECTIVE FUNCTION IS SCALED AND REDUCED.
* TEST VALUES GMAX AND UMAX ARE COMPUTED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RU  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RO  GN(NF)  TRANSFORMED GRADIENT OF THE OBJECTIVE FUNCTION.
*  RI  EPS7  TOLERANCE FOR LINEAR INDEPENDENCE OF CONSTRAINTS.
*  RO  UMAX  MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
*  RO  GMAX  NORM OF THE TRANSFORMED GRADIENT.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  II  IOLD  INDEX OF THE REMOVED CONSTRAINT.
*  IA  KOLD  AUXILIARY VARIABLE.
*
* SUBPROGRAMS USED :
*  S   MXDRMM  PREMULTIPLICATION OF A VECTOR BY A ROWWISE STORED DENSE
*         RECTANGULAR MATRIX.
*  S   MXDPRB  BACK SUBSTITUTION.
*  S   MXVCOP  COPYING OF A VECTOR.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*  RF  MXVMAX  L-INFINITY NORM OF A VECTOR.
*  S   MXVMUL  DIAGONAL PREMULTIPLICATION OF A VECTOR.
*
      SUBROUTINE PYTRBG(NF,N,NC,IX,IC,ICA,CG,CR,CZ,G,GN,UMAX,GMAX,
     & KBF,KBC,IOLD,KOLD)
      INTEGER NF,N,NC,IX(*),IC(*),ICA(*),KBF,KBC,IOLD,KOLD
      DOUBLE PRECISION CG(*),CR(*),CZ(*),G(*),GN(*),UMAX,GMAX
      DOUBLE PRECISION TEMP,MXVMAX,MXVDOT
      INTEGER NCA,NCZ,I,J,K,KC
      IOLD=0
      KOLD=0
      UMAX=0.0D 0
      GMAX=0.0D 0
      IF (KBC.GT.0) THEN
      IF (NF.GT.N) THEN
      NCA=NF-N
      NCZ=N*NF
      CALL MXVCOP(NF,G,GN)
      DO 1 J=1,NCA
      K=ICA(J)
      IF (K.GT.0) THEN
      CZ(NCZ+J)=MXVDOT(NF,CG((K-1)*NF+1),GN)
      ELSE
      I=-K
      CZ(NCZ+J)=GN(I)
      ENDIF
    1 CONTINUE
      CALL MXDPRB(NCA,CR,CZ(NCZ+1),0)
      DO 2 J=1,NCA
      TEMP=CZ(NCZ+J)
      KC=ICA(J)
      IF (KC.GT.0) THEN
      K=IC(KC)
      ELSE
      I=-KC
      K=IX(I)
      ENDIF
      IF (K.LE.-5) THEN
      ELSE IF ((K.EQ.-1.OR.K.EQ.-3).AND.UMAX+TEMP.GE.0.0D 0) THEN
      ELSE IF ((K.EQ.-2.OR.K.EQ.-4).AND.UMAX-TEMP.GE.0.0D 0) THEN
      ELSE
      IOLD=J
      UMAX=ABS(TEMP)
      ENDIF
    2 CONTINUE
      ENDIF
      IF (N.GT.0) THEN
      CALL MXDRMM(NF,N,CZ,G,GN)
      GMAX=MXVMAX(N,GN)
      ENDIF
      ELSE IF (KBF.GT.0) THEN
      J=0
      IOLD=0
      KOLD=0
      DO 3 I=1,NF
      TEMP=G(I)
      K=IX(I)
      IF (K.GE.0) THEN
      J=J+1
      GN(J)=TEMP
      GMAX=MAX(GMAX,ABS(TEMP))
      ELSE IF (K.LE.-5) THEN
      ELSE IF ((K.EQ.-1.OR.K.EQ.-3).AND.UMAX+TEMP.GE.0.0D 0) THEN
      ELSE IF ((K.EQ.-2.OR.K.EQ.-4).AND.UMAX-TEMP.GE.0.0D 0) THEN
      ELSE
      IOLD=I
      KOLD=J+1
      UMAX=ABS(TEMP)
      ENDIF
    3 CONTINUE
      N=J
      ELSE
      DO 4 I=1,NF
      TEMP=G(I)
      GMAX=MAX(GMAX,ABS(TEMP))
    4 CONTINUE
      N=NF
      ENDIF
      RETURN
      END
* SUBROUTINE PYTRBH               ALL SYSTEMS                98/12/01
C PORTABILITY : ALL SYSTEMS
C 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* HESSIAN MATRIX OF THE OBJECTIVE FUNCTION OR ITS APPROXIMATION IS
* SCALED AND REDUCED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  RI  CR(NF*(NF+1)/2)  TRIANGULAR DECOMPOSITION OF KERNEL OF THE
*         ORTHOGONAL PROJECTION.
*  RI  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RI  H(NF*(NF+1)/2)  HESSIAN MATRIX OR ITS APPROXIMATION.
*  RA  S(NF)  AUXILIARY VECTOR.
*  II  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
*
* SUBPROGRAMS USED :
*  S   MXDSMM  MATRIX VECTOR PRODUCT.
*  S   MXVCOP  COPYING OF A VECTOR.
*  RF  MXVDOT  DOT PRODUCT OF TWO VECTORS.
*
      SUBROUTINE PYTRBH(NF,N,CR,CZ,H,S,LD,ITERS)
      INTEGER NF,N,LD,ITERS
      REAL*8 CR(*),CZ(*),H(*),S(*)
      REAL*8 TEMP,MXVDOT
      INTEGER NCA,NCR,ICZ,JCZ,I,J,K
      IF (LD.NE.2.OR.ITERS.EQ.0) RETURN
      IF (N.LE.0) RETURN
      NCA=NF-N
      NCR=NCA*(NCA+1)/2
      K=NCR
      JCZ=1
      DO 4 J=1,N
      CALL MXDSMM(NF,H,CZ(JCZ),S)
      ICZ=1
      DO 3 I=1,J
      K=K+1
      CR(K)=MXVDOT(NF,CZ(ICZ),S)
      ICZ=ICZ+NF
    3 CONTINUE
      JCZ=JCZ+NF
    4 CONTINUE
      CALL MXVCOP(N*(N+1)/2,CR(NCR+1),H)
      RETURN
      END
* SUBROUTINE PYTRBS               ALL SYSTEMS                98/12/01
* PORTABILITY : ALL SYSTEMS
* 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* SCALED AND REDUCED DIRECTION VECTOR IS BACK TRANSFORMED.
* VECTORS X,G AND VALUES F,P ARE SAVED.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  IU  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC  NUMBER OF LINEARIZED CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS.
*  RO  XO(NF)  SAVED VECTOR OF VARIABLES.
*  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
*  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RO  GO(NF)  SAVED GRADIENT OF THE OBJECTIVE FUNCTION.
*  RI  CF(NF)  VECTOR CONTAINING VALUES OF THE CONSTRAINT FUNCYIONS.
*  RO  CFD(NF)  VECTOR CONTAINING INCREMENTS OF THE CONSTRAINT
*         FUNCTIONS.
*  II  IC(NC)  VECTOR CONTAINING TYPES OF CONSTRAINTS.
*  RI  CL(NC)  VECTOR CONTAINING LOWER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CU(NC)  VECTOR CONTAINING UPPER BOUNDS FOR CONSTRAINT FUNCTIONS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  CZ(NF*NF)  MATRIX WHOSE COLUMNS ARE BASIC VECTORS FROM THE
*         CURRENT REDUCED SUBSPACE.
*  RI  SN(NF)  TRANSFORMED DIRECTION VECTOR.
*  RO  S(NF)  DIRECTION VECTOR.
*  RO  RO  SAVED VALUE OF THE STEPSIZE PARAMETER.
*  RO  FP  PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
*  RU  FO  SAVED VALUE OF THE OBJECTIVE FUNCTION.
*  RI  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  PO  SAVED VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  P  VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RU  RMAX  MAXIMUM VALUE OF THE STEPSIZE PARAMETER.
*  II  KBF  SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
*         KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
*  II  KBC  SPECIFICATION OF LINEAR CONSTRAINTS. KBC=0-NO LINEAR
*         CONSTRAINTS. KBC=1-ONE SIDED LINEAR CONSTRAINTS. KBC=2=TWO
*         SIDED LINEAR CONSTRAINTS.
*  IO  KREM  INDICATION OF LINEARLY DEPENDENT GRADIENTS.
*  IO  INEW  INDEX OF THE NEW ACTIVE FUNCTION.
*
* SUBPROGRAMS USED :
*  S   PLMAXS  DETERMINATION OF THE MAXIMUM STEPSIZE USING SIMPLE
*         BOUNDS.
*  S   PLMAXL  DETERMINATION OF THE MAXIMUM STEPSIZE USING LINEAR
*         CONSTRAINTS.
*  S   MXDCMM  MATRIX VECTOR PRODUCT.
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*
      SUBROUTINE PYTRBS(NF,N,NC,X,IX,XO,XL,XU,G,GO,CF,CFD,IC,CL,CU,CG,
     & CZ,SN,S,RO,FP,FO,F,PO,P,RMAX,KBF,KBC,KREM,INEW)
      INTEGER NF,N,NC,IX(*),IC(*),KBF,KBC,KREM,INEW
      DOUBLE PRECISION X(*),XO(*),XL(*),XU(*),G(*),GO(*),CF(*),CFD(*),
     & CL(*),CU(*),CG(*),CZ(*),SN(*),S(*),RO,FP,FO,F,PO,P,RMAX
      INTEGER I,K
      FP = FO
      RO = 0.0D 0
      FO = F
      PO = P
      CALL MXVCOP(NF,X,XO)
      CALL MXVCOP(NF,G,GO)
      IF (KBC.GT.0) THEN
      IF (N.GT.0) THEN
      CALL MXDCMM(NF,N,CZ,SN,S)
      INEW=0
      CALL PLMAXL(NF,NC,CF,CFD,IC,CL,CU,CG,S,RMAX,KBC,KREM,INEW)
      CALL PLMAXS(NF,X,IX,XL,XU,S,RMAX,KBF,KREM,INEW)
      ELSE
      CALL MXVSET(NF,0.0D 0,S)
      ENDIF
      ELSE IF (KBF.GT.0) THEN
      K=N+1
      DO 1 I=NF,1,-1
      IF (IX(I).LT.0) THEN
      S(I)=0.0D 0
      ELSE
      K=K-1
      S(I)=SN(K)
      ENDIF
    1 CONTINUE
      INEW=0
      CALL PLMAXS(NF,X,IX,XL,XU,S,RMAX,KBF,KREM,INEW)
      ENDIF
      RETURN
      END
* SUBROUTINE PYTRFD             ALL SYSTEMS                   90/12/01
* PORTABILITY : ALL SYSTEMS
* 90/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* PREPARATION OF VARIABLE METRIC UPDATE.
*
* PARAMETERS :
*  II  NF  DECLARED NUMBER OF VARIABLES.
*  II  NC  NUMBER OF CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RU  XO(NF)  SAVED VECTOR OF VARIABLES.
*  II  IAA(NF+1)  VECTOR CONTAINING INDICES OF ACTIVE FUNCTIONS.
*  RI  AG(NF*NA)  MATRIX WHOSE COLUMNS ARE GRADIENTS OF THE LINEAR
*          APPROXIMATED FUNCTIONS.
*  RI  AZ(NF+1)  VECTOR OF LAGRANGE MULTIPLIERS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RI  G(NF)  GRADIENT OF THE LAGRANGIAN FUNCTION.
*  RU  GO(NF)  SAVED GRADIENT OF THE LAGRANGIAN FUNCTION.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IU  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  RU  R  VALUE OF THE STEPSIZE PARAMETER.
*  RU  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  SAVED VALUE OF THE OBJECTIVE FUNCTION.
*  RU  P  VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RU  PO  SAVED VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  DMAX  RELATIVE STEPSIZE.
*  IO  ITERS  TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT
*         LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY
*         STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE.
*         ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE.
*         ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION.
*         ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL
*         DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION.
*
* SUBPROGRAMS USED :
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*  S   MXVSAV  DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
*         SUBSTRACTED ONE.
*
      SUBROUTINE PYTRFD(NF,NC,X,XO,IAA,AG,AZ,CG,G,GO,N,KD,LD,R,F,FO,P,
     +                  PO,DMAX,ITERS)
*     .. Scalar Arguments ..
      DOUBLE PRECISION DMAX,F,FO,P,PO,R
      INTEGER          ITERS,KD,LD,N,NC,NF
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION AG(*),AZ(*),CG(*),G(*),GO(*),X(*),XO(*)
      INTEGER          IAA(*)
*     ..
*     .. Local Scalars ..
      INTEGER          I,J,L
*     ..
*     .. External Subroutines ..
      EXTERNAL         MXVDIF,MXVDIR,MXVSAV,MXVSET
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC        ABS,MAX
*     ..
      CALL MXVSET(NF,0.0D0,G)
      DO 10 J = 1,NF - N
          L = IAA(J)
          IF (L.GT.NC) THEN
              L = L - NC
              CALL MXVDIR(NF,-AZ(J),AG((L-1)*NF+1),G,G)

          ELSE IF (L.GT.0) THEN
              CALL MXVDIR(NF,-AZ(J),CG((L-1)*NF+1),G,G)

          ELSE
              L = -L
              G(L) = G(L) - AZ(J)
          END IF

   10 CONTINUE
      IF (ITERS.GT.0) THEN
          CALL MXVDIF(NF,X,XO,XO)
          CALL MXVDIF(NF,G,GO,GO)
          PO = R*PO
          P = R*P

      ELSE
          R = 0.0D0
          F = FO
          P = PO
          CALL MXVSAV(NF,X,XO)
          CALL MXVSAV(NF,G,GO)
          LD = KD
      END IF

      DMAX = 0.0D0
      DO 20 I = 1,NF
          DMAX = MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D0))
   20 CONTINUE
      N = NF
      RETURN

      END
* SUBROUTINE PYTRND             ALL SYSTEMS                   91/12/01
C PORTABILITY : ALL SYSTEMS
C 91/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* DUAL RANGE SPACE QUADRATIC PROGRAMMING METHOD FOR MINIMAX
* APPROXIMATION.
*
* PARAMETERS :
*  II  NF DECLARED NUMBER OF VARIABLES.
*  II  N  ACTUAL NUMBER OF VARIABLES.
*  II  NC NUMBER OF CONSTRAINTS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RI  XN(NF)  VECTOR OF SCALING FACTORS.
*  RO  XO(NF)  SAVED VECTOR OF VARIABLES.
*  II  ICA(NF)  VECTOR CONTAINING INDICES OF ACTIVE CONSTRAINTS.
*  RI  CG(NF*NC)  MATRIX WHOSE COLUMNS ARE NORMALS OF THE LINEAR
*         CONSTRAINTS.
*  RO  CZ(NF)  VECTOR OF LAGRANGE MULTIPLIERS.
*  RO  CZS(NF)  SAVED VECTOR OF LAGRANGE MULTIPLIERS.
*  RI  G(NF)  GRADIENT OF THE LAGRANGIAN FUNCTION.
*  RI  GO(NF)  SAVED GRADIENT OF THE LAGRANGIAN FUNCTION.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RO  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
*  RO  P  NEW VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  CMAX  VALUE OF THE CONSTRAINT VIOLATION.
*  RO  CMAXO  SAVED VALUE OF THE CONSTRAINT VIOLATION.
*  RO  DMAX  MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
*
* COMMON DATA :
*  II  NORMF  SCALING SPECIFICATION. NORMF=0-NO SCALING PERFORMED.
*         NORMF=1-SCALING FACTORS ARE DETERMINED AUTOMATICALLY.
*         NORMF=2-SCALING FACTORS ARE SUPPLIED BY USER.
*  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
*         ITERS=0 FOR ZERO STEP.
*
* SUBPROGRAMS USED :
*  S   MXVCOP  COPYING OF A VECTOR.
*  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  S   MXVSET  INITIATION OF A VECTOR.
*  S   MXVSAV  DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
*         SUBSTRACTED ONE.
*
      SUBROUTINE PYTRND(NF,N,X,XO,ICA,CG,CZ,G,GO,R,F,FO,P,PO,CMAX,CMAXO,
     & DMAX,KD,LD,ITERS)
      INTEGER NF,N,KD,LD,ITERS
      INTEGER ICA(*)
      REAL*8 X(*),XO(*),CG(*),CZ(*),G(*),GO(*),R,F,FO,P,PO,CMAX,CMAXO,
     & DMAX
      INTEGER I,J,L
      DO 2 J=1,NF-N
      L=ICA(J)
      IF (L.GT.0) THEN
      CALL MXVDIR(NF,-CZ(J),CG((L-1)*NF+1),G,G)
      ELSE
      L=-L
      G(L)=G(L)-CZ(J)
      ENDIF
    2 CONTINUE
      IF (ITERS.GT.0) THEN
      CALL MXVDIF(NF,X,XO,XO)
      CALL MXVDIF(NF,G,GO,GO)
      PO=R*PO
      P=R*P
      ELSE
      F = FO
      P = PO
      CMAX = CMAXO
      CALL MXVSAV(NF,X,XO)
      CALL MXVSAV(NF,G,GO)
      LD = KD
      ENDIF
      DMAX = 0.0D0
      DO 1 I=1,NF
      DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D0))
    1 CONTINUE
      N=NF
      RETURN
      END
* SUBROUTINE PYTRUD             ALL SYSTEMS                   98/12/01
* PORTABILITY : ALL SYSTEMS
* 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED
* AND SCALED. TEST VALUE DMAX IS DETERMINED.
*
* PARAMETERS :
*  II  NF DECLARED NUMBER OF VARIABLES.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RI  G(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RU  GO(NF)  GRADIENTS DIFFERENCE.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RO  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
*  RO  P  NEW VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  DMAX  MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IO  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
*         ITERS=0 FOR ZERO STEP.
*
* SUBPROGRAMS USED :
*  S   PYSET1  DEGREE DEFINITION OF THE COMPUTED DERIVATIVES.
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  S   MXVSAV  DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
*         SUBSTRACTED ONE.
*
      SUBROUTINE PYTRUD(NF,X,XO,G,GO,R,F,FO,P,PO,DMAX,KD,LD,ITERS)
      INTEGER NF,KD,LD,ITERS
      REAL*8 X(*),XO(*),G(*),GO(*),R,F,FO,P,PO,DMAX
      INTEGER I
      IF (ITERS.GT.0) THEN
      CALL MXVDIF(NF,X,XO,XO)
      CALL MXVDIF(NF,G,GO,GO)
      PO=R*PO
      P=R*P
      ELSE
      F = FO
      P = PO
      CALL MXVSAV(NF,X,XO)
      CALL MXVSAV(NF,G,GO)
      LD=KD
      ENDIF
      DMAX = 0.0D 0
      DO 1 I=1,NF
      DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0))
    1 CONTINUE
      RETURN
      END
* SUBROUTINE PYTRUF             ALL SYSTEMS                   98/12/01
* PORTABILITY : ALL SYSTEMS
* 98/12/01 LU : ORIGINAL VERSION
*
* PURPOSE :
* VECTORS OF VARIABLES DIFFERENCE AND RIGHT HAND SIDES DIFFERENCE ARE
* COMPUTED AND SCALED. TEST VALUE DMAX IS DETERMINED.
*
* PARAMETERS :
*  II  NF DECLARED NUMBER OF VARIABLES.
*  II  NA NUMBER OF APPROXIMATED FUNCTIONS.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RI  AF(NA)  VECTOR OF RIGHT HAND SIDES.
*  RI  AFO(NA)  VECTOR OF RIGHT HAND SIDES DIFFERENCE.
*  RO  R  VALUE OF THE STEPSIZE PARAMETER.
*  RO  F  NEW VALUE OF THE OBJECTIVE FUNCTION.
*  RI  FO  OLD VALUE OF THE OBJECTIVE FUNCTION.
*  RO  P  NEW VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RI  PO  OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
*  RO  DMAX  MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
*  II  KD  DEGREE OF REQUIRED DERVATIVES.
*  IO  LD  DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES.
*  II  ITERS  TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
*         ITERS=0 FOR ZERO STEP.
*
* SUBPROGRAMS USED :
*  S   PYSET1  DEGREE DEFINITION OF THE COMPUTED DERIVATIVES.
*  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
*  S   MXVSAV  DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
*         SUBSTRACTED ONE.
*
      SUBROUTINE PYTRUF(NF,NA,X,XO,AF,AFO,R,F,FO,P,PO,DMAX,KD,LD,
     & ITERS)
      INTEGER NF,NA,KD,LD,ITERS
      REAL*8 X(*),XO(*),AF(*),AFO(*),R,F,FO,P,PO,DMAX
      INTEGER I
      IF (ITERS.GT.0) THEN
      CALL MXVDIF(NF,X,XO,XO)
      CALL MXVDIF(NA,AF,AFO,AFO)
      PO=R*PO
      P=R*P
      ELSE
      R = 0.0D 0
      F = FO
      P = PO
      CALL MXVSAV(NF,X,XO)
      CALL MXVSAV(NA,AF,AFO)
      LD=KD
      ENDIF
      DMAX = 0.0D 0
      DO 1 I=1,NF
      DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0))
    1 CONTINUE
      RETURN
      END


C ##############################################################################

Generated by  Doxygen 1.6.0   Back to index