luksan/plip.for

   1 ************************************************************************
   2 * SUBROUTINE PLIPU              ALL SYSTEMS                   97/01/22
   3 * PURPOSE :
   4 * EASY TO USE SUBROUTINE FOR LARGE-SCALE UNCONSTRAINED MINIMIZATION.
   5 *
   6 * PARAMETERS :
   7 *  II  NF  NUMBER OF VARIABLES.
   8 *  RI  X(NF)  VECTOR OF VARIABLES.
   9 *  II  IPAR(7)  INTEGER PAREMETERS:
  10 *      IPAR(1)  MAXIMUM NUMBER OF ITERATIONS.
  11 *      IPAR(2)  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
  12 *      IPAR(3)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
  13 *      IPAR(4)  ESTIMATION INDICATOR. IPAR(4)=0-MINIMUM IS NOT
  14 *         ESTIMATED. IPAR(4)=1-MINIMUM IS ESTIMATED BY THE VALUE
  15 *         RPAR(6).
  16 *      IPAR(5)  METHOD USED. IPAR(5)=1-RANK-ONE METHOD.
  17 *         IPAR(5)=2-RANK-TWO METHOD.
  18 *      IPAR(6)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
  19 *      IPAR(7)  MAXIMUM NUMBER OF VARIABLE METRIC UPDATES.
  20 *  RI  RPAR(9)  REAL PARAMETERS:
  21 *      RPAR(1)  MAXIMUM STEPSIZE.
  22 *      RPAR(2)  TOLERANCE FOR THE CHANGE OF VARIABLES.
  23 *      RPAR(3)  TOLERANCE FOR THE CHANGE OF FUNCTION VALUES.
  24 *      RPAR(4)  TOLERANCE FOR THE FUNCTION FALUE.
  25 *      RPAR(5)  TOLERANCE FOR THE GRADIENT NORM.
  26 *      RPAR(6)  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
  27 *      RPAR(7)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
  28 *      RPAR(8)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
  29 *      RPAR(9)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
  30 *  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
  31 *  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
  32 *  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
  33 *         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
  34 *         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
  35 *         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
  36 *         RESULTS.
  37 *  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
  38 *         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
  39 *                   MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS.
  40 *         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
  41 *                   MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS.
  42 *         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
  43 *         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
  44 *         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
  45 *                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
  46 *         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
  47 *         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
  48 *
  49 * VARIABLES IN COMMON /STAT/ (STATISTICS) :
  50 *  IO  NRES  NUMBER OF RESTARTS.
  51 *  IO  NDEC  NUMBER OF MATRIX DECOMPOSITIONS.
  52 *  IO  NIN  NUMBER OF INNER ITERATIONS.
  53 *  IO  NIT  NUMBER OF ITERATIONS.
  54 *  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
  55 *  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
  56 *  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
  57 *
  58 * SUBPROGRAMS USED :
  59 *  S   PLIP  LIMITED MEMORY SHIFTED VARIABLE METRIC METHOD IN THE
  60 *            PRODUCT FORM.
  61 *
  62 * EXTERNAL SUBROUTINES :
  63 *  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION.
  64 *         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER
  65 *         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS
  66 *         THE VALUE OF THE OBJECTIVE FUNCTION.
  67 *  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION.
  68 *         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER
  69 *         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF)
  70 *         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
  71 *
  72       SUBROUTINE PLIPU(NF,X,IPAR,RPAR,F,GMAX,IPRNT,ITERM)
  73       INTEGER NF,IPAR(7),IPRNT,ITERM
  74       DOUBLE PRECISION X(*),RPAR(9),F,GMAX
  75       INTEGER MF,NB,LGF,LS,LXO,LGO,LSO,LXM,LXR,LGR
  76       INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
  77       COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
  78       DOUBLE PRECISION RA(:)
  79       ALLOCATABLE RA
  80       MF=IPAR(7)
  81       IF (MF.LE.0) MF=10
  82       ALLOCATE (RA(5*NF+NF*MF+2*MF))
  83       NB=0
  84 *
  85 *     POINTERS FOR AUXILIARY ARRAYS
  86 *
  87       LGF=1
  88       LS=LGF+NF
  89       LXO=LS+NF
  90       LGO=LXO+NF
  91       LSO=LGO+NF
  92       LXM=LSO+NF
  93       LXR=LXM+NF*MF
  94       LGR=LXR+MF
  95       CALL PLIP(NF,NB,X,IPAR,RA,RA,RA(LGF),RA(LS),RA(LXO),RA(LGO),
  96      & RA(LSO),RA(LXM),RA(LXR),RA(LGR),RPAR(1),RPAR(2),RPAR(3),RPAR(4),
  97      & RPAR(5),RPAR(6),GMAX,F,IPAR(1),IPAR(2),IPAR(4),IPAR(5),MF,IPRNT,
  98      & ITERM)
  99       DEALLOCATE (RA)
 100       RETURN
 101       END
 102 ************************************************************************
 103 * SUBROUTINE PLIPS              ALL SYSTEMS                   97/01/22
 104 * PURPOSE :
 105 * EASY TO USE SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION.
 106 *
 107 * PARAMETERS :
 108 *  II  NF  NUMBER OF VARIABLES.
 109 *  RI  X(NF)  VECTOR OF VARIABLES.
 110 *  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE
 111 *         X(I) IS UNBOUNDED. IX(I)=1-LOWER BOUND XL(I).LE.X(I).
 112 *         IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND
 113 *         XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED.
 114 *  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
 115 *  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
 116 *  II  IPAR(7)  INTEGER PAREMETERS:
 117 *      IPAR(1)  MAXIMUM NUMBER OF ITERATIONS.
 118 *      IPAR(2)  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
 119 *      IPAR(3)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
 120 *      IPAR(4)  ESTIMATION INDICATOR. IPAR(4)=0-MINIMUM IS NOT
 121 *         ESTIMATED. IPAR(4)=1-MINIMUM IS ESTIMATED BY THE VALUE
 122 *         RPAR(6).
 123 *      IPAR(5)  METHOD USED. IPAR(5)=1-RANK-ONE METHOD.
 124 *         IPAR(5)=2-RANK-TWO METHOD.
 125 *      IPAR(6)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
 126 *      IPAR(7)  MAXIMUM NUMBER OF VARIABLE METRIC UPDATES.
 127 *  RI  RPAR(9)  REAL PARAMETERS:
 128 *      RPAR(1)  MAXIMUM STEPSIZE.
 129 *      RPAR(2)  TOLERANCE FOR THE CHANGE OF VARIABLES.
 130 *      RPAR(3)  TOLERANCE FOR THE CHANGE OF FUNCTION VALUES.
 131 *      RPAR(4)  TOLERANCE FOR THE FUNCTION FALUE.
 132 *      RPAR(5)  TOLERANCE FOR THE GRADIENT NORM.
 133 *      RPAR(6)  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
 134 *      RPAR(7)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
 135 *      RPAR(8)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
 136 *      RPAR(9)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PLIP.
 137 *  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
 138 *  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
 139 *  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
 140 *         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
 141 *         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
 142 *         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
 143 *         RESULTS.
 144 *  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
 145 *         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
 146 *                   MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS.
 147 *         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
 148 *                   MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS.
 149 *         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
 150 *         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
 151 *         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
 152 *                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
 153 *         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
 154 *         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
 155 *
 156 * VARIABLES IN COMMON /STAT/ (STATISTICS) :
 157 *  IO  NRES  NUMBER OF RESTARTS.
 158 *  IO  NDEC  NUMBER OF MATRIX DECOMPOSITIONS.
 159 *  IO  NIN  NUMBER OF INNER ITERATIONS.
 160 *  IO  NIT  NUMBER OF ITERATIONS.
 161 *  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
 162 *  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
 163 *  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
 164 *
 165 * SUBPROGRAMS USED :
 166 *  S   PLIP  LIMITED MEMORY SHIFTED VARIABLE METRIC METHOD IN THE
 167 *            PRODUCT FORM.
 168 *
 169 * EXTERNAL SUBROUTINES :
 170 *  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION.
 171 *         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER
 172 *         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS
 173 *         THE VALUE OF THE OBJECTIVE FUNCTION.
 174 *  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION.
 175 *         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER
 176 *         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF)
 177 *         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
 178 *
 179       SUBROUTINE PLIPS(NF,X,IX,XL,XU,IPAR,RPAR,F,GMAX,IPRNT,ITERM)
 180       INTEGER NF,IX(*),IPAR(7),IPRNT,ITERM
 181       DOUBLE PRECISION X(*),XL(*),XU(*),RPAR(9),F,GMAX
 182       INTEGER MF,NB,LGF,LS,LXO,LGO,LSO,LXM,LXR,LGR
 183       INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
 184       COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
 185       DOUBLE PRECISION RA(:)
 186       ALLOCATABLE RA
 187       MF=IPAR(7)
 188       IF (MF.LE.0) MF=10
 189       ALLOCATE (RA(5*NF+NF*MF+2*MF))
 190       NB=1
 191 *
 192 *     POINTERS FOR AUXILIARY ARRAYS
 193 *
 194       LGF=1
 195       LS=LGF+NF
 196       LXO=LS+NF
 197       LGO=LXO+NF
 198       LSO=LGO+NF
 199       LXM=LSO+NF
 200       LXR=LXM+NF*MF
 201       LGR=LXR+MF
 202       CALL PLIP(NF,NB,X,IX,XL,XU,RA(LGF),RA(LS),RA(LXO),RA(LGO),
 203      & RA(LSO),RA(LXM),RA(LXR),RA(LGR),RPAR(1),RPAR(2),RPAR(3),RPAR(4),
 204      & RPAR(5),RPAR(6),GMAX,F,IPAR(1),IPAR(2),IPAR(4),IPAR(5),MF,IPRNT,
 205      & ITERM)
 206       DEALLOCATE (RA)
 207       RETURN
 208       END
 209 ************************************************************************
 210 * SUBROUTINE PLIP               ALL SYSTEMS                   01/09/22
 211 * PURPOSE :
 212 * GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION THAT
 213 * USE THE SHIFTED LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE
 214 * PRODUCT FORM UPDATES.
 215 *
 216 * PARAMETERS :
 217 *  II  NF  NUMBER OF VARIABLES.
 218 *  II  NB  CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED.
 219 *         NB>0-SIMPLE BOUNDS ACCEPTED.
 220 *  RI  X(NF)  VECTOR OF VARIABLES.
 221 *  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE
 222 *         X(I) IS UNBOUNDED. IX(I)=1-LOVER BOUND XL(I).LE.X(I).
 223 *         IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND
 224 *         XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED.
 225 *  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
 226 *  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
 227 *  RA  GF(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
 228 *  RO  S(NF)  DIRECTION VECTOR.
 229 *  RU  XO(NF)  VECTORS OF VARIABLES DIFFERENCE.
 230 *  RI  GO(NF)  GRADIENTS DIFFERENCE.
 231 *  RA  SO(NF)  AUXILIARY VECTOR.
 232 *  RA  XM(NF*MF)  AUXILIARY VECTOR.
 233 *  RA  XR(MF)  AUXILIARY VECTOR.
 234 *  RA  GR(MF)  AUXILIARY VECTOR.
 235 *  RI  XMAX  MAXIMUM STEPSIZE.
 236 *  RI  TOLX  TOLERANCE FOR CHANGE OF VARIABLES.
 237 *  RI  TOLF  TOLERANCE FOR CHANGE OF FUNCTION VALUES.
 238 *  RI  TOLB  TOLERANCE FOR THE FUNCTION VALUE.
 239 *  RI  TOLG  TOLERANCE FOR THE GRADIENT NORM.
 240 *  RI  FMIN  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
 241 *  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
 242 *  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
 243 *  II  MIT  MAXIMUM NUMBER OF ITERATIONS.
 244 *  II  MFV  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
 245 *  II  IEST  ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED.
 246 *         IEST=1-MINIMUM IS ESTIMATED BY THE VALUE FMIN.
 247 *  II  MET  METHOD USED. MET=1-RANK-ONE METHOD. MET=2-RANK-TWO
 248 *         METHOD.
 249 *  II  MF  NUMBER OF LIMITED MEMORY STEPS.
 250 *  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
 251 *         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
 252 *         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
 253 *         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
 254 *         RESULTS.
 255 *  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
 256 *         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
 257 *                   MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS.
 258 *         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
 259 *                   MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS.
 260 *         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
 261 *         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
 262 *         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
 263 *                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
 264 *         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
 265 *         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
 266 *
 267 * VARIABLES IN COMMON /STAT/ (STATISTICS) :
 268 *  IO  NRES  NUMBER OF RESTARTS.
 269 *  IO  NDEC  NUMBER OF MATRIX DECOMPOSITION.
 270 *  IO  NIN  NUMBER OF INNER ITERATIONS.
 271 *  IO  NIT  NUMBER OF ITERATIONS.
 272 *  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
 273 *  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
 274 *  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
 275 *
 276 * SUBPROGRAMS USED :
 277 *  S   PCBS04  ELIMINATION OF BOX CONSTRAINT VIOLATIONS.
 278 *  S   PS1L01  STEPSIZE SELECTION USING LINE SEARCH.
 279 *  S   PULSP3  SHIFTED VARIABLE METRIC UPDATE.
 280 *  S   PULVP3  SHIFTED LIMITED-MEMORY VARIABLE METRIC UPDATE.
 281 *  S   PYADC0  ADDITION OF A BOX CONSTRAINT.
 282 *  S   PYFUT1  TEST ON TERMINATION.
 283 *  S   PYRMC0  DELETION OF A BOX CONSTRAINT.
 284 *  S   PYTRCD  COMPUTATION OF PROJECTED DIFFERENCES FOR THE VARIABLE METRIC
 285 *         UPDATE.
 286 *  S   PYTRCG  COMPUTATION OF THE PROJECTED GRADIENT.
 287 *  S   PYTRCS  COMPUTATION OF THE PROJECTED DIRECTION VECTOR.
 288 *  S   MXDRMM  MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR
 289 *         MATRIX A BY A VECTOR X.
 290 *  S   MXDCMD  MULTIPLICATION OF A COLUMNWISE STORED DENSE RECTANGULAR
 291 *         MATRIX A BY A VECTOR X AND ADDITION OF THE SCALED VECTOR
 292 *         ALF*Y.
 293 *  S   MXUCOP  COPYING OF A VECTOR.
 294 *  S   MXUDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
 295 *  RF  MXUDOT  DOT PRODUCT OF TWO VECTORS.
 296 *  S   MXUNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN.
 297 *  S   MXUZER  VECTOR ELEMENTS CORRESPONDING TO ACTIVE BOUNDS ARE SET
 298 *         TO ZERO.
 299 *  S   MXVCOP  COPYING OF A VECTOR.
 300 *
 301 * EXTERNAL SUBROUTINES :
 302 *  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION.
 303 *         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER
 304 *         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS
 305 *         THE VALUE OF THE OBJECTIVE FUNCTION.
 306 *  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION.
 307 *         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER
 308 *         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF)
 309 *         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
 310 *
 311 * METHOD :
 312 * HYBRID METHOD WITH SPARSE MARWIL UPDATES FOR SPARSE LEAST SQUARES
 313 * PROBLEMS.
 314 *
 315       SUBROUTINE PLIP(NF,NB,X,IX,XL,XU,GF,S,XO,GO,SO,XM,XR,GR,XMAX,
 316      & TOLX,TOLF,TOLB,TOLG,FMIN,GMAX,F,MIT,MFV,IEST,MET,MF,IPRNT,ITERM)
 317       INTEGER NF,NB,IX(*),MIT,MFV,IEST,MET,MF,IPRNT,ITERM
 318       DOUBLE PRECISION X(*),XL(*),XU(*),GF(*),S(*),XO(*),GO(*),SO(*),
 319      & XM(*),XR(*),GR(*),XMAX,TOLX,TOLF,TOLG,TOLB,FMIN,GMAX,F
 320       INTEGER ITERD,ITERS,ITERH,KD,LD,NTESX,NTESF,MTESX,MTESF,MRED,KIT,
 321      & IREST,KBF,MEC,MES,MES1,MES2,MES3,MAXST,ISYS,ITES,INITS,KTERS,
 322      & IRES1,IRES2,NRED,INEW,IOLD,I,NN,N,MFG,META,MET3
 323       DOUBLE PRECISION R,RO,RP,FO,FP,P,PO,PP,GNORM,SNORM,RMIN,RMAX,
 324      & UMAX,FMAX,DMAX,ETA0,ETA9,EPS8,EPS9,ALF1,ALF2,PAR1,PAR2,PAR,TOLD,
 325      & TOLS,TOLP
 326       DOUBLE PRECISION MXUDOT
 327       INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
 328       COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
 329       IF (ABS(IPRNT).GT.1) WRITE(6,'(1X,''ENTRY TO PLIP :'')')
 330 *
 331 *     INITIATION
 332 *
 333       KBF=0
 334       IF (NB.GT.0) KBF=2
 335       NRES=0
 336       NDEC=0
 337       NIN=0
 338       NIT=0
 339       NFV=0
 340       NFG=0
 341       NFH=0
 342       ISYS=0
 343       ITES=1
 344       MTESX=2
 345       MTESF=2
 346       INITS=2
 347       ITERM=0
 348       ITERD=0
 349       ITERS=2
 350       ITERH=0
 351       KTERS=3
 352       IREST=0
 353       IRES1=999
 354       IRES2=0
 355       MRED=10
 356       META=1
 357       MET3=4
 358       MEC=4
 359       MES=4
 360       MES1=2
 361       MES2=2
 362       MES3=2
 363       ETA0=1.0D-15
 364       ETA9=1.0D 120
 365       EPS8=1.00D 0
 366       EPS9=1.00D-8
 367       ALF1=1.0D-10
 368       ALF2=1.0D 10
 369       RMAX=ETA9
 370       DMAX=ETA9
 371       FMAX=1.0D 20
 372       IF (IEST.LE.0) FMIN=-1.0D 60
 373       IF (IEST.GT.0) IEST=1
 374       IF (XMAX.LE.0.0D 0) XMAX=1.0D 16
 375       IF (TOLX.LE.0.0D 0) TOLX=1.0D-16
 376       IF (TOLF.LE.0.0D 0) TOLF=1.0D-14
 377       IF (TOLG.LE.0.0D 0) TOLG=1.0D-6
 378       IF (TOLB.LE.0.0D 0) TOLB=FMIN+1.0D-16
 379       TOLD=1.0D-4
 380       TOLS=1.0D-4
 381       TOLP=0.9D 0
 382       IF (MET.LE.0) MET=2
 383       IF (MIT.LE.0) MIT=9000
 384       IF (MFV.LE.0) MFV=9000
 385       MFG=MFV
 386       KD= 1
 387       LD=-1
 388       KIT=-(IRES1*NF+IRES2)
 389       FO=FMIN
 390 *
 391 *     INITIAL OPERATIONS WITH SIMPLE BOUNDS
 392 *
 393       IF (KBF.GT.0) THEN
 394       DO 2 I = 1,NF
 395       IF ((IX(I).EQ.3.OR.IX(I).EQ.4) .AND. XU(I).LE.XL(I)) THEN
 396       XU(I) = XL(I)
 397       IX(I) = 5
 398       ELSE IF (IX(I).EQ.5 .OR. IX(I).EQ.6) THEN
 399       XL(I) = X(I)
 400       XU(I) = X(I)
 401       IX(I) = 5
 402       END IF
 403     2 CONTINUE
 404       CALL PCBS04(NF,X,IX,XL,XU,EPS9,KBF)
 405       CALL PYADC0(NF,N,X,IX,XL,XU,INEW)
 406       END IF
 407       IF (ITERM.NE.0) GO TO 11190
 408       CALL OBJ(NF,X,F)
 409       NFV=NFV+1
 410       CALL DOBJ(NF,X,GF)
 411       NFG=NFG+1
 412 11120 CONTINUE
 413       CALL PYTRCG(NF,NF,IX,GF,UMAX,GMAX,KBF,IOLD)
 414       IF (ABS(IPRNT).GT.1)
 415      & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X,
 416      & ''F='', G16.9,2X,''G='',E10.3)') NIT,NFV,NFG,F,GMAX
 417       CALL PYFUT1(NF,F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLB,TOLG,KD,
 418      & NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF,ITES,
 419      & IRES1,IRES2,IREST,ITERS,ITERM)
 420       IF (ITERM.NE.0) GO TO 11190
 421       IF (KBF.GT.0.AND.RMAX.GT.0.0D 0) THEN
 422       CALL PYRMC0(NF,N,IX,GF,EPS8,UMAX,GMAX,RMAX,IOLD,IREST)
 423       END IF
 424 11130 CONTINUE
 425       IF (IREST.GT.0) THEN
 426       NN=0
 427       PAR=1.0D 0
 428       LD=MIN(LD,1)
 429       IF (KIT.LT.NIT) THEN
 430         NRES=NRES+1
 431         KIT = NIT
 432       ELSE
 433         ITERM=-10
 434         IF (ITERS.LT.0) ITERM=ITERS-5
 435       END IF
 436       END IF
 437       IF (ITERM.NE.0) GO TO 11190
 438 *
 439 *     DIRECTION DETERMINATION
 440 *
 441       GNORM=SQRT(MXUDOT(NF,GF,GF,IX,KBF))
 442 *
 443 *     NEWTON LIKE STEP
 444 *
 445       CALL MXUNEG(NF,GF,S,IX,KBF)
 446       CALL MXDRMM(NF,NN,XM,S,GR)
 447       CALL MXDCMD(NF,NN,XM,GR,PAR,S,S)
 448       CALL MXUZER(NF,S,IX,KBF)
 449       ITERD=1
 450       SNORM=SQRT(MXUDOT(NF,S,S,IX,KBF))
 451 *
 452 *     TEST ON DESCENT DIRECTION AND PREPARATION OF LINE SEARCH
 453 *
 454       IF (KD.GT.0) P=MXUDOT(NF,GF,S,IX,KBF)
 455       IF (ITERD.LT.0) THEN
 456         ITERM=ITERD
 457       ELSE
 458 *
 459 *     TEST ON DESCENT DIRECTION
 460 *
 461       IF (SNORM.LE.0.0D 0) THEN
 462         IREST=MAX(IREST,1)
 463       ELSE IF (P+TOLD*GNORM*SNORM.LE.0.0D 0) THEN
 464         IREST=0
 465       ELSE
 466 *
 467 *     UNIFORM DESCENT CRITERION
 468 *
 469       IREST=MAX(IREST,1)
 470       END IF
 471       IF (IREST.EQ.0) THEN
 472 *
 473 *     PREPARATION OF LINE SEARCH
 474 *
 475         NRED = 0
 476         RMIN=ALF1*GNORM/SNORM
 477         RMAX=MIN(ALF2*GNORM/SNORM,XMAX/SNORM)
 478       END IF
 479       END IF
 480       IF (ITERM.NE.0) GO TO 11190
 481       IF (IREST.NE.0) GO TO 11130
 482       CALL PYTRCS(NF,X,IX,XO,XL,XU,GF,GO,S,RO,FP,FO,F,PO,P,RMAX,ETA9,
 483      & KBF)
 484       IF (RMAX.EQ.0.0D 0) GO TO 11175
 485 11170 CONTINUE
 486       CALL PS1L01(R,RP,F,FO,FP,P,PO,PP,FMIN,FMAX,RMIN,RMAX,
 487      & TOLS,TOLP,PAR1,PAR2,KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,
 488      & INITS,ITERS,KTERS,MES,ISYS)
 489       IF (ISYS.EQ.0) GO TO 11174
 490       CALL MXUDIR(NF,R,S,XO,X,IX,KBF)
 491       CALL PCBS04(NF,X,IX,XL,XU,EPS9,KBF)
 492       CALL OBJ(NF,X,F)
 493       NFV=NFV+1
 494       CALL DOBJ(NF,X,GF)
 495       NFG=NFG+1
 496       P=MXUDOT(NF,GF,S,IX,KBF)
 497       GO TO 11170
 498 11174 CONTINUE
 499       IF (ITERS.LE.0) THEN
 500       R=0.0D 0
 501       F=FO
 502       P=PO
 503       CALL MXVCOP(NF,XO,X)
 504       CALL MXVCOP(NF,GO,GF)
 505       IREST=MAX(IREST,1)
 506       LD=KD
 507       GO TO 11130
 508       END IF
 509       CALL MXUNEG(NF,GO,S,IX,KBF)
 510       CALL PYTRCD(NF,X,IX,XO,GF,GO,R,F,FO,P,PO,DMAX,KBF,KD,LD,ITERS)
 511       CALL MXUCOP(NF,GF,SO,IX,KBF)
 512       IF (NN.LT.MF) THEN
 513       CALL PULSP3(NF,NN,MF,XM,GR,XO,GO,R,PO,PAR,ITERH,MET3)
 514       ELSE
 515       CALL PULVP3(NF,NN,XM,XR,GR,S,SO,XO,GO,R,PO,PAR,ITERH,MEC,MET3,
 516      & MET)
 517       END IF
 518 11175 CONTINUE
 519       IF (ITERH.NE.0) IREST=MAX(IREST,1)
 520       IF (KBF.GT.0) CALL PYADC0(NF,N,X,IX,XL,XU,INEW)
 521       GO TO 11120
 522 11190 CONTINUE
 523       IF (IPRNT.GT.1.OR.IPRNT.LT.0)
 524      & WRITE(6,'(1X,''EXIT FROM PLIP :'')')
 525       IF (IPRNT.NE.0)
 526      & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X,
 527      & ''F='', G16.9,2X,''G='',E10.3,2X,''ITERM='',I3)') NIT,NFV,NFG,
 528      & F,GMAX,ITERM
 529       IF (IPRNT.LT.0)
 530      & WRITE (6,'(1X,''X='',5(G14.7,1X):/(3X,5(G14.7,1X)))')
 531      & (X(I),I=1,NF)
 532       RETURN
 533       END