************************************************************************
* SUBROUTINE PLIPU              ALL SYSTEMS                   97/01/22
* PURPOSE :
* EASY TO USE SUBROUTINE FOR LARGE-SCALE UNCONSTRAINED MINIMIZATION.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  RI  X(NF)  VECTOR OF VARIABLES.
*  II  IPAR(7)  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 PLIP.
*      IPAR(4)  ESTIMATION INDICATOR. IPAR(4)=0-MINIMUM IS NOT
*         ESTIMATED. IPAR(4)=1-MINIMUM IS ESTIMATED BY THE VALUE
*         RPAR(6).
*      IPAR(5)  METHOD USED. IPAR(5)=1-RANK-ONE METHOD.
*         IPAR(5)=2-RANK-TWO METHOD.
*      IPAR(6)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
*      IPAR(7)  MAXIMUM NUMBER OF VARIABLE METRIC UPDATES.
*  RI  RPAR(9)  REAL PARAMETERS:
*      RPAR(1)  MAXIMUM STEPSIZE.
*      RPAR(2)  TOLERANCE FOR THE CHANGE OF VARIABLES.
*      RPAR(3)  TOLERANCE FOR THE CHANGE OF FUNCTION VALUES.
*      RPAR(4)  TOLERANCE FOR THE FUNCTION FALUE.
*      RPAR(5)  TOLERANCE FOR THE GRADIENT NORM.
*      RPAR(6)  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
*      RPAR(7)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
*      RPAR(8)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
*      RPAR(9)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
*  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) SUBSEQUENT ITERATIONS.
*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
*                   MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS.
*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
*         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
*                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
*
* VARIABLES IN COMMON /STAT/ (STATISTICS) :
*  IO  NRES  NUMBER OF RESTARTS.
*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITIONS.
*  IO  NIN  NUMBER OF INNER ITERATIONS.
*  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   PLIP  LIMITED MEMORY SHIFTED VARIABLE METRIC METHOD IN THE
*            PRODUCT FORM.
*
* 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 THE 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 THE VECTOR OF VARIABLES AND GF(NF)
*         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
*
      SUBROUTINE PLIPU(NF,X,IPAR,RPAR,F,GMAX,IPRNT,ITERM)
      INTEGER NF,IPAR(7),IPRNT,ITERM
      DOUBLE PRECISION X(*),RPAR(9),F,GMAX
      INTEGER MF,NB,LGF,LS,LXO,LGO,LSO,LXM,LXR,LGR
      INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      DOUBLE PRECISION RA(:)
      ALLOCATABLE RA
      MF=IPAR(7)
      IF (MF.LE.0) MF=10
      ALLOCATE (RA(5*NF+NF*MF+2*MF))
      NB=0
*
*     POINTERS FOR AUXILIARY ARRAYS
*
      LGF=1
      LS=LGF+NF
      LXO=LS+NF
      LGO=LXO+NF
      LSO=LGO+NF
      LXM=LSO+NF
      LXR=LXM+NF*MF
      LGR=LXR+MF
      CALL PLIP(NF,NB,X,IPAR,RA,RA,RA(LGF),RA(LS),RA(LXO),RA(LGO),
     & RA(LSO),RA(LXM),RA(LXR),RA(LGR),RPAR(1),RPAR(2),RPAR(3),RPAR(4),
     & RPAR(5),RPAR(6),GMAX,F,IPAR(1),IPAR(2),IPAR(4),IPAR(5),MF,IPRNT,
     & ITERM)
      DEALLOCATE (RA)
      RETURN
      END
