* 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   PLMINN  DETERMINATION OF THE NEW ACTIVE LINEAR CONSTRAINT.
*  S   PLADR1  ADDITION OF A NEW ACTIVE CONSTRAINT.
*  S   PLRMR0  CONSTRAIN DELETION.
*  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
      DOUBLE PRECISION X(*),XL(*),XU(*),CF(*),CFD(*),CL(*),CU(*),CG(*),
     & CR(*),CZ(*),G(*),GO(*),H(*),S(*),ETA2,ETA9,EPS7,EPS9,UMAX,GMAX
      DOUBLE PRECISION 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 /STAT/ 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,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