************************************************************************
* 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.
*  II  IPAR(6)  INTEGER PAREMETERS:
*      IPAR(1)  MAXIMUM NUMBER OF ITERATIONS.
*      IPAR(2)  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
*      IPAR(3)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PSQP.
*      IPAR(4)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PSQP.
*      IPAR(5)  VARIABLE METRIC UPDATE USED. IPAR(5)=1-THE BFGS UPDATE.
*         IPAR(5)-THE HOSHINO UPDATE.
*      IPAR(6)  CORRECTION OF THE VARIABLE METRIC UPDATE IF A NEGATIVE
*         CURVATURE OCCURS. IPAR(6)=1-NO CORRECTION. IPAR(6)=2-POWELL'S
*         CORRECTION.
*  RI  RPAR(5)  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)  PENALTY COEFFICIENT.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE OF THE LAGRANGIAN FUNCTION.
*  RO  CMAX  MAXIMUM CONSTRAINT VIOLATION.
*  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.
*  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
*         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
*                   MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS.
*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
*                   MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS.
*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
*         IF ITERM=-6, THEN THE TERMINATION CRITERION HAS NOT BEEN
*         SATISFIED, BUT THE POINT OBTAINED IF USUALLY ACCEPTABLE.
*
* VARIABLES IN COMMON /STAT/ (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,IPAR,RPAR,
     +                 F,GMAX,CMAX,IPRNT,ITERM)
*
*     POINTERS FOR AUXILIARY ARRAYS
*
      DOUBLE PRECISION F,CMAX,GMAX
      INTEGER          IPRNT,ITERM,NB,NC,NF
      DOUBLE PRECISION CF(*),CL(*),CU(*),RPAR(5),X(*),XL(*),XU(*)
      INTEGER          IC(*),IPAR(6),IX(*)
      INTEGER          NADD,NDEC,NFG,NFH,NFV,NIT,NREM,NRES
      INTEGER          LCFD,LCFO,LCG,LCP,LCR,LCZ,LG,LGC,LGF,LGO,LH,LIA,
     +                 LS,LXO
      COMMON           /STAT/NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      INTEGER IA(:)
      DOUBLE PRECISION RA(:)
      ALLOCATABLE IA,RA
      ALLOCATE (IA(NF),RA((NF+NC+8)*NF+3*NC+1))
      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),CMAX,GMAX,F,IPAR(1),IPAR(2),IPAR(5),
     +          IPAR(6),IPRNT,ITERM)
      DEALLOCATE(IA,RA)
      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(NC)  VECTOR CONTAINING SAVED VALUES OF THE CONSTRAINT
*         FUNCTIONS.
*  RA  CFD(NC)  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  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  MIT  MAXIMUN NUMBER OF ITERATIONS.
*  II  MFV  MAXIMUN NUMBER OF FUNCTION EVALUATIONS.
*  II  MET  VARIABLE METRIC UPDATE USED. MET=1-THE BFGS UPDATE.
*         MET=2-THE HOSHINO UPDATE.
*  II  MEC  CORRECTION IF THE NEGATIVE CURVATURE OCCURS.
*         MEC=1-CORRECTION SUPPRESSED. MEC=2-POWELL'S CORRECTION.
*  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  IOUT  PRINT OUTPUT UNIT NUMBER.
*  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
*         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
*                   MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS.
*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
*                   MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS.
*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
*         IF ITERM=-6, THEN THE TERMINATION CRITERION HAS NOT BEEN
*         SATISFIED, BUT THE POINT OBTAINED IF USUALLY ACCEPTABLE.
*
* VARIABLES IN COMMON /STAT/ (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,
     +                RPF,CMAX,GMAX,F,MIT,MFV,MET,MEC,IPRNT,IOUT,ITERM)
      DOUBLE PRECISION F,CMAX,GMAX,RPF,TOLC,TOLD,TOLG,TOLS,TOLX,XMAX
      INTEGER          IPRNT,ITERM,MET,MET1,MEC,MES,MFV,MIT,NB,NC,NF
      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(*)
      INTEGER          NADD,NDEC,NFG,NFH,NFV,NIT,NREM,NRES
      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,K,
     +                 NTESX,IEST,INITS,KTERS,MAXST,ISYS,MFP,NRED,IPOM,
     +                 LDS
      DOUBLE PRECISION MXVDOT, MXVMAX
      COMMON           /STAT/NRES,NDEC,NREM,NADD,NIT,NFV,NFG,NFH
      IF (ABS(IPRNT).GT.1) WRITE (IOUT,'(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+16
      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
      TOLD = 1.0D-8
      TOLS = 1.0D-4
      IF (RPF.LE.0.0D0) RPF = 1.0D-4
      IF (MET.LE.0) MET = 1
      MET1 = 2
      IF (MEC.LE.0) MEC = 2
      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
      IF (ABS(IPRNT).GT.1)
     + WRITE (IOUT,'(1X,''NIT='',I4,2X,''NFV='',I4,2X,''NFG='',I4,2X,
     + ''F='',G12.6,2X,''C='',E7.1,2X,''G='',E7.1)')
     +  NIT,NFV,NFG,F,CMAX,GMAX
*
*     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 = 11
          GO TO 80
      END IF
      IF (NFV.GE.MFV) THEN
          ITERM = 12
          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)
      IF (ISYS.EQ.0) GO TO 11174
C      GO TO (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)
      GO TO 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,MET1,MEC)
C      IF (MER.GT.0.AND.ITERH.GT.0) IREST=1
*
*     END OF THE ITERATION
*
      GO TO 40

   80 IF (IPRNT.GT.1 .OR. IPRNT.LT.0)
     + WRITE (IOUT,'(1X,''EXIT FROM PSQP :'')')
      IF (IPRNT.NE.0)
     + WRITE (IOUT,'(1X,''NIT='',I4,2X,''NFV='',I4,2X,''NFG='',I4,2X,
     + ''F='',G12.6,2X,''C='',E7.1,2X,''G='',E7.1,2X,''ITERM='',I3)')
     + NIT,NFV,NFG,F,CMAX,GMAX,ITERM
      IF (IPRNT.LT.0)
     + WRITE (IOUT,'(1X,''X='',5(G14.7,1X):/(3X,5(G14.7,1X)))')
     + (X(I),I=1,NF)
      RETURN
      END