************************************************************************
* SUBROUTINE PLIPS              ALL SYSTEMS                   97/01/22
* PURPOSE :
* EASY TO USE SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  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-LOWER 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.
*  II  IPAR(7)  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 PLIP.
*      IPAR(4)  ESTIMATION INDICATOR. IPAR(4)=0-MINIMUM IS NOT
*         ESTIMATED. IPAR(4)=1-MINIMUM IS ESTIMATED BY THE VALUE
*         RPAR(6).
*      IPAR(5)  METHOD USED. IPAR(5)=1-RANK-ONE METHOD.
*         IPAR(5)=2-RANK-TWO METHOD.
*      IPAR(6)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
*      IPAR(7)  MAXIMUM NUMBER OF VARIABLE METRIC UPDATES.
*  RI  RPAR(9)  REAL PARAMETERS:
*      RPAR(1)  MAXIMUM STEPSIZE.
*      RPAR(2)  TOLERANCE FOR THE CHANGE OF VARIABLES.
*      RPAR(3)  TOLERANCE FOR THE CHANGE OF FUNCTION VALUES.
*      RPAR(4)  TOLERANCE FOR THE FUNCTION FALUE.
*      RPAR(5)  TOLERANCE FOR THE GRADIENT NORM.
*      RPAR(6)  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
*      RPAR(7)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
*      RPAR(8)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
*      RPAR(9)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
*  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) SUBSEQUENT ITERATIONS.
*         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
*                   MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS.
*         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
*         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
*         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
*                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
*
* VARIABLES IN COMMON /STAT/ (STATISTICS) :
*  IO  NRES  NUMBER OF RESTARTS.
*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITIONS.
*  IO  NIN  NUMBER OF INNER ITERATIONS.
*  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   PLIP  LIMITED MEMORY SHIFTED VARIABLE METRIC METHOD IN THE
*            PRODUCT FORM.
*
* 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 THE 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 THE VECTOR OF VARIABLES AND GF(NF)
*         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
*
      SUBROUTINE PLIPS(NF,X,IX,XL,XU,IPAR,RPAR,F,GMAX,IPRNT,ITERM)
      INTEGER NF,IX(*),IPAR(7),IPRNT,ITERM
      DOUBLE PRECISION X(*),XL(*),XU(*),RPAR(9),F,GMAX
      INTEGER MF,NB,LGF,LS,LXO,LGO,LSO,LXM,LXR,LGR
      INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      DOUBLE PRECISION RA(:)
      ALLOCATABLE RA
      MF=IPAR(7)
      IF (MF.LE.0) MF=10
      ALLOCATE (RA(5*NF+NF*MF+2*MF))
      NB=1
*
*     POINTERS FOR AUXILIARY ARRAYS
*
      LGF=1
      LS=LGF+NF
      LXO=LS+NF
      LGO=LXO+NF
      LSO=LGO+NF
      LXM=LSO+NF
      LXR=LXM+NF*MF
      LGR=LXR+MF
      CALL PLIP(NF,NB,X,IX,XL,XU,RA(LGF),RA(LS),RA(LXO),RA(LGO),
     & RA(LSO),RA(LXM),RA(LXR),RA(LGR),RPAR(1),RPAR(2),RPAR(3),RPAR(4),
     & RPAR(5),RPAR(6),GMAX,F,IPAR(1),IPAR(2),IPAR(4),IPAR(5),MF,IPRNT,
     & ITERM)
      DEALLOCATE (RA)
      RETURN
      END
