src/algs/luksan/pnet.for

   1 ************************************************************************
   2 * SUBROUTINE PNETU              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)  MAXIMUM NUMBER OF GRADIENT EVALUATIONS.
  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)  CHOICE OF DIRECTION VECTORS AFTER RESTARTS.
  17 *         IPAR(5)=1-THE NEWTON DIRECTIONS ARE USED. IPAR(5)=2-THE
  18 *         STEEPEST DESCENT DIRECTIONS ARE USED.
  19 *      IPAR(6)  CHOICE OF PRECONDITIONING STRATEGY.
  20 *         IPAR(6)=1-PRECONDITIONING IS NOT USED.
  21 *         IPAR(6)=2-PRECONDITIONING BY THE LIMITED MEMORY BFGS METHOD
  22 *         IS USED.
  23 *      IPAR(7)  THE NUMBER OF LIMITED-MEMORY VARIABLE METRIC UPDATES
  24 *         IN EACH ITERATION (THEY USE 2*MF STORED VECTORS).
  25 *  RI  RPAR(9)  REAL PARAMETERS:
  26 *      RPAR(1)  MAXIMUM STEPSIZE.
  27 *      RPAR(2)  TOLERANCE FOR THE CHANGE OF VARIABLES.
  28 *      RPAR(3)  TOLERANCE FOR THE CHANGE OF FUNCTION VALUES.
  29 *      RPAR(4)  TOLERANCE FOR THE FUNCTION FALUE.
  30 *      RPAR(5)  TOLERANCE FOR THE GRADIENT NORM.
  31 *      RPAR(6)  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
  32 *      RPAR(7)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PNET.
  33 *      RPAR(8)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PNET.
  34 *      RPAR(9)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PNET.
  35 *  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
  36 *  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
  37 *  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
  38 *         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
  39 *         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
  40 *         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
  41 *         RESULTS.
  42 *  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
  43 *         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
  44 *                   MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS.
  45 *         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
  46 *                   MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS.
  47 *         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
  48 *         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
  49 *         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
  50 *                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
  51 *         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
  52 *         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
  53 *
  54 * VARIABLES IN COMMON /STAT/ (STATISTICS) :
  55 *  IO  NRES  NUMBER OF RESTARTS.
  56 *  IO  NDEC  NUMBER OF MATRIX DECOMPOSITIONS.
  57 *  IO  NIN  NUMBER OF INNER ITERATIONS.
  58 *  IO  NIT  NUMBER OF ITERATIONS.
  59 *  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
  60 *  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
  61 *  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
  62 *
  63 * SUBPROGRAMS USED :
  64 *  S   PNET  LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG
  65 *         RECURRENCES.
  66 *
  67 * EXTERNAL SUBROUTINES :
  68 *  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION.
  69 *         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER
  70 *         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS
  71 *         THE VALUE OF THE OBJECTIVE FUNCTION.
  72 *  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION.
  73 *         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER
  74 *         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF)
  75 *         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
  76 *
  77       SUBROUTINE PNETU(NF,X,IPAR,RPAR,F,GMAX,IPRNT,ITERM)
  78       INTEGER NF,IPAR(7),IPRNT,ITERM
  79       DOUBLE PRECISION X(*),RPAR(9),F,GMAX
  80       INTEGER MF,NB,LGF,LGN,LS,LXO,LGO,LXS,LGS,LXM,LGM,LU1,LU2
  81       INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
  82       COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
  83       DOUBLE PRECISION RA(:)
  84       ALLOCATABLE RA
  85       MF=IPAR(7)
  86       IF (MF.LE.0) MF=10
  87       ALLOCATE (RA(8*NF+2*NF*MF+2*MF))
  88       NB=0
  89 *
  90 *     POINTERS FOR AUXILIARY ARRAYS
  91 *
  92       LGF=1
  93       LGN=LGF+NF
  94       LS=LGN+NF
  95       LXO=LS+NF
  96       LGO=LXO+NF
  97       LXS=LGO+NF
  98       LGS=LXS+NF
  99       LXM=LGS+NF
 100       LGM=LXM+NF*MF
 101       LU1=LGM+NF*MF
 102       LU2=LU1+MF
 103       CALL PNET(NF,NB,X,IPAR,RA,RA,RA(LGF),RA(LGN),RA(LS),RA(LXO),
 104      & RA(LGO),RA(LXS),RA(LGS),RA(LXM),RA(LGM),RA(LU1),RA(LU2),RPAR(1),
 105      & RPAR(2),RPAR(3),RPAR(4),RPAR(5),RPAR(6),GMAX,F,IPAR(1),IPAR(2),
 106      & IPAR(3),IPAR(4),IPAR(5),IPAR(6),MF,IPRNT,ITERM)
 107       DEALLOCATE (RA)
 108       RETURN
 109       END
 110 ************************************************************************
 111 * SUBROUTINE PNETS              ALL SYSTEMS                   97/01/22
 112 * PURPOSE :
 113 * EASY TO USE SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION.
 114 *
 115 * PARAMETERS :
 116 *  II  NF  NUMBER OF VARIABLES.
 117 *  RI  X(NF)  VECTOR OF VARIABLES.
 118 *  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE
 119 *         X(I) IS UNBOUNDED. IX(I)=1-LOWER BOUND XL(I).LE.X(I).
 120 *         IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND
 121 *         XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED.
 122 *  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
 123 *  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
 124 *  II  IPAR(7)  INTEGER PAREMETERS:
 125 *      IPAR(1)  MAXIMUM NUMBER OF ITERATIONS.
 126 *      IPAR(2)  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
 127 *      IPAR(3)  MAXIMUM NUMBER OF GRADIENT EVALUATIONS.
 128 *      IPAR(4)  ESTIMATION INDICATOR. IPAR(4)=0-MINIMUM IS NOT
 129 *         ESTIMATED. IPAR(4)=1-MINIMUM IS ESTIMATED BY THE VALUE
 130 *         RPAR(6).
 131 *      IPAR(5)  CHOICE OF DIRECTION VECTORS AFTER RESTARTS.
 132 *         IPAR(5)=1-THE NEWTON DIRECTIONS ARE USED. IPAR(5)=2-THE
 133 *         STEEPEST DESCENT DIRECTIONS ARE USED.
 134 *      IPAR(6)  CHOICE OF PRECONDITIONING STRATEGY.
 135 *         IPAR(6)=1-PRECONDITIONING IS NOT USED.
 136 *         IPAR(6)=2-PRECONDITIONING BY THE LIMITED MEMORY BFGS METHOD
 137 *         IS USED.
 138 *      IPAR(7)  THE NUMBER OF LIMITED-MEMORY VARIABLE METRIC UPDATES
 139 *         IN EACH ITERATION (THEY USE 2*MF STORED VECTORS).
 140 *  RI  RPAR(9)  REAL PARAMETERS:
 141 *      RPAR(1)  MAXIMUM STEPSIZE.
 142 *      RPAR(2)  TOLERANCE FOR THE CHANGE OF VARIABLES.
 143 *      RPAR(3)  TOLERANCE FOR THE CHANGE OF FUNCTION VALUES.
 144 *      RPAR(4)  TOLERANCE FOR THE FUNCTION FALUE.
 145 *      RPAR(5)  TOLERANCE FOR THE GRADIENT NORM.
 146 *      RPAR(6)  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
 147 *      RPAR(7)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PNET.
 148 *      RPAR(8)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PNET.
 149 *      RPAR(9)  THIS PARAMETER IS NOT USED IN THE SUBROUTINE PNET.
 150 *  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
 151 *  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
 152 *  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
 153 *         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
 154 *         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
 155 *         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
 156 *         RESULTS.
 157 *  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
 158 *         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
 159 *                   MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS.
 160 *         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
 161 *                   MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS.
 162 *         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
 163 *         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
 164 *         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
 165 *                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
 166 *         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
 167 *         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
 168 *
 169 * VARIABLES IN COMMON /STAT/ (STATISTICS) :
 170 *  IO  NRES  NUMBER OF RESTARTS.
 171 *  IO  NDEC  NUMBER OF MATRIX DECOMPOSITIONS.
 172 *  IO  NIN  NUMBER OF INNER ITERATIONS.
 173 *  IO  NIT  NUMBER OF ITERATIONS.
 174 *  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
 175 *  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
 176 *  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
 177 *
 178 * SUBPROGRAMS USED :
 179 *  S   PNET  LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG
 180 *         RECURRENCES.
 181 *
 182 * EXTERNAL SUBROUTINES :
 183 *  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION.
 184 *         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER
 185 *         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS
 186 *         THE VALUE OF THE OBJECTIVE FUNCTION.
 187 *  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION.
 188 *         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER
 189 *         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF)
 190 *         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
 191 *
 192       SUBROUTINE PNETS(NF,X,IX,XL,XU,IPAR,RPAR,F,GMAX,IPRNT,ITERM)
 193       INTEGER NF,IX(*),IPAR(7),IPRNT,ITERM
 194       DOUBLE PRECISION X(*),XL(*),XU(*),RPAR(9),F,GMAX
 195       INTEGER MF,NB,LGF,LGN,LS,LXO,LGO,LXS,LGS,LXM,LGM,LU1,LU2
 196       INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
 197       COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
 198       DOUBLE PRECISION RA(:)
 199       ALLOCATABLE RA
 200       MF=IPAR(7)
 201       IF (MF.LE.0) MF=10
 202       ALLOCATE (RA(8*NF+2*NF*MF+2*MF))
 203       NB=1
 204 *
 205 *     POINTERS FOR AUXILIARY ARRAYS
 206 *
 207       LGF=1
 208       LGN=LGF+NF
 209       LS=LGN+NF
 210       LXO=LS+NF
 211       LGO=LXO+NF
 212       LXS=LGO+NF
 213       LGS=LXS+NF
 214       LXM=LGS+NF
 215       LGM=LXM+NF*MF
 216       LU1=LGM+NF*MF
 217       LU2=LU1+MF
 218       CALL PNET(NF,NB,X,IX,XL,XU,RA(LGF),RA(LGN),RA(LS),RA(LXO),
 219      & RA(LGO),RA(LXS),RA(LGS),RA(LXM),RA(LGM),RA(LU1),RA(LU2),RPAR(1),
 220      & RPAR(2),RPAR(3),RPAR(4),RPAR(5),RPAR(6),GMAX,F,IPAR(1),IPAR(2),
 221      & IPAR(3),IPAR(4),IPAR(5),IPAR(6),MF,IPRNT,ITERM)
 222       DEALLOCATE (RA)
 223       RETURN
 224       END
 225 ************************************************************************
 226 * SUBROUTINE PNET               ALL SYSTEMS                   01/09/22
 227 * PURPOSE :
 228 * GENERAL SUBROUTINE FOR LARGE-SCALE BOX CONSTRAINED MINIMIZATION THAT
 229 * USE THE LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG
 230 * RECURRENCES.
 231 *
 232 * PARAMETERS :
 233 *  II  NF  NUMBER OF VARIABLES.
 234 *  II  NB  CHOICE OF SIMPLE BOUNDS. NB=0-SIMPLE BOUNDS SUPPRESSED.
 235 *         NB>0-SIMPLE BOUNDS ACCEPTED.
 236 *  RI  X(NF)  VECTOR OF VARIABLES.
 237 *  II  IX(NF)  VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE
 238 *         X(I) IS UNBOUNDED. IX(I)=1-LOVER BOUND XL(I).LE.X(I).
 239 *         IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND
 240 *         XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED.
 241 *  RI  XL(NF)  VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
 242 *  RI  XU(NF)  VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
 243 *  RO  GF(NF)  GRADIENT OF THE OBJECTIVE FUNCTION.
 244 *  RA  GN(NF)  OLD GRADIENT OF THE OBJECTIVE FUNCTION.
 245 *  RO  S(NF)  DIRECTION VECTOR.
 246 *  RA  XO(NF)  ARRAY CONTAINING INCREMENTS OF VARIABLES.
 247 *  RA  GO(NF)  ARRAY CONTAINING INCREMENTS OF GRADIENTS.
 248 *  RA  XS(NF)  AUXILIARY VECTOR.
 249 *  RA  GS(NF)  AUXILIARY VECTOR.
 250 *  RA  XM(NF*MF)  ARRAY CONTAINING INCREMENTS OF VARIABLES.
 251 *  RA  GM(NF*MF)  ARRAY CONTAINING INCREMENTS OF GRADIENTS.
 252 *  RA  U1(MF)  AUXILIARY VECTOR.
 253 *  RA  U2(MF)  AUXILIARY VECTOR.
 254 *  RI  XMAX  MAXIMUM STEPSIZE.
 255 *  RI  TOLX  TOLERANCE FOR CHANGE OF VARIABLES.
 256 *  RI  TOLF  TOLERANCE FOR CHANGE OF FUNCTION VALUES.
 257 *  RI  TOLB  TOLERANCE FOR THE FUNCTION VALUE.
 258 *  RI  TOLG  TOLERANCE FOR THE GRADIENT NORM.
 259 *  RI  FMIN  ESTIMATION OF THE MINIMUM FUNCTION VALUE.
 260 *  RO  GMAX  MAXIMUM PARTIAL DERIVATIVE.
 261 *  RO  F  VALUE OF THE OBJECTIVE FUNCTION.
 262 *  II  MIT  MAXIMUM NUMBER OF ITERATIONS.
 263 *  II  MFV  MAXIMUM NUMBER OF FUNCTION EVALUATIONS.
 264 *  II  MFG  MAXIMUM NUMBER OF GRADIENT EVALUATIONS.
 265 *  II  IEST  ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED.
 266 *         IEST=1-MINIMUM IS ESTIMATED BY THE VALUE FMIN.
 267 *  II  MOS1  CHOICE OF RESTARTS AFTER A CONSTRAINT CHANGE.
 268 *         MOS1=1-RESTARTS ARE SUPPRESSED. MOS1=2-RESTARTS WITH
 269 *         STEEPEST DESCENT DIRECTIONS ARE USED.
 270 *  II  MOS1  CHOICE OF DIRECTION VECTORS AFTER RESTARTS. MOS1=1-THE
 271 *         NEWTON DIRECTIONS ARE USED. MOS1=2-THE STEEPEST DESCENT
 272 *         DIRECTIONS ARE USED.
 273 *  II  MOS2  CHOICE OF PRECONDITIONING STRATEGY. MOS2=1-PRECONDITIONING
 274 *         IS NOT USED. MOS2=2-PRECONDITIONING BY THE LIMITED MEMORY
 275 *         BFGS METHOD IS USED.
 276 *  II  MF  THE NUMBER OF LIMITED-MEMORY VARIABLE METRIC UPDATES
 277 *         IN EACH ITERATION (THEY USE 2*MF STORED VECTORS).
 278 *  II  IPRNT  PRINT SPECIFICATION. IPRNT=0-NO PRINT.
 279 *         ABS(IPRNT)=1-PRINT OF FINAL RESULTS.
 280 *         ABS(IPRNT)=2-PRINT OF FINAL RESULTS AND ITERATIONS.
 281 *         IPRNT>0-BASIC FINAL RESULTS. IPRNT<0-EXTENDED FINAL
 282 *         RESULTS.
 283 *  IO  ITERM  VARIABLE THAT INDICATES THE CAUSE OF TERMINATION.
 284 *         ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN
 285 *                   MTESX (USUALLY TWO) SUBSEQUEBT ITERATIONS.
 286 *         ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN
 287 *                   MTESF (USUALLY TWO) SUBSEQUEBT ITERATIONS.
 288 *         ITERM=3-IF F IS LESS THAN OR EQUAL TO TOLB.
 289 *         ITERM=4-IF GMAX IS LESS THAN OR EQUAL TO TOLG.
 290 *         ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED,
 291 *                   BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE.
 292 *         ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV.
 293 *         ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED.
 294 *
 295 * VARIABLES IN COMMON /STAT/ (STATISTICS) :
 296 *  IO  NRES  NUMBER OF RESTARTS.
 297 *  IO  NDEC  NUMBER OF MATRIX DECOMPOSITION.
 298 *  IO  NIN  NUMBER OF INNER ITERATIONS.
 299 *  IO  NIT  NUMBER OF ITERATIONS.
 300 *  IO  NFV  NUMBER OF FUNCTION EVALUATIONS.
 301 *  IO  NFG  NUMBER OF GRADIENT EVALUATIONS.
 302 *  IO  NFH  NUMBER OF HESSIAN EVALUATIONS.
 303 *
 304 * SUBPROGRAMS USED :
 305 *  S   PCBS04  ELIMINATION OF BOX CONSTRAINT VIOLATIONS.
 306 *  S   PS1L01  STEPSIZE SELECTION USING LINE SEARCH.
 307 *  S   PYADC0  ADDITION OF A BOX CONSTRAINT.
 308 *  S   PYFUT1  TEST ON TERMINATION.
 309 *  S   PYRMC0  DELETION OF A BOX CONSTRAINT.
 310 *  S   PYTRCD  COMPUTATION OF PROJECTED DIFFERENCES FOR THE VARIABLE METRIC
 311 *         UPDATE.
 312 *  S   PYTRCG  COMPUTATION OF THE PROJECTED GRADIENT.
 313 *  S   PYTRCS  COMPUTATION OF THE PROJECTED DIRECTION VECTOR.
 314 *  S   MXDRCB BACKWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION
 315 *         OF THE VECTOR X BY AN IMPLICIT BFGS UPDATE.
 316 *  S   MXDRCF FORWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION
 317 *         OF THE VECTOR X BY AN IMPLICIT BFGS UPDATE.
 318 *  S   MXDRSU SHIFT OF COLUMNS OF THE RECTANGULAR MATRICES A AND B.
 319 *         SHIFT OF ELEMENTS OF THE VECTOR U. THESE SHIFTS ARE USED IN
 320 *         THE LIMITED MEMORY BFGS METHOD.
 321 *  S   MXUDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
 322 *  RF  MXUDOT  DOT PRODUCT OF TWO VECTORS.
 323 *  S   MXVNEG  COPYING OF A VECTOR WITH CHANGE OF THE SIGN.
 324 *  S   MXVCOP  COPYING OF A VECTOR.
 325 *  S   MXVSCL  SCALING OF A VECTOR.
 326 *  S   MXVSET  INITIATINON OF A VECTOR.
 327 *  S   MXVDIF  DIFFERENCE OF TWO VECTORS.
 328 *
 329 * EXTERNAL SUBROUTINES :
 330 *  SE  OBJ  COMPUTATION OF THE VALUE OF THE OBJECTIVE FUNCTION.
 331 *         CALLING SEQUENCE: CALL OBJ(NF,X,FF) WHERE NF IS THE NUMBER
 332 *         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND FF IS
 333 *         THE VALUE OF THE OBJECTIVE FUNCTION.
 334 *  SE  DOBJ  COMPUTATION OF THE GRADIENT OF THE OBJECTIVE FUNCTION.
 335 *         CALLING SEQUENCE: CALL DOBJ(NF,X,GF) WHERE NF IS THE NUMBER
 336 *         OF VARIABLES, X(NF) IS THE VECTOR OF VARIABLES AND GF(NF)
 337 *         IS THE GRADIENT OF THE OBJECTIVE FUNCTION.
 338 *
 339 * METHOD :
 340 * LIMITED MEMORY VARIABLE METRIC METHOD BASED ON THE STRANG
 341 * RECURRENCES.
 342 *
 343       SUBROUTINE PNET(NF,NB,X,IX,XL,XU,GF,GN,S,XO,GO,XS,GS,XM,GM,U1,U2,
 344      & XMAX,TOLX,TOLF,TOLB,TOLG,FMIN,GMAX,F,MIT,MFV,MFG,IEST,MOS1,MOS2,
 345      & MF,IPRNT,ITERM)
 346       INTEGER NF,NB,IX(*),MIT,MFV,MFG,IEST,MOS1,MOS2,MF,IPRNT,ITERM
 347       DOUBLE PRECISION X(*),XL(*),XU(*),GF(*),GN(*),S(*),XO(*),GO(*),
 348      & XS(*),GS(*),XM(*),GM(*),U1(*),U2(*),XMAX,TOLX,TOLF,TOLG,TOLB,
 349      & FMIN,GMAX,F
 350       INTEGER ITERD,ITERS,KD,LD,NTESX,NTESF,MTESX,MTESF,MRED,KIT,
 351      & IREST,KBF,MES,MES1,MES2,MES3,MAXST,ISYS,ITES,INITS,KTERS,
 352      & IRES1,IRES2,INEW,IOLD,I,N,NRED,MX,MMX
 353       DOUBLE PRECISION R,RO,RP,FO,FP,P,PO,PP,GNORM,SNORM,RMIN,RMAX,
 354      & UMAX,FMAX,DMAX,ETA0,ETA9,EPS8,EPS9,ALF,ALF1,ALF2,RHO,RHO1,RHO2,
 355      & PAR,PAR1,PAR2,A,B,TOLD,TOLS,TOLP,EPS
 356       DOUBLE PRECISION MXUDOT
 357       INTEGER NRES,NDEC,NIN,NIT,NFV,NFG,NFH
 358       COMMON /STAT/ NRES,NDEC,NIN,NIT,NFV,NFG,NFH
 359       IF (ABS(IPRNT).GT.1) WRITE(6,'(1X,''ENTRY TO PNET :'')')
 360 *
 361 *     INITIATION
 362 *
 363       KBF=0
 364       IF (NB.GT.0) KBF=2
 365       NRES=0
 366       NDEC=0
 367       NIN=0
 368       NIT=0
 369       NFV=0
 370       NFG=0
 371       NFH=0
 372       ISYS=0
 373       ITES=1
 374       MTESX=2
 375       MTESF=2
 376       INITS=2
 377       ITERM=0
 378       ITERD=0
 379       ITERS=2
 380       KTERS=3
 381       IREST=0
 382       IRES1=999
 383       IRES2=0
 384       MRED=10
 385       MES=4
 386       MES1=2
 387       MES2=2
 388       MES3=2
 389       EPS=0.80D 0
 390       ETA0=1.0D-15
 391       ETA9=1.0D 120
 392       EPS8=1.0D 0
 393       EPS9=1.0D-8
 394       ALF1=1.0D-10
 395       ALF2=1.0D 10
 396       RMAX=ETA9
 397       DMAX=ETA9
 398       FMAX=1.0D 20
 399       IF (IEST.LE.0) FMIN=-1.0D 60
 400       IF (IEST.GT.0) IEST=1
 401       IF (XMAX.LE.0.0D 0) XMAX=1.0D 16
 402       IF (TOLX.LE.0.0D 0) TOLX=1.0D-16
 403       IF (TOLF.LE.0.0D 0) TOLF=1.0D-14
 404       IF (TOLG.LE.0.0D 0) TOLG=1.0D-6
 405       IF (TOLB.LE.0.0D 0) TOLB=FMIN+1.0D-16
 406       TOLD=1.0D-4
 407       TOLS=1.0D-4
 408       TOLP=0.9D 0
 409       IF (MIT.LE.0) MIT=5000
 410       IF (MFV.LE.0) MFV=5000
 411       IF (MFG.LE.0) MFG=30000
 412       IF (MOS1.LE.0) MOS1=1
 413       IF (MOS2.LE.0) MOS2=1
 414       KD= 1
 415       LD=-1
 416       KIT=-(IRES1*NF+IRES2)
 417       FO=FMIN
 418 *
 419 *     INITIAL OPERATIONS WITH SIMPLE BOUNDS
 420 *
 421       IF (KBF.GT.0) THEN
 422       DO 2 I = 1,NF
 423       IF ((IX(I).EQ.3.OR.IX(I).EQ.4) .AND. XU(I).LE.XL(I)) THEN
 424       XU(I) = XL(I)
 425       IX(I) = 5
 426       ELSE IF (IX(I).EQ.5 .OR. IX(I).EQ.6) THEN
 427       XL(I) = X(I)
 428       XU(I) = X(I)
 429       IX(I) = 5
 430       END IF
 431     2 CONTINUE
 432       CALL PCBS04(NF,X,IX,XL,XU,EPS9,KBF)
 433       CALL PYADC0(NF,N,X,IX,XL,XU,INEW)
 434       END IF
 435       CALL OBJ(NF,X,F)
 436       NFV=NFV+1
 437       CALL DOBJ(NF,X,GF)
 438       NFG=NFG+1
 439       LD=KD
 440 11020 CONTINUE
 441       CALL PYTRCG(NF,NF,IX,GF,UMAX,GMAX,KBF,IOLD)
 442       CALL MXVCOP(NF,GF,GN)
 443       IF (ABS(IPRNT).GT.1)
 444      & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X,
 445      & ''F='', G16.9,2X,''G='',E10.3)') NIT,NFV,NFG,F,GMAX
 446       CALL PYFUT1(NF,F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLB,TOLG,KD,NIT,KIT,
 447      & MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF,ITES,IRES1,IRES2,
 448      & IREST,ITERS,ITERM)
 449       IF (ITERM.NE.0) GO TO 11080
 450       IF (KBF.GT.0) THEN
 451       CALL PYRMC0(NF,N,IX,GN,EPS8,UMAX,GMAX,RMAX,IOLD,IREST)
 452       IF (UMAX.GT.EPS8*GMAX) IREST=MAX(IREST,1)
 453       END IF
 454       CALL MXVCOP(NF,X,XO)
 455 11040 CONTINUE
 456 *
 457 *     DIRECTION DETERMINATION
 458 *
 459       IF (IREST.NE.0) THEN
 460       IF (KIT.LT.NIT) THEN
 461       MX=0
 462       NRES=NRES+1
 463       KIT = NIT
 464       ELSE
 465       ITERM=-10
 466       IF (ITERS.LT.0) ITERM=ITERS-5
 467       GO TO 11080
 468       END IF
 469       IF (MOS1.GT.1) THEN
 470       CALL MXVNEG(NF,GN,S)
 471       GNORM=SQRT(MXUDOT(NF,GN,GN,IX,KBF))
 472       SNORM=GNORM
 473       GO TO 12560
 474       END IF
 475       END IF
 476       RHO1=MXUDOT(NF,GN,GN,IX,KBF)
 477       GNORM=SQRT(RHO1)
 478       PAR=MIN(EPS,SQRT(GNORM))
 479       IF (PAR.GT.1.0D 1*1.0D-3) THEN
 480       PAR=MIN(PAR,1.0D 0/DBLE(NIT))
 481       END IF
 482       PAR=PAR*PAR
 483 *
 484 *     CG INITIATION
 485 *
 486       RHO=RHO1
 487       SNORM=0.0D 0
 488       CALL MXVSET(NF,0.0D 0,S)
 489       CALL MXVNEG(NF,GN,GS)
 490       CALL MXVCOP(NF,GS,XS)
 491       IF (MOS2.GT.1) THEN
 492       IF (MX.EQ.0) THEN
 493       B=0.0D 0
 494       ELSE
 495       B=MXUDOT(NF,XM,GM,IX,KBF)
 496       ENDIF
 497       IF (B.GT.0.0D 0) THEN
 498       U1(1)=1.0D 0/B
 499       CALL MXDRCB(NF,MX,XM,GM,U1,U2,XS,IX,KBF)
 500       A=MXUDOT(NF,GM,GM,IX,KBF)
 501       IF (A.GT.0.0D 0) CALL MXVSCL(NF,B/A,XS,XS)
 502       CALL MXDRCF(NF,MX,XM,GM,U1,U2,XS,IX,KBF)
 503       END IF
 504       END IF
 505       RHO=MXUDOT(NF,GS,XS,IX,KBF)
 506 C      SIG=RHO
 507       MMX=NF+3
 508       NRED=0
 509 12520 CONTINUE
 510       NRED=NRED+1
 511       IF (NRED.GT.MMX) GO TO 12550
 512       FO=F
 513       PP=SQRT(ETA0/MXUDOT(NF,XS,XS,IX,KBF))
 514       LD=0
 515       CALL MXUDIR(NF,PP,XS,XO,X,IX,KBF)
 516       CALL DOBJ(NF,X,GF)
 517       NFG=NFG+1
 518       LD=KD
 519       CALL MXVDIF(NF,GF,GN,GO)
 520       F=FO
 521       CALL MXVSCL(NF,1.0D 0/PP,GO,GO)
 522       ALF=MXUDOT(NF,XS,GO,IX,KBF)
 523       IF (ALF.LE.1.0D 0/ETA9) THEN
 524 C      IF (ALF.LE.1.0D-8*SIG) THEN
 525 *
 526 *     CG FAILS (THE MATRIX IS NOT POSITIVE DEFINITE)
 527 *
 528       IF (NRED.EQ.1) THEN
 529       CALL MXVNEG(NF,GN,S)
 530       SNORM=GNORM
 531       END IF
 532       ITERD=0
 533       GO TO 12560
 534       ELSE
 535       ITERD=2
 536       END IF
 537 *
 538 *     CG STEP
 539 *
 540       ALF=RHO/ALF
 541       CALL MXUDIR(NF, ALF,XS,S,S,IX,KBF)
 542       CALL MXUDIR(NF,-ALF,GO,GS,GS,IX,KBF)
 543       RHO2=MXUDOT(NF,GS,GS,IX,KBF)
 544       SNORM=SQRT(MXUDOT(NF,S,S,IX,KBF))
 545       IF (RHO2.LE.PAR*RHO1) GO TO 12560
 546       IF (NRED.GE.MMX) GO TO 12550
 547       IF (MOS2.GT.1) THEN
 548       IF (B.GT.0.0D 0) THEN
 549       CALL MXVCOP(NF,GS,GO)
 550       CALL MXDRCB(NF,MX,XM,GM,U1,U2,GO,IX,KBF)
 551       IF (A.GT.0.0D 0) CALL MXVSCL(NF,B/A,GO,GO)
 552       CALL MXDRCF(NF,MX,XM,GM,U1,U2,GO,IX,KBF)
 553       RHO2=MXUDOT(NF,GS,GO,IX,KBF)
 554       ALF=RHO2/RHO
 555       CALL MXUDIR(NF,ALF,XS,GO,XS,IX,KBF)
 556       ELSE
 557       ALF=RHO2/RHO
 558       CALL MXUDIR(NF,ALF,XS,GS,XS,IX,KBF)
 559       END IF
 560       ELSE
 561       ALF=RHO2/RHO
 562       CALL MXUDIR(NF,ALF,XS,GS,XS,IX,KBF)
 563       END IF
 564       RHO=RHO2
 565 C      SIG=RHO2+ALF*ALF*SIG
 566       GO TO 12520
 567 12550 CONTINUE
 568 *
 569 *     AN INEXACT SOLUTION IS OBTAINED
 570 *
 571 12560 CONTINUE
 572 *
 573 *     ------------------------------
 574 *     END OF DIRECTION DETERMINATION
 575 *     ------------------------------
 576 *
 577       CALL MXVCOP(NF,XO,X)
 578       CALL MXVCOP(NF,GN,GF)
 579       IF (KD.GT.0) P=MXUDOT(NF,GN,S,IX,KBF)
 580       IF (ITERD.LT.0) THEN
 581         ITERM=ITERD
 582       ELSE
 583 *
 584 *     TEST ON DESCENT DIRECTION
 585 *
 586       IF (SNORM.LE.0.0D 0) THEN
 587         IREST=MAX(IREST,1)
 588       ELSE IF (P+TOLD*GNORM*SNORM.LE.0.0D 0) THEN
 589         IREST=0
 590       ELSE
 591 *
 592 *     UNIFORM DESCENT CRITERION
 593 *
 594       IREST=MAX(IREST,1)
 595       END IF
 596       IF (IREST.EQ.0) THEN
 597 *
 598 *     PREPARATION OF LINE SEARCH
 599 *
 600         NRED = 0
 601         RMIN=ALF1*GNORM/SNORM
 602         RMAX=MIN(ALF2*GNORM/SNORM,XMAX/SNORM)
 603       END IF
 604       END IF
 605       LD=KD
 606       IF (ITERM.NE.0) GO TO 11080
 607       IF (IREST.NE.0) GO TO 11040
 608       CALL PYTRCS(NF,X,IX,XO,XL,XU,GF,GO,S,RO,FP,FO,F,PO,P,RMAX,ETA9,
 609      & KBF)
 610       IF (RMAX.EQ.0.0D 0) GO TO 11075
 611 11060 CONTINUE
 612       CALL PS1L01(R,RP,F,FO,FP,P,PO,PP,FMIN,FMAX,RMIN,RMAX,TOLS,TOLP,
 613      & PAR1,PAR2,KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS,
 614      & MES,ISYS)
 615       IF (ISYS.EQ.0) GO TO 11064
 616       CALL MXUDIR(NF,R,S,XO,X,IX,KBF)
 617       CALL PCBS04(NF,X,IX,XL,XU,EPS9,KBF)
 618       CALL OBJ(NF,X,F)
 619       NFV=NFV+1
 620       CALL DOBJ(NF,X,GF)
 621       NFG=NFG+1
 622       LD=KD
 623       P=MXUDOT(NF,GF,S,IX,KBF)
 624       GO TO 11060
 625 11064 CONTINUE
 626       IF (ITERS.LE.0) THEN
 627       R=0.0D 0
 628       F=FO
 629       P=PO
 630       CALL MXVCOP(NF,XO,X)
 631       CALL MXVCOP(NF,GO,GF)
 632       IREST=MAX(IREST,1)
 633       LD=KD
 634       GO TO 11040
 635       END IF
 636       CALL PYTRCD(NF,X,IX,XO,GF,GO,R,F,FO,P,PO,DMAX,KBF,KD,LD,ITERS)
 637       IF (MOS2.GT.1) THEN
 638       MX=MIN(MX+1,MF)
 639       CALL MXDRSU(NF,MX,XM,GM,U1)
 640       CALL MXVCOP(NF,XO,XM)
 641       CALL MXVCOP(NF,GO,GM)
 642       END IF
 643 11075 CONTINUE
 644       IF (KBF.GT.0) THEN
 645       CALL PYADC0(NF,N,X,IX,XL,XU,INEW)
 646       IF (INEW.GT.0) IREST=MAX(IREST,1)
 647       END IF
 648       GO TO 11020
 649 11080 CONTINUE
 650       IF (IPRNT.GT.1.OR.IPRNT.LT.0)
 651      & WRITE(6,'(1X,''EXIT FROM PNET :'')')
 652       IF (IPRNT.NE.0)
 653      & WRITE (6,'(1X,''NIT='',I5,2X,''NFV='',I5,2X,''NFG='',I5,2X,
 654      & ''F='', G16.9,2X,''G='',E10.3,2X,''ITERM='',I3)') NIT,NFV,NFG,
 655      & F,GMAX,ITERM
 656       IF (IPRNT.LT.0)
 657      & WRITE (6,'(1X,''X='',5(G14.7,1X):/(3X,5(G14.7,1X)))')
 658      & (X(I),I=1,NF)
 659       RETURN
 660       END