************************************************************************
* SUBROUTINE PLIP               ALL SYSTEMS                   01/09/22
* PURPOSE :
* GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION THAT
* USE THE SHIFTED LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE
* PRODUCT FORM UPDATES.
*
* PARAMETERS :
*  II  NF  NUMBER OF VARIABLES.
*  II  NB  CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED.
*         NB>0-SIMPLE BOUNDS ACCEPTED.
*  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.
*  RA  GF(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
*  RO  S(NF)  DIRECTION VECTOR.
*  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
*  RI  GO(NF)  GRADIENTS DIFFERENCE.
*  RA  SO(NF)  AUXILIARY VECTOR.
*  RA  XM(NF*MF)  AUXILIARY VECTOR.
*  RA  XR(MF)  AUXILIARY VECTOR.
*  RA  GR(MF)  AUXILIARY VECTOR.
*  RI  XMAX  MAXIMUM STEPSIZE.
*  RI  TOLX  TOLERANCE FOR CHANGE OF VARIABLES.
*  RI  TOLF  TOLERANCE FOR CHANGE OF FUNCTION VALUES.
*  RI  TOLB  TOLERANCE FOR THE FUNCTION VALUE.
*  RI  TOLG  TOLERANCE FOR THE GRADIENT NORM.
*  RI  FMIN  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
*  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
*  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
*  II  MIT  MAXIMUM NUMBER OF ITERATIONS.
*  II  MFV  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
*  II  IEST  ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED.
*         IEST=1-MINIMUM IS ESTIMATED BY THE VALUE FMIN.
*  II  MET  METHOD USED. MET=1-RANK-ONE METHOD. MET=2-RANK-TWO
*         METHOD.
*  II  MF  NUMBER OF LIMITED MEMORY STEPS.
*  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=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
*                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
*         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
*         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
*
* VARIABLES IN COMMON /STAT/ (STATISTICS) :
*  IO  NRES  NUMBER OF RESTARTS.
*  IO  NDEC  NUMBER OF MATRIX DECOMPOSITION.
*  IO  NIN  NUMBER OF INNER ITERATIONS.
*  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   PCBS04  ELIMINATION OF BOX CONSTRAINT VIOLATIONS.
*  S   PS1L01  STEPSIZE SELECTION USING LINE SEARCH.
*  S   PULSP3  SHIFTED VARIABLE METRIC UPDATE.
*  S   PULVP3  SHIFTED LIMITED-MEMORY VARIABLE METRIC UPDATE.
*  S   PYADC0  ADDITION OF A BOX CONSTRAINT.
*  S   PYFUT1  TEST ON TERMINATION.
*  S   PYRMC0  DELETION OF A BOX CONSTRAINT.
*  S   PYTRCD  COMPUTATION OF PROJECTED DIFFERENCES FOR THE VARIABLE METRIC
*         UPDATE.
*  S   PYTRCG  COMPUTATION OF THE PROJECTED GRADIENT.
*  S   PYTRCS  COMPUTATION OF THE PROJECTED DIRECTION VECTOR.
*  S   MXDRMM  MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR
*         MATRIX A BY A VECTOR X.
*  S   MXDCMD  MULTIPLICATION OF A COLUMNWISE STORED DENSE RECTANGULAR
*         MATRIX A BY A VECTOR X AND ADDITION OF THE SCALED VECTOR
*         ALF*Y.
*  S   MXUCOP  COPYING OF A VECTOR.
*  S   MXUDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
*  RF  MXUDOT  DOT PRODUCT OF TWO VECTORS.
*  S   MXUNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN.
*  S   MXUZER  VECTOR ELEMENTS CORRESPONDING TO ACTIVE BOUNDS ARE SET
*         TO ZERO.
*  S   MXVCOP  COPYING 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 THE 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 THE VECTOR OF VARIABLES AND GF(NF)
*         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
*
* METHOD :
* HYBRID METHOD WITH SPARSE MARWIL UPDATES FOR SPARSE LEAST SQUARES
* PROBLEMS.
*
      SUBROUTINE PLIP(NF,NB,X,IX,XL,XU,GF,S,XO,GO,SO,XM,XR,GR,XMAX,
     & TOLX,TOLF,TOLB,TOLG,FMIN,GMAX,F,MIT,MFV,IEST,MET,MF,IPRNT,ITERM)
      INTEGER NF,NB,IX(*),MIT,MFV,IEST,MET,MF,IPRNT,ITERM
      DOUBLE PRECISION X(*),XL(*),XU(*),GF(*),S(*),XO(*),GO(*),SO(*),
     & XM(*),XR(*),GR(*),XMAX,TOLX,TOLF,TOLG,TOLB,FMIN,GMAX,F
      INTEGER ITERD,ITERS,ITERH,KD,LD,NTESX,NTESF,MTESX,MTESF,MRED,KIT,
     & IREST,KBF,MEC,MES,MES1,MES2,MES3,MAXST,ISYS,ITES,INITS,KTERS,
     & IRES1,IRES2,NRED,INEW,IOLD,I,NN,N,MFG,META,MET3
      DOUBLE PRECISION R,RO,RP,FO,FP,P,PO,PP,GNORM,SNORM,RMIN,RMAX,
     & UMAX,FMAX,DMAX,ETA0,ETA9,EPS8,EPS9,ALF1,ALF2,PAR1,PAR2,PAR,TOLD,
     & TOLS,TOLP
      DOUBLE PRECISION MXUDOT
      INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
      IF (ABS(IPRNT).GT.1) WRITE(6,'(1X,''ENTRY TO PLIP :'')')
*
*     INITIATION
*
      KBF=0
      IF (NB.GT.0) KBF=2
      NRES=0
      NDEC=0
      NIN=0
      NIT=0
      NFV=0
      NFG=0
      NFH=0
      ISYS=0
      ITES=1
      MTESX=2
      MTESF=2
      INITS=2
      ITERM=0
      ITERD=0
      ITERS=2
      ITERH=0
      KTERS=3
      IREST=0
      IRES1=999
      IRES2=0
      MRED=10
      META=1
      MET3=4
      MEC=4
      MES=4
      MES1=2
      MES2=2
      MES3=2
      ETA0=1.0D-15
      ETA9=1.0D 120
      EPS8=1.00D 0
      EPS9=1.00D-8
      ALF1=1.0D-10
      ALF2=1.0D 10
      RMAX=ETA9
      DMAX=ETA9
      FMAX=1.0D 20
      IF (IEST.LE.0) FMIN=-1.0D 60
      IF (IEST.GT.0) IEST=1
      IF (XMAX.LE.0.0D 0) XMAX=1.0D 16
      IF (TOLX.LE.0.0D 0) TOLX=1.0D-16
      IF (TOLF.LE.0.0D 0) TOLF=1.0D-14
      IF (TOLG.LE.0.0D 0) TOLG=1.0D-6
      IF (TOLB.LE.0.0D 0) TOLB=FMIN+1.0D-16
      TOLD=1.0D-4
      TOLS=1.0D-4
      TOLP=0.9D 0
      IF (MET.LE.0) MET=2
      IF (MIT.LE.0) MIT=9000
      IF (MFV.LE.0) MFV=9000
      MFG=MFV
      KD= 1
      LD=-1
      KIT=-(IRES1*NF+IRES2)
      FO=FMIN
*
*     INITIAL OPERATIONS WITH SIMPLE BOUNDS
*
      IF (KBF.GT.0) THEN
      DO 2 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
    2 CONTINUE
      CALL PCBS04(NF,X,IX,XL,XU,EPS9,KBF)
      CALL PYADC0(NF,N,X,IX,XL,XU,INEW)
      END IF
      IF (ITERM.NE.0) GO TO 11190
      CALL OBJ(NF,X,F)
      NFV=NFV+1
      CALL DOBJ(NF,X,GF)
      NFG=NFG+1
11120 CONTINUE
      CALL PYTRCG(NF,NF,IX,GF,UMAX,GMAX,KBF,IOLD)
      IF (ABS(IPRNT).GT.1)
     & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X,
     & ''F='', G16.9,2X,''G='',E10.3)') NIT,NFV,NFG,F,GMAX
      CALL PYFUT1(NF,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)
      IF (ITERM.NE.0) GO TO 11190
      IF (KBF.GT.0.AND.RMAX.GT.0.0D 0) THEN
      CALL PYRMC0(NF,N,IX,GF,EPS8,UMAX,GMAX,RMAX,IOLD,IREST)
      END IF
11130 CONTINUE
      IF (IREST.GT.0) THEN
      NN=0
      PAR=1.0D 0
      LD=MIN(LD,1)
      IF (KIT.LT.NIT) THEN
        NRES=NRES+1
        KIT = NIT
      ELSE
        ITERM=-10
        IF (ITERS.LT.0) ITERM=ITERS-5
      END IF
      END IF
      IF (ITERM.NE.0) GO TO 11190
*
*     DIRECTION DETERMINATION
*
      GNORM=SQRT(MXUDOT(NF,GF,GF,IX,KBF))
*
*     NEWTON LIKE STEP
*
      CALL MXUNEG(NF,GF,S,IX,KBF)
      CALL MXDRMM(NF,NN,XM,S,GR)
      CALL MXDCMD(NF,NN,XM,GR,PAR,S,S)
      CALL MXUZER(NF,S,IX,KBF)
      ITERD=1
      SNORM=SQRT(MXUDOT(NF,S,S,IX,KBF))
*
*     TEST ON DESCENT DIRECTION AND PREPARATION OF LINE SEARCH
*
      IF (KD.GT.0) P=MXUDOT(NF,GF,S,IX,KBF)
      IF (ITERD.LT.0) THEN
        ITERM=ITERD
      ELSE
*
*     TEST ON DESCENT DIRECTION
*
      IF (SNORM.LE.0.0D 0) THEN
        IREST=MAX(IREST,1)
      ELSE IF (P+TOLD*GNORM*SNORM.LE.0.0D 0) THEN
        IREST=0
      ELSE
*
*     UNIFORM DESCENT CRITERION
*
      IREST=MAX(IREST,1)
      END IF
      IF (IREST.EQ.0) THEN
*
*     PREPARATION OF LINE SEARCH
*
        NRED = 0
        RMIN=ALF1*GNORM/SNORM
        RMAX=MIN(ALF2*GNORM/SNORM,XMAX/SNORM)
      END IF
      END IF
      IF (ITERM.NE.0) GO TO 11190
      IF (IREST.NE.0) GO TO 11130
      CALL PYTRCS(NF,X,IX,XO,XL,XU,GF,GO,S,RO,FP,FO,F,PO,P,RMAX,ETA9,
     & KBF)
      IF (RMAX.EQ.0.0D 0) GO TO 11175
11170 CONTINUE
      CALL PS1L01(R,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)
      IF (ISYS.EQ.0) GO TO 11174
      CALL MXUDIR(NF,R,S,XO,X,IX,KBF)
      CALL PCBS04(NF,X,IX,XL,XU,EPS9,KBF)
      CALL OBJ(NF,X,F)
      NFV=NFV+1
      CALL DOBJ(NF,X,GF)
      NFG=NFG+1
      P=MXUDOT(NF,GF,S,IX,KBF)
      GO TO 11170
11174 CONTINUE
      IF (ITERS.LE.0) THEN
      R=0.0D 0
      F=FO
      P=PO
      CALL MXVCOP(NF,XO,X)
      CALL MXVCOP(NF,GO,GF)
      IREST=MAX(IREST,1)
      LD=KD
      GO TO 11130
      END IF
      CALL MXUNEG(NF,GO,S,IX,KBF)
      CALL PYTRCD(NF,X,IX,XO,GF,GO,R,F,FO,P,PO,DMAX,KBF,KD,LD,ITERS)
      CALL MXUCOP(NF,GF,SO,IX,KBF)
      IF (NN.LT.MF) THEN
      CALL PULSP3(NF,NN,MF,XM,GR,XO,GO,R,PO,PAR,ITERH,MET3)
      ELSE
      CALL PULVP3(NF,NN,XM,XR,GR,S,SO,XO,GO,R,PO,PAR,ITERH,MEC,MET3,
     & MET)
      END IF
11175 CONTINUE
      IF (ITERH.NE.0) IREST=MAX(IREST,1)
      IF (KBF.GT.0) CALL PYADC0(NF,N,X,IX,XL,XU,INEW)
      GO TO 11120
11190 CONTINUE
      IF (IPRNT.GT.1.OR.IPRNT.LT.0)
     & WRITE(6,'(1X,''EXIT FROM PLIP :'')')
      IF (IPRNT.NE.0)
     & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X,
     & ''F='', G16.9,2X,''G='',E10.3,2X,''ITERM='',I3)') NIT,NFV,NFG,
     & F,GMAX,ITERM
      IF (IPRNT.LT.0)
     & WRITE (6,'(1X,''X='',5(G14.7,1X):/(3X,5(G14.7,1X)))')
     & (X(I),I=1,NF)
      RETURN
      END