1* SUBROUTINE PA0GS3 ALL SYSTEMS 91/12/01 2* PURPOSE : 3* NUMERICAL COMPUTATION OF THE GRADIENT OF THE APPROXIMATED 4* FUNCTION. 5* 6* PARAMETERS : 7* II N NUMBER OF VARIABLES. 8* II KA INDEX OF THE APPROXIMATED FUNCTION. 9* RI X(N) VECTOR OF VARIABLES. 10* RO FA VALUE OF THE APPROXIMATED FUNCTION. 11* RA GA(N) GRADIENT OF THE APPROXIMATED FUNCTION. 12* II IAG(N+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 13* II JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 14* RI ETA1 PRECISION OF THE COMPUTED FUNCTION VALUES. 15* IU NAV NUMBER OF APPROXIMATED FUNCTION EVALUATIONS. 16* 17* SUBPROGRAMS USED : 18* SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION. 19* 20 SUBROUTINE PA0GS3(N,KA,X,FA,GA,IAG,JAG,ETA1,NAV) 21 DOUBLE PRECISION ETA1,FA 22 INTEGER KA,N,NAV 23 DOUBLE PRECISION GA(*),X(*) 24 INTEGER IAG(*),JAG(*) 25 DOUBLE PRECISION ETA,FTEMP,XSTEP,XTEMP 26 INTEGER IVAR,KVAR 27 ETA = SQRT(ETA1) 28 FTEMP = FA 29 DO 10 KVAR = IAG(KA),IAG(KA+1) - 1 30 IVAR = JAG(KVAR) 31* 32* STEP SELECTION 33* 34 XSTEP = ETA*MAX(ABS(X(IVAR)),1.0D0)*SIGN(1.0D0,X(IVAR)) 35 XTEMP = X(IVAR) 36 X(IVAR) = X(IVAR) + XSTEP 37 XSTEP = X(IVAR) - XTEMP 38 NAV = NAV + 1 39 CALL FUN(N,KA,X,FA) 40* 41* NUMERICAL DIFFERENTIATION 42* 43 GA(IVAR) = (FA-FTEMP)/XSTEP 44 X(IVAR) = XTEMP 45 10 CONTINUE 46 FA = FTEMP 47 RETURN 48 END 49* SUBROUTINE PA0HS3 ALL SYSTEMS 99/12/01 50* PURPOSE : 51* NUMERICAL COMPUTATION OF THE HESSIAN MATRIX OF THE APPROXIMATED 52* FUNCTION USING ITS VALUES. 53* 54* PARAMETERS : 55* II NF NUMBER OF VARIABLES. 56* II KA INDEX OF THE SELECTED FUNCTION. 57* RI X(NF) VECTOR OF VARIABLES. 58* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. 59* RO HA(M) HESSIAN MATRIX OF THE APPROXIMATED FUNCTION. 60* RA GO(NF) AUXILIARY VECTOR. 61* RA GS(NF) AUXILIARY VECTOR. 62* RI IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 63* RI JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 64* RI FA VALUE OF THE SELECTED FUNCTION. 65* RI ETA1 PRECISION OF THE COMPUTED VALUES. 66* II KBF TYPE OF BOUNDS. KBF=0-BOUNDS ARE NOT USED. KBF=1-ONE SIDED 67* BOUNDS. KBF=1-TWO SIDED BOUNDS. 68* IO NAV NUMBER OF APPROXIMATED FUNTION VALUES. 69* 70* SUBPROGRAMS USED : 71* SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION. 72* 73 SUBROUTINE PA0HS3(NF,KA,X,IX,HA,GO,GS,IAG,JAG,FA,ETA1,KBF,NAV) 74 INTEGER NF,KA,IX(*),IAG(*),JAG(*),KBF,NAV 75 DOUBLE PRECISION X(*),HA(*),GO(*),GS(*),FA,ETA1 76 DOUBLE PRECISION XTEMPI,XTEMPJ,FTEMP,ETA 77 INTEGER I,J,IJ 78 INTEGER IVAR,JVAR,KVAR,LVAR,MVAR 79 ETA=ETA1**(1.0D 0/3.0D 0) 80 FTEMP=FA 81 MVAR=IAG(KA)-1 82 DO 4 KVAR=MVAR+1,IAG(KA+1)-1 83 IVAR=ABS(JAG(KVAR)) 84 IF (KBF.GT.0) THEN 85 IF (IX(IVAR).LE.-5) GO TO 4 86 END IF 87* 88* STEP SELECTION 89* 90 XTEMPI=X(IVAR) 91 IF (XTEMPI.GE.0.0D 0) THEN 92 GO(IVAR)= ETA*MAX(ABS(XTEMPI),1.0D 0) 93 ELSE 94 GO(IVAR)=-ETA*MAX(ABS(XTEMPI),1.0D 0) 95 END IF 96 X(IVAR)=X(IVAR)+GO(IVAR) 97 GO(IVAR)=X(IVAR)-XTEMPI 98 CALL FUN(NF,KA,X,FA) 99 NAV=NAV+1 100 GS(IVAR)=FA 101 X(IVAR)=XTEMPI 102 4 CONTINUE 103* 104* NUMERICAL DIFFERENTIATION 105* 106 DO 10 KVAR=MVAR+1,IAG(KA+1)-1 107 IVAR=ABS(JAG(KVAR)) 108 IF (KBF.GT.0) THEN 109 IF (IX(IVAR).LE.-5) GO TO 10 110 END IF 111 XTEMPI=X(IVAR) 112 X(IVAR)=XTEMPI+GO(IVAR) 113 DO 9 LVAR=KVAR,IAG(KA+1)-1 114 JVAR=ABS(JAG(LVAR)) 115 IF (KBF.GT.0) THEN 116 IF (IX(JVAR).LE.-5) GO TO 9 117 END IF 118 XTEMPJ=X(JVAR) 119 X(JVAR)=X(JVAR)+GO(JVAR) 120 CALL FUN(NF,KA,X,FA) 121 NAV=NAV+1 122 I=KVAR-MVAR 123 J=LVAR-MVAR 124 IJ=MAX(I,J)*(MAX(I,J)-1)/2+MIN(I,J) 125 HA(IJ)=((FTEMP-GS(IVAR))+(FA-GS(JVAR)))/(GO(IVAR)*GO(JVAR)) 126 X(JVAR)=XTEMPJ 127 9 CONTINUE 128 X(IVAR)=XTEMPI 129 10 CONTINUE 130 FA=FTEMP 131 RETURN 132 END 133* SUBROUTINE PA0SQ3 ALL SYSTEMS 92/12/01 134* PURPOSE : 135* COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE FUNCTION 136* WHICH IS DEFINED AS A SUM OF SQUARES. 137* 138* PARAMETERS: 139* II N NUMBER OF VARIABLES. 140* RI X(N) VECTOR OF VARIABLES. 141* RO F VALUE OF THE OBJECTIVE FUNCTION. 142* RO AF(N) VALUES OF THE APPROXIMATED FUNCTIONS. 143* RA GA(N) GRADIENT OF THE APPROXIMATED FUNCTION. 144* RI AG(IAG(N+1)-1) SPARSE RECTANGULAR MATRIX WHICH IS USED FOR THE 145* DIRECTION VECTOR DETERMINATION. 146* II IAG(N+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 147* II JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 148* RI G(N) GRADIENT OF THE OBJECTIVE FUNCTION. 149* RI ETA1 PRECISION OF THE COMPUTED FUNCTION VALUES. 150* II KD DEGREE OF REQUIRED DERIVATIVES. 151* IU LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. 152* IO NFV NUMBER OF FUNCTION EVALUATIONS. 153* IO NFG NUMBER OF GRADIENT EVALUATIONS. 154* II IDER DEGREE OF ANALYTICALLY COMPUTED DERIVATIVES (0 OR 1). 155* 156* SUBPROGRAMS USED : 157* SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION. 158* S PA0GS3 NUMERICAL DIFFERENTIATION. 159* S MXVSET INITIATION OF A VECTOR. 160* 161 SUBROUTINE PA0SQ3(N,X,F,AF,GA,AG,IAG,JAG,G,ETA1,KD,LD,NFV,NFG, 162 & IDER) 163 DOUBLE PRECISION ETA1,F 164 INTEGER IDER,KD,LD,N,NFV,NFG 165 DOUBLE PRECISION AF(*),AG(*),G(*),GA(*),X(*) 166 INTEGER IAG(*),JAG(*) 167 DOUBLE PRECISION FA 168 INTEGER J,JP,K,KA,L,NAV 169 IF (KD.LE.LD) RETURN 170 IF (KD.GE.0 .AND. LD.LT.0) THEN 171 F = 0.0D0 172 NFV=NFV+1 173 END IF 174 IF (KD.GE.1 .AND. LD.LT.1) THEN 175 CALL MXVSET(N,0.0D0,G) 176 IF (IDER.GT.0) NFG=NFG+1 177 END IF 178 NAV=0 179 DO 30 KA = 1,N 180 IF (KD.LT.0) GO TO 30 181 IF (LD.GE.0) THEN 182 FA = AF(KA) 183 ELSE 184 CALL FUN(N,KA,X,FA) 185 AF(KA) = FA 186 END IF 187 IF (LD.GE.0) GO TO 10 188 F = F + FA*FA 189 10 IF (KD.LT.1) GO TO 30 190 IF (IDER.EQ.0) THEN 191 CALL PA0GS3(N,KA,X,FA,GA,IAG,JAG,ETA1,NAV) 192 ELSE 193 CALL DFUN(N,KA,X,GA) 194 END IF 195 K = IAG(KA) 196 L = IAG(KA+1) - K 197 DO 20 J = 1,L 198 JP = JAG(K) 199 G(JP) = G(JP) + FA*GA(JP) 200 AG(K) = GA(JP) 201 K = K + 1 202 20 CONTINUE 203 30 CONTINUE 204 IF (KD.GE.0 .AND. LD.LT.0) F = 0.5D0*F 205 IF (IDER.EQ.0) NFV=NFV+NAV/N 206 LD = KD 207 RETURN 208 END 209* SUBROUTINE PA1HS3 ALL SYSTEMS 99/12/01 210* PURPOSE : 211* NUMERICAL COMPUTATION OF THE HESSIAN MATRIX OF THE APPROXIMATED 212* FUNCTION USING ITS GRADIENTS. 213* 214* PARAMETERS : 215* II NF NUMBER OF VARIABLES. 216* II KA INDEX OF THE SELECTED FUNCTION. 217* RI X(NF) VECTOR OF VARIABLES. 218* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. 219* RO HA(M) HESSIAN MATRIX OF THE APPROXIMATED FUNCTION. 220* RI GA(NF) GRADIENT OF THE APPROXIMATED FUNCTION. 221* RA GO(NF) AUXILIARY VECTOR. 222* RI IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 223* RI JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 224* RI FA VALUE OF THE SELECTED FUNCTION. 225* RI ETA1 PRECISION OF THE COMPUTED VALUES. 226* II KBF TYPE OF BOUNDS. KBF=0-BOUNDS ARE NOT USED. KBF=1-ONE SIDED 227* BOUNDS. KBF=2-TWO SIDED BOUNDS. 228* IO NAG NUMBER OF APPROXIMATED FUNTION GRADIENTS. 229* 230* SUBPROGRAMS USED : 231* SE DFUN COMPUTATION OF THE GRADIENT OF THE APPROXIMATED FUNCTION. 232* 233 SUBROUTINE PA1HS3(NF,KA,X,IX,HA,GA,GO,IAG,JAG,FA,ETA1,KBF,NAG) 234 INTEGER NF,KA,IX(*),IAG(*),JAG(*),KBF,NAG 235 DOUBLE PRECISION X(*),HA(*),GA(*),GO(*),FA,ETA1 236 DOUBLE PRECISION XSTEP,XTEMP,FTEMP,ETA 237 INTEGER I,J,IJ 238 INTEGER IVAR,JVAR,KVAR,LVAR,MVAR 239 ETA=SQRT(ETA1) 240 FTEMP=FA 241 MVAR=IAG(KA)-1 242 DO 5 KVAR=MVAR+1,IAG(KA+1)-1 243 IVAR=ABS(JAG(KVAR)) 244 IF (KBF.GT.0) THEN 245 IF (IX(IVAR).LE.-5) GO TO 5 246 END IF 247* 248* STEP SELECTION 249* 250 XTEMP=X(IVAR) 251 IF (XTEMP.GE.0.0D 0) THEN 252 XSTEP= ETA*MAX(ABS(XTEMP),1.0D 0) 253 ELSE 254 XSTEP=-ETA*MAX(ABS(XTEMP),1.0D 0) 255 END IF 256 X(IVAR)=XTEMP+XSTEP 257 XSTEP=X(IVAR)-XTEMP 258 CALL DFUN(NF,KA,X,GA) 259 NAG=NAG+1 260* 261* NUMERICAL DIFFERENTIATION 262* 263 DO 4 LVAR=MVAR+1,IAG(KA+1)-1 264 JVAR=ABS(JAG(LVAR)) 265 IF (KBF.GT.0) THEN 266 IF (IX(JVAR).LE.-5) GO TO 4 267 END IF 268 I=KVAR-MVAR 269 J=LVAR-MVAR 270 IJ=MAX(I,J)*(MAX(I,J)-1)/2+MIN(I,J) 271 IF (LVAR .GE. KVAR) THEN 272 HA(IJ)=(GA(JVAR)-GO(JVAR))/XSTEP 273 ELSE 274 HA(IJ)=0.5D 0*(HA(IJ)+(GA(JVAR)-GO(JVAR))/XSTEP) 275 END IF 276 4 CONTINUE 277 X(IVAR)=XTEMP 278 5 CONTINUE 279 FA=FTEMP 280 RETURN 281 END 282* SUBROUTINE PA1SF3 ALL SYSTEMS 97/12/01 283* PURPOSE : 284* COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE FUNCTION 285* WHICH IS DEFINED AS A SUM OF SQUARES. 286* 287* PARAMETERS: 288* II NF NUMBER OF VARIABLES. 289* II NA NUMBER OF APPROXIMATED FUNCTIONS. 290* RI X(NF) VECTOR OF VARIABLES. 291* RU GA(NF) GRADIENT OF THE APPROXIMATED FUNCTION. 292* RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 293* RO AG(MA) SPARSE JACOBIAN MATRIX. 294* RI IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 295* RI JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 296* RO F VALUE OF THE OBJECTIVE FUNCTION. 297* RO AF(NA) VECTOR CONTAINING VALUES OF THE APPROXIMATED 298* FUNCTIONS. 299* II KD DEGREE OF REQUIRED DERIVATIVES. 300* IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. 301* II ISNA SAVING SPECIFICATION. ISNA=0-NO SAVING. ISNA=1-FUNCTION 302* VALUES ARE SAVED. ISNA=2-FUNCTION VALUES AND GRADIENTS ARE 303* SAVED. 304* IU NFV NUMBER OF OBJECTIVE FUNCTION VALUES COMPUTED. 305* IU NFG NUMBER OF OBJECTIVE FUNCTION GRADIENTS COMPUTED. 306* 307* SUBPROGRAMS USED : 308* SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION. 309* SE DFUN COMPUTATION OF THE GRADIENT OF THE APPROXIMATED FUNCTION. 310* S MXVSET INITIATION OF A VECTOR. 311* 312 SUBROUTINE PA1SF3(NF,NA,X,GA,G,AG,IAG,JAG,F,AF,KD,LD,ISNA, 313 & NFV,NFG) 314 INTEGER NF,NA,IAG(*),JAG(*),KD,LD,ISNA,NFV,NFG 315 DOUBLE PRECISION X(*),GA(*),G(*),AG(*),F,AF(*) 316 INTEGER J,JP,K,L,KA 317 DOUBLE PRECISION FA 318 IF (KD.LE.LD) RETURN 319 IF (KD.GE.0.AND.LD.LT.0) THEN 320 F=0.0D 0 321 NFV=NFV+1 322 END IF 323 IF (KD.GE.1.AND.LD.LT.1) THEN 324 CALL MXVSET(NF,0.0D 0,G) 325 NFG=NFG+1 326 END IF 327 DO 5 KA=1,NA 328 IF (KD.LT.0) GO TO 5 329 IF (LD.LT.0) THEN 330 CALL FUN(NF,KA,X,FA) 331 F=F+FA 332 AF(KA)=FA 333 ELSE 334 FA=AF(KA) 335 END IF 336 IF (KD.LT.1) GO TO 5 337 IF (LD.LT.1) THEN 338 CALL DFUN(NF,KA,X,GA) 339 K=IAG(KA) 340 L=IAG(KA+1)-K 341 DO 4 J=1,L 342 JP=ABS(JAG(K)) 343 G(JP)=G(JP)+GA(JP) 344 IF (ISNA.GT.1) AG(K)=GA(JP) 345 K=K+1 346 4 CONTINUE 347 END IF 348 5 CONTINUE 349 LD=KD 350 RETURN 351 END 352* SUBROUTINE PA2SF4 ALL SYSTEMS 97/12/01 353* PURPOSE : 354* COMPUTATION OF THE VALUE AND THE GRADIENT AND THE HESSIAN MATRIX 355* OF THE OBJECTIVE FUNCTION WHICH IS DEFINED AS A SUM OF SQUARES. 356* 357* PARAMETERS: 358* II NF NUMBER OF VARIABLES. 359* II NA NUMBER OF APPROXIMATED FUNCTIONS. 360* RI X(NF) VECTOR OF VARIABLES. 361* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. 362* RU GA(NF) GRADIENT OF THE APPROXIMATED FUNCTION. 363* RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 364* RA GO(NF) AUXILIARY VECTOR. 365* RU HA(MB) HESSIAN MATRIX OF THE APPROXIMATED FUNCTION. 366* RO H(M) SPARSE HESSIAN MATRIX OF THE OBJECTIVE FUNCTION. 367* II IH(NF+1) POINTERS OF THE DIAGONAL ELEMENTS OF H. 368* II JH(M) INDICES OF THE NONZERO ELEMENTS OF H. 369* RI IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 370* RI JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 371* RO AF(NA) VECTOR CONTAINING VALUES OF THE APPROXIMATED 372* FUNCTIONS. 373* RO F VALUE OF THE OBJECTIVE FUNCTION. 374* RI ETA1 PRECISION OF THE COMPUTED FUNCTION VALUES. 375* II KBF TYPE OF BOUNDS. KBF=0-BOUNDS ARE NOT USED. KBF=1-ONE SIDED 376* BOUNDS. KBF=2-TWO SIDED BOUNDS. 377* II KD DEGREE OF REQUIRED DERVATIVES. 378* IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. 379* IU NFV NUMBER OF OBJECTIVE FUNCTION VALUES COMPUTED. 380* IU NFG NUMBER OF OBJECTIVE FUNCTION GRADIENTS COMPUTED. 381* IU IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. 382* 383* SUBPROGRAMS USED : 384* SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION. 385* SE DFUN COMPUTATION OF THE GRADIENT OF THE APPROXIMATED FUNCTION. 386* S MXVSET INITIATION OF A VECTOR. 387* S PA1HS3 NUMERICAL COMPUTATION OF THE PARTIAL HESSIAN MATRIX. 388* S PASSH2 ADDITION OF THE PARTIAL HESSIAN MATRIX TO THE SPARSE 389* HESSIAN MATRIX. 390* 391 SUBROUTINE PA2SF4(NF,NA,X,IX,GA,G,GO,HA,H,IH,JH,IAG,JAG,AF,F, 392 & ETA1,KBF,KD,LD,NFV,NFG,IDECF) 393 INTEGER NF,NA,IX(*),IH(*),JH(*),IAG(*),JAG(*),KBF,KD,LD,NFV,NFG, 394 & IDECF 395 DOUBLE PRECISION X(*),GA(*),G(*),GO(*),HA(*),H(*),AF(*),F,ETA1 396 DOUBLE PRECISION FA 397 INTEGER J,JP,K,KA,L,NAG 398 IF (KD.LE.LD) RETURN 399 IF (KD.GE.0.AND.LD.LT.0) THEN 400 F=0.0D 0 401 NFV=NFV+1 402 END IF 403 IF (KD.GE.1.AND.LD.LT.1) THEN 404 CALL MXVSET(NF,0.0D 0,G) 405 NFG=NFG+1 406 END IF 407 IF (KD.GE.2.AND.LD.LT.2) CALL MXVSET(IH(NF+1)-1,0.0D 0,H) 408 NAG=0 409 DO 9 KA=1,NA 410 IF (KD.LT.0) GO TO 9 411 IF (LD.LT.0) THEN 412 CALL FUN(NF,KA,X,FA) 413 F=F+FA 414 AF(KA)=FA 415 ELSE 416 FA=AF(KA) 417 END IF 418 IF (KD.LT.1) GO TO 9 419 CALL DFUN(NF,KA,X,GA) 420 IF (LD.LT.1) THEN 421 K=IAG(KA) 422 L=IAG(KA+1)-K 423 DO 1 J=1,L 424 JP=ABS(JAG(K)) 425 G(JP)=G(JP)+GA(JP) 426 K=K+1 427 1 CONTINUE 428 END IF 429 IF (KD.LT.2) GO TO 9 430 IDECF=0 431 CALL PA1HS3(NF,KA,X,IX,HA,GO,GA,IAG,JAG,FA,ETA1,KBF,NAG) 432 CALL PASSH2(H,IH,JH,HA,IAG,JAG,KA,1.0D 0) 433 9 CONTINUE 434 NFG=NFG+NAG/NA 435 LD=KD 436 RETURN 437 END 438* SUBROUTINE PA2SQ4 ALL SYSTEMS 97/12/01 439* PURPOSE : 440* COMPUTATION OF THE VALUE AND THE GRADIENT AND THE HESSIAN MATRIX 441* OF THE OBJECTIVE FUNCTION WHICH IS DEFINED AS A SUM OF SQUARES. 442* 443* PARAMETERS: 444* II NF NUMBER OF VARIABLES. 445* II NA NUMBER OF APPROXIMATED FUNCTIONS. 446* RI X(NF) VECTOR OF VARIABLES. 447* RU GA(NF) GRADIENT OF THE APPROXIMATED FUNCTION. 448* RO AG(MA) SPARSE JACOBIAN MATRIX. 449* RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 450* RO H(M) SPARSE HESSIAN MATRIX OF THE OBJECTIVE FUNCTION. 451* II IH(NF+1) POINTERS OF THE DIAGONAL ELEMENTS OF H. 452* II JH(M) INDICES OF THE NONZERO ELEMENTS OF H. 453* RI IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 454* RI JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 455* RI AF(NA) VECTOR CONTAINING VALUES OF THE APPROXIMATED 456* FUNCTIONS. 457* RI F VALUE OF THE OBJECTIVE FUNCTION. 458* RI ETA1 PRECISION OF THE COMPUTED FUNCTION VALUES. 459* II KD DEGREE OF REQUIRED DERIVATIVES. 460* IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. 461* II ISNA SAVING SPECIFICATION. ISNA=0-NO SAVING. ISNA=1-FUNCTION 462* VALUES ARE SAVED. ISNA=2-FUNCTION VALUES AND GRADIENTS ARE 463* SAVED. 464* IU NFV NUMBER OF OBJECTIVE FUNCTION VALUES COMPUTED. 465* IU NFG NUMBER OF OBJECTIVE FUNCTION GRADIENTS COMPUTED. 466* II IDER DEGREE OF ANALYTICALLY COMPUTED DERIVATIVES (0 OR 1). 467* IU IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. 468* 469* SUBPROGRAMS USED : 470* SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION. 471* SE DFUN COMPUTATION OF THE GRADIENT OF THE APPROXIMATED FUNCTION. 472* S MXVSET INITIATION OF A VECTOR. 473* S PASSH1 ADDITION OF THE PARTIAL GAUSS-NEWTON TERM TO THE SPARSE 474* HESSIAN MATRIX. 475* 476 SUBROUTINE PA2SQ4(NF,NA,X,GA,AG,G,H,IH,JH,IAG,JAG,AF,F,ETA1,KD, 477 & LD,ISNA,NFV,NFG,IDER,IDECF) 478 INTEGER NF,NA,IH(*),JH(*),IAG(*),JAG(*),KD,LD,ISNA,NFV,NFG,IDER, 479 & IDECF 480 DOUBLE PRECISION X(*),GA(*),AG(*),G(*),H(*),AF(*),F,ETA1 481 INTEGER J,JP,K,KA,L,NAV 482 DOUBLE PRECISION FA 483 IF (KD.LE.LD) RETURN 484 IF (KD.GE.0.AND.LD.LT.0) THEN 485 F=0.0D 0 486 NFV=NFV+1 487 END IF 488 IF (KD.GE.1.AND.LD.LT.1) THEN 489 CALL MXVSET(NF,0.0D 0,G) 490 IF (IDER.GT.0) NFG=NFG+1 491 END IF 492 IF (KD.GE.2.AND.LD.LT.2) CALL MXVSET(IH(NF+1)-1,0.0D 0,H) 493 NAV=0 494 DO 3 KA=1,NA 495 IF (KD.LT.0) GO TO 3 496 IF (LD.LT.0) THEN 497 CALL FUN(NF,KA,X,FA) 498 F=F+FA*FA 499 AF(KA)=FA 500 ELSE 501 FA=AF(KA) 502 END IF 503 IF (KD.LT.1) GO TO 3 504 IF (IDER.EQ.0) THEN 505 CALL PA0GS3(NF,KA,X,FA,GA,IAG,JAG,ETA1,NAV) 506 ELSE 507 CALL DFUN(NF,KA,X,GA) 508 END IF 509 IF (LD.GE.1) GO TO 2 510 K=IAG(KA) 511 L=IAG(KA+1)-K 512 DO 1 J=1,L 513 JP=ABS(JAG(K)) 514 G(JP)=G(JP)+FA*GA(JP) 515 IF (ISNA.GT.1) AG(K)=GA(JP) 516 K=K+1 517 1 CONTINUE 518 2 IF (KD.LT.2) GO TO 3 519 IDECF=0 520 CALL PASSH1(H,IH,JH,IAG,JAG,GA,KA,1.0D 0) 521 3 CONTINUE 522 IF (KD.GE.0.AND.LD.LT.0) F=5.0D-1*F 523 IF (IDER.EQ.0) NFV=NFV+NAV/NA 524 LD=KD 525 RETURN 526 END 527* SUBROUTINE PA2SQ8 ALL SYSTEMS 97/12/01 528* PURPOSE : 529* COMPUTATION OF THE VALUE AND THE GRADIENT AND THE HESSIAN MATRIX 530* OF THE OBJECTIVE FUNCTION WHICH IS DEFINED AS A SUM OF SQUARES. 531* 532* PARAMETERS: 533* II NF NUMBER OF VARIABLES. 534* II NA NUMBER OF APPROXIMATED FUNCTIONS. 535* RI X(NF) VECTOR OF VARIABLES. 536* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. 537* RU GA(NF) GRADIENT OF THE APPROXIMATED FUNCTION. 538* RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 539* RA GO(NF) AUXILIARY VECTOR. 540* RA GS(NF) AUXILIARY VECTOR. 541* RU HA(ME) HESSIAN MATRIX OF THE APPROXIMATED FUNCTION. 542* RO H(M) SPARSE HESSIAN MATRIX OF THE OBJECTIVE FUNCTION. 543* II IH(NF+1) POINTERS OF THE DIAGONAL ELEMENTS OF H. 544* II JH(M) INDICES OF THE NONZERO ELEMENTS OF H. 545* RI IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 546* RI JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 547* RO AF(NA) VECTOR CONTAINING VALUES OF THE APPROXIMATED 548* FUNCTIONS. 549* RO F VALUE OF THE OBJECTIVE FUNCTION. 550* RI ETA1 PRECISION OF THE COMPUTED FUNCTION VALUES. 551* II KBF TYPE OF BOUNDS. KBF=0-BOUNDS ARE NOT USED. KBF=1-ONE SIDED 552* BOUNDS. KBF=2-TWO SIDED BOUNDS. 553* II KD DEGREE OF REQUIRED DERIVATIVES. 554* IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. 555* IU NFV NUMBER OF OBJECTIVE FUNCTION VALUES COMPUTED. 556* IU NFG NUMBER OF OBJECTIVE FUNCTION GRADIENTS COMPUTED. 557* II IPOM1 CORRECTION OPTION. IPOM1=0-THE NEWTON CORRECTION IS USED. 558* IPOM1=1-CORRECTION IS NOT USED. 559* II IDER DEGREE OF ANALYTICALLY COMPUTED DERIVATIVES (0 OR 1). 560* IU IDECF DECOMPOSITION INDICATOR. IDECF=0-NO DECOMPOSITION. 561* 562* SUBPROGRAMS USED : 563* SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION. 564* SE DFUN COMPUTATION OF THE GRADIENT OF THE APPROXIMATED FUNCTION. 565* S MXVSET INITIATION OF A VECTOR. 566* S PA0HS3 NUMERICAL COMPUTATION OF THE PARTIAL HESSIAN MATRIX. 567* S PA1HS3 NUMERICAL COMPUTATION OF THE PARTIAL HESSIAN MATRIX. 568* S PASSH1 ADDITION OF THE PARTIAL GAUSS-NEWTON TERM TO THE SPARSE 569* HESSIAN MATRIX. 570* S PASSH2 ADDITION OF THE PARTIAL HESSIAN MATRIX TO THE SPARSE 571* HESSIAN MATRIX. 572* 573 SUBROUTINE PA2SQ8(NF,NA,X,IX,GA,G,GO,GS,HA,H,IH,JH,IAG,JAG,AF,F, 574 & ETA1,KBF,KD,LD,NFV,NFG,IPOM1,IDER,IDECF) 575 INTEGER NF,NA,IX(*),IH(*),JH(*),IAG(*),JAG(*),KBF,KD,LD,NFV,NFG, 576 & IPOM1,IDER,IDECF 577 DOUBLE PRECISION X(*),GA(*),G(*),GO(*),GS(*),HA(*),H(*),AF(*),F, 578 & ETA1 579 INTEGER J,JP,K,KA,L,NAV,NAG 580 DOUBLE PRECISION FA 581 IF (KD.LE.LD) RETURN 582 IF (KD.GE.0.AND.LD.LT.0) THEN 583 F=0.0D 0 584 NFV=NFV+1 585 END IF 586 IF (KD.GE.1.AND.LD.LT.1) THEN 587 CALL MXVSET(NF,0.0D 0,G) 588 IF (IDER.GT.0) NFG=NFG+1 589 END IF 590 IF (KD.GE.2.AND.LD.LT.2) CALL MXVSET(IH(NF+1)-1,0.0D 0,H) 591 NAV=0 592 NAG=0 593 DO 9 KA=1,NA 594 IF (KD.LT.0) GO TO 9 595 IF (LD.LT.0) THEN 596 CALL FUN(NF,KA,X,FA) 597 F=F+FA*FA 598 AF(KA)=FA 599 ELSE 600 FA=AF(KA) 601 END IF 602 IF (KD.LT.1) GO TO 9 603 IF (IDER.EQ.0) THEN 604 CALL PA0GS3(NF,KA,X,FA,GA,IAG,JAG,ETA1,NAV) 605 ELSE 606 CALL DFUN(NF,KA,X,GA) 607 END IF 608 IF (LD.LT.1) THEN 609 K=IAG(KA) 610 L=IAG(KA+1)-K 611 DO 1 J=1,L 612 JP=ABS(JAG(K)) 613 G(JP)=G(JP)+FA*GA(JP) 614 K=K+1 615 1 CONTINUE 616 END IF 617 IF (KD.LT.2) GO TO 9 618 IDECF=0 619 IF (IPOM1.EQ.0) THEN 620 IF (IDER.EQ.0) THEN 621 CALL PA0HS3(NF,KA,X,IX,HA,GO,GS,IAG,JAG,FA,ETA1,KBF,NAV) 622 ELSE 623 CALL PA1HS3(NF,KA,X,IX,HA,GO,GA,IAG,JAG,FA,ETA1,KBF,NAG) 624 END IF 625 END IF 626 CALL PASSH1(H,IH,JH,IAG,JAG,GA,KA,1.0D 0) 627 IF (IPOM1.EQ.0) CALL PASSH2(H,IH,JH,HA,IAG,JAG,KA,FA) 628 9 CONTINUE 629 IF (KD.GE.0.AND.LD.LT.0) F=5.0D-1*F 630 IF (IDER.EQ.0) NFV=NFV+NAV/NA 631 IF (IDER.GT.0) NFG=NFG+NAG/NA 632 LD=KD 633 RETURN 634 END 635* SUBROUTINE PALNG3 ALL SYSTEMS 97/12/01 636* PURPOSE : 637* COMPUTATION OF THE GRADIENT OF THE LINEAR APPROXIMATED FUNCTION. 638* 639* PARAMETERS : 640* RO AG(MA) SPARSE JACOBIAN MATRIX. 641* RI IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 642* RI JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 643* RO GA(NF) GRADIENT OF THE APPROXIMATED FUNCTION. 644* II KA INDEX OF THE SELECTED FUNCTION. 645* 646 SUBROUTINE PALNG3(AG,IAG,JAG,GA,KA) 647 DOUBLE PRECISION AG(*),GA(*) 648 INTEGER IAG(*),JAG(*),KA 649 INTEGER J,JP,K,L 650 K=IAG(KA) 651 L=IAG(KA+1)-K 652 DO 2 J=1,L 653 JP=ABS(JAG(K)) 654 GA(JP)=AG(K) 655 K=K+1 656 2 CONTINUE 657 RETURN 658 END 659* SUBROUTINE PASED3 ALL SYSTEMS 07/12/01 660* PURPOSE : 661* COMPRESSED SPARSE STRUCTURE OF THE JACOBIAN MATRIX IS COMPUTED FROM 662* THE COORDINATE FORM. 663* 664* PARAMETERS : 665* II NA NUMBER OF APPROXIMATED FUNCTIONS. 666* II MA NUMBER OF NONZERO ELEMENTS IN THE SPARSE JACOBIAN MATRIX. 667* IU IAG(MA+NA) ON INPUT ROW INDICES OF NONZERO ELEMENTS IN THE FIELD AG. 668* ON OUTPUT POSITIONS OF THE FIRST ROW ELEMENTS IN THE FIELD AG. 669* II JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 670* IO IER ERROR MESAGE. IER=0-THE STANDARD INPUT DATA ARE CORRECT. 671* IER=1-ERROR IN THE ARRAY IAG. IER=2-ERROR IN THE ARRAY JAG. 672* 673 SUBROUTINE PASED3(NA,MA,IAG,JAG,IER) 674 INTEGER NA,MA,IAG(*),JAG(*),IER 675 INTEGER I,J,K,L,KA 676 IER=0 677 CALL MXVSR7(MA,IAG,JAG) 678 IF (IAG(1).LT.1.OR.IAG(MA).GT.NA) THEN 679 IER=1 680 RETURN 681 END IF 682 CALL MXVINS(NA,0,IAG(MA+1)) 683 DO 1 J=1,MA 684 IAG(IAG(J)+MA)=IAG(IAG(J)+MA)+1 685 1 CONTINUE 686 IAG(1)=1 687 DO 2 KA=1,NA 688 IAG(KA+1)=IAG(KA)+IAG(KA+MA) 689 2 CONTINUE 690 I=0 691 DO 4 KA=1,NA 692 K=IAG(KA) 693 L=IAG(KA+1)-K 694 IF (L.GT.0) THEN 695 CALL MXVSRT(L,JAG(K)) 696 IF (JAG(K).LT.1.OR.JAG(K+L-1).GT.NA) THEN 697 IER=2 698 RETURN 699 END IF 700 END IF 701 IAG(KA)=IAG(KA)-I 702 DO 3 J=1,L 703 IF (J.GT.1.AND.JAG(K).EQ.JAG(K-1)) THEN 704 I=I+1 705 ELSE 706 JAG(K-I)=JAG(K) 707 END IF 708 K=K+1 709 3 CONTINUE 710 4 CONTINUE 711 IAG(NA+1)=IAG(NA+1)-I 712 MA=IAG(NA+1)-1 713 RETURN 714 END 715* SUBROUTINE PASSH1 ALL SYSTEMS 98/12/01 716* PURPOSE : 717* COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE SPARSE JACOBIAN 718* MATRIX. 719* 720* PARAMETERS : 721* RU H(M) NONZERO ELEMENTS OF THE SPARSE HESSIAN MATRIX. 722* II IH(NF+1) POINTERS OF THE DIAGONAL ELEMENTS OF H. 723* II JH(M) COLUMN INDICES OF THE NONZERO ELEMENTS OF H. 724* II IAG(NA+1) POSITIONS OF THE FIRST ROWS ELEMENTS IN THE SPARSE 725* JACOBIAN STRUCTURE. 726* II JAG(MA) COLUMN INDICES OF ELEMENTS IN THE SPARSE JACOBIAN 727* STRUCTURE. 728* RI GA(NF) GRADIENT OF THE SELECTED FUNCTION. 729* II KA INDEX OF THE SELECTED FUNCTION (ROW OF THE SPARSE JACOBIAN 730* MATRIX). 731* RI FACTOR SCALING FACTOR. 732* 733 SUBROUTINE PASSH1(H,IH,JH,IAG,JAG,GA,KA,FACTOR) 734 INTEGER IH(*),JH(*),IAG(*),JAG(*),KA 735 DOUBLE PRECISION H(*),GA(*),FACTOR 736 DOUBLE PRECISION TEMP 737 INTEGER I,J,JF,JA,K,LA 738 LA=IAG(KA+1)-1 739 DO 6 K=IAG(KA),LA 740 I=ABS(JAG(K)) 741 TEMP=FACTOR*GA(I) 742 JF=IH(I) 743 DO 5 JA=K,LA 744 J=ABS(JAG(JA)) 745 2 IF (ABS(JH(JF)).LT.J) THEN 746 JF=JF+1 747 GO TO 2 748 END IF 749 H(JF)=H(JF)+TEMP*GA(J) 750 5 CONTINUE 751 6 CONTINUE 752 RETURN 753 END 754* SUBROUTINE PASSH2 ALL SYSTEMS 98/12/01 755* PURPOSE : 756* COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE SPARSE JACOBIAN 757* MATRIX. 758* 759* PARAMETERS : 760* RU H(M) NONZERO ELEMENTS OF THE SPARSE HESSIAN MATRIX. 761* II IH(NF+1) POINTERS OF THE DIAGONAL ELEMENTS OF H. 762* II JH(M) COLUMN INDICES OF THE NONZERO ELEMENTS OF H. 763* II HA(ME) PACKED HESSIAN MATRIX OF THE SELECTED FUNCTION. 764* II IAG(NA+1) POSITIONS OF THE FIRST ROWS ELEMENTS IN THE SPARSE 765* JACOBIAN STRUCTURE. 766* II JAG(MA) COLUMN INDICES OF ELEMENTS IN THE SPARSE JACOBIAN 767* STRUCTURE. 768* II KA INDEX OF THE SELECTED FUNCTION (ROW OF THE SPARSE JACOBIAN 769* MATRIX). 770* RI FACTOR SCALING FACTOR. 771* 772 SUBROUTINE PASSH2(H,IH,JH,HA,IAG,JAG,KA,FACTOR) 773 INTEGER IH(*),JH(*),IAG(*),JAG(*),KA 774 DOUBLE PRECISION H(*),HA(*),FACTOR 775 INTEGER I,II,IA,J,JJ,JA,JF,K,KK,L 776 KK=0 777 II=IAG(KA) 778 L=IAG(KA+1)-II 779 DO 6 IA=1,L 780 KK=KK+IA 781 I=ABS(JAG(II)) 782 JF=IH(I) 783 JJ=II 784 K=KK 785 DO 4 JA=IA,L 786 J=ABS(JAG(JJ)) 787 2 IF (ABS(JH(JF)).LT.J) THEN 788 JF=JF+1 789 GO TO 2 790 END IF 791 H(JF)=H(JF)+FACTOR*HA(K) 792 K=K+JA 793 JJ=JJ+1 794 4 CONTINUE 795 II=II+1 796 6 CONTINUE 797 RETURN 798 END 799* SUBROUTINE PASSH3 ALL SYSTEMS 98/12/01 800* PURPOSE : 801* COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE SPARSE JACOBIAN 802* MATRIX. 803* 804* PARAMETERS : 805* RU H(M) NONZERO ELEMENTS OF THE SPARSE HESSIAN MATRIX. 806* II IH(NF+1) POINTERS OF THE DIAGONAL ELEMENTS OF H. 807* II JH(M) COLUMN INDICES OF THE NONZERO ELEMENTS OF H. 808* II IAG(NA+1) POSITIONS OF THE FIRST ROWS ELEMENTS IN THE SPARSE 809* JACOBIAN STRUCTURE. 810* II JAG(MA) COLUMN INDICES OF ELEMENTS IN THE SPARSE JACOBIAN 811* STRUCTURE. 812* RI GA(NF) GRADIENT OF THE SELECTED FUNCTION. 813* II KA INDEX OF THE SELECTED FUNCTION (ROW OF THE SPARSE JACOBIAN 814* MATRIX). 815* RI FACTOR SCALING FACTOR. 816* 817 SUBROUTINE PASSH3(H,IH,JH,IAG,JAG,GA,KA,FACTOR) 818 INTEGER IH(*),JH(*),IAG(*),JAG(*),KA 819 DOUBLE PRECISION H(*),GA(*),FACTOR 820 DOUBLE PRECISION TEMP 821 INTEGER I,J,JF,JA,K,LA 822 LA=IAG(KA+1)-1 823 DO 6 K=IAG(KA),LA 824 I=ABS(JAG(K)) 825 IF (I.LE.0) GO TO 6 826 TEMP=FACTOR*GA(I) 827 JF=IH(I) 828 DO 5 JA=K,LA 829 J=ABS(JAG(JA)) 830 IF (J.LE.0) GO TO 5 831 2 IF (ABS(JH(JF)).LT.J) THEN 832 JF=JF+1 833 GO TO 2 834 END IF 835 H(JF)=H(JF)+TEMP*GA(J) 836 5 CONTINUE 837 6 CONTINUE 838 RETURN 839 END 840* SUBROUTINE PCBS04 ALL SYSTEMS 98/12/01 841* PURPOSE : 842* INITIATION OF THE VECTOR CONTAINING TYPES OF CONSTRAINTS. 843* 844* PARAMETERS : 845* II NF NUMBER OF VARIABLES. 846* RI X(NF) VECTOR OF VARIABLES. 847* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. 848* RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. 849* RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. 850* RI EPS9 TOLERANCE FOR ACTIVE CONSTRAINTS. 851* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. 852* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. 853* 854 SUBROUTINE PCBS04(NF,X,IX,XL,XU,EPS9,KBF) 855 INTEGER NF,IX(*),KBF 856 DOUBLE PRECISION X(*),XL(*),XU(*),EPS9 857 DOUBLE PRECISION TEMP 858 INTEGER I,IXI 859 IF (KBF.GT.0) THEN 860 DO 1 I=1,NF 861 TEMP=1.0D 0 862 IXI=ABS(IX(I)) 863 IF ((IXI.EQ.1.OR.IXI.EQ.3.OR.IXI.EQ.4).AND.X(I).LE.XL(I)+ 864 & EPS9*MAX(ABS(XL(I)),TEMP)) X(I)=XL(I) 865 IF ((IXI.EQ.2.OR.IXI.EQ.3.OR.IXI.EQ.4).AND.X(I).GE.XU(I)- 866 & EPS9*MAX(ABS(XU(I)),TEMP)) X(I)=XU(I) 867 1 CONTINUE 868 END IF 869 RETURN 870 END 871* SUBROUTINE PDSGM1 ALL SYSTEMS 01/09/22 872* PURPOSE : 873* COMPUTATION OF A TRUST-REGION STEP BY THE DOG-LEG METHOD WITH DIRECT 874* MATRIX DECOMPOSITIONS. 875* 876* PARAMETERS : 877* II NF NUMBER OF VARIABLES. 878* II MMAX MAXIMUM DIMENSION OF THE SPARSE TABLEAU. 879* II MH POINTER OBTAINED BY THE SUBROUTINE MXSPCC. 880* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE 881* X(I) IS UNBOUNDED. IX(I)=1-LOWER BOUND XL(I).LE.X(I). 882* IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND 883* XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED. 884* RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 885* RA H(MMAX) NONZERO ELEMENTS OF THE APPROXIMATION OF THE SPARSE 886* HESSIAN MATRIX TOGETHER WITH AN ADDITIONAL SPACE USED FOR 887* THE NUMERICAL DIFFERENTIATION. 888* II IH(NF+1) POINTERS OF DIAGONAL ELEMENTS OF THE MATRIX H. 889* IU JH(MMAX) INDICES OF NONZERO ELEMENTS OF THE MATRIX H 890* TOGETHER WITH AN ADDITIONAL SPACE USED FOR THE NUMERICAL 891* DIFFERENTIATION. 892* RO S(NF) DIRECTION VECTOR. 893* RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. 894* RI GO(NF) GRADIENTS DIFFERENCE. 895* RA XS(NF) AUXILIARY VECTOR. 896* II PSL(NF+1) POINTER VECTOR OF THE COMPACT FORM OF THE TRIANGULAR 897* FACTOR OF THE HESSIAN APPROXIMATION. 898* IA PERM(NF) PERMUTATION VECTOR. 899* IA WN11(NF+1) AUXILIARY VECTOR. 900* IA WN12(NF+1) AUXILIARY VECTOR. 901* RI XMAX MAXIMUM STEPSIZE. 902* RU XDEL TRUST REGION RADIUS. 903* RO GNORM NORM OF THE GRADIENT VECTOR. 904* RO SNORM NORM OF THE DIRECTION VECTOR. 905* RI FMIN ESTIMATION OF THE MINIMUM FUNCTION VALUE. 906* RI F VALUE OF THE OBJECTIVE FUNCTION. 907* RO P VALUE OF THE DIRECTIONAL DERIVATIVE. 908* RO PP VALUE OF THE QUADRATIC TERM. 909* RI ETA2 TOLERANCE FOR POSITIVE DEFINITENESS. 910* RI ALF2 TOLERANCE FOR THE GRADIENT NORM. 911* II KD ORDER OF COMPUTED DERIVATIVES. 912* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. 913* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. 914* II IEST ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED. 915* IEST=1-MINIMUM IS ESTIMATED BY THE VALUE FMIN. 916* IU IDEC DECOMPOSITION INDICATOR. IDEC=0-NO DECOMPOSITION. 917* IU NDEC NUMBER OF MATRIX DECOMPOSITIONS. 918* II ITERD CAUSE OF TERMINATION IN THE DIRECTION DETERMINATION. 919* ITERD<0-BAD DECOMPOSITION. ITERD=0-DESCENT DIRECTION. 920* ITERD=1-NEWTON LIKE STEP. ITERD=2-INEXACT NEWTON LIKE STEP. 921* ITERD=3-BOUNDARY STEP. ITERD=4-DIRECTION WITH THE NEGATIVE 922* CURVATURE. ITERD=5-MARQUARDT STEP. 923* IO ITERM VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. 924* ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN 925* MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS. 926* ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN 927* MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS. 928* ITERM=3-IF F WAS LESS THAN OR EQUAL TO TOLB. 929* ITERM=4-IF GMAX WAS LESS THAN OR EQUAL TO TOLG. 930* ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, 931* BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. 932* ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. 933* ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. 934* VALUES ITERM<=-40 DETECT A LACK OF SPACE. 935* 936* SUBPROGRAMS USED : 937* S PNSTEP COMPUTATION OF THE BOUNDARY STEP. 938* S MXSPCB BACK SUBSTITUTION USING THE SPARSE DECOMPOSITION 939* OBTAINED BY MXSPCF. 940* S MXSPCD COMPUTATION OF A DIRECTION OF NEGATIVE CURVATURE USING 941* THE SPARSE DECOMPOSITION OBTAINED BY MXSPCF. 942* S MXSPCF GILL-MURRAY DECOMPOSITION OD A SPARSE SYMMETRIC MATRIX. 943* S MXSPCM MATRIX-VECTOR PRODUCT USING THE SPARSE DECOMPOSITION 944* OBTAINED BY MXSPCF. 945* RF MXSPCQ GENERALIZED DOT PRODUCT USING THE SPARSE DECOMPOSITION 946* OBTAINED BY MXSPCF. 947* S MXSPCT COPYING A SPARSE SYMMETRIC MATRIX INTO THE PERMUTED 948* FACTORIZED COMPACT SCHEME. 949* RF MXSSMQ COMPUTATION OF THE SPARSE QUADRATIC TERM. 950* S MXUCOP COPYING OF A VECTOR. 951* S MXUDIF DIFFERENCE OF TWO VECTORS. 952* S MXUDIR VECTOR AUGMENTED BY THE SCALED VECTOR. 953* RF MXUDOT DOT PRODUCT OF TWO VECTORS. 954* S MXUNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. 955* S MXVSBP INVERSE PERMUTATION OF A VECTOR 956* S MXVSCL SCALING OF A VECTOR. 957* S MXVSET INITIATION OF A VECTOR. 958* S MXVSFP PERMUTATION OF A VECTOR. 959* 960* METHOD : 961* J.E.DENNIS, H.H.W.MEI: AN UNCONSTRAINED OPTIMIZATION ALGORITHM WHICH 962* USES FUNCTION AND GRADIENT VALUES. REPORT NO. TR-75-246, DEPT. OF 963* COMPUTER SCIENCE, CORNELL UNIVERSITY 1975. 964* 965 SUBROUTINE PDSGM1(NF,MMAX,MH,IX,G,H,IH,JH,S,XO,GO,XS,PSL,PERM, 966 & WN11,WN12,XMAX,XDEL,GNORM,SNORM,FMIN,F,P,PP,ETA2,ALF2,KD,KBF, 967 & IEST,IDEC,NDEC,ITERD,ITERM) 968 INTEGER NF,MMAX,MH,IX(*),IH(*),JH(*),PSL(*),PERM(*),WN11(*), 969 & WN12(*),KD,KBF,IEST,IDEC,NDEC,ITERD,ITERM 970 DOUBLE PRECISION G(*),H(*),S(*),XO(*),GO(*),XS(*),XMAX,XDEL, 971 & GNORM,SNORM,FMIN,F,P,PP,ETA2,ALF2 972 INTEGER MM,INF,MODE 973 DOUBLE PRECISION B1,B2,B3,D3,S1,S2 974 DOUBLE PRECISION MXSSMQ,MXSPCQ,MXUDOT 975 SAVE INF 976* 977* DIRECTION DETERMINATION 978* 979 IF (IDEC.LT.0) IDEC=0 980 IF (IDEC.EQ.0) THEN 981 ELSE IF (IDEC.EQ.1) THEN 982 ELSE 983 ITERD=-1 984 GO TO 13130 985 END IF 986 MM=IH(NF+1)-1 987 B2=MXUDOT(NF,G,G,IX,KBF) 988 GNORM=SQRT(B2) 989 MODE=1 990 IF (ALF2*GNORM.LE.XDEL) THEN 991 MODE=2 992 IF (IDEC.EQ.0) THEN 993 CALL MXSPCT(NF,MM,MH,MMAX,H,JH,PSL,ITERM) 994 IF (ITERM.NE.0) GO TO 13130 995* 996* SPARSE GILL-MURRAY DECOMPOSITION 997* 998 S1=ETA2 999 CALL MXSPCF(NF,H(MM+1),PSL,JH(MM+1),WN11,WN12,XS,INF,S1,S2) 1000 NDEC=NDEC+1 1001 IDEC=1 1002 END IF 1003 IF (INF.GT.0) THEN 1004 CALL MXSPCD(NF,H(MM+1),PSL,JH(MM+1),S,INF) 1005 CALL MXVSBP(NF,PERM,S,XS) 1006* 1007* DIRECTION OF NEGATIVE CURVATURE 1008* 1009 SNORM=SQRT(MXUDOT(NF,S,S,IX,KBF)) 1010 IF (SNORM*SNORM*GNORM+S1*XDEL.LE.0.0D 0) THEN 1011 CALL MXVSCL(NF,XDEL/SNORM,S,S) 1012 SNORM=XDEL 1013 ITERD=4 1014 GO TO 13120 1015 END IF 1016 ELSE IF (GNORM.LE.0.0D 0) THEN 1017* 1018* ZERO DIRECTION 1019* 1020 SNORM=0.0D 0 1021 CALL MXVSET(NF,0.0D 0,S) 1022 GO TO 13120 1023 END IF 1024 END IF 1025 IF (IDEC.EQ.0) THEN 1026 B1=MXSSMQ(NF,H,IH,JH,G,G) 1027 ELSE 1028 CALL MXUCOP(NF,G,GO,IX,KBF) 1029 CALL MXVSFP(NF,PERM,GO,XS) 1030 CALL MXSPCM(NF,H(MM+1),PSL,JH(MM+1),GO,XS,1) 1031 B1=MXSPCQ(NF,H(MM+1),PSL,GO) 1032 END IF 1033 IF (XDEL.LE.0.0D 0) THEN 1034* 1035* INITIAL TRUST REGION BOUND 1036* 1037 IF (B1.LE.0.0D 0) THEN 1038 XDEL=GNORM 1039 ELSE 1040 XDEL=(B2/B1)*GNORM 1041 END IF 1042 IF (IEST.EQ.1) XDEL=MIN(XDEL,4.0D 0*(F-FMIN)/GNORM) 1043 XDEL=MIN(XDEL,XMAX) 1044 END IF 1045 IF (B1.LE.0.0D 0.OR.B2*GNORM.GE.B1*XDEL) THEN 1046* 1047* SCALED STEEPEST DESCENT DIRECTION IS ACCEPTED 1048* 1049 CALL MXVSCL(NF,-XDEL/GNORM,G,S) 1050 SNORM=XDEL 1051 ITERD=3 1052 GO TO 13120 1053 END IF 1054 IF (IDEC.EQ.0) THEN 1055 CALL MXSPCT(NF,MM,MH,MMAX,H,JH,PSL,ITERM) 1056 IF (ITERM.NE.0) THEN 1057 GO TO 13130 1058 END IF 1059* 1060* SPARSE GILL-MURRAY DECOMPOSITION 1061* 1062 S1=ETA2 1063 CALL MXSPCF(NF,H(MM+1),PSL,JH(MM+1),WN11,WN12,XS,INF,S1,S2) 1064 NDEC=NDEC+1 1065 IDEC=1 1066 END IF 1067* 1068* COMPUTATION OF THE NEWTON DIRECTION 1069* 1070 CALL MXUCOP(NF,G,GO,IX,KBF) 1071 CALL MXVSFP(NF,PERM,GO,XS) 1072 CALL MXSPCB(NF,H(MM+1),PSL,JH(MM+1),GO,0) 1073 CALL MXVSBP(NF,PERM,GO,XS) 1074 D3=SQRT(MXUDOT(NF,GO,GO,IX,KBF)) 1075* 1076* COMPUTATION OF THE STEEPEST DESCENT DIRECTION 1077* 1078 B2=B2/B1 1079 SNORM=B2*GNORM 1080 CALL MXVSCL(NF,-B2,G,S) 1081 CALL MXUNEG(NF,GO,GO,IX,KBF) 1082 CALL MXUDIF(NF,GO,S,XO,IX,KBF) 1083 B1=MXUDOT(NF,S,XO,IX,KBF) 1084 B2=MXUDOT(NF,XO,XO,IX,KBF) 1085 IF (B2.LE.1.0D-8*XDEL*XDEL) THEN 1086* 1087* NEWTON AND THE STEEPEST DESCENT DIRECTION ARE 1088* APPROXIMATELY EQUAL 1089* 1090 CALL MXUCOP(NF,GO,S,IX,KBF) 1091 SNORM=D3 1092 ITERD=1 1093 ELSE IF (B1.LE.0.0D 0) THEN 1094* 1095* BOUNDARY STEP WITH NEGATIVE INCREMENT 1096* 1097 CALL PNSTEP(XDEL,SNORM,-B1,B2,B3) 1098 CALL MXUDIR(NF,-B3,XO,S,S,IX,KBF) 1099 SNORM=XDEL 1100 ITERD=3 1101 ELSE IF (D3.LE.XDEL) THEN 1102* 1103* NEWTON DIRECTION IS ACCEPTED 1104* 1105 CALL MXUCOP(NF,GO,S,IX,KBF) 1106 SNORM=D3 1107 ITERD=1 1108 ELSE 1109* 1110* DOUBLE DOGLEG STRATEGY 1111* 1112 D3=XDEL/D3 1113 B3=MXUDOT(NF,S,GO,IX,KBF) 1114 D3=MAX(D3,SNORM*SNORM/B3) 1115 CALL MXUDIR(NF,-D3,GO,S,XO,IX,KBF) 1116 B1=SNORM*SNORM-D3*B3 1117 B2=MXUDOT(NF,XO,XO,IX,KBF) 1118 CALL PNSTEP(XDEL,SNORM,-B1,B2,B3) 1119 CALL MXUDIR(NF,-B3,XO,S,S,IX,KBF) 1120 SNORM=XDEL 1121 ITERD=3 1122 END IF 112313120 CONTINUE 1124 IF (IDEC.EQ.0) THEN 1125 PP=MXSSMQ(NF,H,IH,JH,S,S)*0.5D 0 1126 ELSE 1127 CALL MXUCOP(NF,S,GO,IX,KBF) 1128 CALL MXVSFP(NF,PERM,GO,XS) 1129 CALL MXSPCM(NF,H(MM+1),PSL,JH(MM+1),GO,XS,1) 1130 PP=MXSPCQ(NF,H(MM+1),PSL,GO)*0.5D 0 1131 IF (ITERD.EQ.1.AND.INF.NE.0) ITERD=2 1132 END IF 113313130 CONTINUE 1134 IF (KD.GT.0) P=MXUDOT(NF,G,S,IX,KBF) 1135 RETURN 1136 END 1137* SUBROUTINE PDSGM4 ALL SYSTEMS 01/09/22 1138* PURPOSE : 1139* COMPUTATION OF A TRUST-REGION STEP BY THE SHIFTED STEIHAUG-TOINT 1140* METHOD WITH CONJUGATE GRADIENT ITERATIONS. 1141* 1142* PARAMETERS : 1143* II NF NUMBER OF VARIABLES. 1144* II MMAX MAXIMUM DIMENSION OF THE SPARSE TABLEAU. 1145* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE 1146* X(I) IS UNBOUNDED. IX(I)=1-LOWER BOUND XL(I).LE.X(I). 1147* IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND 1148* XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED. 1149* RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 1150* RA H(MMAX) NONZERO ELEMENTS OF THE APPROXIMATION OF THE SPARSE 1151* HESSIAN MATRIX TOGETHER WITH AN ADDITIONAL SPACE USED FOR 1152* THE NUMERICAL DIFFERENTIATION. 1153* II IH(NF+1) POINTERS OF DIAGONAL ELEMENTS OF THE MATRIX H. 1154* IU JH(MMAX) INDICES OF NONZERO ELEMENTS OF THE MATRIX H 1155* TOGETHER WITH AN ADDITIONAL SPACE USED FOR THE NUMERICAL 1156* DIFFERENTIATION. 1157* RO S(NF) DIRECTION VECTOR. 1158* RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. 1159* RI GO(NF) GRADIENTS DIFFERENCE. 1160* RA XS(NF) AUXILIARY VECTOR. 1161* RA GS(NF) AUXILIARY VECTOR. 1162* IA IW(NF+1) AUXILIARY VECTOR. 1163* RI XMAX MAXIMUM STEPSIZE. 1164* RU XDEL TRUST REGION RADIUS. 1165* RO GNORM NORM OF THE GRADIENT VECTOR. 1166* RO GNORMO OLD NORM OF THE GRADIENT VECTOR. 1167* RO SNORM NORM OF THE DIRECTION VECTOR. 1168* RI FMIN ESTIMATION OF THE MINIMUM FUNCTION VALUE. 1169* RI F VALUE OF THE OBJECTIVE FUNCTION. 1170* RO P VALUE OF THE DIRECTIONAL DERIVATIVE. 1171* RO PP VALUE OF THE QUADRATIC TERM. 1172* RI ETA0 MACHINE PRECISION. 1173* RI ETA2 TOLERANCE FOR POSITIVE DEFINITENESS. 1174* RI DEL1 LOWER TOLERANCE FOR THE TRUST-REGION RADIUS. 1175* II KD ORDER OF COMPUTED DERIVATIVES. 1176* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. 1177* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. 1178* II MOS1 NUMBER OF LANCZOS STEPS IN THE SHIFTED STEIHAUG-TOINT 1179* METHOD. 1180* II MOS2 TYPE OF PRECONDITIONING. MOS2=1-PRECONDITIONING IS NOT 1181* USED. MOS2=2-PRECONDITIONING BY THE INCOMPLETE GILL-MURRAY 1182* DECOMPOSITION. MOS2=3-PRECONDITIONING BY THE INCOMPLETE 1183* GILL-MURRAY DECOMPOSITION WITH A PRELIMINARY SOLUTION OF 1184* THE PRECONDITIONED SYSTEM WHICH IS USED IF IT SATISFIES 1185* THE TERMINATION CRITERION. 1186* II MOS3 PRECONDITIONING IN ILL-CONTITIONED AND INDEFINITE CASES. 1187* MOS3=0-PRECONDITIONING IN BOTH THESE CASES IS SUPPRESSED. 1188* MOS3=1-PRECONDITIONING IN ILL-CONDITIONED CASE IS SUPPRESSED. 1189* MOS3=2-PRECONDITIONING IS ALWAYS USED. 1190* II IEST ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED. 1191* IEST=1-MINIMUM IS ESTIMATED BY THE VALUE FMIN. 1192* IU IDEC DECOMPOSITION INDICATOR. IDEC=0-NO DECOMPOSITION. 1193* IU NDEC NUMBER OF MATRIX DECOMPOSITIONS. 1194* II NIT NUMBER OF OUTER ITERATIONS. 1195* IU NIN NUMBER OF INNER CONJUGATE GRADIENT ITERATIONS. 1196* II ITERD CAUSE OF TERMINATION IN THE DIRECTION DETERMINATION. 1197* ITERD<0-BAD DECOMPOSITION. ITERD=0-DESCENT DIRECTION. 1198* ITERD=1-NEWTON LIKE STEP. ITERD=2-INEXACT NEWTON LIKE STEP. 1199* ITERD=3-BOUNDARY STEP. ITERD=4-DIRECTION WITH THE NEGATIVE 1200* CURVATURE. ITERD=5-MARQUARDT STEP. 1201* IO ITERM VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. 1202* ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN 1203* MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS. 1204* ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN 1205* MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS. 1206* ITERM=3-IF F WAS LESS THAN OR EQUAL TO TOLB. 1207* ITERM=4-IF GMAX WAS LESS THAN OR EQUAL TO TOLG. 1208* ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, 1209* BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. 1210* ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. 1211* ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. 1212* VALUES ITERM<=-40 DETECT A LACK OF SPACE. 1213* 1214* SUBPROGRAMS USED : 1215* S PNSTEP COMPUTATION OF THE BOUNDARY STEP. 1216* S MXSPTB BACK SUBSTITUTION AFTER THE GILL-MURRAY DECOMPOSITION. 1217* S MXSPTF INCOMPLETE GILL-MURRAY DECOMPOSITION. 1218* S MXSSDA SPARSE SYMMETRIC MATRIX IS AUGMENTED BY THE SCALED UNIT 1219* MATRIX. 1220* S MXSSMD MATRIX-VECTOR PRODUCT FOLLOWED BY THE ADDITION OF A 1221* SCALED VECTOR. 1222* S MXSSMM MATRIX-VECTOR PRODUCT. 1223* RF MXSSMQ COMPUTATION OF THE SPARSE QUADRATIC TERM. 1224* S MXTPGB BACK SUBSTITUTION FOR A DECOMPOSED TRIDIAGONAL MATRIX. 1225* S MXTPGF CHOLESKI DECOMPOSITION OF A TRIDIAGONAL MATRIX. 1226* S MXUDIR VECTOR AUGMENTED BY THE SCALED VECTOR. 1227* RF MXUDEL NORM OF VECTOR DIFFERENCE. 1228* RF MXUDOT DOT PRODUCT OF TWO VECTORS. 1229* RF MXUNOR EUCLIDEAN NORM OF A VECTOR. 1230* S MXVCOP COPYING OF A VECTOR. 1231* S MXVCOR CORRECTION OF A VECTOR (ZERO ELEMENTS ARE REPLACED BY 1232* THE NONZERO NUMBER). 1233* RF MXVDOT DOT PRODUCT OF TWO VECTORS. 1234* S MXVNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. 1235* S MXVSCL SCALING OF A VECTOR. 1236* S MXVSET INITIATION OF A VECTOR. 1237* S MXVSUM SUM OF TWO VECTORS. 1238* RF MXVVDP GENERALIZED DOT PRODUCT. 1239* 1240* METHOD : 1241* L.LUKSAN, C.MATONOHA, J.VLCEK: A SHIFTED STEIHAUG-TOINT METHOD FOR 1242* COMPUTING TRUST-REGION STEP. REPORT NO. V-914, INST. OF COMPUTER 1243* SCIENCE, CZECH ACADEMY OF SCIENCES, 2004. 1244* 1245 SUBROUTINE PDSGM4(NF,MMAX,IX,G,H,IH,JH,S,XO,GO,XS,GS,IW,XMAX, 1246 & XDEL,GNORM,GNORMO,SNORM,FMIN,F,P,PP,ETA0,ETA2,DEL1,KD,KBF, 1247 & MOS1,MOS2,MOS3,IEST,IDEC,NDEC,NIT,NIN,ITERD,ITERM) 1248 INTEGER NF,MMAX,IX(*),IH(*),JH(*),IW(*),KD,KBF,MOS1,MOS2,MOS3, 1249 & IEST,IDEC,NDEC,NIT,NIN,ITERD,ITERM 1250 DOUBLE PRECISION G(*),H(*),S(*),XO(*),GO(*),XS(*),GS(*),XMAX, 1251 & XDEL,GNORM,GNORMO,SNORM,FMIN,F,P,PP,ETA0,ETA2,DEL1 1252 INTEGER NOS1,NOS2,NRED,I,M,INF 1253 DOUBLE PRECISION T,EL,EU,PAR,ALF,EPS,RHO,RHO1,RHO2,SIG,TAU 1254 DOUBLE PRECISION MXSSMQ,MXUDOT,MXUDEL,MXUNOR,MXVDOT,MXVVDP 1255 SAVE EPS 1256* 1257* DIRECTION DETERMINATION 1258* 1259 IF (NIT.LE.1) THEN 1260 EPS=0.9D 0 1261 GNORMO=1.0D 60 1262 END IF 1263 IF (IDEC.LT.0) IDEC=0 1264 IF (IDEC.NE.0.AND.IDEC.NE.1) THEN 1265 ITERD=-1 1266 GO TO 13180 1267 END IF 1268 GNORM=SQRT(MXUDOT(NF,G,G,IX,KBF)) 1269 IF (GNORM.GE.1.0D 3*GNORMO) EPS=1.0D-6 1270 GNORMO=GNORM 1271 RHO1=MXUDOT(NF,G,G,IX,KBF) 1272 IF (XDEL.LE.0.0D 0) THEN 1273* 1274* INITIAL TRUST REGION BOUND 1275* 1276 RHO2=MXSSMQ(NF,H,IH,JH,G,G) 1277 IF (RHO2.LE.0.0D 0) THEN 1278 XDEL=GNORM 1279 ELSE 1280 XDEL=(GNORM*GNORM/RHO2)*GNORM 1281 END IF 1282 IF (IEST.EQ.1) XDEL=MIN(XDEL,4.0D 0*(F-FMIN)/GNORM) 1283 XDEL=MIN(XDEL,XMAX) 1284 END IF 1285 PAR=MIN(EPS,SQRT(GNORM)) 1286 IF (PAR.GT.1.0D-2) THEN 1287 PAR=MIN(PAR,1.0D 0/DBLE(NIT)) 1288 END IF 1289 PAR=PAR*PAR 1290 NOS1=MIN(NF,MOS1) 1291 IF (NOS1.LE.1) THEN 1292 T=0.0D 0 1293 ELSE 1294* 1295* INCOMPLETE LANCZOS TRIDIAGONALIZATION 1296* 1297 INF=0 1298 CALL MXVCOP(NF,G,XS) 1299 CALL MXVSET(NF,0.0D 0,GS) 1300 CALL MXVSCL(NF,1.0D 0/MXUNOR(NF,XS,IX,KBF),XS,XS) 1301 DO 13111 NRED=1,NOS1 1302 IF (NRED.GT.1) THEN 1303 DO 13112 I=1,NF 1304 EL=XS(I) 1305 XS(I)=GS(I)/EU 1306 GS(I)=-EU*EL 130713112 CONTINUE 1308 END IF 1309 CALL MXSSMD(NF,H,IH,JH,XS,1.0D 0,GS,GS) 1310 EL=MXUDOT(NF,XS,GS,IX,KBF) 1311 CALL MXUDIR(NF,-EL,XS,GS,GS,IX,KBF) 1312 EU=MXUNOR(NF,GS,IX,KBF) 1313 IF (EU.LE.0.0D 0) THEN 1314 INF=NRED 1315 GO TO 13116 1316 END IF 1317 XO(NRED)=EL 1318 GO(NRED)=EU 131913111 CONTINUE 132013116 CONTINUE 1321 CALL MXVCOR(NOS1,ETA0,XO) 1322 T=0.0D 0 1323 RHO2=DEL1*XDEL 1324 DO 13117 NRED=1,10 1325 T=MIN(T,1.0D 5) 1326 IF (T.GE.1.0D 5) GO TO 13118 1327* 1328* SOLUTION OF THE TRIDIAGONAL SYSTEM 1329* 1330 ALF=ETA0 1331 CALL MXVSET(NOS1,T,XS) 1332 CALL MXVSUM(NOS1,XO,XS,XS) 1333 CALL MXVCOP(NOS1,GO,GS) 1334 CALL MXTPGF(NOS1,XS,GS,INF,ALF,TAU) 1335 CALL MXVSET(NOS1,0.0D 0,S) 1336 S(1)=GNORM 1337 CALL MXTPGB(NOS1,XS,GS,S,0) 1338 RHO=MXVDOT(NOS1,S,S) 1339 IF (RHO.LE.XDEL**2) GO TO 13118 1340 CALL MXTPGB(NOS1,XS,GS,S,1) 1341* 1342* MARQUARDT PARAMETER T IS COMPUTED USING THE ONE-DIMENSIONAL 1343* NEWTON METHOD 1344* 1345 T=T+(RHO/MXVVDP(NOS1,XS,S,S))*((SQRT(RHO)-RHO2)/RHO2) 134613117 CONTINUE 1347 END IF 134813118 CONTINUE 1349 CALL MXVNEG(NF,G,XO) 1350 NOS2=MOS2-1 1351 IF (NOS2.GT.0) THEN 1352* 1353* INCOMPLETE GILL-MURRAY DECOMPOSITION 1354* 1355 ALF=ETA2 1356 M=IH(NF+1)-1 1357 IF (2*M.GE.MMAX) THEN 1358 ITERM=-48 1359 GO TO 13180 1360 END IF 1361 CALL MXVCOP(M,H,H(M+1)) 1362 IF (T.GT.0.0D 0) CALL MXSSDA(NF,H(M+1),IH,T) 1363 CALL MXSPTF(NF,H(M+1),IH,JH,IW,INF,ALF,SIG) 1364 IF (INF+10.LT.0) THEN 1365 ITERM=-48 1366 GO TO 13180 1367 END IF 1368 IF (MOS3.EQ.0) THEN 1369 IF (INF.NE.0) NOS2=0 1370 ELSE IF (MOS3.EQ.1) THEN 1371 IF (INF.LT.0) NOS2=0 1372 END IF 1373 NDEC=NDEC+1 1374 IDEC=1 1375 IF (NOS2.GT.1) THEN 1376* 1377* PRELIMINARY INEXACT SOLUTION 1378* 1379 CALL MXSPTB(NF,H(M+1),IH,JH,XO,0) 1380 SNORM=SQRT(MXUDOT(NF,XO,XO,IX,KBF)) 1381 IF (SNORM.LE.XDEL*1.0D 5) THEN 1382 CALL MXVCOP(NF,XO,S) 1383 IF (SNORM.LE.XDEL) THEN 1384 ITERD=2 1385 ELSE 1386 CALL MXVSCL(NF,XDEL/SNORM,S,S) 1387 SNORM=XDEL 1388 ITERD=3 1389 END IF 1390 CALL MXSSMD(NF,H,IH,JH,S,1.0D 0,G,GO) 1391 IF (MXUDOT(NF,GO,GO,IX,KBF).LE.1.0D-2*PAR*RHO1) GO TO 13180 1392 END IF 1393 END IF 1394 END IF 1395* 1396* CG INITIATION 1397* 1398 RHO=RHO1 1399 SNORM=0.0D 0 1400 CALL MXVSET(NF,0.0D 0,S) 1401 CALL MXVNEG(NF,G,XS) 1402 IF (NOS2.EQ.0) THEN 1403 ELSE IF (NOS2.EQ.1) THEN 1404 CALL MXSPTB(NF,H(M+1),IH,JH,XO,0) 1405 RHO=MXUDOT(NF,XS,XO,IX,KBF) 1406 ELSE 1407 RHO=MXUDOT(NF,XS,XO,IX,KBF) 1408 END IF 1409 DO 13120 NRED=1,NF+3 1410 IF (T.GT.0.0D 0) THEN 1411 CALL MXSSMD(NF,H,IH,JH,XO,T,XO,GO) 1412 ELSE 1413 CALL MXSSMM(NF,H,IH,JH,XO,GO) 1414 END IF 1415 ALF=MXUDOT(NF,XO,GO,IX,KBF) 1416 IF (ALF.LE.0.0D 0) GO TO 13160 1417 ALF=RHO/ALF 1418 RHO2=SQRT(MXUDEL(NF,ALF,XO,S,IX,KBF)) 1419 IF (RHO2.GE.XDEL) GO TO 13160 1420* 1421* CG STEP 1422* 1423 CALL MXUDIR(NF, ALF,XO,S,S,IX,KBF) 1424 CALL MXUDIR(NF,-ALF,GO,XS,XS,IX,KBF) 1425 NIN=NIN+1 1426 SNORM=RHO2 1427 RHO2=MXUDOT(NF,XS,XS,IX,KBF) 1428 IF (RHO2.LE.PAR*RHO1) GO TO 13150 1429 IF (NRED.GE.NF+3) GO TO 13150 1430 IF (NOS2.NE.0) THEN 1431 CALL MXVCOP(NF,XS,GO) 1432 CALL MXSPTB(NF,H(M+1),IH,JH,GO,0) 1433 RHO2=MXUDOT(NF,XS,GO,IX,KBF) 1434 ALF=RHO2/RHO 1435 CALL MXUDIR(NF,ALF,XO,GO,XO,IX,KBF) 1436 ELSE 1437 ALF=RHO2/RHO 1438 CALL MXUDIR(NF,ALF,XO,XS,XO,IX,KBF) 1439 END IF 1440 RHO=RHO2 144113120 CONTINUE 1442* 1443* AN INEXACT SOLUTION IS OBTAINED 1444* 144513150 CONTINUE 1446 SNORM=SQRT(MXUDOT(NF,S,S,IX,KBF)) 1447 ITERD=2 1448 GO TO 13180 1449* 1450* BOUNDARY STEP IS COMPUTED 1451* 145213160 CONTINUE 1453 RHO1=MXUDOT(NF,XO,S,IX,KBF) 1454 RHO2=MXUDOT(NF,XO,XO,IX,KBF) 1455 CALL PNSTEP(XDEL,SNORM,RHO1,RHO2,ALF) 1456 CALL MXUDIR(NF,ALF,XO,S,S,IX,KBF) 1457 SNORM=XDEL 1458 ITERD=3 1459 NRED=-NRED 146013180 CONTINUE 1461 PP=MXSSMQ(NF,H,IH,JH,S,S)*0.5D 0 1462 IF (KD.GT.0) P=MXUDOT(NF,G,S,IX,KBF) 1463 RETURN 1464 END 1465* SUBROUTINE PDSGM7 ALL SYSTEMS 01/09/22 1466* PURPOSE : 1467* COMPUTATION OF A TRUST-REGION STEP BY THE MORE-SORENSEN METHOD WITH 1468* DIRECT MATRIX DECOMPOSITIONS. 1469* 1470* PARAMETERS : 1471* II NF NUMBER OF VARIABLES. 1472* II MMAX MAXIMUM DIMENSION OF THE SPARSE TABLEAU. 1473* II MH POINTER OBTAINED BY THE SUBROUTINE MXSPCC. 1474* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE 1475* X(I) IS UNBOUNDED. IX(I)=1-LOWER BOUND XL(I).LE.X(I). 1476* IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND 1477* XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED. 1478* RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 1479* RA H(MMAX) NONZERO ELEMENTS OF THE APPROXIMATION OF THE SPARSE 1480* HESSIAN MATRIX TOGETHER WITH AN ADDITIONAL SPACE USED FOR 1481* THE NUMERICAL DIFFERENTIATION. 1482* II IH(NF+1) POINTERS OF DIAGONAL ELEMENTS OF THE MATRIX H. 1483* IU JH(MMAX) INDICES OF NONZERO ELEMENTS OF THE MATRIX H 1484* TOGETHER WITH AN ADDITIONAL SPACE USED FOR THE NUMERICAL 1485* DIFFERENTIATION. 1486* RO S(NF) DIRECTION VECTOR. 1487* RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. 1488* RI GO(NF) GRADIENTS DIFFERENCE. 1489* II PSL(NF+1) POINTER VECTOR OF THE COMPACT FORM OF THE TRIANGULAR 1490* FACTOR OF THE HESSIAN APPROXIMATION. 1491* IA PERM(NF) PERMUTATION VECTOR. 1492* IA WN11(NF+1) AUXILIARY VECTOR. 1493* IA WN12(NF+1) AUXILIARY VECTOR. 1494* RI XMAX MAXIMUM STEPSIZE. 1495* RI XDEL TRUST REGION RADIUS. 1496* RO XDELO OLD TRUST REGION RADIUS. 1497* RO GNORM NORM OF THE GRADIENT VECTOR. 1498* RO SNORM NORM OF THE DIRECTION VECTOR. 1499* RI FMIN ESTIMATION OF THE MINIMUM FUNCTION VALUE. 1500* RO F VALUE OF THE OBJECTIVE FUNCTION. 1501* RO P VALUE OF THE DIRECTIONAL DERIVATIVE. 1502* RO PP VALUE OF THE QUADRATIC TERM. 1503* RI ETA2 TOLERANCE FOR POSITIVE DEFINITENESS. 1504* RI DEL1 LOWER TOLERANCE FOR THE TRUST-REGION RADIUS. 1505* RI DEL2 UPPER TOLERANCE FOR THE TRUST-REGION RADIUS. 1506* II KD ORDER OF COMPUTED DERIVATIVES. 1507* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. 1508* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. 1509* II IEST ESTIMATION INDICATOR. IEST=0-MINIMUM IS NOT ESTIMATED. 1510* IEST=1-MINIMUM IS ESTIMATED BY THE VALUE FMIN. 1511* II IDIR TRUST-REGION CHANGE INDICATOR. 1512* IU IDEC DECOMPOSITION INDICATOR. IDEC=0-NO DECOMPOSITION. 1513* IU NDEC NUMBER OF MATRIX DECOMPOSITIONS. 1514* II ITERD CAUSE OF TERMINATION IN THE DIRECTION DETERMINATION. 1515* ITERD<0-BAD DECOMPOSITION. ITERD=0-DESCENT DIRECTION. 1516* ITERD=1-NEWTON LIKE STEP. ITERD=2-INEXACT NEWTON LIKE STEP. 1517* ITERD=3-BOUNDARY STEP. ITERD=4-DIRECTION WITH THE NEGATIVE 1518* CURVATURE. ITERD=5-MARQUARDT STEP. 1519* IO ITERM VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. 1520* ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN 1521* MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS. 1522* ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN 1523* MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS. 1524* ITERM=3-IF F WAS LESS THAN OR EQUAL TO TOLB. 1525* ITERM=4-IF GMAX WAS LESS THAN OR EQUAL TO TOLG. 1526* ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, 1527* BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. 1528* ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. 1529* ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. 1530* VALUES ITERM<=-40 DETECT A LACK OF SPACE. 1531* 1532* SUBPROGRAMS USED : 1533* S PNSTEP COMPUTATION OF THE BOUNDARY STEP. 1534* S MXSPCA ADDITION OF THE LEVENBERG-MARQUARDT TERM TO THE SPARSE 1535* SYMMETRIC MATRIX. 1536* S MXSPCB BACK SUBSTITUTION USING THE SPARSE DECOMPOSITION 1537* OBTAINED BY MXSPCF. 1538* S MXSPCD COMPUTATION OF A DIRECTION OF NEGATIVE CURVATURE USING 1539* THE SPARSE DECOMPOSITION OBTAINED BY MXSPCF. 1540* S MXSPCF GILL-MURRAY DECOMPOSITION OD A SPARSE SYMMETRIC MATRIX. 1541* S MXSPCN ESTIMATION OF THE MINIMUM EIGENVALUE AND THE 1542* CORRESPONDING EIGENVECTOR OF A SYMMETRIC MATRIX USING THE 1543* SPARSE DECOMPOSITION OBTAINED BY MXSPCF. 1544* RF MXSPCP GENERALIZED DOT PRODUCT USING THE SPARSE DECOMPOSITION 1545* OBTAINED BY MXSPCF. 1546* S MXSPCT COPYING A SPARSE SYMMETRIC MATRIX INTO THE PERMUTED 1547* FACTORIZED COMPACT SCHEME. 1548* RF MXSSDL DETERMINATION OF A MINIMUM DIAGONAL ELEMENT OF A SPARSE 1549* SYMMETRIC MATRIX. 1550* S MXSSMG GERSHGORIN BOUNDS FOR EIGENVALUES OF A SPARSE SYMMETRIC 1551* MATRIX 1552* RF MXSSMQ COMPUTATION OF THE SPARSE QUADRATIC TERM. 1553* S MXUCOP COPYING OF A VECTOR. 1554* S MXUDIR VECTOR AUGMENTED BY THE SCALED VECTOR. 1555* RF MXUDOT DOT PRODUCT OF TWO VECTORS. 1556* S MXUNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. 1557* S MXVSBP INVERSE PERMUTATION OF A VECTOR 1558* S MXVSFP PERMUTATION OF A VECTOR. 1559* 1560* METHOD : 1561* J.J.MORE, D.C.SORENSEN: COMPUTING A TRUST REGION STEP. REPORT NO. 1562* ANL-81-83, ARGONNE NATIONAL LAB. 1981. 1563* 1564 SUBROUTINE PDSGM7(NF,MMAX,MH,IX,G,H,IH,JH,S,XO,GO,PSL,PERM,WN11, 1565 & WN12,XMAX,XDEL,XDELO,GNORM,SNORM,FMIN,F,P,PP,ETA2,DEL1,DEL2,KD, 1566 & KBF,IEST,IDIR,IDEC,NDEC,ITERD,ITERM) 1567 INTEGER NF,MMAX,MH,IX(*),IH(*),JH(*),PSL(*),PERM(*),WN11(*), 1568 & WN12(*),KD,KBF,IEST,IDIR,IDEC,NDEC,ITERD,ITERM 1569 DOUBLE PRECISION G(*),H(*),S(*),XO(*),GO(*),XMAX,XDEL,XDELO, 1570 & GNORM,SNORM,FMIN,F,P,PP,ETA2,DEL1,DEL2 1571 INTEGER NRED,MM,INF,MODE 1572 DOUBLE PRECISION T,TL,TU,E,EL,EU,ALF,RHO,RHO1,RHO2,CON 1573 DOUBLE PRECISION MXSSMQ,MXSPCP,MXSSDL,MXUDOT 1574 SAVE T,TL,TU,E,EL,EU 1575* 1576* DIRECTION DETERMINATION 1577* 1578 IF (IDEC.LT.0) IDEC=0 1579 IF (IDEC.NE.0) THEN 1580 ITERD=-1 1581 GO TO 13250 1582 END IF 1583 MM=IH(NF+1)-1 1584 GNORM=SQRT(MXUDOT(NF,G,G,IX,KBF)) 1585 IF (XDEL.LE.0.0D 0) THEN 1586* 1587* INITIAL TRUST REGION BOUND 1588* 1589 RHO1=MXSSMQ(NF,H,IH,JH,G,G) 1590 RHO2=GNORM*GNORM 1591 IF (RHO1.LE.0.0D 0) THEN 1592 XDEL=GNORM 1593 ELSE 1594 XDEL=(RHO2/RHO1)*GNORM 1595 END IF 1596 IF (IEST.EQ.1) XDEL=MIN(XDEL,4.0D 0*(F-FMIN)/GNORM) 1597 XDEL=MIN(XDEL,XMAX) 1598 END IF 1599* 1600* INITIAL BOUNDS FOR THE PARAMETER T 1601* 1602 NRED=0 1603 IF (IDIR.LE.0) THEN 1604 T=0.0D 0 1605 E=-MXSSDL(NF,H,IH,JH,INF) 1606 CALL MXSSMG(NF,H,IH,JH,EL,EU,S) 1607 TL=GNORM/XDEL-EU 1608 TU=GNORM/XDEL-EL 1609 ELSE IF (IDIR.EQ.1) THEN 1610 T=T*XDELO/XDEL 1611 TL=MAX(TL,GNORM/XDEL-EU) 1612 TU=GNORM/XDEL-EL 1613 ELSE IF (IDIR.EQ.2) THEN 1614 T=T*XDELO/XDEL 1615 TL=GNORM/XDEL-EU 1616 TU=MIN(TU,GNORM/XDEL-EL) 1617 END IF 1618 TL=MAX(TL,0.0D 0,E) 1619 TU=MAX(TL,TU) 1620 T=MAX(T,TL) 1621 T=MIN(T,TU) 162213220 CONTINUE 1623 TL=MAX(TL,E) 1624 IF (T.LE.E.AND.NRED.NE.0) THEN 1625* 1626* THE PARAMETER T IS SHIFTED 1627* 1628 T=SQRT(TL*TU) 1629 T=MAX(T,TL+0.1D 0*(TU-TL)) 1630 T=MIN(T,TL+0.9D 0*(TU-TL)) 1631 END IF 1632 ALF=ETA2 1633 CALL MXSPCT(NF,MM,MH,MMAX,H,JH,PSL,ITERM) 1634 IF (ITERM.NE.0) THEN 1635 GO TO 13250 1636 END IF 1637* 1638* SPARSE GILL-MURRAY DECOMPOSITION 1639* 1640 CALL MXSPCA(NF,MM,MH,H,IH,JH,T) 1641 CALL MXSPCF(NF,H(MM+1),PSL,JH(MM+1),WN11,WN12,GO,INF,ALF,RHO) 1642 NDEC=NDEC+1 1643 IF (INF.GT.0) THEN 1644* 1645* NEW ESTIMATION E IS COMPUTED (THE MATRIX IS NOT POSITIVE DEFINITE) 1646* 1647 IF (E.GE.TU) THEN 1648 ITERD=-2 1649 GO TO 13250 1650 ELSE 1651 MODE=2 1652 CALL MXSPCD(NF,H(MM+1),PSL,JH(MM+1),S,INF) 1653 CALL MXVSBP(NF,PERM,S,GO) 1654 E=MAX(E,T-ALF/MXUDOT(NF,S,S,IX,KBF)) 1655 NRED=NRED+1 1656 GO TO 13220 1657 END IF 1658 ELSE 1659* 1660* STEP S IS COMPUTED 1661* 1662 CALL MXUNEG(NF,G,S,IX,KBF) 1663 CALL MXVSFP(NF,PERM,S,GO) 1664 CALL MXSPCB(NF,H(MM+1),PSL,JH(MM+1),S,0) 1665 CALL MXVSBP(NF,PERM,S,GO) 1666 SNORM=SQRT(MXUDOT(NF,S,S,IX,KBF)) 1667 MODE=1 1668 END IF 1669 IF (TU-TL.LE.1.0D-8) THEN 1670* 1671* INTERVAL IS TOO SMALL 1672* 1673 IF (T.NE.0.0D 0) THEN 1674 ITERD=5 1675 ELSE 1676 ITERD=1 1677 END IF 1678 GO TO 13240 1679 ELSE IF (NRED.GE.20) THEN 1680* 1681* MAXIMUM NUMBER OF OLC REDUCTIONS 1682* 1683 ITERD=6 1684 GO TO 13240 1685 ELSE IF (SNORM.GT.DEL2*XDEL) THEN 1686* 1687* STEP IS TOO LARGE 1688* 1689 TL=MAX(TL,T) 1690 GO TO 13230 1691 ELSE IF (SNORM.LT.DEL1*XDEL) THEN 1692 IF (T.NE.0.0D 0) THEN 1693* 1694* STEP IS TOO SMAL 1695* 1696 TU=MIN(TU,T) 1697 ELSE 1698* 1699* STEP IS ACCEPTABLE 1700* 1701 ITERD=1 1702 GO TO 13240 1703 END IF 1704 ELSE 1705 ITERD=3 1706 GO TO 13240 1707 END IF 1708* 1709* TRYING TO USE BOUNDARY STEP 1710* 1711 CALL MXSPCN(NF,H(MM+1),PSL,JH(MM+1),XO,RHO,1) 1712 CALL MXVSBP(NF,PERM,XO,GO) 1713 RHO1=MXUDOT(NF,XO,S,IX,KBF) 1714 RHO2=MXUDOT(NF,XO,XO,IX,KBF) 1715 CALL PNSTEP(XDEL,SNORM,ABS(RHO1),RHO2,ALF) 1716 CON=(1.0D 0-DEL1)*(1.0D 0+DEL1) 1717 IF (ALF*ALF*RHO.LE.CON*(T*XDEL*XDEL-MXUDOT(NF,G,S,IX,KBF))) THEN 1718 IF (RHO1.LT.0.0D 0) ALF=-ALF 1719 CALL MXUDIR(NF,ALF,XO,S,S,IX,KBF) 1720 SNORM=XDEL 1721 ITERD=3 1722 GO TO 13240 1723 ELSE 1724 E=MAX(E,T-RHO) 1725 END IF 172613230 CONTINUE 1727 IF (GNORM.LE.0.0D 0) THEN 1728 T=E 1729 ELSE 1730* 1731* NEW T IS COMPUTED USING ONE STEP OF THE NEWTON METHOD FOR 1732* NONLINEAR EQUATION 1733* 1734 CALL MXUCOP(NF,S,XO,IX,KBF) 1735 CALL MXVSFP(NF,PERM,XO,GO) 1736 CALL MXSPCB(NF,H(MM+1),PSL,JH(MM+1),XO,1) 1737 T=T+(SNORM*SNORM/MXSPCP(NF,H(MM+1),PSL,XO))*(SNORM-XDEL)/XDEL 1738 CALL MXVSBP(NF,PERM,XO,GO) 1739 END IF 1740 NRED=NRED+1 1741 GO TO 13220 174213240 CONTINUE 1743 PP=MXSSMQ(NF,H,IH,JH,S,S)*0.5D 0 174413250 CONTINUE 1745 IF (KD.GT.0) P=MXUDOT(NF,G,S,IX,KBF) 1746 RETURN 1747 END 1748* SUBROUTINE PDSLM1 ALL SYSTEMS 01/09/22 1749* PURPOSE : 1750* DIRECTION DETERMINATION FOR LINE SEARCH USING DIRECT MATRIX 1751* DECOMPOSITIONS. 1752* 1753* PARAMETERS : 1754* II NF NUMBER OF VARIABLES. 1755* II MMAX MAXIMUM DIMENSION OF THE SPARSE TABLEAU. 1756* II MH POINTER OBTAINED BY THE SUBROUTINE MXSPCC. 1757* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE 1758* X(I) IS UNBOUNDED. IX(I)=1-LOWER BOUND XL(I).LE.X(I). 1759* IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND 1760* XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED. 1761* RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 1762* RA H(MMAX) NONZERO ELEMENTS OF THE APPROXIMATION OF THE SPARSE 1763* HESSIAN MATRIX TOGETHER WITH AN ADDITIONAL SPACE USED FOR 1764* THE NUMERICAL DIFFERENTIATION. 1765* II IH(NF+1) POINTERS OF DIAGONAL ELEMENTS OF THE MATRIX H. 1766* IU JH(MMAX) INDICES OF NONZERO ELEMENTS OF THE MATRIX H 1767* TOGETHER WITH AN ADDITIONAL SPACE USED FOR THE NUMERICAL 1768* DIFFERENTIATION. 1769* RO S(NF) DIRECTION VECTOR. 1770* RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. 1771* II PSL(NF+1) POINTER VECTOR OF THE COMPACT FORM OF THE TRIANGULAR 1772* FACTOR OF THE HESSIAN APPROXIMATION. 1773* IA PERM(NF) PERMUTATION VECTOR. 1774* IA WN11(NF+1) AUXILIARY VECTOR. 1775* IA WN12(NF+1) AUXILIARY VECTOR. 1776* RO GNORM NORM OF THE GRADIENT VECTOR. 1777* RO SNORM NORM OF THE DIRECTION VECTOR. 1778* RI ETA2 TOLERANCE FOR POSITIVE DEFINITENESS. 1779* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. 1780* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. 1781* IU IDEC DECOMPOSITION INDICATOR. IDEC=0-NO DECOMPOSITION. 1782* IU NDEC NUMBER OF MATRIX DECOMPOSITIONS. 1783* II ITERD CAUSE OF TERMINATION IN THE DIRECTION DETERMINATION. 1784* ITERD<0-BAD DECOMPOSITION. ITERD=0-DESCENT DIRECTION. 1785* ITERD=1-NEWTON LIKE STEP. ITERD=2-INEXACT NEWTON LIKE STEP. 1786* ITERD=3-BOUNDARY STEP. ITERD=4-DIRECTION WITH THE NEGATIVE 1787* CURVATURE. ITERD=5-MARQUARDT STEP. 1788* IO ITERM VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. 1789* ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN 1790* MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS. 1791* ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN 1792* MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS. 1793* ITERM=3-IF F WAS LESS THAN OR EQUAL TO TOLB. 1794* ITERM=4-IF GMAX WAS LESS THAN OR EQUAL TO TOLG. 1795* ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, 1796* BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. 1797* ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. 1798* ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. 1799* VALUES ITERM<=-40 DETECT A LACK OF SPACE. 1800* 1801* SUBPROGRAMS USED : 1802* S MXSPCB BACK SUBSTITUTION USING THE SPARSE DECOMPOSITION 1803* OBTAINED BY MXSPCF. 1804* S MXSPCF GILL-MURRAY DECOMPOSITION OD A SPARSE SYMMETRIC MATRIX. 1805* S MXSPCT COPYING A SPARSE SYMMETRIC MATRIX INTO THE PERMUTED 1806* FACTORIZED COMPACT SCHEME. 1807* RF MXUDOT DOT PRODUCT OF TWO VECTORS. 1808* S MXUNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. 1809* S MXVSBP INVERSE PERMUTATION OF A VECTOR 1810* S MXVSFP PERMUTATION OF A VECTOR. 1811* 1812 SUBROUTINE PDSLM1(NF,MMAX,MH,IX,G,H,IH,JH,S,XO,PSL,PERM,WN11, 1813 & WN12,GNORM,SNORM,ETA2,KBF,IDEC,NDEC,ITERD,ITERM) 1814 INTEGER NF,MMAX,MH,IX(*),IH(*),JH(*),PSL(*),PERM(*),WN11(*), 1815 & WN12(*),KBF,IDEC,NDEC,ITERD,ITERM 1816 DOUBLE PRECISION G(*),H(*),S(*),XO(*),GNORM,SNORM,ETA2 1817 INTEGER MM,INF 1818 DOUBLE PRECISION ALF,BET 1819 DOUBLE PRECISION MXUDOT 1820* 1821* DIRECTION DETERMINATION 1822* 1823 IF (IDEC.LT.0) IDEC=0 1824 MM=IH(NF+1)-1 1825 IF (IDEC.EQ.0) THEN 1826 CALL MXSPCT(NF,MM,MH,MMAX,H,JH,PSL,ITERM) 1827 IF (ITERM.NE.0) RETURN 1828* 1829* SPARSE GILL-MURRAY DECOMPOSITION 1830* 1831 ALF=ETA2 1832 CALL MXSPCF(NF,H(MM+1),PSL,JH(MM+1),WN11,WN12,XO,INF,ALF,BET) 1833 NDEC=NDEC+1 1834 IDEC=1 1835 ELSE IF (IDEC.EQ.1) THEN 1836 ELSE 1837 ITERD=-1 1838 RETURN 1839 END IF 1840 GNORM=SQRT(MXUDOT(NF,G,G,IX,KBF)) 1841* 1842* NEWTON LIKE STEP 1843* 1844 CALL MXUNEG(NF,G,S,IX,KBF) 1845 CALL MXVSFP(NF,PERM,S,XO) 1846 CALL MXSPCB(NF,H(MM+1),PSL,JH(MM+1),S,0) 1847 CALL MXVSBP(NF,PERM,S,XO) 1848 ITERD=1 1849 SNORM=SQRT(MXUDOT(NF,S,S,IX,KBF)) 1850 RETURN 1851 END 1852* SUBROUTINE PDSLM3 ALL SYSTEMS 01/09/22 1853* PURPOSE : 1854* DIRECTION DETERMINATION FOR LINE SEARCH USING CONJUGATE GRADIENT 1855* ITERATIONS. 1856* 1857* PARAMETERS : 1858* II NF NUMBER OF VARIABLES. 1859* II M NUMBER OF NONZERO ELEMENTS IN THE HESSIAN MATRIX. 1860* II MMAX MAXIMUM DIMENSION OF THE SPARSE TABLEAU. 1861* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. IX(I)=0-VARIABLE 1862* X(I) IS UNBOUNDED. IX(I)=1-LOWER BOUND XL(I).LE.X(I). 1863* IX(I)=2-UPPER BOUND X(I).LE.XU(I). IX(I)=3-TWO SIDE BOUND 1864* XL(I).LE.X(I).LE.XU(I). IX(I)=5-VARIABLE X(I) IS FIXED. 1865* RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 1866* RA H(MMAX) NONZERO ELEMENTS OF THE APPROXIMATION OF THE SPARSE 1867* HESSIAN MATRIX TOGETHER WITH AN ADDITIONAL SPACE USED FOR 1868* THE NUMERICAL DIFFERENTIATION. 1869* II IH(NF+1) POINTERS OF DIAGONAL ELEMENTS OF THE MATRIX H. 1870* IU JH(MMAX) INDICES OF NONZERO ELEMENTS OF THE MATRIX H 1871* TOGETHER WITH AN ADDITIONAL SPACE USED FOR THE NUMERICAL 1872* DIFFERENTIATION. 1873* RO S(NF) DIRECTION VECTOR. 1874* RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. 1875* RI GO(NF) GRADIENTS DIFFERENCE. 1876* RA XS(NF) AUXILIARY VECTOR. 1877* RA IW(NF+1) AUXILIARY VECTOR. 1878* RO GNORM NORM OF THE GRADIENT VECTOR. 1879* RO SNORM NORM OF THE DIRECTION VECTOR. 1880* RI ETA2 TOLERANCE FOR POSITIVE DEFINITENESS. 1881* RI ETA9 MAXIMUM FOR REAL NUMBERS. 1882* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. 1883* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. 1884* II MOS2 TYPE OF PRECONDITIONING. MOS2=1-PRECONDITIONING IS NOT 1885* USED. MOS2=2-PRECONDITIONING BY THE INCOMPLETE GILL-MURRAY 1886* DECOMPOSITION. MOS2=3-PRECONDITIONING BY THE INCOMPLETE 1887* GILL-MURRAY DECOMPOSITION WITH A PRELIMINARY SOLUTION OF 1888* THE PRECONDITIONED SYSTEM WHICH IS USED IF IT SATISFIES 1889* THE TERMINATION CRITERION. 1890* IU IDEC DECOMPOSITION INDICATOR. IDEC=0-NO DECOMPOSITION. 1891* IU NDEC NUMBER OF MATRIX DECOMPOSITIONS. 1892* II NIT NUMBER OF OUTER ITERATIONS. 1893* IU NIN NUMBER OF INNER CONJUGATE GRADIENT ITERATIONS. 1894* II ITERD CAUSE OF TERMINATION IN THE DIRECTION DETERMINATION. 1895* ITERD<0-BAD DECOMPOSITION. ITERD=0-DESCENT DIRECTION. 1896* ITERD=1-NEWTON LIKE STEP. ITERD=2-INEXACT NEWTON LIKE STEP. 1897* ITERD=3-BOUNDARY STEP. ITERD=4-DIRECTION WITH THE NEGATIVE 1898* CURVATURE. ITERD=5-MARQUARDT STEP. 1899* IO ITERM VARIABLE THAT INDICATES THE CAUSE OF TERMINATION. 1900* ITERM=1-IF ABS(X-XO) WAS LESS THAN OR EQUAL TO TOLX IN 1901* MTESX (USUALLY TWO) SUBSEQUENT ITERATIONS. 1902* ITERM=2-IF ABS(F-FO) WAS LESS THAN OR EQUAL TO TOLF IN 1903* MTESF (USUALLY TWO) SUBSEQUENT ITERATIONS. 1904* ITERM=3-IF F WAS LESS THAN OR EQUAL TO TOLB. 1905* ITERM=4-IF GMAX WAS LESS THAN OR EQUAL TO TOLG. 1906* ITERM=6-IF THE TERMINATION CRITERION WAS NOT SATISFIED, 1907* BUT THE SOLUTION OBTAINED IS PROBABLY ACCEPTABLE. 1908* ITERM=11-IF NIT EXCEEDED MIT. ITERM=12-IF NFV EXCEEDED MFV. 1909* ITERM=13-IF NFG EXCEEDED MFG. ITERM<0-IF THE METHOD FAILED. 1910* VALUES ITERM<=-40 DETECT A LACK OF SPACE. 1911* 1912* SUBPROGRAMS USED : 1913* S MXSPTB BACK SUBSTITUTION AFTER THE GILL-MURRAY DECOMPOSITION. 1914* S MXSPTF INCOMPLETE GILL-MURRAY DECOMPOSITION. 1915* S MXSSMD MATRIX-VECTOR PRODUCT FOLLOWED BY THE ADDITION OF A 1916* SCALED VECTOR. 1917* S MXSSMM MATRIX-VECTOR PRODUCT. 1918* S MXUDIR VECTOR AUGMENTED BY THE SCALED VECTOR. 1919* RF MXUDOT DOT PRODUCT OF TWO VECTORS. 1920* S MXVCOP COPYING OF A VECTOR. 1921* S MXVNEG COPYING OF A VECTOR WITH CHANGE OF THE SIGN. 1922* S MXVSET INITIATION OF A VECTOR. 1923* 1924 SUBROUTINE PDSLM3(NF,M,MMAX,IX,G,H,IH,JH,S,XO,GO,XS,IW,GNORM, 1925 & SNORM,ETA2,ETA9,KBF,MOS2,IDEC,NDEC,NIT,NIN,ITERD,ITERM) 1926 INTEGER NF,M,MMAX,IX(*),IH(*),JH(*),IW(*),KBF,MOS2,IDEC,NDEC, 1927 & NIT,NIN,ITERD,ITERM 1928 DOUBLE PRECISION G(*),H(*),S(*),XO(*),GO(*),XS(*),GNORM,SNORM, 1929 & ETA2,ETA9 1930 INTEGER NOS2,NRED,MMX,INF 1931 DOUBLE PRECISION PAR,ALF,EPS,RHO,RHO1,RHO2,SIG 1932 DOUBLE PRECISION MXUDOT 1933 SAVE EPS 1934* 1935* DIRECTION DETERMINATION 1936* 1937 IF (NIT.LE.1) THEN 1938 EPS=0.9D 0 1939 END IF 1940 NOS2=MOS2-1 1941 IF (IDEC.LT.0) IDEC=0 1942 IF (IDEC.NE.0.AND.IDEC.NE.1) THEN 1943 ITERD=-1 1944 RETURN 1945 ELSE IF (IDEC.EQ.0) THEN 1946 IF (MOS2.GT.1) THEN 1947* 1948* INCOMPLETE GILL-MURRAY DECOMPOSITION 1949* 1950 ALF=ETA2 1951 IF (2*M.GE.MMAX) THEN 1952 ITERM=-48 1953 RETURN 1954 END IF 1955 CALL MXVCOP(M,H,H(M+1)) 1956 CALL MXSPTF(NF,H(M+1),IH,JH,IW,INF,ALF,SIG) 1957 IF (INF+10.LT.0) THEN 1958 ITERM=-48 1959 RETURN 1960 END IF 1961 IF (INF.NE.0) NOS2=0 1962 NDEC=NDEC+1 1963 IDEC=1 1964 END IF 1965 END IF 1966 RHO1=MXUDOT(NF,G,G,IX,KBF) 1967 GNORM=SQRT(RHO1) 1968 PAR=MIN(EPS,SQRT(GNORM)) 1969 IF (PAR.GT.1.0D-2) THEN 1970 PAR=MIN(PAR,1.0D 0/DBLE(NIT)) 1971 END IF 1972 PAR=PAR*PAR 1973 IF (MOS2.GT.2) THEN 1974* 1975* PRELIMINARY INEXACT SOLUTION 1976* 1977 CALL MXVNEG(NF,G,XO) 1978 IF (NOS2.NE.0) THEN 1979 CALL MXSPTB(NF,H(M+1),IH,JH,XO,0) 1980 CALL MXVCOP(NF,XO,S) 1981 CALL MXSSMD(NF,H,IH,JH,S,1.0D 0,G,GO) 1982 IF (MXUDOT(NF,GO,GO,IX,KBF).LE.1.0D-2*PAR*RHO1) THEN 1983 SNORM=SQRT(MXUDOT(NF,S,S,IX,KBF)) 1984 ITERD=2 1985 RETURN 1986 END IF 1987 END IF 1988 END IF 1989* 1990* CG INITIATION 1991* 1992 RHO=RHO1 1993 SNORM=0.0D 0 1994 CALL MXVSET(NF,0.0D 0,S) 1995 CALL MXVNEG(NF,G,XS) 1996 IF (NOS2.EQ.0) THEN 1997 CALL MXVNEG(NF,G,XO) 1998 ELSE IF (MOS2.GT.2) THEN 1999 RHO=MXUDOT(NF,XS,XO,IX,KBF) 2000 ELSE 2001 CALL MXVNEG(NF,G,XO) 2002 CALL MXSPTB(NF,H(M+1),IH,JH,XO,0) 2003 RHO=MXUDOT(NF,XS,XO,IX,KBF) 2004 END IF 2005C SIG=RHO 2006 MMX=NF+3 2007 DO 10 NRED=1,MMX 2008 CALL MXSSMM(NF,H,IH,JH,XO,GO) 2009 ALF=MXUDOT(NF,XO,GO,IX,KBF) 2010 IF (ALF.LE.1.0D 0/ETA9) THEN 2011C IF (ALF.LE.1.0D-8*SIG) THEN 2012* 2013* CG FAILS (THE MATRIX IS NOT POSITIVE DEFINITE) 2014* 2015 IF (NRED.EQ.1) THEN 2016 CALL MXVNEG(NF,G,S) 2017 SNORM=GNORM 2018 END IF 2019 ITERD=0 2020 RETURN 2021 ELSE 2022 ITERD=2 2023 END IF 2024* 2025* CG STEP 2026* 2027 ALF=RHO/ALF 2028 CALL MXUDIR(NF, ALF,XO,S,S,IX,KBF) 2029 CALL MXUDIR(NF,-ALF,GO,XS,XS,IX,KBF) 2030 NIN=NIN+1 2031 RHO2=MXUDOT(NF,XS,XS,IX,KBF) 2032 SNORM=SQRT(MXUDOT(NF,S,S,IX,KBF)) 2033 IF (RHO2.LE.PAR*RHO1) RETURN 2034 IF (NRED.GE.MMX) RETURN 2035 IF (NOS2.NE.0) THEN 2036 CALL MXVCOP(NF,XS,GO) 2037 CALL MXSPTB(NF,H(M+1),IH,JH,GO,0) 2038 RHO2=MXUDOT(NF,XS,GO,IX,KBF) 2039 ALF=RHO2/RHO 2040 CALL MXUDIR(NF,ALF,XO,GO,XO,IX,KBF) 2041 ELSE 2042 ALF=RHO2/RHO 2043 CALL MXUDIR(NF,ALF,XO,XS,XO,IX,KBF) 2044 END IF 2045 RHO=RHO2 2046C SIG=RHO2+ALF*ALF*SIG 2047 10 CONTINUE 2048 RETURN 2049 END 2050* SUBROUTINE PF1HS2 ALL SYSTEMS 99/12/01 2051* PURPOSE : 2052* NUMERICAL COMPUTATION OF THE HESSIAN MATRIX OF THE MODEL FUNCTION 2053* USING ITS GRADIENTS - SPARSE VERSION USING DIRECT COLOURING METHOD. 2054* 2055* PARAMETERS : 2056* II NF NUMBER OF VARIABLES. 2057* II ML SIZE OF THE COMPACT FACTOR. 2058* II M NUMBER OF NONZERO ELEMENTS OF THE SPARSE HESSIAN MATRIX. 2059* RI X(NF) VECTOR OF VARIABLES. 2060* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. 2061* RA XO(NF) AUXILIARY VECTOR. 2062* RO HF(M) HESSIAN MATRIX OF THE MODEL FUNCTION. 2063* IU IH(NF+1) POINTER VECTOR OF SPARSE HESSIAN MATRIX. 2064* IU JH(M) INDEX VECTOR OF THE HESSIAN MATRIX. 2065* RI GF(NF) GRADIENT OF THE MODEL FUNCTION. 2066* RA GO(NF) AUXILIARY VECTOR. 2067* II COL(NF) VECTOR DISCERNING GROUPS OF THE HESSIAN COLUMN OF THE 2068* SAME COLOUR. 2069* IA WN11(NF+1) AUXILIARY VECTOR. 2070* IA WN12(NF+1) AUXILIARY VECTOR. 2071* RA XS(NF) AUXILIARY VECTOR USED FOR STEP SIZES. 2072* RI FF VALUE OF THE MODEL FUNCTION. 2073* RI ETA1 PRECISION OF THE COMPUTED VALUES. 2074* II KBF TYPE OF BOUNDS. KBF=0-BOUNDS ARE NOT USED. KBF=1-ONE SIDED 2075* BOUNDS. KBF=2-TWO SIDED BOUNDS. 2076* IU ITERM TERMINATION INDICATOR. 2077* IU ISYS CONTROL PARAMETER. 2078* 2079* SUBPROGRAMS USED : 2080* S MXSTG1 WIDTHEN THE STRUCTURE. 2081* S MXSTL2 SHRINK THE STRUCTURE. 2082* S MXVCOP COPYING OF A VECTOR. 2083* S MXVSET INITIATION OF A VECTOR. 2084* 2085 SUBROUTINE PF1HS2(NF,ML,M,X,IX,XO,HF,IH,JH,GF,GO,COL,WN11, 2086 & WN12,XS,FF,ETA1,KBF,ITERM,ISYS) 2087 INTEGER NF,ML,M,IX(*),IH(*),JH(*),COL(*),WN11(*), 2088 & WN12(*),KBF,ITERM,ISYS 2089 DOUBLE PRECISION X(*),XO(*),HF(*),GF(*),GO(*),XS(*), 2090 & FF,ETA1 2091 DOUBLE PRECISION XTEMP,FTEMP,ETA 2092 INTEGER I,J,J1,K,K1,L,MX,MM,IVAR,JVAR 2093 SAVE MX,MM,IVAR,JVAR 2094 SAVE XTEMP,FTEMP,ETA 2095 IF (ITERM.NE.0) GO TO 12 2096 IF (ISYS.EQ.1) GO TO 3 2097 MM=IH(NF+1)-1 2098 IF (3*MM-NF+ML.GE.M) THEN 2099 ITERM=-45 2100 ISYS=0 2101 RETURN 2102 END IF 2103 ETA=SQRT(ETA1) 2104 FTEMP=FF 2105 CALL MXVCOP(NF,X,XO) 2106* 2107* WIDTHEN THE STRUCTURE 2108* 2109 K=2*MM-NF 2110 DO 50 I=ML+MM,1,-1 2111 JH(K+I)=JH(MM+I) 2112 50 CONTINUE 2113 CALL MXSTG1(NF,MX,IH,JH,WN12,WN11) 2114 CALL MXVSET(K,0.0D 0,HF) 2115 IVAR=1 2116 2 CONTINUE 2117 IF (IVAR.GT.NF) GO TO 870 2118 DO 200 J=IVAR,NF 2119 IF (COL(J).GE.1) THEN 2120 GO TO 200 2121 ELSE 2122 JVAR=J 2123 GO TO 300 2124 END IF 2125 200 CONTINUE 2126 300 CONTINUE 2127 DO 400 J=IVAR,JVAR 2128 L=ABS(COL(J)) 2129 IF (KBF.GT.0) THEN 2130 IF (IX(L).LE.-7) GO TO 400 2131 END IF 2132* 2133* STEP SELECTION 2134* 2135 XS(L)=ETA*MAX(ABS(X(L)),1.0D 0)*SIGN(1.0D 0,X(L)) 2136 XTEMP=X(L) 2137 X(L)=XTEMP+XS(L) 2138 XS(L)=X(L)-XTEMP 2139 400 CONTINUE 2140 ISYS=1 2141 RETURN 2142 3 CONTINUE 2143* 2144* NUMERICAL DIFFERENTIATION 2145* 2146* 2147* SET AUXILIARY VECTOR DISCERNING THE SINGLETONS IN A GROUP TO ZERO 2148* 2149 DO 450 J1=1,NF 2150 WN11(J1)=0 2151 450 CONTINUE 2152* 2153* DISCERN SINGLETONS OF THE GROUP OF THE SAME COLOR. 2154* 2155 DO 600 J1=IVAR,JVAR 2156 L=ABS(COL(J1)) 2157 DO 500 K=IH(L),IH(L+1)-1 2158 K1=ABS(JH(K)) 2159 IF (WN11(K1).EQ.0) THEN 2160 WN11(K1)=J1 2161 ELSE 2162 WN11(K1)=-1 2163 END IF 2164 500 CONTINUE 2165 600 CONTINUE 2166* 2167* NUMERICAL VALUES COMPUTATION 2168* 2169 DO 800 J1=IVAR,JVAR 2170 L=ABS(COL(J1)) 2171 DO 700 K=IH(L),IH(L+1)-1 2172 K1=ABS(JH(K)) 2173 IF (WN11(K1).GT.0) THEN 2174 HF(K)=(GF(K1)-GO(K1))/XS(L) 2175 END IF 2176 700 CONTINUE 2177 800 CONTINUE 2178* 2179* SET THE ORIGINAL VALUE OF X FOR THE COMPONENTS OF THE ACTUAL COLOR. 2180* 2181 DO 850 J=IVAR,JVAR 2182 L=ABS(COL(J)) 2183 X(L)=XO(L) 2184 850 CONTINUE 2185 IVAR=JVAR+1 2186 GO TO 2 2187 870 CONTINUE 2188* 2189* MOVE THE ELEMENTS OF THE HESSIAN APPROXIMATION INTO THE UPPER 2190* TRIANGULAR PART 2191* 2192 DO 900 I=1,NF 2193 WN11(I)=WN12(I)+1 2194 900 CONTINUE 2195 DO 1100 I=1,NF 2196 IVAR=IH(I) 2197 JVAR=WN12(I)-1 2198 DO 1000 J=IVAR,JVAR 2199 K=ABS(JH(J)) 2200 L=WN11(K) 2201 IF (HF(L).EQ.0) THEN 2202 HF(L)=HF(J) 2203 ELSE IF (HF(L).NE.0.AND.HF(J).NE.0) THEN 2204 HF(L)=0.5D 0*(HF(J)+HF(L)) 2205 END IF 2206 WN11(K)=WN11(K)+1 2207 1000 CONTINUE 2208 1100 CONTINUE 2209 FF=FTEMP 2210* 2211* SHRINK THE STRUCTURE 2212* 2213 CALL MXSTL2(NF,MX,HF,IH,JH,WN12) 2214 K=2*MM-NF 2215 DO 1200 I=1,ML+MM 2216 JH(MM+I)=JH(K+I) 2217 1200 CONTINUE 2218* 2219* RETRIEVE VALUES 2220* 2221 CALL MXVCOP(NF,XO,X) 2222 12 CONTINUE 2223 ISYS=0 2224 RETURN 2225 END 2226* SUBROUTINE PFSEB4 ALL SYSTEMS 98/12/01 2227* PURPOSE : 2228* COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE PARTITIONED HESSIAN 2229* MATRIX. 2230* 2231* PARAMETERS : 2232* II NC NUMBER OF CONSTRAINTS. 2233* RU B(M) ELEMENTS OF THE SPARSE MATRIX B. 2234* IO IH(NF+1) POINTERS OF THE DIAGONAL ELEMENTS OF B. 2235* IO JH(M) INDICES OF THE NONZERO ELEMENTS OF B. 2236* II CH(MB) ELEMENTS OF THE PARTITIONED MATRIX H. 2237* II ICG(NC+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 2238* II JCG(MC) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 2239* II ICA(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. 2240* RI CZL(NC) VECTOR CONTAINING LOWER MULTIPLIERS FOR CONSTRAINTS. 2241* RI CZU(NC) VECTOR CONTAINING UPPER MULTIPLIERS FOR CONSTRAINTS. 2242* II JOB SUBJECTS OF UPDATES. JOB=0-CONSTRAINT FUNCTIONS. 2243* JOB=1-CONSTRAINT FUNCTIONS MULTIPLIED BY SIGNS OF THE 2244* LAGRANGIAN MULTIPLIERS. JOB-2-ACTIVE TERMS OF THE LAGRANGIAN 2245* FUNCTION. JOB-3-ALL TERMS OF THE LAGRANGIAN FUNCTION. 2246* 2247 SUBROUTINE PFSEB4(NC,B,IH,JH,CH,ICG,JCG,ICA,CZL,CZU,JOB) 2248 INTEGER NC,IH(*),JH(*),ICG(*),JCG(*),ICA(*),JOB 2249 DOUBLE PRECISION B(*),CH(*),CZL(*),CZU(*) 2250 INTEGER I,II,IC,J,JJ,JC,JF,K,KK,L,LL,KC 2251 DOUBLE PRECISION TEMP 2252 KK=0 2253 DO 7 KC=1,NC 2254 IF (JOB.LE.1) THEN 2255 LL=ABS(ICA(KC)) 2256 IF (LL.EQ.3.OR.LL.EQ.4) THEN 2257 TEMP= CZU(KC)-CZL(KC) 2258 ELSE IF (LL.EQ.1) THEN 2259 TEMP=-CZL(KC) 2260 ELSE IF (LL.EQ.2) THEN 2261 TEMP= CZU(KC) 2262 ELSE IF (LL.EQ.5) THEN 2263 TEMP= CZL(KC) 2264 END IF 2265 IF (JOB.EQ.1) TEMP=ABS(TEMP) 2266 ELSE IF (JOB.EQ.2) THEN 2267 IF (ICA(KC).GE.0) GO TO 7 2268 TEMP=1.0D 0 2269 ELSE 2270 TEMP=1.0D 0 2271 END IF 2272 II=ICG(KC) 2273 L=ICG(KC+1)-II 2274 DO 6 IC=1,L 2275 KK=KK+IC 2276 I=JCG(II) 2277 IF (I.LE.0) GO TO 5 2278 JF=IH(I) 2279 JJ=II 2280 K=KK 2281 DO 4 JC=IC,L 2282 J=JCG(JJ) 2283 IF (J.LE.0) GO TO 3 2284 2 IF (JH(JF).LT.J) THEN 2285 JF=JF+1 2286 GO TO 2 2287 END IF 2288 B(JF)=B(JF)+TEMP*CH(K) 2289 3 K=K+JC 2290 JJ=JJ+1 2291 4 CONTINUE 2292 5 II=II+1 2293 6 CONTINUE 2294 7 CONTINUE 2295 RETURN 2296 END 2297* SUBROUTINE PFSEB5 ALL SYSTEMS 06/12/01 2298* PURPOSE : 2299* COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE PARTITIONED HESSIAN 2300* MATRIX. 2301* 2302* PARAMETERS : 2303* II NC NUMBER OF CONSTRAINTS. 2304* RU B(M) ELEMENTS OF THE SPARSE MATRIX B. 2305* IO IH(NF+1) POINTERS OF THE DIAGONAL ELEMENTS OF B. 2306* IO JH(M) INDICES OF THE NONZERO ELEMENTS OF B. 2307* II CH(MB) ELEMENTS OF THE PARTITIONED MATRIX H. 2308* II ICG(NC+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 2309* II JCG(MC) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 2310* RI CZ(NC) VECTOR CONTAINING LAGRANGE MULTIPLIERS FOR CONSTRAINTS. 2311* II JOB SUBJECTS OF UPDATES. JOB=0-CONSTRAINT FUNCTIONS. 2312* JOB=1-CONSTRAINT FUNCTIONS MULTIPLIED BY SIGNS OF THE 2313* LAGRANGIAN MULTIPLIERS. JOB-2-ACTIVE TERMS OF THE LAGRANGIAN 2314* FUNCTION. JOB-3-ALL TERMS OF THE LAGRANGIAN FUNCTION. 2315* 2316 SUBROUTINE PFSEB5(NC,B,IH,JH,CH,ICG,JCG,CZ,JOB) 2317 INTEGER NC,IH(*),JH(*),ICG(*),JCG(*),JOB 2318 DOUBLE PRECISION B(*),CH(*),CZ(*) 2319 INTEGER I,II,IC,J,JJ,JC,JF,K,KK,L,KC 2320 DOUBLE PRECISION TEMP 2321 KK=0 2322 DO 7 KC=1,NC 2323 IF (JOB.EQ.0) THEN 2324 TEMP=CZ(KC) 2325 ELSE IF (JOB.EQ.1) THEN 2326 TEMP=ABS(CZ(KC)) 2327 ELSE 2328 TEMP=1.0D 0 2329 END IF 2330 II=ICG(KC) 2331 L=ICG(KC+1)-II 2332 DO 6 IC=1,L 2333 KK=KK+IC 2334 I=JCG(II) 2335 IF (I.LE.0) GO TO 5 2336 JF=IH(I) 2337 JJ=II 2338 K=KK 2339 DO 4 JC=IC,L 2340 J=JCG(JJ) 2341 IF (J.LE.0) GO TO 3 2342 2 IF (JH(JF).LT.J) THEN 2343 JF=JF+1 2344 GO TO 2 2345 END IF 2346 B(JF)=B(JF)+TEMP*CH(K) 2347 3 K=K+JC 2348 JJ=JJ+1 2349 4 CONTINUE 2350 5 II=II+1 2351 6 CONTINUE 2352 7 CONTINUE 2353 RETURN 2354 END 2355* SUBROUTINE PFSED3 ALL SYSTEMS 07/12/01 2356* PURPOSE : 2357* COMPRESSED SPARSE STRUCTURE OF THE HESSIAN MATRIX IS COMPUTED FROM 2358* THE COORDINATE FORM. 2359* 2360* PARAMETERS : 2361* II NF NUMBER OF VARIABLES. 2362* II M NUMBER OF NONZERO ELEMENTS IN THE UPPER PART OF THE SPARSE 2363* HESSIAN MATRIX. 2364* IU IH(M+NF) ON INPUT ROW INDICES OF NONZERO ELEMENTS IN THE FIELD 2365* H. ON OUTPUT POSITIONS OF DIAGONAL ELEMENTS IN THE FIELD H. 2366* II JH(M+NF) COLUMN INDICES OF NONZERO ELEMENTS IN THE FIELD H. 2367* IO IER ERROR MESAGE. IER=0-THE STANDARD INPUT DATA ARE CORRECT. 2368* IER=1-ERROR IN THE ARRAY IH. IER=2-ERROR IN THE ARRAY JH. 2369* 2370 SUBROUTINE PFSED3(NF,M,IH,JH,IER) 2371 INTEGER NF,M,IH(*),JH(*),IER 2372 INTEGER I,J,K,L,LL 2373 IER=0 2374 DO 1 J=1,M 2375 IF (IH(J).GT.JH(J)) THEN 2376 K=IH(J) 2377 IH(J)=JH(J) 2378 JH(J)=K 2379 END IF 2380 1 CONTINUE 2381 DO 2 I=1,NF 2382 IH(M+I)=I 2383 JH(M+I)=I 2384 2 CONTINUE 2385 CALL MXVSR7(M+NF,IH,JH) 2386 IF (IH(1).LT.1.OR.IH(M+NF).GT.NF) THEN 2387 IER=1 2388 RETURN 2389 END IF 2390 K=1 2391 DO 3 J=1,M+NF 2392 IF (IH(J).EQ.K) THEN 2393 IH(K)=J 2394 K=K+1 2395 END IF 2396 3 CONTINUE 2397 IH(K)=J 2398 LL=0 2399 DO 5 I=1,NF 2400 K=IH(I) 2401 L=IH(I+1)-K 2402 IF (L.GT.0) THEN 2403 CALL MXVSRT(L,JH(K)) 2404 IF (JH(K).LT.1.OR.JH(K+L-1).GT.NF) THEN 2405 IER=2 2406 RETURN 2407 END IF 2408 END IF 2409 IH(I)=IH(I)-LL 2410 DO 4 J=1,L 2411 IF (J.GT.1.AND.JH(K).EQ.JH(K-1)) THEN 2412 LL=LL+1 2413 ELSE 2414 JH(K-LL)=JH(K) 2415 END IF 2416 K=K+1 2417 4 CONTINUE 2418 5 CONTINUE 2419 IH(NF+1)=IH(NF+1)-LL 2420 M=IH(NF+1)-1 2421 RETURN 2422 END 2423* SUBROUTINE PFSET2 ALL SYSTEMS 97/12/01 2424* PURPOSE : 2425* COMPUTATION OF THE NUMBER OF NONZERO ELEMENTS OF THE SPARSE 2426* HESSIAN MATRIX STORED IN THE BLOCKED FORM. 2427* 2428* PARAMETERS : 2429* II NA NUMBER OF APPROXIMATED FUNCTIONS. 2430* II MB NUMBER OF NONZERO ELEMENTS OF THE PARTITIONED HESSIAN MATRIX 2431* II MC MAXIMUM NUMBER OF ELEMENTS OF THE PARTIAL HESSIAN MATRIX. 2432* RI IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE SPARSE 2433* JACOBIAN MATRIX. 2434* 2435 SUBROUTINE PFSET2(NA,MB,MC,IAG) 2436 INTEGER NA,MB,MC,IAG(*) 2437 INTEGER K,L,KA 2438 MB=0 2439 MC=0 2440 DO 1 KA=1,NA 2441 K=IAG(KA) 2442 L=IAG(KA+1)-K 2443 MB=MB+L*(L+1)/2 2444 MC=MAX(MC,L*(L+1)/2) 2445 1 CONTINUE 2446 RETURN 2447 END 2448* SUBROUTINE PFSET3 ALL SYSTEMS 97/12/01 2449* PURPOSE : 2450* COMPUTATION OF THE SPARSE STRUCTURE OF THE HESSIAN MATRIX FROM THE 2451* SPARSE STRUCTURE OF THE JACOBIAN MATRIX. 2452* 2453* PARAMETERS : 2454* II NF NUMBER OF VARIABLES. 2455* II NA NUMBER OF APPROXIMATED FUNCTIONS. 2456* IO M NUMBER OF NONZERO ELEMENTS OF THE HESSIAN MATRIX. 2457* II MMAX DECLARED LENGHT OF THE ARRAYS H AND JH. 2458* IO IH(NF+1) POINTERS OF THE DIAGONAL ELEMENTS OF H. 2459* IO JH(M) INDICES OF THE NONZERO ELEMENTS OF H. 2460* II IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 2461* II JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 2462* IU ITERM TERMINATION INDICATOR. 2463* 2464 SUBROUTINE PFSET3(NF,NA,M,MMAX,IH,JH,IAG,JAG,ITERM) 2465 INTEGER NF,NA,M,MMAX,IH(*),JH(*),IAG(*),JAG(*),ITERM 2466 INTEGER I,J,JF,JA,K,LF,LA,KA 2467 M=IH(NF+1)-1 2468 IF (M.GT.MMAX) THEN 2469 ITERM=-40 2470 RETURN 2471 END IF 2472 DO 7 KA=1,NA 2473 LA=IAG(KA+1)-1 2474 DO 6 K=IAG(KA),LA 2475 I=JAG(K) 2476 JF=IH(I) 2477 LF=IH(I+1)-1 2478 DO 5 JA=K,LA 2479 J=JAG(JA) 2480 2 IF (JH(JF).LT.J.AND.JF.LE.LF) THEN 2481 JF=JF+1 2482 IF (JF.LE.LF) GO TO 2 2483 END IF 2484 IF (JH(JF).GT.J .OR.JF.GT.LF) THEN 2485 DO 3 J=I+1,NF+1 2486 IH(J)=IH(J)+1 2487 3 CONTINUE 2488 DO 4 J=M,JF,-1 2489 JH(J+1)=JH(J) 2490 4 CONTINUE 2491 JH(JF)=JAG(JA) 2492 JF=JF+1 2493 LF=LF+1 2494 M=M+1 2495 IF (M.GT.MMAX) THEN 2496 ITERM=-40 2497 RETURN 2498 END IF 2499 END IF 2500 5 CONTINUE 2501 6 CONTINUE 2502 7 CONTINUE 2503 RETURN 2504 END 2505* SUBROUTINE PFSET4 ALL SYSTEMS 98/12/01 2506* PURPOSE : 2507* COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE PARTITIONED HESSIAN 2508* MATRIX. 2509* 2510* PARAMETERS : 2511* II NA NUMBER OF APPROXIMATED FUNCTIONS. 2512* RU B(M) ELEMENTS OF THE SPARSE MATRIX B. 2513* IO IH(NF+1) POINTERS OF THE DIAGONAL ELEMENTS OF B. 2514* IO JH(M) INDICES OF THE NONZERO ELEMENTS OF B. 2515* II AH(MB) ELEMENTS OF THE PARTITIONED MATRIX H. 2516* II IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 2517* II JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 2518* 2519 SUBROUTINE PFSET4(NA,B,IH,JH,AH,IAG,JAG) 2520 INTEGER NA,IH(*),JH(*),IAG(*),JAG(*) 2521 DOUBLE PRECISION B(*),AH(*) 2522 INTEGER I,II,IA,J,JJ,JA,JF,K,KK,L,KA 2523 KK=0 2524 DO 7 KA=1,NA 2525 II=IAG(KA) 2526 L=IAG(KA+1)-II 2527 DO 6 IA=1,L 2528 KK=KK+IA 2529 I=JAG(II) 2530 IF (I.LE.0) GO TO 5 2531 JF=IH(I) 2532 JJ=II 2533 K=KK 2534 DO 4 JA=IA,L 2535 J=JAG(JJ) 2536 IF (J.LE.0) GO TO 3 2537 2 IF (JH(JF).LT.J) THEN 2538 JF=JF+1 2539 GO TO 2 2540 END IF 2541 B(JF)=B(JF)+AH(K) 2542 3 K=K+JA 2543 JJ=JJ+1 2544 4 CONTINUE 2545 5 II=II+1 2546 6 CONTINUE 2547 7 CONTINUE 2548 RETURN 2549 END 2550* FUNCTION PNFUZ1 ALL SYSTEMS 01/09/22 2551* PURPOSE : 2552* COMPUTATION OF LOWER AND UPPER LAGRANGE MULTIPLIERS. 2553* 2554* PARAMETERS : 2555* RO Z SLACK VARIABLE IN THE NONLINEAR PROGRAMMING FORMULATION OF 2556* A MINIMAX PROBLEM. 2557* II NA NUMBER OF APPROXIMATED FUNCTIONS. 2558* RI RPF3 BARRIER PARAMETER. 2559* RO AF(NA) VECTOR CONTAINING VALUES OF THE APPROXIMATED 2560* FUNCTIONS. 2561* RA AZL(NA) VECTOR OF LOWER LAGRANGE MULTIPLIERS. 2562* RA AZU(NA) VECTOR OF UPPER LAGRANGE MULTIPLIERS. 2563* II IEXT TYPE OF MINIMAX. IEXT<0-MINIMIZATION OF THE MAXIMUM 2564* PARTIAL VALUE. IEXT-0-MINIMIZATION OF THE MAXIMUM PARTIAL 2565* ABSOLUTE VALUE. IEXT>0-MAXIMIZATION OF THE MINIMUM PARTIAL 2566* VALUE. 2567* 2568 FUNCTION PNFUZ1(Z,NA,RPF3,AF,AZL,AZU,IEXT) 2569 INTEGER NA,IEXT 2570 DOUBLE PRECISION Z,RPF3,AF(*),AZL(*),AZU(*),PNFUZ1 2571 INTEGER KA 2572 PNFUZ1=1.0D 0 2573 DO 1 KA=1,NA 2574 IF (IEXT.LE.0) THEN 2575 AZU(KA)=RPF3/(Z-AF(KA)) 2576 PNFUZ1=PNFUZ1-AZU(KA) 2577 END IF 2578 IF (IEXT.GE.0) THEN 2579 AZL(KA)=RPF3/(Z+AF(KA)) 2580 PNFUZ1=PNFUZ1-AZL(KA) 2581 END IF 2582 1 CONTINUE 2583 RETURN 2584 END 2585* SUBROUTINE PNINT1 ALL SYSTEMS 91/12/01 2586* PURPOSE : 2587* EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH WITH DIRECTIONAL 2588* DERIVATIVES. 2589* 2590* PARAMETERS : 2591* RI RL LOWER VALUE OF THE STEPSIZE PARAMETER. 2592* RI RU UPPER VALUE OF THE STEPSIZE PARAMETER. 2593* RI FL VALUE OF THE OBJECTIVE FUNCTION FOR R=RL. 2594* RI FU VALUE OF THE OBJECTIVE FUNCTION FOR R=RU. 2595* RI PL DIRECTIONAL DERIVATIVE FOR R=RL. 2596* RI PU DIRECTIONAL DERIVATIVE FOR R=RU. 2597* RO R VALUE OF THE STEPSIZE PARAMETER OBTAINED. 2598* II MODE MODE OF LINE SEARCH. 2599* II MTYP METHOD SELECTION. MTYP=1-BISECTION. MTYP=2-QUADRATIC 2600* INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE). 2601* MTYP=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL 2602* DERIVATIVES). MTYP=4-CUBIC INTERPOLATION. MTYP=5-CONIC 2603* INTERPOLATION. 2604* IO MERR ERROR INDICATOR. MERR=0 FOR NORMAL RETURN. 2605* 2606* METHOD : 2607* EXTRAPOLATION OR INTERPOLATION WITH STANDARD MODEL FUNCTIONS. 2608* 2609 SUBROUTINE PNINT1(RL,RU,FL,FU,PL,PU,R,MODE,MTYP,MERR) 2610 DOUBLE PRECISION RL, RU, FL, FU, PL, PU, R 2611 INTEGER MODE,MTYP,MERR,NTYP 2612 DOUBLE PRECISION A,B,C,D,DIS,DEN 2613 DOUBLE PRECISION C1L,C1U,C2L,C2U,C3L 2614 PARAMETER (C1L=1.1D 0,C1U=1.0D 3,C2L=1.0D-2,C2U=0.9D 0, 2615 & C3L=0.1D 0) 2616 MERR=0 2617 IF (MODE.LE.0) RETURN 2618 IF (PL.GE.0.0D 0) THEN 2619 MERR=2 2620 RETURN 2621 ELSE IF (RU.LE.RL) THEN 2622 MERR=3 2623 RETURN 2624 END IF 2625 DO 1 NTYP=MTYP,1,-1 2626 IF (NTYP.EQ.1) THEN 2627* 2628* BISECTION 2629* 2630 IF (MODE.EQ.1) THEN 2631 R=4.0D 0*RU 2632 RETURN 2633 ELSE 2634 R=0.5D 0*(RL+RU) 2635 RETURN 2636 END IF 2637 ELSE IF (NTYP.EQ.MTYP) THEN 2638 A = (FU-FL)/(PL*(RU-RL)) 2639 B = PU/PL 2640 END IF 2641 IF (NTYP.EQ.2) THEN 2642* 2643* QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH ONE DIRECTIONAL 2644* DERIVATIVE 2645* 2646 DEN = 2.0D 0*(1.0D 0-A) 2647 ELSE IF (NTYP.EQ.3) THEN 2648* 2649* QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH TWO DIRECTIONAL 2650* DERIVATIVES 2651* 2652 DEN = 1.0D 0 - B 2653 ELSE IF (NTYP.EQ.4) THEN 2654* 2655* CUBIC EXTRAPOLATION OR INTERPOLATION 2656* 2657 C = B - 2.0D 0*A + 1.0D 0 2658 D = B - 3.0D 0*A + 2.0D 0 2659 DIS = D*D - 3.0D 0*C 2660 IF (DIS.LT.0.0D 0) GO TO 1 2661 DEN = D + SQRT(DIS) 2662 ELSE IF (NTYP.EQ.5) THEN 2663* 2664* CONIC EXTRAPOLATION OR INTERPOLATION 2665* 2666 DIS = A*A - B 2667 IF (DIS.LT.0.0D 0) GO TO 1 2668 DEN = A + SQRT(DIS) 2669 IF (DEN.LE.0.0D 0) GO TO 1 2670 DEN = 1.0D 0 - B*(1.0D 0/DEN)**3 2671 END IF 2672 IF (MODE.EQ.1.AND.DEN.GT.0.0D 0.AND.DEN.LT.1.0D 0) THEN 2673* 2674* EXTRAPOLATION ACCEPTED 2675* 2676 R = RL + (RU-RL)/DEN 2677 R = MAX(R,C1L*RU) 2678 R = MIN(R,C1U*RU) 2679 RETURN 2680 ELSE IF (MODE.EQ.2.AND.DEN.GT.1.0D 0) THEN 2681* 2682* INTERPOLATION ACCEPTED 2683* 2684 R = RL + (RU-RL)/DEN 2685 IF (RL.EQ.0.0D 0) THEN 2686 R = MAX(R,RL+C2L*(RU-RL)) 2687 ELSE 2688 R = MAX(R,RL+C3L*(RU-RL)) 2689 END IF 2690 R = MIN(R,RL+C2U*(RU-RL)) 2691 RETURN 2692 END IF 2693 1 CONTINUE 2694 END 2695* SUBROUTINE PNINT3 ALL SYSTEMS 91/12/01 2696* PURPOSE : 2697* EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH WITHOUT DIRECTIONAL 2698* DERIVATIVES. 2699* 2700* PARAMETERS : 2701* RI RO INITIAL VALUE OF THE STEPSIZE PARAMETER. 2702* RI RL LOWER VALUE OF THE STEPSIZE PARAMETER. 2703* RI RU UPPER VALUE OF THE STEPSIZE PARAMETER. 2704* RI RI INNER VALUE OF THE STEPSIZE PARAMETER. 2705* RI FO VALUE OF THE OBJECTIVE FUNCTION FOR R=RO. 2706* RI FL VALUE OF THE OBJECTIVE FUNCTION FOR R=RL. 2707* RI FU VALUE OF THE OBJECTIVE FUNCTION FOR R=RU. 2708* RI FI VALUE OF THE OBJECTIVE FUNCTION FOR R=RI. 2709* RO PO INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE. 2710* RO R VALUE OF THE STEPSIZE PARAMETER OBTAINED. 2711* II MODE MODE OF LINE SEARCH. 2712* II MTYP METHOD SELECTION. MTYP=1-BISECTION. MTYP=2-TWO POINT 2713* QUADRATIC INTERPOLATION. MTYP=2-THREE POINT QUADRATIC 2714* INTERPOLATION. 2715* IO MERR ERROR INDICATOR. MERR=0 FOR NORMAL RETURN. 2716* 2717* METHOD : 2718* EXTRAPOLATION OR INTERPOLATION WITH STANDARD MODEL FUNCTIONS. 2719* 2720 SUBROUTINE PNINT3(RO,RL,RU,RI,FO,FL,FU,FI,PO,R,MODE,MTYP,MERR) 2721 DOUBLE PRECISION RO,RL,RU,RI,FO,FL,FU,FI,PO,R 2722 INTEGER MODE,MTYP,MERR,NTYP 2723 DOUBLE PRECISION AL,AU,AI,DEN,DIS 2724 LOGICAL L1,L2 2725 DOUBLE PRECISION ZERO,HALF,ONE,TWO,THREE,FOUR,C1L,C1U,C2L,C2U,C3L 2726 PARAMETER(ZERO=0.0D 0,HALF=0.5D 0,ONE=1.0D 0,TWO=2.0D 0, 2727 & THREE=3.0D 0,FOUR=4.0D 0,C1L=1.1D 0,C1U=1.0D 3, 2728 & C2L=1.0D-2,C2U=0.9D 0,C3L=1.0D-1) 2729 MERR = 0 2730 IF (MODE .LE. 0) RETURN 2731 IF (PO .GE. ZERO) THEN 2732 MERR = 2 2733 RETURN 2734 ELSE IF (RU .LE. RL) THEN 2735 MERR = 3 2736 RETURN 2737 END IF 2738 L1 = RL .LE. RO 2739 L2 = RI .LE. RL 2740 DO 1 NTYP = MTYP, 1, -1 2741 IF (NTYP .EQ. 1) THEN 2742* 2743* BISECTION 2744* 2745 IF (MODE .EQ. 1) THEN 2746 R = TWO * RU 2747 RETURN 2748 ELSE IF (RI-RL.LE.RU-RI) THEN 2749 R=HALF*(RI+RU) 2750 RETURN 2751 ELSE 2752 R=HALF*(RL+RI) 2753 RETURN 2754 END IF 2755 ELSE IF (NTYP.EQ.MTYP.AND.L1) THEN 2756 IF (.NOT.L2) AI=(FI-FO)/(RI*PO) 2757 AU=(FU-FO)/(RU*PO) 2758 END IF 2759 IF (L1.AND.(NTYP.EQ.2.OR.L2)) THEN 2760* 2761* TWO POINT QUADRATIC EXTRAPOLATION OR INTERPOLATION 2762* 2763 IF (AU.GE.ONE) GO TO 1 2764 R=HALF*RU/(ONE-AU) 2765 ELSE IF (.NOT.L1.OR..NOT.L2.AND.NTYP.EQ.3) THEN 2766* 2767* THREE POINT QUADRATIC EXTRAPOLATION OR INTERPOLATION 2768* 2769 AL=(FI-FL)/(RI-RL) 2770 AU=(FU-FI)/(RU-RI) 2771 DEN=AU-AL 2772 IF (DEN.LE.ZERO) GO TO 1 2773 R=RI-HALF*(AU*(RI-RL)+AL*(RU-RI))/DEN 2774 ELSE IF (L1.AND..NOT.L2.AND.NTYP.EQ.4) THEN 2775* 2776* THREE POINT CUBIC EXTRAPOLATION OR INTERPOLATION 2777* 2778 DIS=(AI-ONE)*(RU/RI) 2779 DEN=(AU-ONE)*(RI/RU)-DIS 2780 DIS=AU+AI-DEN-TWO*(ONE+DIS) 2781 DIS=DEN*DEN-THREE*DIS 2782 IF (DIS.LT.ZERO) GO TO 1 2783 DEN=DEN+SQRT(DIS) 2784 IF (DEN.EQ.ZERO) GO TO 1 2785 R=(RU-RI)/DEN 2786 ELSE 2787 GO TO 1 2788 END IF 2789 IF (MODE .EQ. 1 .AND. R .GT. RU) THEN 2790* 2791* EXTRAPOLATION ACCEPTED 2792* 2793 R = MAX( R, C1L*RU) 2794 R = MIN( R, C1U*RU) 2795 RETURN 2796 ELSE IF (MODE .EQ. 2 .AND. R .GT. RL .AND. R .LT. RU) THEN 2797* 2798* INTERPOLATION ACCEPTED 2799* 2800 IF (RI.EQ.ZERO.AND.NTYP.NE.4) THEN 2801 R = MAX( R, RL + C2L*(RU-RL)) 2802 ELSE 2803 R = MAX( R, RL + C3L*(RU-RL)) 2804 END IF 2805 R = MIN( R, RL + C2U*(RU-RL)) 2806 IF (R.EQ.RI) GO TO 1 2807 RETURN 2808 END IF 2809 1 CONTINUE 2810 END 2811* SUBROUTINE PNNEQ1 ALL SYSTEMS 92/12/01 2812* PURPOSE : 2813* SOLUTION OF A SINGLE NONLINEAR EQUATION. 2814* 2815* PARAMETERS : 2816* RI AA LEFT ENDPOINT OF THE INTERVAL. 2817* RI BB RIGHT ENDPOINT OF THE INTERVAL. 2818* RO X COMPUTED SOLUTION POINT. 2819* RO F COMPUTED VALUE OF THE NONLINEAR FUNCTION. 2820* RF FUN EXTERNAL FUNCTION. 2821* RI EPSX REQUIRED PRECISION FOR THE SOLUTION POINT. 2822* RI EPSF REQUIRED PRECISION FOR THE NONLINEAR FUNCTION. 2823* IO IC NUMBER OF ITERATIONS. 2824* IO IE ERROR SPECIFICATION. 2825* IU ISYS CONTROL PARAMETER. 2826* 2827* METHOD : 2828* D.LEE: THREE NEW RAPIDLY CONVERGENT ALGORITHMS FOR FINDING A ZERO 2829* OF A FUNCTION, SIAM J. SCI. STAT. COMPUT. 6 (1985) 193-208. 2830* 2831 SUBROUTINE PNNEQ1(AA,BB,X,F,EPSX,EPSF,IC,IE,ISYS) 2832 DOUBLE PRECISION AA,BB,X,F,EPSX,EPSF 2833 INTEGER IC,IE,ISYS 2834 INTEGER ITER,ITMAX,K,L 2835 DOUBLE PRECISION FA,FB,X1,X2,X3,F1,F2,F3,R,R1,RA,RB,D,D1,A,B,C,Z, 2836 & W,FW,GW,DEL,DDL,F21,F32 2837 DOUBLE PRECISION ZERO,ONE,TWO,THREE,FOUR,HALF,CON 2838 PARAMETER (ZERO=0.0D 0,ONE=1.0D 0,TWO=2.0D 0,THREE=3.0D 0, 2839 & FOUR=4.0D 0,HALF=0.5D 0,CON=0.1D 0) 2840 SAVE A,B,C,FA,FB,X1,X2,X3,F1,F2,F3,R,D,FW 2841 SAVE L,ITER,ITMAX 2842 GO TO (1,2,3,4,6) ISYS+1 2843 1 IE=0 2844 ITMAX=IC 2845 IF (ITMAX.LE.0) ITMAX=100 2846 X=AA 2847 ISYS=1 2848 IC=1 2849 RETURN 2850 2 CONTINUE 2851 IF (ABS(F).LE.EPSF) GO TO 7 2852 FA=F 2853 X=BB 2854 ISYS=2 2855 IC=2 2856 RETURN 2857 3 CONTINUE 2858 IF (ABS(F).LE.EPSF) GO TO 7 2859 FB=F 2860 IF (FA*FB.GT.0.0D 0) THEN 2861 X=AA 2862 F=FA 2863 IE=-2 2864 GO TO 7 2865 END IF 2866 X1=AA 2867 F1=FA 2868 X=HALF*(AA+BB) 2869 ISYS=3 2870 IC=3 2871 RETURN 2872 4 CONTINUE 2873 X2=X 2874 F2=F 2875 IF (F1*F2.GT.0.0D 0) THEN 2876 X3=X1 2877 F3=F1 2878 X1=BB 2879 F1=FB 2880 ELSE 2881 X3=BB 2882 F3=FB 2883 END IF 2884 L=0 2885 D=0.0D 0 2886 R=0.0D 0 2887 ITER=1 2888 5 CONTINUE 2889 D1=D 2890 R1=R 2891 D=ABS(X1-X2) 2892 IF (ABS(F1).LT.ABS(F2)) THEN 2893 X=X1 2894 F=F1 2895 ELSE 2896 X=X2 2897 F=F2 2898 END IF 2899 DEL=EPSX*(ABS(X)+ONE) 2900 IF (ABS(F).LE.EPSF.OR.D.LE.TWO*DEL) GO TO 7 2901 Z=X1+HALF*(X2-X1) 2902 DDL=MAX(CON*D,DEL) 2903 IF (THREE*D.LE.TWO*D1) THEN 2904 K=0 2905 ELSE 2906 K=1 2907 END IF 2908 IF (X2.EQ.X1) THEN 2909 F21=0.0D 0 2910 ELSE 2911 F21=(F2-F1)/(X2-X1) 2912 ENDIF 2913 IF (X3.EQ.X2) THEN 2914 F32=0.0D 0 2915 ELSE 2916 F32=(F3-F2)/(X3-X2) 2917 ENDIF 2918 A=(F32-F21)/(X3-X1) 2919 B=A*(X2+X1)-F21 2920 C=F2-(A*X2-B)*X2 2921 IF (ABS(A).LE.1.0D-10) THEN 2922 R=(F2*X1-F1*X2)/(F2-F1) 2923 ELSE 2924 R=B*B-FOUR*A*C 2925 IF (R.LT.0.0D 0) THEN 2926 R=(F2*X1-F1*X2)/(F2-F1) 2927 ELSE 2928 R=SQRT(R) 2929 RA=HALF*(B+R)/A 2930 RB=HALF*(B-R)/A 2931 IF (ABS(RA-Z).LE.ABS(RB-Z)) THEN 2932 R=RA 2933 ELSE 2934 R=RB 2935 END IF 2936 IF (R.LE.MIN(X1,X2).OR.R.GE.MAX(X1,X2)) THEN 2937 R=(F2*X1-F1*X2)/(F2-F1) 2938 END IF 2939 END IF 2940 END IF 2941 IF (L.GE.2) THEN 2942 W=R 2943 IF (ABS(W-X).LT.DEL) W=X+DEL*SIGN(ONE,Z-X) 2944 ELSE IF (K.EQ.1.OR.ABS(R-X).GE.ABS(Z-X)) THEN 2945 W=Z 2946 ELSE 2947 W=R+HALF*ABS(R-R1)*SIGN(ONE,R-X) 2948 IF (ABS(W-X).LT.DDL) W=X+DDL*SIGN(ONE,Z-X) 2949 IF (ABS(W-X).GE.ABS(Z-X)) W=Z 2950 END IF 2951 X=W 2952 FW=F 2953 ISYS=4 2954 IC=IC+1 2955 RETURN 2956 6 CONTINUE 2957 GW=(A*X-B)*X+C 2958 IF (ABS(F-GW).LE.1.0D-1*ABS(FW).OR.ABS(FW).LE.1.0D-3* 2959 *MAX(ABS(F1),ABS(F2)).AND.L.GE.2) THEN 2960 L=L+1 2961 ELSE 2962 L=0 2963 END IF 2964 IF (F*SIGN(ONE,F1).GE.0.0D 0) THEN 2965 IF (D.LE.ABS(X3-X)) THEN 2966 X3=X1 2967 F3=F1 2968 X1=X2 2969 F1=F2 2970 X2=X 2971 F2=F 2972 ELSE 2973 X1=X 2974 F1=F 2975 END IF 2976 ELSE 2977 X3=X2 2978 F3=F2 2979 X2=X 2980 F2=F 2981 END IF 2982 ITER=ITER+1 2983 IF (ITER.LE.ITMAX) GO TO 5 2984 IE=-1 2985 7 ISYS=0 2986 RETURN 2987 END 2988* SUBROUTINE PNSTEP ALL SYSTEMS 89/12/01 2989* PURPOSE : 2990* DETERMINATION OF A SCALING FACTOR FOR THE BOUNDARY STEP. 2991* 2992* PARAMETERS : 2993* RI DEL MAXIMUM STEPSIZE. 2994* RI A INPUT PARAMETER. 2995* RI B INPUT PARAMETER. 2996* RI C INPUT PARAMETER. 2997* RO ALF SCALING FACTOR FOR THE BOUNDARY STEP SUCH THAT 2998* A**2+2*B*ALF+C*ALF**2=DEL**2. 2999* 3000 SUBROUTINE PNSTEP(DEL,A,B,C,ALF) 3001 DOUBLE PRECISION DEL, A, B, C, ALF 3002 DOUBLE PRECISION DEN, DIS 3003 ALF = 0.0D 0 3004 DEN = (DEL+A) * (DEL-A) 3005 IF (DEN .LE. 0.0D 0) RETURN 3006 DIS = B*B + C*DEN 3007 IF (B .GE. 0.0D 0) THEN 3008 ALF = DEN / (SQRT(DIS) + B) 3009 ELSE 3010 ALF = (SQRT(DIS) - B) / C 3011 END IF 3012 RETURN 3013 END 3014* SUBROUTINE PNSTP4 ALL SYSTEMS 99/12/01 3015* PURPOSE : 3016* STEPSIZE SELECTION USING POLYHEDRAL APPROXIMATION 3017* FOR DESCENT STEP IN NONCONVEX VARIABLE METRIC METHOD. 3018* 3019* PARAMETERS : 3020* II N ACTUAL NUMBER OF VARIABLES. 3021* II MA DECLARED NUMBER OF LINEAR APPROXIMATED FUNCTIONS 3022* II MAL CURRENT NUMBER OF LINEAR APPROXIMATED FUNCTIONS. 3023* RU X(N) VECTOR OF VARIABLES. 3024* RI AF(4*MA) VECTOR OF BUNDLE FUNCTIONS VALUES. 3025* RI AG(N*MA) MATRIX WHOSE COLUMNS ARE BUNDLE SUBGRADIENTS. 3026* RI AY(N*MA) MATRIX WHOSE COLUMNS ARE VARIABLE VECTORS. 3027* RI S(N) DIRECTION VECTOR. 3028* RI F VALUE OF THE OBJECTIVE FUNCTION. 3029* RI DF DIRECTIONAL DERIVATIVE. 3030* RO T VALUE OF THE STEPSIZE PARAMETER. 3031* RO TB BUNDLE PARAMETER FOR MATRIX SCALING. 3032* RI ETA5 DISTANCE MEASURE PARAMETER. 3033* RI ETA9 MAXIMUM FOR REAL NUMBERS. 3034* RI MOS3 LOCALITY MEASURE PARAMETER. 3035* 3036 SUBROUTINE PNSTP4(N,MA,MAL,X,AF,AG,AY,S,F,DF,T,TB,ETA5,ETA9,MOS3) 3037 DOUBLE PRECISION DF,ETA5,ETA9,F,T,TB 3038 INTEGER MA,MAL,MOS3,N 3039 DOUBLE PRECISION AF(*),AG(*),AY(*),S(*),X(*) 3040 DOUBLE PRECISION ALF,ALFL,ALFR,BET,BETL,BETR,DX,Q,R,W 3041 INTEGER I,J,JN,K,L,LQ 3042 W = DF*T* (1.0D0-T*0.5D0) 3043* 3044* INITIAL CHOICE OF POSSIBLY ACTIVE LINES 3045* 3046 K = 0 3047 L = -1 3048 JN = 0 3049 TB = SQRT(ETA9) 3050 BETR = -ETA9 3051 DO 20 J = 1,MAL - 1 3052 R = 0.0D0 3053 BET = 0.0D0 3054 ALFL = AF(J) - F 3055 DO 10 I = 1,N 3056 DX = X(I) - AY(JN+I) 3057 Q = AG(JN+I) 3058 R = R + DX*DX 3059 ALFL = ALFL + DX*Q 3060 BET = BET + S(I)*Q 3061 10 CONTINUE 3062 IF (MOS3.NE.2) R = R** (DBLE(MOS3)*0.5D0) 3063 ALF = MAX(ABS(ALFL),ETA5*R) 3064 R = 1.0D0 - BET/DF 3065 IF (R*R+ (ALF+ALF)/DF.GT.1.0D-6) THEN 3066 K = K + 1 3067 AF(MA+K) = ALF 3068 AF(MA+MA+K) = BET 3069 R = T*BET - ALF 3070 IF (R.GT.W) THEN 3071 W = R 3072 L = K 3073 END IF 3074 END IF 3075 IF (BET.GT.0.0D0) TB = MIN(TB,ALF/ (BET-DF)) 3076 BETR = MAX(BETR,BET-ALF) 3077 JN = JN + N 3078 20 CONTINUE 3079 LQ = -1 3080 IF (BETR.LE.DF*0.5D0) RETURN 3081 LQ = 1 3082 IF (L.LT.0) RETURN 3083 BETR = AF(MA+MA+L) 3084 IF (BETR.LE.0.0D0) THEN 3085 IF (T.LT.1.0D0 .OR. BETR.EQ.0.0D0) RETURN 3086 LQ = 2 3087 END IF 3088 ALFR = AF(MA+L) 3089* 3090* ITERATION LOOP 3091* 3092 30 IF (LQ.GE.1) THEN 3093 Q = 1.0D0 - BETR/DF 3094 R = Q + SQRT(Q*Q+ (ALFR+ALFR)/DF) 3095 IF (BETR.GE.0.0D0) R = - (ALFR+ALFR)/ (DF*R) 3096 R = MIN(1.95D0,MAX(0.0D0,R)) 3097 ELSE 3098 IF (ABS(BETR-BETL)+ABS(ALFR-ALFL).LT.-1.0D-4*DF) RETURN 3099 R = (ALFR-ALFL)/ (BETR-BETL) 3100 END IF 3101 IF (ABS(T-R).LT.1.0D-4) RETURN 3102 T = R 3103 AF(MA+L) = -1.0D0 3104 W = T*BETR - ALFR 3105 L = -1 3106 DO 40 J = 1,K 3107 ALF = AF(MA+J) 3108 IF (ALF.LT.0.0D0) GO TO 40 3109 BET = AF(MA+MA+J) 3110 R = T*BET - ALF 3111 IF (R.GT.W) THEN 3112 W = R 3113 L = J 3114 END IF 3115 40 CONTINUE 3116 IF (L.LT.0) RETURN 3117 BET = AF(MA+MA+L) 3118 IF (BET.EQ.0.0D0) RETURN 3119* 3120* NEW INTERVAL SELECTION 3121* 3122 ALF = AF(MA+L) 3123 IF (BET.LT.0.0D0) THEN 3124 IF (LQ.EQ.2) THEN 3125 ALFR = ALF 3126 BETR = BET 3127 ELSE 3128 ALFL = ALF 3129 BETL = BET 3130 LQ = 0 3131 END IF 3132 ELSE 3133 IF (LQ.EQ.2) THEN 3134 ALFL = ALFR 3135 BETL = BETR 3136 LQ = 0 3137 END IF 3138 ALFR = ALF 3139 BETR = BET 3140 END IF 3141 GO TO 30 3142 END 3143* SUBROUTINE PNSTP5 ALL SYSTEMS 99/12/01 3144* PURPOSE : 3145* STEPSIZE SELECTION USING POLYHEDRAL APPROXIMATION 3146* FOR NULL STEP IN NONCONVEX VARIABLE METRIC METHOD. 3147* 3148* PARAMETERS : 3149* II N ACTUAL NUMBER OF VARIABLES. 3150* II MA DECLARED NUMBER OF LINEAR APPROXIMATED FUNCTIONS. 3151* II MAL CURRENT NUMBER OF LINEAR APPROXIMATED FUNCTIONS. 3152* RU X(N) VECTOR OF VARIABLES. 3153* RI AF(4*MA) VECTOR OF BUNDLE FUNCTIONS VALUES. 3154* RI AG(N*MA) MATRIX WHOSE COLUMNS ARE BUNDLE SUBGRADIENTS. 3155* RI AY(N*MA) MATRIX WHOSE COLUMNS ARE VARIABLE VECTORS. 3156* RI S(N) DIRECTION VECTOR. 3157* RI F VALUE OF THE OBJECTIVE FUNCTION. 3158* RI DF DIRECTIONAL DERIVATIVE. 3159* RO T VALUE OF THE STEPSIZE PARAMETER. 3160* RO TB BUNDLE PARAMETER FOR MATRIX SCALING. 3161* RI ETA5 DISTANCE MEASURE PARAMETER. 3162* RI ETA9 MAXIMUM FOR REAL NUMBERS. 3163* RI MOS3 LOCALITY MEASURE PARAMETER. 3164* 3165 SUBROUTINE PNSTP5(N,MA,MAL,X,AF,AG,AY,S,F,DF,T,TB,ETA5,ETA9,MOS3) 3166 DOUBLE PRECISION DF,ETA5,ETA9,F,T,TB 3167 INTEGER MA,MAL,MOS3,N 3168 DOUBLE PRECISION AF(*),AG(*),AY(*),S(*),X(*) 3169 DOUBLE PRECISION ALF,ALFL,ALFR,BET,BETL,BETR,DX,Q,R,W 3170 INTEGER I,J,JN,K,L 3171 W = DF*T 3172* 3173* INITIAL CHOICE OF POSSIBLY ACTIVE PARABOLAS 3174* 3175 K = 0 3176 L = -1 3177 JN = 0 3178 TB = SQRT(ETA9) 3179 BETR = -ETA9 3180 DO 20 J = 1,MAL - 1 3181 BET = 0.0D0 3182 R = 0.0D0 3183 ALFL = AF(J) - F 3184 DO 10 I = 1,N 3185 DX = X(I) - AY(JN+I) 3186 R = R + DX*DX 3187 Q = AG(JN+I) 3188 ALFL = ALFL + DX*Q 3189 BET = BET + S(I)*Q 3190 10 CONTINUE 3191 IF (MOS3.NE.2) R = R** (DBLE(MOS3)*0.5D0) 3192 ALF = MAX(ABS(ALFL),ETA5*R) 3193 IF (BET+BET.GT.DF) TB = MIN(TB,ALF/ (BET-DF)) 3194 BETR = MAX(BETR,BET-ALF) 3195 IF (ALF.LT.BET-DF) THEN 3196 K = K + 1 3197 R = T*BET - ALF 3198 AF(MA+K) = ALF 3199 AF(MA+MA+K) = BET 3200 IF (R.GT.W) THEN 3201 W = R 3202 L = K 3203 END IF 3204 END IF 3205 JN = JN + N 3206 20 CONTINUE 3207 IF (L.LT.0) RETURN 3208 BETR = AF(MA+MA+L) 3209 ALFR = AF(MA+L) 3210 ALF = ALFR 3211 BET = BETR 3212 ALFL = 0.0D0 3213 BETL = DF 3214* 3215* ITERATION LOOP 3216* 3217 30 W = BET/DF 3218 IF (ABS(BETR-BETL)+ABS(ALFR-ALFL).LT.-1.0D-4*DF) RETURN 3219 IF (BETR-BETL.EQ.0.0D0) STOP 11 3220 R = (ALFR-ALFL)/ (BETR-BETL) 3221 IF (ABS(T-W).LT.ABS(T-R)) R = W 3222 Q = T 3223 T = R 3224 IF (ABS(T-Q).LT.1.0D-3) RETURN 3225 AF(MA+L) = -1.0D0 3226 W = T*BET - ALF 3227 L = -1 3228 DO 40 J = 1,K 3229 ALF = AF(MA+J) 3230 IF (ALF.LT.0.0D0) GO TO 40 3231 BET = AF(MA+MA+J) 3232 R = T*BET - ALF 3233 IF (R.GT.W) THEN 3234 W = R 3235 L = J 3236 END IF 3237 40 CONTINUE 3238 IF (L.LT.0) RETURN 3239 BET = AF(MA+MA+L) 3240 Q = BET - T*DF 3241 IF (Q.EQ.0.0D0) RETURN 3242* 3243* NEW INTERVAL SELECTION 3244* 3245 ALF = AF(MA+L) 3246 IF (Q.LT.0.0D0) THEN 3247 ALFL = ALF 3248 BETL = BET 3249 ELSE 3250 ALFR = ALF 3251 BETR = BET 3252 END IF 3253 GO TO 30 3254 END 3255* SUBROUTINE PP0BA1 ALL SYSTEMS 05/12/01 3256* PURPOSE : 3257* EVALUATION OF THE BARRIER FUNCTION FOR THE SUM OF ABSOLUTE VALUES. 3258* 3259* PARAMETERS : 3260* II NA NUMBER OF APPROXIMATED FUNCTIONS. 3261* RI AS(NA) SUM OF ABSOLUTE VALUE SLACK VARIABLES. 3262* RI RPF3 BARRIER COEFFICIENT. 3263* RO F VALUE OF THE BARRIER FUNCTION. 3264* 3265 SUBROUTINE PP0BA1(NA,AS,RPF3,F) 3266 INTEGER NA 3267 DOUBLE PRECISION AS(*),RPF3,F 3268 INTEGER KA 3269 F=-DBLE(NA)*RPF3*LOG(2.0D 0*RPF3) 3270 DO 1 KA=1,NA 3271 F=F+AS(KA)-RPF3*LOG(AS(KA)) 3272 1 CONTINUE 3273 RETURN 3274 END 3275* SUBROUTINE PP0BX1 ALL SYSTEMS 05/12/01 3276* PURPOSE : 3277* EVALUATION OF THE BARRIER FUNCTION FOR THE MINIMAX OPTIMIZATION. 3278* 3279* PARAMETERS : 3280* II NA NUMBER OF APPROXIMATED FUNCTIONS. 3281* RI Z MINIMAX SLACK VARIABLE. 3282* RI AF(NA) VECTOR CONTAINING VALUES OF APPROXIMATED FUNCTIONS. 3283* RO F VALUE OF THE BARRIERY FUNCTION. 3284* RI FF VALUE OF THE THE OBJECTIVE FUNCTION. 3285* RI PAR PARAMETER OF THE BEN-TAL BARRIER FUNCTION. 3286* RI RPF3 BARRIER COEFFICIENT. 3287* II MEP MERIT FUNCTION USED. MEP=1-LOGARITHMIC BARIER FUNCTION. 3288* MEP=2-BEN-TAL BARRIER FUNCTION. MEP=3-COMPOSITE BARRIER 3289* FUNCTION. 3290* II IEXT KIND OF THE MINIMAX APPROXIMATION. IEXT=0-CHEBYSHEV 3291* APPROXIMATION. IEXT=-1-MINIMAX. IEXT=+1-MAXIMIN. 3292* 3293 SUBROUTINE PP0BX1(NA,Z,AF,F,FF,PAR,RPF3,MEP,IEXT) 3294 INTEGER NA,MEP,IEXT 3295 DOUBLE PRECISION Z,AF(*),PAR,RPF3,F,FF 3296 DOUBLE PRECISION FA 3297 INTEGER KA 3298 IF (Z.LE.FF) THEN 3299 F=1.0D 60 3300 ELSE 3301 F=Z 3302 IF (MEP.EQ.1) THEN 3303 DO 11 KA=1,NA 3304 FA=AF(KA) 3305 IF (IEXT.LE.0) THEN 3306 F=F-RPF3*LOG(Z-FA) 3307 END IF 3308 IF (IEXT.GE.0) THEN 3309 F=F-RPF3*LOG(Z+FA) 3310 END IF 3311 11 CONTINUE 3312 ELSE IF (MEP.EQ.2) THEN 3313 DO 21 KA=1,NA 3314 FA=AF(KA) 3315 IF (IEXT.LE.0) THEN 3316 IF (Z-FA.LE.PAR) THEN 3317 F=F-RPF3*LOG(Z-FA) 3318 ELSE 3319 F=F+(2.0D 0-0.5D 0*PAR/(Z-FA))*RPF3*PAR/(Z-FA) 3320 END IF 3321 END IF 3322 IF (IEXT.GE.0) THEN 3323 IF (Z+FA.LE.PAR) THEN 3324 F=F-RPF3*LOG(Z+FA) 3325 ELSE 3326 F=F+(2.0D 0-0.5D 0*PAR/(Z+FA))*RPF3*PAR/(Z+FA) 3327 END IF 3328 END IF 3329 21 CONTINUE 3330 ELSE IF (MEP.EQ.3) THEN 3331 DO 31 KA=1,NA 3332 FA=AF(KA) 3333 IF (IEXT.LE.0) THEN 3334 F=F+RPF3*LOG(1.0D 0/(Z-FA)+1.0D 0) 3335 END IF 3336 IF (IEXT.GE.0) THEN 3337 F=F+RPF3*LOG(1.0D 0/(Z+FA)+1.0D 0) 3338 END IF 3339 31 CONTINUE 3340 ELSE IF (MEP.EQ.4) THEN 3341 DO 41 KA=1,NA 3342 FA=AF(KA) 3343 IF (IEXT.LE.0) THEN 3344 F=F+RPF3*RPF3/(Z-FA) 3345 END IF 3346 IF (IEXT.GE.0) THEN 3347 F=F+RPF3*RPF3/(Z+FA) 3348 END IF 3349 41 CONTINUE 3350 END IF 3351 END IF 3352 RETURN 3353 END 3354* SUBROUTINE PP1MX3 ALL SYSTEMS 05/12/01 3355* PURPOSE : 3356* COMPUTATION OF THE VALUE AND THE GRADIENT OF THE LAGRANGIAN FUNCTION 3357* FOR THE MINIMAX OPTIMIZATION. 3358* 3359* PARAMETERS: 3360* II NF NUMBER OF VARIABLES. 3361* II NA NUMBER OF APPROXIMATED FUNCTIONS. 3362* RI X(NF) VECTOR OF VARIABLES. 3363* RI GA(NF) GRADIENT OF THE APPROXIMATED FUNCTION. 3364* RI AG(IAG(N+1)-1) SPARSE RECTANGULAR MATRIX WHICH IS USED FOR THE 3365* DIRECTION VECTOR DETERMINATION. 3366* II IAG(N+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 3367* II JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 3368* RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 3369* RI AZL(NA) LOWER LAGRANGE MULTIPLIERS. 3370* RI AZU(NA) UPPER LAGRANGE MULTIPLIERS. 3371* RI FA VALUE OF THE SELECTED FUNCTION. 3372* RI AF(NA) VALUES OF THE APPROXIMATED FUNCTIONS. 3373* RO F VALUE OF THE OBJECTIVE FUNCTION. 3374* II KD DEGREE OF REQUIRED DERIVATIVES. 3375* IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. 3376* IU NFV NUMBER OF OBJECTIVE FUNCTION VALUES COMPUTED. 3377* IU NFG NUMBER OF OBJECTIVE FUNCTION GRADIENTS COMPUTED. 3378* II ISNA INDICATOR FOR STORING ELEMENTAL FUNCTION VALUES AND 3379* GRADIENTS. ISNA=0-STORING SUPPRESSED. ISNA=1-STORING 3380* ELEMENTAL FUNCTION VALUES. ISNA=2-STORING ELEMENTAL 3381* FUNCTION VALUES AND GRADIENTS. 3382* II IEXT TYPE OF MINIMAX. IEXT=0-MINIMIZATION OF THE MAXIMUM VALUE. 3383* IEXT=1-MINIMIZATION OF THE MAXIMUM ABSOLUTE VALUE. 3384* 3385* SUBPROGRAMS USED : 3386* SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION. 3387* SE DFUN COMPUTATION OF THE GRADIENT OF THE APPROXIMATED FUNCTION. 3388* S MXVSET INITIATION OF A VECTOR. 3389* 3390 SUBROUTINE PP1MX3(NF,NA,X,GA,AG,IAG,JAG,G,AZL,AZU,FA,AF, 3391 & F,KD,LD,NFV,NFG,ISNA,IEXT) 3392 INTEGER NF,NA,IAG(*),JAG(*),KD,LD,NFV,NFG,ISNA,IEXT 3393 DOUBLE PRECISION X(*),GA(*),AG(*),G(*),AZL(*),AZU(*),FA,AF(*),F 3394 INTEGER J,JP,K,KA,L 3395 IF (KD.LE.LD) RETURN 3396 IF (KD.GE.0.AND.LD.LT.0) THEN 3397 NFV=NFV+1 3398 END IF 3399 IF (KD.GE.1.AND.LD.LT.1) THEN 3400 CALL MXVSET(NF,0.0D 0,G) 3401 NFG=NFG+1 3402 END IF 3403 DO 3 KA=1,NA 3404 IF (LD.GE.0) GO TO 1 3405 CALL FUN(NF,KA,X,FA) 3406 IF (ISNA.GE.1) AF(KA)=FA 3407 IF (IEXT.EQ.0) THEN 3408 IF (KA.EQ.1) F=ABS(FA) 3409 F=MAX(F,ABS(FA)) 3410 ELSE IF (IEXT.LT.0) THEN 3411 IF (KA.EQ.1) F= FA 3412 F=MAX(F, FA) 3413 ELSE IF (IEXT.GT.0) THEN 3414 IF (KA.EQ.1) F=-FA 3415 F=MAX(F,-FA) 3416 END IF 3417 1 IF (KD.LT.1) GO TO 3 3418 IF (LD.GE.1) GO TO 3 3419 CALL DFUN(NF,KA,X,GA) 3420 K=IAG(KA) 3421 L=IAG(KA+1)-K 3422 DO 2 J=1,L 3423 JP=ABS(JAG(K)) 3424 IF (IEXT.EQ.0) THEN 3425 G(JP)=G(JP)+(AZU(KA)-AZL(KA))*GA(JP) 3426 ELSE IF (IEXT.LT.0) THEN 3427 G(JP)=G(JP)+AZU(KA)*GA(JP) 3428 ELSE IF (IEXT.GT.0) THEN 3429 G(JP)=G(JP)-AZL(KA)*GA(JP) 3430 END IF 3431 IF (ISNA.GE.2) AG(K)=GA(JP) 3432 K=K+1 3433 2 CONTINUE 3434 3 CONTINUE 3435 RETURN 3436 END 3437* SUBROUTINE PP1SA3 ALL SYSTEMS 05/12/01 3438* PURPOSE : 3439* COMPUTATION OF THE VALUE AND THE GRADIENT OF THE LAGRANGIAN FUNCTION 3440* FOR THE SUM OF ABSOLUTE VALUES. 3441* 3442* PARAMETERS: 3443* II NF NUMBER OF VARIABLES. 3444* II NA NUMBER OF APPROXIMATED FUNCTIONS. 3445* RI X(NF) VECTOR OF VARIABLES. 3446* RI GA(NF) GRADIENT OF THE APPROXIMATED FUNCTION. 3447* RI AG(IAG(N+1)-1) SPARSE RECTANGULAR MATRIX WHICH IS USED FOR THE 3448* DIRECTION VECTOR DETERMINATION. 3449* II IAG(N+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 3450* II JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 3451* RO G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 3452* RI AZ(NA) VECTOR OF LAGRANGE MULTIPLIERS. 3453* RI FA VALUE OF THE SELECTED FUNCTION. 3454* RI AF(NA) VALUES OF THE APPROXIMATED FUNCTIONS. 3455* RO F VALUE OF THE OBJECTIVE FUNCTION. 3456* II KD DEGREE OF REQUIRED DERIVATIVES. 3457* IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES. 3458* IU NFV NUMBER OF OBJECTIVE FUNCTION VALUES COMPUTED. 3459* IU NFG NUMBER OF OBJECTIVE FUNCTION GRADIENTS COMPUTED. 3460* II ISNA INDICATOR FOR STORING ELEMENTAL FUNCTION VALUES AND 3461* GRADIENTS. ISNA=0-STORING SUPPRESSED. ISNA=1-STORING 3462* ELEMENTAL FUNCTION VALUES. ISNA=2-STORING ELEMENTAL 3463* FUNCTION VALUES AND GRADIENTS. 3464* 3465* SUBPROGRAMS USED : 3466* SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION. 3467* SE DFUN COMPUTATION OF THE GRADIENT OF THE APPROXIMATED FUNCTION. 3468* S MXVSET INITIATION OF A VECTOR. 3469* 3470 SUBROUTINE PP1SA3(NF,NA,X,GA,AG,IAG,JAG,G,AZ,FA,AF,F,KD,LD,NFV, 3471 & NFG,ISNA) 3472 INTEGER NF,NA,IAG(*),JAG(*),KD,LD,NFV,NFG,ISNA 3473 DOUBLE PRECISION X(*),GA(*),AG(*),G(*),AZ(*),FA,AF(*),F 3474 INTEGER J,JP,K,KA,L 3475 IF (KD.LE.LD) RETURN 3476 IF (KD.GE.0.AND.LD.LT.0) THEN 3477 F=0.0D 0 3478 NFV=NFV+1 3479 END IF 3480 IF (KD.GE.1.AND.LD.LT.1) THEN 3481 CALL MXVSET(NF,0.0D 0,G) 3482 NFG=NFG+1 3483 END IF 3484 DO 3 KA=1,NA 3485 IF (LD.GE.0) GO TO 1 3486 CALL FUN(NF,KA,X,FA) 3487 IF (ISNA.GE.1) AF(KA)=FA 3488 F=F+ABS(FA) 3489 1 IF (KD.LT.1) GO TO 3 3490 IF (LD.GE.1) GO TO 3 3491 CALL DFUN(NF,KA,X,GA) 3492 K=IAG(KA) 3493 L=IAG(KA+1)-K 3494 DO 2 J=1,L 3495 JP=ABS(JAG(K)) 3496 G(JP)=G(JP)+AZ(KA)*GA(JP) 3497 IF (ISNA.GE.2) AG(K)=GA(JP) 3498 K=K+1 3499 2 CONTINUE 3500 3 CONTINUE 3501 RETURN 3502 END 3503* SUBROUTINE PPLAG1 ALL SYSTEMS 05/12/01 3504* PURPOSE : 3505* COMPUTATION OF THE LAGRANGE MULTIPLIERS FOR THE SUM OF ABSOLUTE 3506* VALUES. 3507* 3508* PARAMETERS : 3509* II NA NUMBER OF APPROXIMATED FUNCTIONS. 3510* RI AF(NA) VECTOR CONTAINING VALUES OF APPROXIMATED FUNCTIONS. 3511* RA AS(NA) AUXILIARY ARRAY. 3512* RO AZ(NA) LAGRANGE MULTIPLIERS. 3513* RI RPF3 BARRIER COEFFICIENT. 3514* 3515 SUBROUTINE PPLAG1(NA,AF,AS,AZ,RPF3) 3516 INTEGER NA 3517 DOUBLE PRECISION AF(*),AS(*),AZ(*),RPF3 3518 DOUBLE PRECISION FA 3519 INTEGER KA 3520 DO 1 KA=1,NA 3521 FA=AF(KA) 3522 AS(KA)=RPF3+SQRT(RPF3**2+FA**2) 3523 AZ(KA)=FA/AS(KA) 3524 1 CONTINUE 3525 RETURN 3526 END 3527* SUBROUTINE PS0G01 ALL SYSTEMS 97/12/01 3528* PURPOSE : 3529* SIMPLE SEARCH WITH TRUST REGION UPDATE. 3530* 3531* PARAMETERS : 3532* RO R VALUE OF THE STEPSIZE PARAMETER. 3533* RO F VALUE OF THE OBJECTIVE FUNCTION. 3534* RI FO INITIAL VALUE OF THE OBJECTIVE FUNCTION. 3535* RI PO INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE. 3536* RI PP QUADRATIC PART OF THE PREDICTED FUNCTION VALUE. 3537* RU XDEL TRUST REGION BOUND. 3538* RO XDELO PREVIOUS TRUST REGION BOUND. 3539* RI XMAX MAXIMUM STEPSIZE. 3540* RI RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER. 3541* RI SNORM EUCLIDEAN NORM OF THE DIRECTION VECTOR. 3542* RI BET1 LOWER BOUND FOR STEPSIZE REDUCTION. 3543* RI BET2 UPPER BOUND FOR STEPSIZE REDUCTION. 3544* RI GAM1 LOWER BOUND FOR STEPSIZE EXPANSION. 3545* RI GAM2 UPPER BOUND FOR STEPSIZE EXPANSION. 3546* RI EPS4 FIRST TOLERANCE FOR RATIO DF/DFPRED. STEP BOUND IS 3547* DECREASED IF DF/DFPRED<EPS4. 3548* RI EPS5 SECOND TOLERANCE FOR RATIO DF/DFPRED. STEP BOUND IS 3549* INCREASED IF IT IS ACTIVE AND DF/DFPRED>EPS5. 3550* II KD DEGREE OF REQUIRED DERIVATIVES. 3551* IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES OF OBJECTIVE 3552* FUNCTION. 3553* IU IDIR INDICATOR FOR DIRECTION DETERMINATION. 3554* IDIR=0-BASIC DETERMINATION. IDIR=1-DETERMINATION 3555* AFTER STEPSIZE REDUCTION. IDIR=2-DETERMINATION AFTER 3556* STEPSIZE EXPANSION. 3557* IO ITERS TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-STEP 3558* BOUND WAS DECREASED. ITERS=2-STEP BOUND WAS UNCHANGED. 3559* ITERS=3-STEP BOUND WAS INCREASED. ITERS=6-FIRST STEPSIZE. 3560* II ITERD CAUSE OF TERMINATION IN THE DIRECTION DETERMINATION. 3561* ITERD<0-BAD DECOMPOSITION. ITERD=0-DESCENT DIRECTION. 3562* ITERD=1-NEWTON LIKE STEP. ITERD=2-INEXACT NEWTON LIKE STEP. 3563* ITERD=3-BOUNDARY STEP. ITERD=4-DIRECTION WITH THE NEGATIVE 3564* CURVATURE. ITERD=5-MARQUARDT STEP. 3565* IO MAXST MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM 3566* STEPSIZE WAS NOT OR WAS REACHED. 3567* IO NRED ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. 3568* II MRED MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. 3569* II KTERS TERMINATION SELECTION. KTERS=1-NORMAL TERMINATION. 3570* KTERS=6-FIRST STEPSIZE. 3571* II MES1 SWITCH FOR EXTRAPOLATION. MES1=1-CONSTANT INCREASING OF 3572* THE INTERVAL. MES1=2-EXTRAPOLATION SPECIFIED BY THE PARAMETER 3573* MES. MES1=3 SUPPRESSED EXTRAPOLATION. 3574* II MES2 SWITCH FOR TERMINATION. MES2=1-NORMAL TERMINATION. 3575* MES2=2-TERMINATION AFTER AT LEAST TWO STEPS (ASYMPTOTICALLY 3576* PERFECT LINE SEARCH). 3577* II MES3 SAFEGUARD AGAINST ROUNDING ERRORS. MES3=0-SAFEGUARD 3578* SUPPRESSED. MES3=1-FIRST LEVEL OF SAFEGUARD. MES3=2-SECOND 3579* LEVEL OF SAFEGUARD. 3580* IU ISYS CONTROL PARAMETER. 3581* 3582* METHOD : 3583* G.A.SCHULTZ, R.B.SCHNABEL, R.H.BYRD: A FAMILY OF TRUST-REGION-BASED 3584* ALGORITHMS FOR UNCONSTRAINED MINIMIZATION WITH STRONG GLOBAL 3585* CONVERGENCE PROPERTIES, SIAM J. NUMER.ANAL. 22 (1985) PP. 47-67. 3586* 3587 SUBROUTINE PS0G01(R,F,FO,PO,PP,XDEL,XDELO,XMAX,RMAX,SNORM,BET1, 3588 & BET2,GAM1,GAM2,EPS4,EPS5,KD,LD,IDIR,ITERS,ITERD,MAXST,NRED,MRED, 3589 & KTERS,MES1,MES2,MES3,ISYS) 3590 INTEGER KD,LD,IDIR,ITERS,ITERD,MAXST,NRED,MRED,KTERS,MES1,MES2, 3591 & MES3,ISYS 3592 DOUBLE PRECISION R,F,FO,PO,PP,XDEL,XDELO,XMAX,RMAX,SNORM,BET1, 3593 & BET2,GAM1,GAM2,EPS4,EPS5 3594 DOUBLE PRECISION DF,DFPRED 3595 INTEGER NRED1,NRED2 3596 SAVE NRED1,NRED2 3597 IF (ISYS.EQ.1) GO TO 2 3598 IF (IDIR.EQ.0) THEN 3599 NRED1=0 3600 NRED2=0 3601 END IF 3602 IDIR=0 3603 XDELO=XDEL 3604* 3605* COMPUTATION OF THE NEW FUNCTION VALUE 3606* 3607 R=MIN(1.0D 0,RMAX) 3608 KD= 0 3609 LD=-1 3610 ISYS=1 3611 RETURN 3612 2 CONTINUE 3613 IF (KTERS.LT.0.OR.KTERS.GT.5) THEN 3614 ITERS=6 3615 ELSE 3616 DF=FO-F 3617 DFPRED=-R*(PO+R*PP) 3618 IF (DF.LT.EPS4*DFPRED) THEN 3619* 3620* STEP IS TOO LARGE, IT HAS TO BE REDUCED 3621* 3622 IF (MES1.EQ.1) THEN 3623 XDEL=BET2*SNORM 3624 ELSE IF (MES1.EQ.2) THEN 3625 XDEL=BET2*MIN(0.5D 0*XDEL,SNORM) 3626 ELSE 3627 XDEL=0.5D 0*PO*SNORM/(PO+DF) 3628 XDEL=MAX(XDEL,BET1*SNORM) 3629 XDEL=MIN(XDEL,BET2*SNORM) 3630 END IF 3631 ITERS=1 3632 IF (MES3.LE.1) THEN 3633 NRED2=NRED2+1 3634 ELSE 3635 IF (ITERD.GT.2) NRED2=NRED2+1 3636 END IF 3637 ELSE IF (DF.LE.EPS5*DFPRED) THEN 3638* 3639* STEP IS SUITABLE 3640* 3641 ITERS=2 3642 ELSE 3643* 3644* STEP IS TOO SMALL, IT HAS TO BE ENLARGED 3645* 3646 IF (MES2.EQ.2) THEN 3647 XDEL=MAX(XDEL,GAM1*SNORM) 3648 ELSE IF (ITERD.GT.2) THEN 3649 XDEL=GAM1*XDEL 3650 END IF 3651 ITERS=3 3652 END IF 3653 XDEL=MIN(XDEL,XMAX,GAM2*SNORM) 3654 IF (FO.LE.F) THEN 3655 IF (NRED1.GE.MRED) THEN 3656 ITERS=-1 3657 ELSE 3658 IDIR=1 3659 ITERS=0 3660 NRED1=NRED1+1 3661 END IF 3662 END IF 3663 END IF 3664 MAXST=0 3665 IF (XDEL.GE.XMAX) MAXST=1 3666 IF (MES3.EQ.0) THEN 3667 NRED=NRED1 3668 ELSE 3669 NRED=NRED2 3670 END IF 3671 ISYS=0 3672 RETURN 3673 END 3674* SUBROUTINE PS0L02 ALL SYSTEMS 97/12/01 3675* PURPOSE : 3676* EXTENDED LINE SEARCH WITHOUT DIRECTIONAL DERIVATIVES. 3677* 3678* PARAMETERS : 3679* RO R VALUE OF THE STEPSIZE PARAMETER. 3680* RO RO INITIAL VALUE OF THE STEPSIZE PARAMETER. 3681* RO RP PREVIOUS VALUE OF THE STEPSIZE PARAMETER. 3682* RO F VALUE OF THE OBJECTIVE FUNCTION. 3683* RI FO INITIAL VALUE OF THE OBJECTIVE FUNCTION. 3684* RO FP PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. 3685* RI PO INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE. 3686* RO PP PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE. 3687* RI FMIN LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. 3688* RI FMAX UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. 3689* RI RMIN MINIMUM VALUE OF THE STEPSIZE PARAMETER 3690* RI RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER 3691* RI TOLS TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE 3692* CHANGE OF THE FUNCTION VALUE). 3693* II KD DEGREE OF REQUIRED DERIVATIVES. 3694* IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES OF OBJECTIVE 3695* II NIT ACTUAL NUMBER OF ITERATIONS. 3696* II KIT NUMBER OF THE ITERATION AFTER LAST RESTART. 3697* IO NRED ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. 3698* II MRED MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. 3699* IO MAXST MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM 3700* STEPSIZE WAS NOT OR WAS REACHED. 3701* II IEST LOWER BOUND SPECIFICATION. IEST=0 OR IEST=1 IF LOWER BOUND 3702* IS NOT OR IS GIVEN. 3703* II INITS CHOICE OF THE INITIAL STEPSIZE. INITS=0-INITIAL STEPSIZE 3704* IS SPECIFIED IN THE CALLING PROGRAM. INITS=1-UNIT INITIAL 3705* STEPSIZE. INITS=2-COMBINED UNIT AND QUADRATICALLY ESTIMATED 3706* INITIAL STEPSIZE. INITS=3-QUADRATICALLY ESTIMATED INITIAL 3707* STEPSIZE. 3708* IO ITERS TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT 3709* LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY 3710* STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE. 3711* ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE. 3712* ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION. 3713* ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL 3714* DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION. 3715* II KTERS TERMINATION SELECTION. KTERS=1-PERFECT LINE SEARCH. 3716* KTERS=2-GOLDSTEIN STEPSIZE. KTERS=3-CURRY STEPSIZE. 3717* KTERS=4-EXTENDED CURRY STEPSIZE. KTERS=5-ARMIJO STEPSIZE. 3718* KTERS=6-FIRST STEPSIZE. 3719* II MES METHOD SELECTION. MES=1-BISECTION. MES=2-QUADRATIC 3720* INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE). 3721* MES=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL 3722* DERIVATIVES). MES=4-CUBIC INTERPOLATION. MES=5-CONIC 3723* INTERPOLATION. 3724* IU ISYS CONTROL PARAMETER. 3725* 3726* SUBPROGRAM USED : 3727* S PNINT3 EXTRAPOLATION OR INTERPOLATION WITHOUT DIRECTIONAL 3728* DERIVATIVES. 3729* 3730* METHOD : 3731* SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH EXTENDED TERMINATION 3732* CRITERIA. 3733* 3734 SUBROUTINE PS0L02(R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX,TOLS, 3735 & KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS,MES,ISYS) 3736 INTEGER KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS,MES, 3737 & ISYS 3738 DOUBLE PRECISION R,RO,RP,F,FO,FP,PO,PP,FMIN,FMAX,RMIN,RMAX,TOLS 3739 DOUBLE PRECISION RL,FL,RU,FU,RI,FI,RTEMP,TOL 3740 INTEGER MTYP,MERR,MODE,INIT1,MES1,MES2 3741 LOGICAL L1,L2,L3,L4,L6,L7 3742 PARAMETER(TOL=1.0D-2) 3743 SAVE MTYP,MODE,MES1,MES2 3744 SAVE RL,FL,RU,FU,RI,FI 3745 IF (ISYS.EQ.1) GO TO 3 3746 MES1=2 3747 MES2=2 3748 ITERS=0 3749 IF (PO.GE.0.0D 0) THEN 3750 R=0.0D 0 3751 ITERS=-2 3752 GO TO 4 3753 END IF 3754 IF (RMAX.LE.0.0D 0) THEN 3755 ITERS= 0 3756 GO TO 4 3757 END IF 3758* 3759* INITIAL STEPSIZE SELECTION 3760* 3761 IF (INITS.GT.0) THEN 3762 RTEMP=FMIN-F 3763 ELSE IF (IEST.EQ.0) THEN 3764 RTEMP=F-FP 3765 ELSE 3766 RTEMP=MAX(F-FP,FMIN-F) 3767 END IF 3768 INIT1=ABS(INITS) 3769 RP=0.0D 0 3770 FP=FO 3771 PP=PO 3772 IF (INIT1.EQ.0) THEN 3773 ELSE IF (INIT1.EQ.1.OR.INITS.GE.1.AND.IEST.EQ.0) THEN 3774 R=1.0D 0 3775 ELSE IF (INIT1.EQ.2) THEN 3776 R=MIN(1.0D 0,4.0D 0*RTEMP/PO) 3777 ELSE IF (INIT1.EQ.3) THEN 3778 R=MIN(1.0D 0, 2.0D 0*RTEMP/PO) 3779 ELSE IF (INIT1.EQ.4) THEN 3780 R=2.0D 0*RTEMP/PO 3781 END IF 3782 RTEMP=R 3783 R=MAX(R,RMIN) 3784 R=MIN(R,RMAX) 3785 MODE=0 3786 RL=0.0D 0 3787 FL=FO 3788 RU=0.0D 0 3789 FU=FO 3790 RI=0.0D 0 3791 FI=FO 3792* 3793* NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION) 3794* 3795 2 CALL PNINT3(RO,RL,RU,RI,FO,FL,FU,FI,PO,R,MODE,MTYP,MERR) 3796 IF (MERR.GT.0) THEN 3797 ITERS=-MERR 3798 GO TO 4 3799 ELSE IF (MODE.EQ.1) THEN 3800 NRED=NRED-1 3801 R=MIN(R,RMAX) 3802 ELSE IF (MODE.EQ.2) THEN 3803 NRED=NRED+1 3804 END IF 3805* 3806* COMPUTATION OF THE NEW FUNCTION VALUE 3807* 3808 KD= 0 3809 LD=-1 3810 ISYS=1 3811 RETURN 3812 3 CONTINUE 3813 IF (ITERS.NE.0) GO TO 4 3814 IF (F.LE.FMIN) THEN 3815 ITERS=7 3816 GO TO 4 3817 ELSE 3818 L1=R.LE.RMIN.AND.NIT.NE.KIT 3819 L2=R.GE.RMAX 3820 L3=F-FO.LE.TOLS*R*PO.OR.F-FMIN.LE.(FO-FMIN)/1.0D 1 3821 L4=F-FO.GE.(1.0D 0-TOLS)*R*PO.OR.MES2.EQ.2.AND.MODE.EQ.2 3822 L6=RU-RL.LE.TOL*RU.AND.MODE.EQ.2 3823 L7=MES2.LE.2.OR.MODE.NE.0 3824 MAXST=0 3825 IF (L2) MAXST=1 3826 END IF 3827* 3828* TEST ON TERMINATION 3829* 3830 IF (L1.AND..NOT.L3) THEN 3831 ITERS=0 3832 GO TO 4 3833 ELSE IF (L2.AND..NOT.F.GE.FU) THEN 3834 ITERS=7 3835 GO TO 4 3836 ELSE IF (L6) THEN 3837 ITERS=1 3838 GO TO 4 3839 ELSE IF (L3.AND.L7.AND.KTERS.EQ.5) THEN 3840 ITERS=5 3841 GO TO 4 3842 ELSE IF (L3.AND.L4.AND.L7.AND.(KTERS.EQ.2.OR.KTERS.EQ.3.OR. 3843 * KTERS.EQ.4)) THEN 3844 ITERS=2 3845 GO TO 4 3846 ELSE IF (KTERS.LT.0.OR.KTERS.EQ.6.AND.L7) THEN 3847 ITERS=6 3848 GO TO 4 3849 ELSE IF (ABS(NRED).GE.MRED) THEN 3850 ITERS=-1 3851 GO TO 4 3852 ELSE 3853 RP=R 3854 FP=F 3855 MODE=MAX(MODE,1) 3856 MTYP=ABS(MES) 3857 IF (F.GE.FMAX) MTYP=1 3858 END IF 3859 IF (MODE.EQ.1) THEN 3860* 3861* INTERVAL CHANGE AFTER EXTRAPOLATION 3862* 3863 RL=RI 3864 FL=FI 3865 RI=RU 3866 FI=FU 3867 RU=R 3868 FU=F 3869 IF (F.GE.FI) THEN 3870 NRED=0 3871 MODE=2 3872 ELSE IF ( MES1 .EQ. 1) THEN 3873 MTYP=1 3874 END IF 3875* 3876* INTERVAL CHANGE AFTER INTERPOLATION 3877* 3878 ELSE IF (R.LE.RI) THEN 3879 IF (F.LE.FI) THEN 3880 RU=RI 3881 FU=FI 3882 RI=R 3883 FI=F 3884 ELSE 3885 RL=R 3886 FL=F 3887 END IF 3888 ELSE 3889 IF (F.LE.FI) THEN 3890 RL=RI 3891 FL=FI 3892 RI=R 3893 FI=F 3894 ELSE 3895 RU=R 3896 FU=F 3897 END IF 3898 END IF 3899 GO TO 2 3900 4 ISYS=0 3901 RETURN 3902 END 3903* SUBROUTINE PS1L01 ALL SYSTEMS 97/12/01 3904* PURPOSE : 3905* STANDARD LINE SEARCH WITH DIRECTIONAL DERIVATIVES. 3906* 3907* PARAMETERS : 3908* RO R VALUE OF THE STEPSIZE PARAMETER. 3909* RO RP PREVIOUS VALUE OF THE STEPSIZE PARAMETER. 3910* RO F VALUE OF THE OBJECTIVE FUNCTION. 3911* RI FO INITIAL VALUE OF THE OBJECTIVE FUNCTION. 3912* RO FP PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. 3913* RO P VALUE OF THE DIRECTIONAL DERIVATIVE. 3914* RI PO INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE. 3915* RO PP PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE. 3916* RI FMIN LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. 3917* RI FMAX UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION. 3918* RI RMIN MINIMUM VALUE OF THE STEPSIZE PARAMETER 3919* RI RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER 3920* RI TOLS TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE 3921* CHANGE OF THE FUNCTION VALUE). 3922* RI TOLP TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE 3923* CHANGE OF THE DIRECTIONAL DERIVATIVE). 3924* RO PAR1 PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC 3925* UPDATES. 3926* RO PAR2 PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC 3927* UPDATES. 3928* II KD DEGREE OF REQUIRED DERIVATIVES. 3929* IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES OF OBJECTIVE 3930* II NIT ACTUAL NUMBER OF ITERATIONS. 3931* II KIT NUMBER OF THE ITERATION AFTER LAST RESTART. 3932* IO NRED ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. 3933* II MRED MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS. 3934* IO MAXST MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM 3935* STEPSIZE WAS NOT OR WAS REACHED. 3936* II IEST LOWER BOUND SPECIFICATION. IEST=0 OR IEST=1 IF LOWER BOUND 3937* IS NOT OR IS GIVEN. 3938* II INITS CHOICE OF THE INITIAL STEPSIZE. INITS=0-INITIAL STEPSIZE 3939* IS SPECIFIED IN THE CALLING PROGRAM. INITS=1-UNIT INITIAL 3940* STEPSIZE. INITS=2-COMBINED UNIT AND QUADRATICALLY ESTIMATED 3941* INITIAL STEPSIZE. INITS=3-QUADRATICALLY ESTIMATED INITIAL 3942* STEPSIZE. 3943* IO ITERS TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT 3944* LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY 3945* STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE. 3946* ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE. 3947* ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION. 3948* ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL 3949* DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION. 3950* II KTERS TERMINATION SELECTION. KTERS=1-PERFECT LINE SEARCH. 3951* KTERS=2-GOLDSTEIN STEPSIZE. KTERS=3-CURRY STEPSIZE. 3952* KTERS=4-EXTENDED CURRY STEPSIZE. KTERS=5-ARMIJO STEPSIZE. 3953* KTERS=6-FIRST STEPSIZE. 3954* II MES METHOD SELECTION. MES=1-BISECTION. MES=2-QUADRATIC 3955* INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE). 3956* MES=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL 3957* DERIVATIVES). MES=4-CUBIC INTERPOLATION. MES=5-CONIC 3958* INTERPOLATION. 3959* IU ISYS CONTROL PARAMETER. 3960* 3961* SUBPROGRAM USED : 3962* S PNINT1 EXTRAPOLATION OR INTERPOLATION WITH DIRECTIONAL 3963* DERIVATIVES. 3964* 3965* METHOD : 3966* SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH STANDARD TERMINATION 3967* CRITERIA. 3968* 3969 SUBROUTINE PS1L01(R,RP,F,FO,FP,P,PO,PP,FMIN,FMAX,RMIN,RMAX, 3970 & TOLS,TOLP,PAR1,PAR2,KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS, 3971 & ITERS,KTERS,MES,ISYS) 3972 INTEGER KD,LD,NIT,KIT,NRED,MRED,MAXST,IEST,INITS,ITERS,KTERS, 3973 & MES,ISYS 3974 DOUBLE PRECISION R,RP,F,FO,FP,P,PO,PP,FMIN,FMAX,RMIN,RMAX, 3975 & TOLS,TOLP,PAR1,PAR2 3976 DOUBLE PRECISION RL,FL,PL,RU,FU,PU,RTEMP 3977 INTEGER MTYP,MERR,MODE,INIT1,MES1,MES2,MES3 3978 LOGICAL L1,L2,L3,L5,L7,M1,M2,M3 3979 DOUBLE PRECISION CON,CON1 3980 PARAMETER (CON=1.0D-2,CON1=1.0D-13) 3981 SAVE MTYP,MODE,MES1,MES2,MES3 3982 SAVE RL,FL,PL,RU,FU,PU 3983 IF (ISYS.EQ.1) GO TO 3 3984 MES1=2 3985 MES2=2 3986 MES3=2 3987 ITERS=0 3988 IF (PO.GE.0.0D 0) THEN 3989 R=0.0D 0 3990 ITERS=-2 3991 GO TO 4 3992 END IF 3993 IF (RMAX.LE.0.0D 0) THEN 3994 ITERS=0 3995 GO TO 4 3996 END IF 3997* 3998* INITIAL STEPSIZE SELECTION 3999* 4000 IF (INITS.GT.0) THEN 4001 RTEMP=FMIN-F 4002 ELSE IF (IEST.EQ.0) THEN 4003 RTEMP=F-FP 4004 ELSE 4005 RTEMP=MAX(F-FP,FMIN-F) 4006 END IF 4007 INIT1=ABS(INITS) 4008 RP=0.0D 0 4009 FP=FO 4010 PP=PO 4011 IF (INIT1.EQ.0) THEN 4012 ELSE IF (INIT1.EQ.1.OR.INITS.GE.1.AND.IEST.EQ.0) THEN 4013 R=1.0D 0 4014 ELSE IF (INIT1.EQ.2) THEN 4015 R=MIN(1.0D 0,4.0D 0*RTEMP/PO) 4016 ELSE IF (INIT1.EQ.3) THEN 4017 R=MIN(1.0D 0, 2.0D 0*RTEMP/PO) 4018 ELSE IF (INIT1.EQ.4) THEN 4019 R=2.0D 0*RTEMP/PO 4020 END IF 4021 R=MAX(R,RMIN) 4022 R=MIN(R,RMAX) 4023 MODE=0 4024 RU=0.0D 0 4025 FU=FO 4026 PU=PO 4027* 4028* NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION) 4029* 4030 2 CALL PNINT1(RL,RU,FL,FU,PL,PU,R,MODE,MTYP,MERR) 4031 IF (MERR.GT.0) THEN 4032 ITERS=-MERR 4033 GO TO 4 4034 ELSE IF (MODE.EQ.1) THEN 4035 NRED=NRED-1 4036 R=MIN(R,RMAX) 4037 ELSE IF (MODE.EQ.2) THEN 4038 NRED=NRED+1 4039 END IF 4040* 4041* COMPUTATION OF THE NEW FUNCTION VALUE AND THE NEW DIRECTIONAL 4042* DERIVATIVE 4043* 4044 KD= 1 4045 LD=-1 4046 ISYS=1 4047 RETURN 4048 3 CONTINUE 4049 IF (MODE.EQ.0) THEN 4050 PAR1=P/PO 4051 PAR2=F-FO 4052 END IF 4053 IF (ITERS.NE.0) GO TO 4 4054 IF (F.LE.FMIN) THEN 4055 ITERS=7 4056 GO TO 4 4057 ELSE 4058 L1=R.LE.RMIN.AND.NIT.NE.KIT 4059 L2=R.GE.RMAX 4060 L3=F-FO.LE.TOLS*R*PO 4061 L5=P.GE.TOLP*PO.OR.MES2.EQ.2.AND.MODE.EQ.2 4062 L7=MES2.LE.2.OR.MODE.NE.0 4063 M1=.FALSE. 4064 M2=.FALSE. 4065 M3=L3 4066 IF (MES3.GE.1) THEN 4067 M1=ABS(P).LE.CON*ABS(PO).AND.FO-F.GE.(CON1/CON)*ABS(FO) 4068 L3=L3.OR.M1 4069 END IF 4070 IF (MES3.GE.2) THEN 4071 M2=ABS(P).LE.0.5D 0*ABS(PO).AND.ABS(FO-F).LE.2.0D 0*CON1*ABS(FO) 4072 L3=L3.OR.M2 4073 END IF 4074 MAXST=0 4075 IF (L2) MAXST=1 4076 END IF 4077* 4078* TEST ON TERMINATION 4079* 4080 IF (L1.AND..NOT.L3) THEN 4081 ITERS=0 4082 GO TO 4 4083 ELSE IF (L2.AND.L3.AND..NOT.L5) THEN 4084 ITERS=7 4085 GO TO 4 4086 ELSE IF (M3.AND.MES1.EQ.3) THEN 4087 ITERS=5 4088 GO TO 4 4089 ELSE IF (L3.AND.L5.AND.L7) THEN 4090 ITERS=4 4091 GO TO 4 4092 ELSE IF (KTERS.LT.0.OR.KTERS.EQ.6.AND.L7) THEN 4093 ITERS=6 4094 GO TO 4 4095 ELSE IF (ABS(NRED).GE.MRED) THEN 4096 ITERS=-1 4097 GO TO 4 4098 ELSE 4099 RP=R 4100 FP=F 4101 PP=P 4102 MODE=MAX(MODE,1) 4103 MTYP=ABS(MES) 4104 IF (F.GE.FMAX) MTYP=1 4105 END IF 4106 IF (MODE.EQ.1) THEN 4107* 4108* INTERVAL CHANGE AFTER EXTRAPOLATION 4109* 4110 RL=RU 4111 FL=FU 4112 PL=PU 4113 RU=R 4114 FU=F 4115 PU=P 4116 IF (.NOT.L3) THEN 4117 NRED=0 4118 MODE=2 4119 ELSE IF ( MES1 .EQ. 1) THEN 4120 MTYP=1 4121 END IF 4122 ELSE 4123* 4124* INTERVAL CHANGE AFTER INTERPOLATION 4125* 4126 IF (.NOT.L3) THEN 4127 RU=R 4128 FU=F 4129 PU=P 4130 ELSE 4131 RL=R 4132 FL=F 4133 PL=P 4134 END IF 4135 END IF 4136 GO TO 2 4137 4 ISYS=0 4138 RETURN 4139 END 4140* SUBROUTINE PS1L18 ALL SYSTEMS 99/12/01 4141* PURPOSE : 4142* SPECIAL LINE SEARCH FOR NONSMOOTH NONCONVEX VARIABLE METRIC METHOD. 4143* 4144* PARAMETERS : 4145* II N ACTUAL NUMBER OF VARIABLES. 4146* II MA DECLARED NUMBER OF LINEAR APPROXIMATED FUNCTIONS 4147* II MAL CURRENT NUMBER OF LINEAR APPROXIMATED FUNCTIONS. 4148* RU X(N) VECTOR OF VARIABLES. 4149* RO G(N) SUBGRADIENT OF THE OBJECTIVE FUNCTION. 4150* RI S(N) DIRECTION VECTOR. 4151* RU U(N) PREVIOUS VECTOR OF VARIABLES. 4152* RI AF(4*MA) VECTOR OF BUNDLE FUNCTIONS VALUES. 4153* RI AG(N*MA) MATRIX WHOSE COLUMNS ARE BUNDLE SUBGRADIENTS. 4154* RI AY(N*MA) MATRIX WHOSE COLUMNS ARE VARIABLE VECTORS. 4155* RO T VALUE OF THE STEPSIZE PARAMETER. 4156* RO TB BUNDLE PARAMETER FOR MATRIX SCALING. 4157* RO FO PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. 4158* RO F VALUE OF THE OBJECTIVE FUNCTION. 4159* RU PO PREVIOUS DIRECTIONAL DERIVATIVE. 4160* RU P DIRECTIONAL DERIVATIVE. 4161* RI TMIN MINIMUM VALUE OF THE STEPSIZE PARAMETER. 4162* RI TMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER. 4163* RI SNORM EUCLIDEAN NORM OF THE DIRECTION VECTOR. 4164* RI WK STOPPING PARAMETER. 4165* RI EPS1 TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE 4166* CHANGE OF THE FUNCTION VALUE). 4167* RI EPS2 TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE 4168* DIRECTIONAL DERIVATIVE). 4169* RI ETA5 DISTANCE MEASURE PARAMETER. 4170* RI ETA9 MAXIMUM FOR REAL NUMBERS. 4171* II KD DEGREE OF REQUIRED DERIVATIVES. 4172* IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES OF OBJECTIVE 4173* II JE EXTRAPOLATION INDICATOR. 4174* RI MOS3 LOCALITY MEASURE PARAMETER. 4175* IO ITERS NULL STEP INDICATOR. ITERS=0-NULL STEP. ITERS=1-DESCENT 4176* STEP. 4177* IU ISYS CONTROL PARAMETER. 4178* 4179* VARIABLES IN COMMON /STAT/ (STATISTICS) : 4180* IO NRES NUMBER OF RESTARTS. 4181* IO NDEC NUMBER OF MATRIX DECOMPOSITIONS. 4182* IO NIN NUMBER OF INNER ITERATIONS. 4183* IO NIT NUMBER OF ITERATIONS. 4184* IO NFV NUMBER OF FUNCTION EVALUATIONS. 4185* IO NFG NUMBER OF GRADIENT EVALUATIONS. 4186* IO NFH NUMBER OF HESSIAN EVALUATIONS. 4187* 4188* SUBPROGRAMS USED : 4189* S PNINT1 EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH 4190* S PNSTP4 STEPSIZE DETERMINATION FOR DESCENT STEPS. 4191* S PNSTP5 STEPSIZE DETERMINATION FOR NULL STEPS. 4192* WITH DIRECTIONAL DERIVATIVES. 4193* S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. 4194* RF MXVDOT DOT PRODUCT OF TWO VECTORS. 4195* 4196* METHOD : 4197* SPECIAL METHOD OF STEP LENGTH DETERMINATION. 4198* 4199 SUBROUTINE PS1L18(N,MA,MAL,X,G,S,U,AF,AG,AY,T,TB,FO,F,PO,P,TMIN, 4200 & TMAX,SNORM,WK,EPS1,EPS2,ETA5,ETA9,KD,LD,JE,MOS3,ITERS,ISYS) 4201 DOUBLE PRECISION EPS1,EPS2,ETA5,ETA9,F,FO,P,PO,SNORM,T,TB,TMAX, 4202 & TMIN,WK 4203 INTEGER ITERS,ISYS,JE,KD,LD,MA,MAL,MOS3,N 4204 DOUBLE PRECISION AF(*),AG(*),AY(*),G(*),S(*),U(*),X(*) 4205 DOUBLE PRECISION BET,FL,FU,PL,PU,TL,TU 4206 INTEGER IER 4207 DOUBLE PRECISION MXVDOT 4208 SAVE FL,FU,PL,PU,TL,TU 4209 IF (ISYS.GT.0) GO TO 25 4210 IF (JE.GT.0) T = DBLE(2-JE/99)*T 4211 IF (JE.LE.0) T = MIN(1.0D0,TMAX) 4212 IF (PO.EQ.0.0D0 .OR. JE.GT.0) GO TO 10 4213 IF (ITERS.EQ.1) THEN 4214 CALL PNSTP4(N,MA,MAL,X,AF,AG,AY,S,F,PO,T,TB,ETA5,ETA9,MOS3) 4215 ELSE 4216 CALL PNSTP5(N,MA,MAL,X,AF,AG,AY,S,F,PO,T,TB,ETA5,ETA9,MOS3) 4217 END IF 4218 10 T = MIN(MAX(T,TMIN),TMAX) 4219 TL = 0.0D0 4220 TU = T 4221 FL = FO 4222 PL = PO 4223* 4224* FUNCTION AND GRADIENT EVALUATION AT A NEW POINT 4225* 4226 20 CALL MXVDIR(N,T,S,U,X) 4227 KD= 1 4228 LD=-1 4229 ISYS=1 4230 RETURN 4231 25 CONTINUE 4232 P = MXVDOT(N,G,S) 4233* 4234* NULL/DESCENT STEP TEST (ITERS=0/1) 4235* 4236 ITERS = 1 4237 IF (F.LE.FO-T* (EPS1+EPS1)*WK) THEN 4238 TL = T 4239 FL = F 4240 PL = P 4241 ELSE 4242 TU = T 4243 FU = F 4244 PU = P 4245 END IF 4246 BET = MAX(ABS(FO-F+P*T),ETA5* (SNORM*T)**MOS3) 4247 IF (F.LE.FO-T*EPS1*WK .AND. (T.GE.TMIN.OR. 4248 & BET.GT.EPS1*WK)) GO TO 40 4249 IF (P-BET.GE.-EPS2*WK .OR. TU-TL.LT.TMIN*1.0D-1) GO TO 30 4250 IF (TL.EQ.0.0D0 .AND. PL.LT.0.0D0) THEN 4251 CALL PNINT1(TL,TU,FL,FU,PL,PU,T,2,2,IER) 4252 ELSE 4253 T = 5.0D-1* (TU+TL) 4254 END IF 4255 GO TO 20 4256 30 ITERS = 0 4257 40 CONTINUE 4258 ISYS=0 4259 RETURN 4260 END 4261* SUBROUTINE PUBBM1 ALL SYSTEMS 97/12/01 4262* PURPOSE : 4263* PARTITIONED VARIABLE METRIC UPDATE. 4264* 4265* PARAMETERS : 4266* II NA NUMBER OF BLOCKS OF THE MATRIX H. 4267* RU AH(MB) APPROXIMATION OF THE PARTITIONED HESSIAN MATRIX. 4268* RI IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 4269* RI JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 4270* RA S(NF) AUXILIARY VECTOR. 4271* RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. 4272* RI AGO(MA) GRADIENTS DIFFERENCE. 4273* RI ETA0 MACHINE PRECISION. 4274* RI ETA9 MAXIMUM MACHINE NUMBER. 4275* IU ICOR SWITCH BETWEEN UPDATES. ICOR=0-THE BFGS UPDATE. 4276* ICOR=1-THE RANK ONE UPDATE. 4277* II NIT ACTUAL NUMBER OF ITERATIONS. 4278* II KIT NUMBER OF THE ITERATION AFTER LAST RESTART. 4279* IO ITERH TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION. 4280* ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS. 4281* II MET METHOD SELECTION. MET=0-NO UPDATE. MET=1-BFGS UPDATE. 4282* MET=2-COMBINATION OF BFGS AND RANK-ONE UPDATES. 4283* II MET1 SELECTION OF SELF SCALING. MET1=1-SELF SCALING SUPPRESSED. 4284* MET1=2 SELF SCALING IN THE FIRST ITERATION AFTER RESTART. 4285* MET1=3-SELF SCALING IN EACH ITERATION. 4286* 4287* SUBPROGRAMS USED : 4288* S MXBSBM MULTIPLICATION OF A PARTITIONED MATRIX BY A VECTOR. 4289* S MXBSBU UPDATE OF A PARTITIONED MATRIX. 4290* S MXDSMS SCALING OF A DENSE SYMMETRIC MATRIX. 4291* S MXWDIR VECTOR AUGMENTED BY THE SCALED VECTOR. 4292* RF MXWDOT DOT PRODUCT OF TWO SPARSE VECTORS. 4293* 4294 SUBROUTINE PUBBM1(NA,AH,IAG,JAG,S,XO,AGO,ETA0,ETA9,ICOR,NIT,KIT, 4295 & ITERH,MET,MET1) 4296 INTEGER NA,IAG(*),JAG(*),ICOR,NIT,KIT,ITERH,MET, 4297 & MET1 4298 DOUBLE PRECISION AH(*),S(*),XO(*),AGO(*),ETA0,ETA9 4299 DOUBLE PRECISION A,B,C,GAM,POM,DEN,MXWDOT 4300 INTEGER K,L,KA,NB,INEG 4301 LOGICAL L1,L3 4302 IF (MET.LE.0) GO TO 22 4303 L1=MET1.GE.3.OR.MET1.EQ.2.AND.NIT.EQ.KIT 4304 L3=.NOT.L1 4305 NB=0 4306 INEG=0 4307 DO 21 KA=1,NA 4308 K=IAG(KA) 4309 L=IAG(KA+1)-K 4310* 4311* DETERMINATION OF THE PARAMETERS B, C 4312* 4313 B=MXWDOT(L,JAG(K),AGO(K),XO,2) 4314 IF (MET.EQ.1) THEN 4315 IF (B.LE.1.0D 0/ETA9) GO TO 20 4316 ELSE 4317 IF (ABS(B).LE.1.0D 0/ETA9) GO TO 20 4318 END IF 4319 A=0.0D 0 4320 CALL MXBSBM(L,AH(NB+1),JAG(K),XO,S,1) 4321 C=MXWDOT(L,JAG(K),XO,S,1) 4322 IF (ABS(C).LE.1.0D 0/ETA9) GO TO 20 4323 IF (L1) THEN 4324* 4325* DETERMINATION OF THE PARAMETER GAM (SELF SCALING) 4326* 4327 GAM=C/B 4328 IF (L3) THEN 4329 GAM=1.0D 0 4330 END IF 4331 ELSE 4332 GAM=1.0D 0 4333 END IF 4334 IF (MET.EQ.1) THEN 4335* 4336* BFGS UPDATE 4337* 4338 POM=0.0D 0 4339 CALL MXBSBU(L,AH(NB+1),JAG(K),GAM/B,AGO(K),2) 4340 CALL MXBSBU(L,AH(NB+1),JAG(K),-1.0D 0/C,S,1) 4341 ELSE 4342 IF (B.LT.0.0D 0) INEG=INEG+1 4343 IF (ICOR.GT.0) THEN 4344* 4345* RANK ONE UPDATE 4346* 4347 DEN=GAM*B-C 4348 IF (ABS(DEN).GT.ETA0*ABS(C)) THEN 4349 POM=GAM*B/DEN 4350 CALL MXWDIR(L,JAG(K),-GAM,AGO(K),S,AGO(K),2) 4351 CALL MXBSBU(L,AH(NB+1),JAG(K), 1.0D 0/DEN,AGO(K),2) 4352 ELSE 4353 GO TO 20 4354 END IF 4355 ELSE IF (B.LT.0.0D 0) THEN 4356 GO TO 20 4357 ELSE 4358* 4359* BFGS UPDATE 4360* 4361 POM=0.0D 0 4362 CALL MXBSBU(L,AH(NB+1),JAG(K),GAM/B,AGO(K),2) 4363 CALL MXBSBU(L,AH(NB+1),JAG(K),-1.0D 0/C,S,1) 4364 END IF 4365 END IF 4366 ITERH=0 4367 IF (GAM.NE.1.0D 0) THEN 4368 CALL MXDSMS(L,AH(NB+1),1.0D 0/GAM) 4369 END IF 4370 20 CONTINUE 4371 NB=NB+L*(L+1)/2 4372 21 CONTINUE 4373 IF (INEG.GE.NA/2) ICOR=1 4374 22 CONTINUE 4375 RETURN 4376 END 4377* SUBROUTINE PUBBM2 ALL SYSTEMS 97/12/01 4378* PURPOSE : 4379* PARTITIONED VARIABLE METRIC UPDATE. 4380* 4381* PARAMETERS : 4382* II NA NUMBER OF BLOCKS OF THE MATRIX H. 4383* RU AH(MB) APPROXIMATION OF THE PARTITIONED HESSIAN MATRIX. 4384* RI IAG(NA+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 4385* RI JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 4386* RA S(NF) AUXILIARY VECTOR. 4387* RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. 4388* RI AGO(MA) GRADIENTS DIFFERENCE. 4389* RI ETA0 MACHINE PRECISION. 4390* RI ETA9 MAXIMUM MACHINE NUMBER. 4391* II NIT ACTUAL NUMBER OF ITERATIONS. 4392* II KIT NUMBER OF THE ITERATION AFTER LAST RESTART. 4393* IO ITERH TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION. 4394* ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS. 4395* II MET VARIABLE METRIC UPDATE. MET=1-THE BFGS UPDATE. MET=2-THE 4396* DFP UPDATE. MET=3-THE HOSHINO UPDATE. MET=4-THE RANK ONE 4397* UPDATE. 4398* II MET1 SELECTION OF SELF SCALING. MET1=1-SELF SCALING SUPPRESSED. 4399* MET1=2 SELF SCALING IN THE FIRST ITERATION AFTER RESTART. 4400* MET1=3-CONTROLLED SELF SCALING. 4401* II MET3 CORRECTION OF THE UPDATE. MET3=1-CORRECTION IS SUPPRESSED. 4402* MET3=2-THE POWELL UPDATE. 4403* 4404* SUBPROGRAMS USED : 4405* S MXBSBM MULTIPLICATION OF A PARTITIONED MATRIX BY A VECTOR. 4406* S MXBSBU UPDATE OF A PARTITIONED MATRIX. 4407* S MXDSMS SCALING OF A DENSE SYMMETRIC MATRIX. 4408* S MXWDIR VECTOR AUGMENTED BY THE SCALED VECTOR. 4409* RF MXWDOT DOT PRODUCT OF TWO SPARSE VECTORS. 4410* 4411 SUBROUTINE PUBBM2(NA,AH,IAG,JAG,S,XO,AGO,ETA0,ETA9,NIT,KIT,ITERH, 4412 & MET,MET1,MET3) 4413 INTEGER NA,IAG(*),JAG(*),NIT,KIT,ITERH,MET,MET1,MET3 4414 DOUBLE PRECISION AH(*),S(*),XO(*),AGO(*),ETA0,ETA9 4415 DOUBLE PRECISION A,B,C,GAM,POM,DEN,DIS,MXWDOT 4416 INTEGER K,L,KA,NB 4417 LOGICAL L1,L3 4418 DOUBLE PRECISION CON,CON1,CON2 4419 PARAMETER (CON=0.1D 0,CON1=0.5D 0,CON2=4.0D 0) 4420 L1=MET1.GE.3.OR.MET1.EQ.2.AND.NIT.EQ.KIT 4421 L3=.NOT.L1 4422 NB=0 4423 DO 21 KA=1,NA 4424 K=IAG(KA) 4425 L=IAG(KA+1)-K 4426* 4427* DETERMINATION OF THE PARAMETERS B, C 4428* 4429 B=MXWDOT(L,JAG(K),AGO(K),XO,2) 4430 IF (MET3.EQ.1) THEN 4431 IF (B.LE.1.0D 0/ETA9) GO TO 20 4432 ELSE 4433 IF (ABS(B).LE.1.0D 0/ETA9) GO TO 20 4434 END IF 4435 A=0.0D 0 4436 CALL MXBSBM(L,AH(NB+1),JAG(K),XO,S,1) 4437 C=MXWDOT(L,JAG(K),XO,S,1) 4438 IF (MET3.EQ.3) THEN 4439 IF (ABS(C).LE.1.0D 0/ETA9) GO TO 20 4440 ELSE 4441 IF (C.LE.1.0D 0/ETA9) GO TO 20 4442 END IF 4443 IF (MET3.EQ.2) THEN 4444 IF (B.LE.0.0D 0) THEN 4445* 4446* POWELL'S CORRECTION 4447* 4448 DIS=(1.0D 0-CON)*C/(C-B) 4449 CALL MXWDIR(L,JAG(K),-1.0D 0,AGO(K),S,AGO(K),2) 4450 CALL MXWDIR(L,JAG(K),-DIS,AGO(K),S,AGO(K),2) 4451 B=C+DIS*(B-C) 4452 END IF 4453 END IF 4454 IF (L1) THEN 4455* 4456* DETERMINATION OF THE PARAMETER GAM (SELF SCALING) 4457* 4458 GAM=C/B 4459 IF (MET1.EQ.3) THEN 4460 IF (NIT.NE.KIT) THEN 4461 L3=GAM.LT.CON1.OR.GAM.GT.CON2 4462 END IF 4463 ELSE IF (MET1.EQ.4) THEN 4464 GAM=MAX(1.0D 0,GAM) 4465 END IF 4466 IF (L3) THEN 4467 GAM=1.0D 0 4468 END IF 4469 ELSE 4470 GAM=1.0D 0 4471 END IF 4472 IF (MET.EQ.1) THEN 4473 GO TO 18 4474 ELSE IF (MET.EQ.2) THEN 4475* 4476* DFP UPDATE 4477* 4478 DEN=GAM*B+C 4479 DIS=GAM+C/B 4480 POM=1.0D 0 4481 CALL MXWDIR(L,JAG(K),-DIS,AGO(K),S,AGO(K),2) 4482 CALL MXBSBU(L,AH(NB+1),JAG(K),1.0D 0/DEN,AGO(K),2) 4483 CALL MXBSBU(L,AH(NB+1),JAG(K),-1.0D 0/DEN,S,1) 4484 GO TO 19 4485 ELSE IF (MET.EQ.3) THEN 4486* 4487* HOSHINO UPDATE 4488* 4489 DEN=GAM*B+C 4490 DIS=0.5D 0*B 4491 POM=GAM*B/DEN 4492 CALL MXBSBU(L,AH(NB+1),JAG(K),GAM/DIS,AGO(K),2) 4493 CALL MXWDIR(L,JAG(K),GAM,AGO(K),S,AGO(K),2) 4494 CALL MXBSBU(L,AH(NB+1),JAG(K),-1.0D 0/DEN,AGO(K),2) 4495 GO TO 19 4496 ELSE IF (MET.EQ.4) THEN 4497* 4498* RANK ONE UPDATE 4499* 4500 DEN=GAM*B-C 4501 IF (MET3.EQ.3) THEN 4502 IF (ABS(DEN).LE.ETA0*ABS(C)) GO TO 18 4503 ELSE 4504 IF (DEN.LE.ETA0*C) GO TO 18 4505 END IF 4506 POM=GAM*B/DEN 4507 CALL MXWDIR(L,JAG(K),-GAM,AGO(K),S,AGO(K),2) 4508 CALL MXBSBU(L,AH(NB+1),JAG(K), 1.0D 0/DEN,AGO(K),2) 4509 GO TO 19 4510 END IF 4511 18 CONTINUE 4512* 4513* BFGS UPDATE 4514* 4515 POM=0.0D 0 4516 CALL MXBSBU(L,AH(NB+1),JAG(K),GAM/B,AGO(K),2) 4517 CALL MXBSBU(L,AH(NB+1),JAG(K),-1.0D 0/C,S,1) 4518 19 CONTINUE 4519 ITERH=0 4520 IF (GAM.NE.1.0D 0) THEN 4521 CALL MXDSMS(L,AH(NB+1),1.0D 0/GAM) 4522 END IF 4523 20 CONTINUE 4524 NB=NB+L*(L+1)/2 4525 21 CONTINUE 4526 RETURN 4527 END 4528* SUBROUTINE PUBVI2 ALL SYSTEMS 04/12/01 4529* PURPOSE : 4530* NONSMOOTH VARIABLE METRIC UPDATE OF THE INVERSE HESSIAN MATRIX. 4531* 4532* PARAMETERS : 4533* II NF ACTUAL NUMBER OF VARIABLES. 4534* II NA NUMBER OF APPROXIMATED FUNCTIONS. 4535* II MA NUMBER OF ELEMENTS IN THE FIELD AG. 4536* II MB NUMBER OF NONZERO ELEMENTS OF THE PARTITIONED HESSIAN MATRIX. 4537* RU AH(MB) NUMERICAL VALUES OF ELEMENTS OF THE PARTITIONED HESSIAN 4538* MATRIX. 4539* II IAG(NA+1) POINTERS OF THE JACOBIAN MATRIX. 4540* RI JAG(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 4541* RI AG(NF) NEW GENERALIZED JACOBIAN MATRIX. 4542* RI AGO(NF) OLD GENERALIZED JACOBIAN MATRIX. 4543* RI XO(N) VECTOR OF VARIABLES DIFFERENCE. 4544* RO S(NF) AUXILIARY VECTOR. 4545* RO U(NF) AUXILIARY VECTOR. 4546* RI ETA9 MAXIMUM MACHINE NUMBER. 4547* II NNK CONSECUTIVE NULL STEPS COUNTER. 4548* II NIT ACTUAL NUMBER OF ITERATIONS. 4549* IO ITERH TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION. 4550* ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS. 4551* 4552* SUBPROGRAMS USED : 4553* S MXBSBM MULTIPLICATION OF A DENSE SYMMETRIC MATRIX BY A VECTOR. 4554* S MXBSBU UPDATE OF A PARTITIONED SYMMETRIC MATRIX. 4555* S MXDSMS SCALING OF A DENSE SYMMETRIC MATRIX. 4556* S MXVDIF DIFFERENCE OF TWO VECTORS. 4557* S MXWDIR VECTOR AUGMENTED BY THE SCALED VECTOR. 4558* RF MXWDOT DOT PRODUCT OF VECTORS. 4559* 4560 SUBROUTINE PUBVI2(NA,AH,IAG,JAG,AG,AGO,XO,S,U,ETA9,NNK,NIT,ITERH) 4561 INTEGER NA,IAG(*),JAG(*),NNK,NIT,ITERH 4562 DOUBLE PRECISION AH(*),AG(*),AGO(*),XO(*),S(*),U(*),ETA9 4563 DOUBLE PRECISION GAM,A,B,C,Q,DEN,POM,MXWDOT 4564 INTEGER KA,K,L,NB,INEG 4565 LOGICAL LB,LR 4566 NB=0 4567 INEG=0 4568 DO 21 KA=1,NA 4569 K=IAG(KA) 4570 L=IAG(KA+1)-K 4571 CALL MXVDIF(L,AG(K),AGO(K),U) 4572* 4573* DETERMINATION OF THE PARAMETERS B, C 4574* 4575 B=MXWDOT(L,JAG(K),U,XO,2) 4576 IF (ABS(B).LE.1.0D 0/ETA9) GO TO 20 4577 A=0.0D 0 4578 CALL MXBSBM(L,AH(NB+1),JAG(K),XO,S,1) 4579 C=MXWDOT(L,JAG(K),XO,S,1) 4580 IF (ABS(C).LE.1.0D 0/ETA9) GO TO 20 4581 GAM=1.0D 0 4582 IF (NIT.EQ.1) THEN 4583 Q=1.0D 0 4584 IF (C.NE.0.0D 0) Q=C/B 4585 IF ((Q-2.5D-1)*(Q-3.0D 0).GT.0.0D 0) GAM=MIN(3.0D 0, 4586 & MAX(2.0D-2,Q)) 4587 END IF 4588 IF (B.LT.0.0D 0) INEG=INEG+1 4589 LB=NNK.EQ.0 4590 LR=NNK.NE.0.AND.C.LT.GAM*B 4591 IF (LB)THEN 4592 IF (B.LT.0.0D 0) GO TO 20 4593* 4594* BFGS UPDATE 4595* 4596 POM=0.0D 0 4597 CALL MXBSBU(L,AH(NB+1),JAG(K),GAM/B,U,2) 4598 CALL MXBSBU(L,AH(NB+1),JAG(K),-1.0D 0/C,S,1) 4599 ITERH=0 4600 IF (GAM.NE.1.0D 0) THEN 4601 CALL MXDSMS(L,AH(NB+1),1.0D 0/GAM) 4602 END IF 4603 ELSE IF (LR) THEN 4604 DEN=GAM*B-C 4605 POM=GAM*B/DEN 4606 CALL MXWDIR(L,JAG(K),-GAM,U,S,U,2) 4607 CALL MXBSBU(L,AH(NB+1),JAG(K), 1.0D 0/DEN,U,2) 4608 END IF 4609 20 CONTINUE 4610 NB=NB+L*(L+1)/2 4611 21 CONTINUE 4612 RETURN 4613 END 4614* SUBROUTINE PULCI3 ALL SYSTEMS 96/12/01 4615* PURPOSE : 4616* LIMITED STORAGE INVERSE COLUMN UPDATE METHODS. 4617* 4618* PARAMETERS : 4619* II N NUMBER OF VARIABLES. 4620* RI A(IAG(N+1)-1) SPARSE RECTANGULAR MATRIX WHICH IS USED FOR THE 4621* DIRECTION VECTOR DETERMINATION. 4622* II IA(N+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD AG. 4623* II JA(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD AG. 4624* IU IP(N) PERMUTATION VECTOR. 4625* IU ID(N) POSITION OF THE DIAGONAL ELEMENTS IN THE FIELD AG. 4626* RU XM(N*MF) SET OF VECTORS FOR INVERSE COLUMN UPDATE. 4627* RU GM(MF) SET OF VALUES FOR INVERSE COLUMN UPDATE. 4628* IU IM(MF) SET OF INDICES FOR INVERSE COLUMN UPDATE. 4629* RA XO(N) AUXILIARY VECTOR. 4630* RI AFO(N) GRADIENTS DIFERENCES. 4631* RO S(N) DIRECTION VECTOR. 4632* II MF NUMBER OF VARIABLE METRIC UPDATES. 4633* II NIT NUMBER OF ITERATIONS. 4634* II KIT NUMBER OF THE ITERATION AFTER LAST RESTART. 4635* IO ITERH TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION. 4636* ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS. 4637* IU IREST RESTART INDICATOR. 4638* 4639* SUBPROGRAMS USED : 4640* S MXLIIM MATRIX MULTIPLICATION FOR LIMITED STORAGE INVERSE 4641* COLUMN UPDATE METHOD. 4642* S MXVDIR VECTOR AUGMENTED BY THE SCALED VECTOR. 4643* RF MXVMX1 DOT PRODUCT OF VECTORS. 4644* 4645* METHOD : 4646* LIMITED STORAGE VARIABLE METRIC METHODS. 4647* 4648 SUBROUTINE PULCI3(N,A,IA,JA,IP,ID,XM,GM,IM,XO,AFO,S,MF,NIT,KIT, 4649 + ITERH,IREST) 4650 INTEGER IREST,ITERH,NIT,KIT,MF,N 4651 DOUBLE PRECISION A(*),AFO(*),GM(*),S(*),XM(*),XO(*) 4652 INTEGER IA(*),ID(*),IM(*),IP(*),JA(*) 4653 DOUBLE PRECISION TEMP 4654 INTEGER II,MA,MM 4655 DOUBLE PRECISION MXVMX1 4656 MA = IA(N+1) - 1 4657 MM = MIN(NIT-KIT,MF) 4658 IF (MM.GE.MF) THEN 4659 ITERH = 1 4660 IREST = 1 4661 ELSE 4662 II = N*MM + 1 4663 CALL MXLIIM(N,MM,A(MA+1),IA,JA,IP,ID,XM,GM,IM,AFO,XM(II),S) 4664 CALL MXVDIR(N,-1.0D0,XM(II),XO,XM(II)) 4665 MM = MM + 1 4666 TEMP = MXVMX1(N,AFO,II) 4667 IF (TEMP.LE.0.0D0) THEN 4668 ITERH = 2 4669 ELSE 4670 IM(MM) = II 4671 GM(MM) = AFO(II) 4672 ITERH = 0 4673 END IF 4674 END IF 4675 RETURN 4676 END 4677* SUBROUTINE PULSP3 ALL SYSTEMS 02/12/01 4678* PURPOSE : 4679* LIMITED STORAGE VARIABLE METRIC UPDATE. 4680* 4681* PARAMETERS : 4682* II N NUMBER OF VARIABLES (NUMBER OF ROWS OF XM). 4683* II M NUMBER OF COLUMNS OF XM. 4684* II MF MAXIMUM NUMBER OF COLUMNS OF XM. 4685* RI XM(N*M) RECTANGULAR MATRIX IN THE PRODUCT FORM SHIFTED BROYDEN 4686* METHOD (STORED COLUMNWISE): H-SIGMA*I=XM*TRANS(XM) 4687* RO GR(M) MATRIX TRANS(XM)*GO. 4688* RU XO(N) VECTORS OF VARIABLES DIFFERENCE XO AND VECTOR XO-TILDE. 4689* RU GO(N) GRADIENT DIFFERENCE GO AND VECTOR XM*TRANS(XM)*GO. 4690* RI R STEPSIZE PARAMETER. 4691* RI PO OLD DIRECTIONAL DERIVATIVE (MULTIPLIED BY R) 4692* RU SIG SCALING PARAMETER (ZETA AND SIGMA). 4693* IO ITERH TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION. 4694* ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS. 4695* II MET3 CHOICE OF SIGMA (1-CONSTANT, 2-QUADRATIC EQUATION). 4696* 4697* SUBPROGRAMS USED : 4698* S MXDRMM MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR 4699* MATRIX BY A VECTOR. 4700* S MXDCMU UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX. 4701* WITH CONTROLLING OF POSITIVE DEFINITENESS. 4702* S MXVDIR VECTOR AUGMENTED BY A SCALED VECTOR. 4703* RF MXVDOT DOT PRODUCT OF VECTORS. 4704* S MXVSCL SCALING OF A VECTOR. 4705* 4706* METHOD : 4707* SHIFTED BFGS METHOD IN THE PRODUCT FORM. 4708* 4709 SUBROUTINE PULSP3(N,M,MF,XM,GR,XO,GO,R,PO,SIG,ITERH,MET3) 4710 INTEGER N,M,MF,ITERH,MET3 4711 DOUBLE PRECISION XM(*),GR(*),XO(*),GO(*),R,PO,SIG 4712 DOUBLE PRECISION DEN,POM,A,B,C,AA,AH,BB,PAR,MXVDOT 4713 IF (M.GE.MF) RETURN 4714 B=MXVDOT(N,XO,GO) 4715 IF (B.LE.0.0D 0) THEN 4716 ITERH=2 4717 GO TO 22 4718 END IF 4719 CALL MXDRMM(N,M,XM,GO,GR) 4720 AH=MXVDOT(N,GO,GO) 4721 AA=MXVDOT(M,GR,GR) 4722 A=AA+AH*SIG 4723 C=-R*PO 4724* 4725* DETERMINATION OF THE PARAMETER SIG (SHIFT) 4726* 4727 PAR=1.0D 0 4728 POM=B/AH 4729 IF (A.GT.0.0D 0) THEN 4730 DEN=MXVDOT(N,XO,XO) 4731 IF (MET3.LE.4) THEN 4732 SIG=SQRT(MAX(0.0D 0,1.0D 0-AA/A))/(1.0D 0+ 4733 & SQRT(MAX(0.0D 0,1.0D 0-B*B/(DEN*AH))))*POM 4734 ELSE 4735 SIG=SQRT(MAX(0.0D 0,SIG*AH/A))/(1.0D 0+ 4736 & SQRT(MAX(0.0D 0,1.0D 0-B*B/(DEN*AH))))*POM 4737 END IF 4738 SIG=MAX(SIG,2.0D-1*POM) 4739 SIG=MIN(SIG,8.0D-1*POM) 4740 ELSE 4741 SIG=2.5D-1*POM 4742 END IF 4743* 4744* COMPUTATION OF SHIFTED XO AND SHIFTED B 4745* 4746 BB=B-AH*SIG 4747 CALL MXVDIR(N,-SIG,GO,XO,XO) 4748* 4749* BFGS-BASED SHIFTED BFGS UPDATE 4750* 4751 POM=1.0D 0 4752 CALL MXDCMU(N,M,XM,-1.0D 0/BB,XO,GR) 4753 CALL MXVSCL(N,SQRT(PAR/BB),XO,XM(N*M+1)) 4754 M=M+1 4755 22 CONTINUE 4756 ITERH=0 4757 RETURN 4758 END 4759* SUBROUTINE PULVP3 ALL SYSTEMS 03/12/01 4760* PURPOSE : 4761* RANK-TWO LIMITED-STORAGE VARIABLE-METRIC METHODS IN THE PRODUCT FORM. 4762* 4763* PARAMETERS : 4764* II N NUMBER OF VARIABLES (NUMBER OF ROWS OF XM). 4765* II M NUMBER OF COLUMNS OF XM. 4766* RI XM(N*M) RECTANGULAR MATRIX IN THE PRODUCT FORM SHIFTED BROYDEN 4767* METHOD (STORED COLUMNWISE): H-SIGMA*I=XM*TRANS(XM) 4768* RO XR(M) VECTOR TRANS(XM)*H**(-1)*XO. 4769* RO GR(M) MATRIX TRANS(XM)*GO. 4770* RA S(N) AUXILIARY VECTORS (H**(-1)*XO AND U). 4771* RA SO(N) AUXILIARY VECTORS ((H-SIGMA*I)*H**(-1)*XO AND V). 4772* RU XO(N) VECTORS OF VARIABLES DIFFERENCE XO AND VECTOR XO-TILDE. 4773* RU GO(N) GRADIENT DIFFERENCE GO AND VECTOR XM*TRANS(XM)*GO. 4774* RI R STEPSIZE PARAMETER. 4775* RI PO OLD DIRECTIONAL DERIVATIVE (MULTIPLIED BY R) 4776* RU SIG SCALING PARAMETER (ZETA AND SIGMA). 4777* IO ITERH TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION. 4778* ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS. 4779* II MET2 CHOICE OF THE CORRECTION PARAMETER (1-THE UNIT VALUE, 4780* 2-THE BALANCING VALUE, 3-THE SQUARE ROOT, 4-THE GEOMETRIC 4781* MEAN). 4782* II MET3 CHOICE OF THE SHIFT PARAMETER (4-THE FIRST FORMULA, 4783* 5-THE SECOND FORMULA). 4784* II MET5 CHOICE OF THE METHOD (1-RANK-ONE METHOD, 2-RANK-TWO 4785* METHOD). 4786* 4787* SUBPROGRAMS USED : 4788* S MXDRMM MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR 4789* MATRIX BY A VECTOR. 4790* S MXDCMU UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX. 4791* WITH CONTROLLING OF POSITIVE DEFINITENESS. RANK-ONE FORMULA. 4792* S MXDCMV UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX. 4793* WITH CONTROLLING OF POSITIVE DEFINITENESS. RANK-TWO FORMULA. 4794* S MXVDIR VECTOR AUGMENTED BY A SCALED VECTOR. 4795* RF MXVDOT DOT PRODUCT OF VECTORS. 4796* S MXVLIN LINEAR COMBINATION OF TWO VECTORS. 4797* S MXVSCL SCALING OF A VECTOR. 4798* 4799* METHOD : 4800* RANK-ONE LIMITED-STORAGE VARIABLE-METRIC METHOD IN THE PRODUCT FORM. 4801* 4802 SUBROUTINE PULVP3(N,M,XM,XR,GR,S,SO,XO,GO,R,PO,SIG,ITERH,MET2, 4803 & MET3,MET5) 4804 INTEGER N,M,ITERH,MET2,MET3,MET5 4805 DOUBLE PRECISION XM(*),XR(*),GR(*),S(*),SO(*),XO(*),GO(*), 4806 & R,PO,SIG 4807 DOUBLE PRECISION MXVDOT 4808 DOUBLE PRECISION DEN,POM,A,B,C,AA,BB,CC,AH,PAR,ZET 4809 ZET=SIG 4810* 4811* COMPUTATION OF B 4812* 4813 B=MXVDOT(N,XO,GO) 4814 IF (B.LE.0.0D 0) THEN 4815 ITERH=2 4816 GO TO 22 4817 END IF 4818* 4819* COMPUTATION OF GR=TRANS(XM)*GO, XR=TRANS(XM)*H**(-1)*XO 4820* AND S=H**(-1)*XO, SO=(H-SIGMA*I)*H**(-1)*XO. COMPUTATION 4821* OF AA=GR*GR, BB=GR*XR, CC=XR*XR. COMPUTATION OF A AND C. 4822* 4823 CALL MXDRMM(N,M,XM,GO,GR) 4824 CALL MXVSCL(N,R,S,S) 4825 CALL MXDRMM(N,M,XM,S,XR) 4826 CALL MXVDIR(N,-SIG,S,XO,SO) 4827 AH=MXVDOT(N,GO,GO) 4828 AA=MXVDOT(M,GR,GR) 4829 BB=MXVDOT(M,GR,XR) 4830 CC=MXVDOT(M,XR,XR) 4831 A=AA+AH*SIG 4832 C=-R*PO 4833* 4834* DETERMINATION OF THE PARAMETER SIG (SHIFT) 4835* 4836 POM=B/AH 4837 IF (A.GT.0.0D 0) THEN 4838 DEN=MXVDOT(N,XO,XO) 4839 IF (MET3.LE.4) THEN 4840 SIG=SQRT(MAX(0.0D 0,1.0D 0-AA/A))/(1.0D 0+ 4841 & SQRT(MAX(0.0D 0,1.0D 0-B*B/(DEN*AH))))*POM 4842 ELSE 4843 SIG=SQRT(MAX(0.0D 0,SIG*AH/A))/(1.0D 0+ 4844 & SQRT(MAX(0.0D 0,1.0D 0-B*B/(DEN*AH))))*POM 4845 END IF 4846 SIG=MAX(SIG,2.0D-1*POM) 4847 SIG=MIN(SIG,8.0D-1*POM) 4848 ELSE 4849 SIG=2.5D-1*POM 4850 END IF 4851* 4852* COMPUTATION OF SHIFTED XO AND SHIFTED B 4853* 4854 B=B-AH*SIG 4855 CALL MXVDIR(N,-SIG,GO,XO,XO) 4856* 4857* COMPUTATION OF THE PARAMETER RHO (CORRECTION) 4858* 4859 IF (MET2.LE.1) THEN 4860 PAR=1.0D 0 4861 ELSE IF (MET2.EQ.2) THEN 4862 PAR=SIG*AH/B 4863 ELSE IF (MET2.EQ.3) THEN 4864 PAR=SQRT(1.0D 0-AA/A) 4865 ELSE IF (MET2.EQ.4) THEN 4866 PAR=SQRT(SQRT(1.0D 0-AA/A)*(SIG*AH/B)) 4867 ELSE 4868 PAR=ZET/(ZET+SIG) 4869 END IF 4870* 4871* COMPUTATION OF THE PARAMETER THETA (BFGS) 4872* 4873 POM=SIGN(SQRT(PAR*B/CC),BB) 4874* 4875* COMPUTATION OF Q AND P 4876* 4877 IF (MET5.EQ.1) THEN 4878* 4879* RANK ONE UPDATE OF XM 4880* 4881 CALL MXVDIR(M,POM,XR,GR,XR) 4882 CALL MXVLIN(N,PAR,XO,POM,SO,S) 4883 CALL MXDCMU(N,M,XM,-1.0D 0/(PAR*B+POM*BB),S,XR) 4884 ELSE 4885* 4886* RANK TWO UPDATE OF XM 4887* 4888 CALL MXVDIR(N,PAR/POM-BB/B,XO,SO,S) 4889 CALL MXDCMV(N,M,XM,-1.0D 0/B,XO,GR,-1.0D 0/CC,S,XR) 4890 END IF 4891 22 CONTINUE 4892 ITERH=0 4893 RETURN 4894 END 4895* SUBROUTINE PUSMM1 ALL SYSTEMS 97/12/01 4896* PURPOSE : 4897* VARIABLE METRIC UPDATE OF A SPARSE SYMMETRIC POSITIVE DEFINITE MATRIX 4898* USING THE MARWIL PROJECTION. 4899* 4900* PARAMETERS : 4901* II NF NUMBER OF VARIABLES. 4902* RU H(M) POSITIVE DEFINITE APPROXIMATION OF THE SPARSE HESSIAN 4903* MATRIX. 4904* II IH(NF) POINTERS OF THE DIAGONAL ELEMENTS OF H. 4905* II JH(M) INDICES OF THE NONZERO ELEMENTS OF H. 4906* RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 4907* RA XS(NF) AUXILIARY VECTOR. 4908* RA S(NF) AUXILIARY VECTOR. 4909* RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. 4910* RI GO(NF) GRADIENTS DIFFERENCE. 4911* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. 4912* RO R VALUE OF THE STEPSIZE PARAMETER. 4913* RI PO INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE. 4914* II NIT ACTUAL NUMBER OF ITERATIONS. 4915* II KIT NUMBER OF THE ITERATION AFTER LAST RESTART. 4916* II MET1 SELECTION OF SELF SCALING. MET1=1-SELF SCALING SUPPRESSED. 4917* MET1=2-INITIAL SELF SCALING. MET1=3-SELF SCALING IN EACH 4918* ITERATION. 4919* II ITERD CAUSE OF TERMINATION IN THE DIRECTION DETERMINATION. 4920* ITERD<0-BAD DECOMPOSITION. ITERD=0-DESCENT DIRECTION. 4921* ITERD=1-NEWTON LIKE STEP. ITERD=2-INEXACT NEWTON LIKE STEP. 4922* ITERD=3-BOUNDARY STEP. ITERD=4-DIRECTION WITH THE NEGATIVE 4923* CURVATURE. ITERD=5-MARQUARDT STEP. 4924* IO ITERH TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION. 4925* ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS. 4926* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. 4927* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. 4928* 4929* SUBPROGRAMS USED : 4930* S MXSSMM MATRIX-VECTOR PRODUCT. 4931* S MXSSMY MARWILL CORRECTION OF A SPARSE SYMMETRIC MATRIX. 4932* S MXUDIF DIFFERENCE OF TWO VECTORS. 4933* S MXUDIR VECTOR AUGMENTED BY THE SCALED VECTOR. 4934* RF MXUDOT DOT PRODUCT OF VECTORS. 4935* S MXVSCL SCALING OF A VECTOR. 4936* 4937 SUBROUTINE PUSMM1(NF,H,IH,JH,G,XS,S,XO,GO,IX,R,PO,NIT,KIT, 4938 & MET1,ITERD,ITERH,KBF) 4939 INTEGER NF,IH(*),JH(*),IX(*),NIT,KIT,MET1,ITERD,ITERH,KBF 4940 DOUBLE PRECISION H(*),G(*),S(*),XO(*),GO(*),XS(*),R,PO 4941 INTEGER MM 4942 DOUBLE PRECISION MXUDOT 4943 DOUBLE PRECISION A,B,C,GAM 4944 LOGICAL L1 4945 MM=IH(NF+1)-1 4946* 4947* DETERMINATION OF THE PARAMETER C AND THE VECTOR S 4948* 4949 A=0.0D 0 4950 L1=MET1.GE.3.OR.MET1.GE.2.AND.NIT.EQ.KIT 4951 IF (ITERD.NE.1) THEN 4952 CALL MXSSMM(NF,H,IH,JH,XO,S) 4953 IF (L1) C=MXUDOT(NF,XO,S,IX,KBF) 4954 ELSE 4955 CALL MXUDIF(NF,GO,G,S,IX,KBF) 4956 CALL MXVSCL(NF,R,S,S) 4957 IF (L1) C=-R*PO 4958 END IF 4959 GAM=1.0D 0 4960 IF (L1) THEN 4961* 4962* SELF SCALING 4963* 4964 B=MXUDOT(NF,XO,GO,IX,KBF) 4965 IF (B.GT.0.0D 0.AND.C.GT.0.0D 0) THEN 4966 GAM=C/B 4967 CALL MXVSCL(MM,1.0D 0/GAM,H,H) 4968 CALL MXVSCL(NF,1.0D 0/GAM,S,S) 4969 END IF 4970 END IF 4971 CALL MXUDIR(NF,-1.0D 0,S,GO,S,IX,KBF) 4972* 4973* RANK-ONE UPDATE PROJECTED USING MXSSMY 4974* 4975 CALL MXSSMY(NF,H,IH,JH,XS,S,XO) 4976 ITERH=0 4977 RETURN 4978 END 4979* SUBROUTINE PUSSD5 ALL SYSTEMS 97/12/01 4980* PURPOSE : 4981* INITIATION OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX 4982* 4983* PARAMETERS : 4984* II NA NUMBER OF APPROXIMATED FUNCTIONS. 4985* RI AF(NA) VECTOR CONTAINING VALUES OF THE APPROXIMATED 4986* FUNCTIONS. 4987* RU AH(MB) POSITIVE DEFINITE APPROXIMATION OF THE PARTITIONED 4988* HESSIAN MATRIX. 4989* II IAG(NA+1) POINTERS OF THE SPARSE JACOBIAN MATRIX. 4990* II JAG(MA) COLUMN INDICES OF THE SPARSE JACOBIAN MATRIX. 4991* RU H(M) POSITIVE DEFINITE APPROXIMATION OF THE SPARSE HESSIAN 4992* MATRIX 4993* II IH(NF+1) POINTERS OF THE DIAGONAL ELEMENTS OF THE SPARSE 4994* HESSIAN MATRIX. 4995* II JH(M) INDICES OF THE NONZERO ELEMENTS OF THE SPARSE HESSIAN 4996* MATRIX IN THE PACKED ROW FORM. 4997* 4998* SUBPROGRAMS USED : 4999* S PASSH2 COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE 5000* PARTITIONED HESSIAN MATRIX. 5001* 5002 SUBROUTINE PUSSD5(NA,AF,AH,IAG,JAG,H,IH,JH) 5003 INTEGER NA,IAG(*),JAG(*),IH(*),JH(*) 5004 DOUBLE PRECISION AF(*),AH(*),H(*) 5005 INTEGER K,KA,L,LL,NB 5006 NB=0 5007 DO 2 KA=1,NA 5008 K=IAG(KA) 5009 L=IAG(KA+1)-K 5010 LL=L*(L+1)/2 5011 CALL PASSH2(H,IH,JH,AH(NB+1),IAG,JAG,KA,AF(KA)) 5012 NB=NB+LL 5013 2 CONTINUE 5014 RETURN 5015 END 5016* SUBROUTINE PYABU1 ALL SYSTEMS 04/12/01 5017* PURPOSE : 5018* SUBGRADIENT AGGREGATION FOR NONSMOOTH VARIABLE METRIC METHOD. 5019* 5020* PARAMETERS : 5021* II NF NUMBER OF VARIABLES. 5022* RI H(M) POSITIVE DEFINITE APPROXIMATION OF THE SPARSE HESSIAN 5023* MATRIX. 5024* IO JH(M) INDICES OF THE NONZERO ELEMENTS OF H. 5025* II PSL(NF+1) POINTER ARRAY OF THE FACTORIZED SPARSE MATRIX 5026* II PERM(NF) PERMUTATION VECTOR 5027* RI G(NF) NEW SUBGRADIENT OF THE OBJECTIVE FUNCTION. 5028* RI GO(NF) OLD SUBGRADIENT OF THE OBJECTIVE FUNCTION. 5029* RU GV(NF) AGGREGATED SUBGRADIENT OF THE OBJECTIVE FUNCTION. 5030* RI S(NF) DIRECTION VECTOR. 5031* RA U(NF) AUXILIARY VECTOR. 5032* RA V(NF) AUXILIARY VECTOR. 5033* RO ALF LINEARIZATION TERM. 5034* RU ALFV AGGREGATED LINEARIZATION TERM. 5035* RI RHO CORRECTION PARAMETER. 5036* II JC CORRECTION INDICATOR. 5037* 5038* SUBPROGRAMS USED : 5039* S MXSPCB BACK SUBSTITUTION USING THE SPARSE DECOMPOSITION 5040* OBTAINED BY MXSPCF. 5041* RF MXVDOT DOT PRODUCT OF TWO VECTORS. 5042* S MXVSBP INVERSE PERMUTATION OF A VECTOR 5043* S MXVSFP PERMUTATION OF A VECTOR. 5044* 5045 SUBROUTINE PYABU1(NF,H,JH,PSL,PERM,G,GO,GV,S,U,V,ALF,ALFV,RHO, 5046 & JC) 5047 INTEGER NF,JH(*),PSL(*),PERM(*),JC 5048 DOUBLE PRECISION H(*),G(*),GO(*),GV(*),S(*),U(*),V(*),ALF,ALFV, 5049 & RHO 5050 DOUBLE PRECISION A,B,ALFM,LAM1,LAM2,PQ,PR,PRQR,QQP,QR,RR,RRP,RRQ, 5051 & W,W1 5052 INTEGER I 5053 DOUBLE PRECISION ZERO,ONE,MXVDOT 5054 PARAMETER (ZERO=0.0D 0,ONE=1.0D 0) 5055 ALFM=ZERO 5056* 5057* General routine - here always input parameter ALFM=0 5058* 5059 RR=ALFV+ALFV 5060 RRP=ALFV-ALFM 5061 RRQ=ALFV-ALF 5062 DO 1 I=1,NF 5063 A=S(I) 5064 U(I)=GO(I)-GV(I) 5065 S(I)=G(I)-GV(I) 5066 RR=RR-A*GV(I) 5067 RRP=RRP+A*U(I) 5068 RRQ=RRQ+A*S(I) 5069 1 CONTINUE 5070 PQ=ZERO 5071 PR=ZERO 5072 QR=ZERO 5073 PRQR=ZERO 5074 QQP=ZERO 5075 IF (JC.GE.1) THEN 5076 DO 2 I=1,NF 5077 PQ=PQ+RHO*(S(I)-U(I))**2 5078 PR=PR+RHO*U(I)**2 5079 QR=QR+RHO*S(I)**2 5080 PRQR=PRQR+RHO*U(I)*S(I) 5081 QQP=QQP+RHO+G(I)*(S(I)-U(I)) 5082 2 CONTINUE 5083 END IF 5084 QQP=QQP+ALF-ALFM 5085 CALL MXVSFP(NF,PERM,U,V) 5086 CALL MXSPCB(NF,H,PSL,JH,U,1) 5087 CALL MXVSFP(NF,PERM,S,V) 5088 CALL MXSPCB(NF,H,PSL,JH,S,1) 5089 DO 4 I=1,NF 5090 W1=ONE/H(PSL(I)+I-1) 5091 PQ=PQ+W1*(S(I)-U(I))**2 5092 PR=PR+W1*U(I)**2 5093 QR=QR+W1*S(I)**2 5094 PRQR=PRQR+W1*U(I)*S(I) 5095 S(I)=W1*(S(I)-U(I)) 5096 4 CONTINUE 5097 CALL MXSPCB(NF,H,PSL,JH,S,-1) 5098 CALL MXVSBP(NF,PERM,S,V) 5099 QQP=QQP+MXVDOT(NF,G,S) 5100 IF (PR.LE.ZERO.OR.QR.LE.ZERO) GO TO 10 5101 A=RRQ/QR 5102 B=PRQR/QR 5103 W=PRQR*B-PR 5104 IF (W.EQ.ZERO) GO TO 10 5105 LAM1=(A*PRQR-RRP)/W 5106 LAM2=A-LAM1*B 5107 IF (LAM1*(LAM1-ONE).LT.ZERO.AND.LAM2*(LAM1+LAM2-ONE).LT.ZERO) 5108 & GO TO 40 5109* 5110* MINIMUM ON THE BOUNDARY 5111* 5112 10 LAM1=ZERO 5113 LAM2=ZERO 5114 IF (ALF.LE.ALFV) LAM2=ONE 5115 IF (QR.GT.ZERO) LAM2=MIN(ONE,MAX(ZERO,RRQ/QR)) 5116 W=(LAM2*QR-RRQ-RRQ)*LAM2 5117 A=ZERO 5118 IF (ALFM.LE.ALFV) A=ONE 5119 IF (PR.GT.ZERO) A=MIN(ONE,MAX(ZERO,RRP/PR)) 5120 B=(A*PR-RRP-RRP)*A 5121 IF (B.LT.W)THEN 5122 W=B 5123 LAM1=A 5124 LAM2=ZERO 5125 END IF 5126 IF (QQP*(QQP-PQ).GE.ZERO) GO TO 40 5127 IF (QR-RRQ-RRQ-QQP*QQP/PQ.GE.W) GO TO 40 5128 LAM1=QQP/PQ 5129 LAM2=ONE-LAM1 5130 40 IF (LAM1.EQ.ZERO.AND.LAM2*(LAM2-ONE).LT.ZERO.AND.RRP-LAM2*PRQR 5131 & .GT.ZERO.AND.PR.GT.ZERO) LAM1=MIN(ONE-LAM2,(RRP-LAM2*PRQR)/PR) 5132 A=ONE-LAM1-LAM2 5133 ALFV=LAM1*ALFM+LAM2*ALF+A*ALFV 5134 DO 5 I=1,NF 5135 GV(I)=LAM1*GO(I)+LAM2*G(I)+A*GV(I) 5136 5 CONTINUE 5137 RETURN 5138 END 5139* SUBROUTINE PYABU2 ALL SYSTEMS 04/12/01 5140* PURPOSE : 5141* SIMPLIFIED AGGREGATION FOR NONSMOOTH VARIABLE METRIC METHOD. 5142* 5143* PARAMETERS : 5144* II NF NUMBER OF VARIABLES. 5145* RI H(M) POSITIVE DEFINITE APPROXIMATION OF THE SPARSE HESSIAN 5146* MATRIX. 5147* IO JH(M) INDICES OF THE NONZERO ELEMENTS OF H. 5148* II PSL(NF+1) POINTER ARRAY OF THE FACTORIZED SPARSE MATRIX 5149* II PERM(NF) PERMUTATION VECTOR 5150* RI G(NF) ACTUAL SUBGRADIENT OF THE OBJECTIVE FUNCTION. 5151* RU GV(NF) AGGREGATED SUBGRADIENT OF THE OBJECTIVE FUNCTION. 5152* RA S(NF) DIRECTION VECTOR. 5153* RA V(NF) AUXILIARY VECTOR. 5154* RO ALF LINEARIZATION TERM. 5155* RU ALFV AGGREGATED LINEARIZATION TERM. 5156* RI RHO CORRECTION PARAMETER. 5157* II JC CORRECTION INDICATOR. 5158* 5159* SUBPROGRAMS USED : 5160* S MXSPCB BACK SUBSTITUTION USING THE SPARSE DECOMPOSITION 5161* OBTAINED BY MXSPCF. 5162* S MXVSFP PERMUTATION OF A VECTOR. 5163* 5164 SUBROUTINE PYABU2(NF,H,JH,PSL,PERM,G,GV,S,V,ALF,ALFV,RHO,JC) 5165 INTEGER NF,JH(*),PSL(*),PERM(NF),JC 5166 DOUBLE PRECISION H(*),G(*),GV(*),S(*),V(*),ALF,ALFV,RHO 5167 DOUBLE PRECISION P,Q,W,LAM 5168 INTEGER I 5169 DOUBLE PRECISION ZERO,ONE 5170 PARAMETER (ZERO=0.0D 0,ONE=1.0D 0) 5171 P=ALFV-ALF 5172 DO 1 I=1,NF 5173 W=S(I) 5174 P=P+W*S(I) 5175 S(I)=G(I)-GV(I) 5176 1 CONTINUE 5177 Q=ZERO 5178 IF (JC.GE.1) THEN 5179 DO 2 I=1,NF 5180 Q=Q+RHO*S(I)**2 5181 2 CONTINUE 5182 END IF 5183 CALL MXVSFP(NF,PERM,S,V) 5184 CALL MXSPCB(NF,H,PSL,JH,S,1) 5185 DO 4 I=1,NF 5186 W=ONE/H(PSL(I)+I-1) 5187 Q=Q+W*S(I)**2 5188 4 CONTINUE 5189 LAM=0.5D 0+SIGN(0.5D 0,P) 5190 IF (Q.GT.ZERO) LAM=MIN(ONE,MAX(ZERO,P/Q)) 5191 P=ONE-LAM 5192 ALFV=LAM*ALF+P*ALFV 5193 DO 5 I=1,NF 5194 GV(I)=LAM*G(I)+P*GV(I) 5195 5 CONTINUE 5196 RETURN 5197 END 5198* SUBROUTINE PYADC0 ALL SYSTEMS 98/12/01 5199* PURPOSE : 5200* NEW SIMPLE BOUNDS ARE ADDED TO THE ACTIVE SET. 5201* 5202* PARAMETERS : 5203* II NF DECLARED NUMBER OF VARIABLES. 5204* II N REDUCED NUMBER OF VARIABLES. 5205* RI X(NF) VECTOR OF VARIABLES. 5206* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. 5207* RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. 5208* RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. 5209* IO INEW NUMBER OF ACTIVE CONSTRAINTS. 5210* 5211 SUBROUTINE PYADC0(NF,N,X,IX,XL,XU,INEW) 5212 INTEGER NF,N,IX(NF),INEW 5213 DOUBLE PRECISION X(*),XL(*),XU(*) 5214 INTEGER I,II,IXI 5215 N=NF 5216 INEW=0 5217 DO 1 I=1,NF 5218 II=IX(I) 5219 IXI=ABS(II) 5220 IF (IXI.GE.5) THEN 5221 IX(I)=-IXI 5222 ELSE IF ((IXI.EQ.1.OR.IXI.EQ.3.OR.IXI.EQ.4).AND.X(I).LE.XL(I)) 5223 & THEN 5224 X(I)=XL(I) 5225 IF (IXI.EQ.4) THEN 5226 IX(I)=-3 5227 ELSE 5228 IX(I)=-IXI 5229 END IF 5230 N=N-1 5231 IF (II.GT.0) INEW=INEW+1 5232 ELSE IF ((IXI.EQ.2.OR.IXI.EQ.3.OR.IXI.EQ.4).AND.X(I).GE.XU(I)) 5233 & THEN 5234 X(I)=XU(I) 5235 IF (IXI.EQ.3) THEN 5236 IX(I)=-4 5237 ELSE 5238 IX(I)=-IXI 5239 END IF 5240 N=N-1 5241 IF (II.GT.0) INEW=INEW+1 5242 END IF 5243 1 CONTINUE 5244 RETURN 5245 END 5246* SUBROUTINE PYBUN1 ALL SYSTEMS 97/12/01 5247* PURPOSE : 5248* BUNDLE UPDATING. 5249* 5250* PARAMETERS : 5251* II N ACTUAL NUMBER OF VARIABLES. 5252* II MB DECLARED NUMBER OF LINEAR APPROXIMATED FUNCTIONS. 5253* II NB CURRENT NUMBER OF LINEAR APPROXIMATED FUNCTIONS. 5254* RU X(N) VECTOR OF VARIABLES. 5255* RO G(N) SUBGRADIENT OF THE OBJECTIVE FUNCTION. 5256* RO F VALUE OF THE OBJECTIVE FUNCTION. 5257* RI AY(N*MB) MATRIX WHOSE COLUMNS ARE VARIABLE VECTORS. 5258* RI AG(N*MB) MATRIX WHOSE COLUMNS ARE BUNDLE SUBGRADIENTS. 5259* RI AF(4*MB) VECTOR OF BUNDLE FUNCTIONS VALUES. 5260* IO ITERS NULL STEP INDICATOR. ITERS=0-NULL STEP. ITERS=1-DESCENT 5261* STEP. 5262* 5263* SUBPROGRAMS USED : 5264* S MXVCOP COPYING OF A VECTOR. 5265* 5266 SUBROUTINE PYBUN1(N,MB,NB,X,G,F,AY,AG,AF,ITERS) 5267 INTEGER N,MB,NB,ITERS 5268 DOUBLE PRECISION X(*),G(*),F,AY(*),AG(*),AF(*) 5269 INTEGER I,IND,K,KN,L 5270 L=0 5271 IF (ITERS.EQ.0) L=1 5272* 5273* BUNDLE REDUCTION 5274* 5275 KN=0 5276 IF (NB.GE.MB) THEN 5277 DO 2 K=1,NB-1 5278 KN=K*N-N 5279 DO 1 I=1,N 5280 IF (G(I).NE.AG(KN+I)) GO TO 2 5281 1 CONTINUE 5282 IND=K 5283 GO TO 3 5284 2 CONTINUE 5285 IND=1 5286 3 DO 4 K=IND,NB-1 5287 AF(K)=AF(K+1) 5288 AF(K+MB*3)=AF(K+1+MB*3) 5289 KN=K*N+1 5290 CALL MXVCOP(N,AG(KN),AG(KN-N)) 5291 CALL MXVCOP(N,AY(KN),AY(KN-N)) 5292 4 CONTINUE 5293 NB=NB-1 5294 END IF 5295* 5296* BUNDLE COMPLETION 5297* 5298 IF (L.GT.0.AND.KN.EQ.0) THEN 5299 AF(NB+1)=AF(NB) 5300 AF(3*MB+NB+1)=AF(3*MB+NB) 5301 KN=NB*N+1 5302 CALL MXVCOP(N,AG(KN-N),AG(KN)) 5303 CALL MXVCOP(N,AY(KN-N),AY(KN)) 5304 END IF 5305 NB=NB+1 5306 KN=NB-L 5307 AF(KN)=F 5308 AF(KN+MB*3)=L 5309 K=(KN-1)*N+1 5310 CALL MXVCOP(N,G,AG(K)) 5311 CALL MXVCOP(N,X,AY(K)) 5312 RETURN 5313 END 5314* SUBROUTINE PYCSER ALL SYSTEMS 98/12/01 5315* PURPOSE : 5316* GROUP OF THE SAME COLOUR FOR THE POWELL-TOINT ALGORITHM FOR SPARSE 5317* HESSIANS APPROXIMATIONS IS CREATED. 5318* 5319* PARAMETERS : 5320* IU IH(MCOLS+1) POINTER VECTOR OF SPARSE HESSIAN MATRIX. 5321* IU JH(M) INDEX VECTOR OF THE HESSIAN MATRIX. 5322* IA WN02(MCOLS) AUXILIARY VECTOR. 5323* RA WN03(MCOLS) AUXILIARY VECTOR. 5324* RI DEG(MCOLS) DEGREES OF THE ADJACENCY GRAPH. 5325* IA WN01(NF) AUXILIARY VECTOR USED FOR INDICES OF THE COLUMNS 5326* THAT HAVE NOT BEEN COLOURED YET. 5327* II COL(NF) VECTOR DISCERNING GROUPS OF THE HESSIAN COLUMN OF THE 5328* SAME COLOUR. 5329* IU NCOL NUMBER OF COLOURS USED SO FAR. 5330* IU CNM NUMBER OF COLUMNS THAT HAVE NOT BEEN COLOURED SO FAR. 5331* 5332 SUBROUTINE PYCSER(JH,IH,WN02,WN03,DEG,WN01,COL,NCOL,CNM) 5333 INTEGER JH(*),IH(*),COL(*) 5334 INTEGER WN01(*),WN02(*) 5335 DOUBLE PRECISION WN03(*),DEG(*) 5336 INTEGER NCOL,CNM,I,J,K,L,IP 5337* 5338* DEFINITION OF THE INCIDENCE ARRAY A 5339* 5340 L=WN01(1) 5341* 5342* ELEMENT WAS MARKED THAT IT IS INSERTED 5343* 5344 DO 100 I=IH(L),IH(L+1)-1 5345 K=JH(I) 5346* 5347* COLUMN OF THIS NUMBER HAS APPEARED IN ONE OF THE PREVIOUS GROUPS 5348* 5349 IF (COL(K).LT.NCOL) GO TO 100 5350 DEG(K)=DEG(K)-1 5351 WN02(K)=NCOL 5352100 CONTINUE 5353* 5354* COLUMN IS INSERTED 5355* 5356 COL(L)=NCOL 5357* 5358* THE CYCLE OF COMPARING COLUMN WITH THE ARRAY A 5359* A2 IS AN HELP ARRAY CONTAINING COLUMNS THAT ARE 5360* BEEING EXAMINED BUT THAT WERE NOT YET ACCEPTED 5361* P IS ITS POINTER 5362* 5363 IF (CNM.EQ.1) GO TO 250 5364 DO 200 I=2,CNM 5365* 5366* TRANSFORMATION OF THE EXAMINED COLUMN I IS 5367* 5368 IP=1 5369 L=WN01(I) 5370 DO 300 J=IH(L),IH(L+1)-1 5371 K=JH(J) 5372 IF (COL(K).LT.NCOL) GO TO 300 5373 IF (WN02(K).GE.NCOL) GO TO 200 5374 WN03(IP)=K 5375 IP=IP+1 5376300 CONTINUE 5377 IF (IP.NE.1) THEN 5378* 5379* COPY OF THE WN03 ARRAY INTO WN02 FOR THE COLUMN WAS ACCEPTED 5380* 5381 DO 400 K=1,IP-1 5382 WN02(INT(WN03(K)))=NCOL 5383 DEG(INT(WN03(K)))=DEG(INT(WN03(K)))-1 5384400 CONTINUE 5385 END IF 5386* 5387* INSERT THE COLUMN INTO THE PROCESSED GROUP 5388* 5389 COL(L)=NCOL 5390* 5391* END OF THE MAIN CYCLE 5392* 5393200 CONTINUE 5394* 5395* JUMP LABEL 5396* 5397250 CONTINUE 5398* 5399* INVP SHIFT 5400* 5401 K=1 5402 DO 500 I=1,CNM 5403 L=WN01(I) 5404 IF (COL(L).EQ.NCOL) THEN 5405 ELSE 5406 WN01(K)=L 5407 K=K+1 5408 END IF 5409500 CONTINUE 5410* 5411* CNM UPDATE 5412* 5413 CNM=K-1 5414 RETURN 5415 END 5416* SUBROUTINE PYFUT1 ALL SYSTEMS 98/12/01 5417* PURPOSE : 5418* TERMINATION CRITERIA AND TEST ON RESTART. 5419* 5420* PARAMETERS : 5421* II N ACTUAL NUMBER OF VARIABLES. 5422* RI F NEW VALUE OF THE OBJECTIVE FUNCTION. 5423* RI FO OLD VALUE OF THE OBJECTIVE FUNCTION. 5424* RI UMAX MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER. 5425* RO GMAX NORM OF THE TRANSFORMED GRADIENT. 5426* RI DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES. 5427* RI TOLX LOWER BOUND FOR STEPLENGTH. 5428* RI TOLF LOWER BOUND FOR FUNCTION DECREASE. 5429* RI TOLB LOWER BOUND FOR FUNCTION VALUE. 5430* RI TOLG LOWER BOUND FOR GRADIENT. 5431* II KD DEGREE OF REQUIRED DERIVATIVES. 5432* IU NIT ACTUAL NUMBER OF ITERATIONS. 5433* II KIT NUMBER OF THE ITERATION AFTER RESTART. 5434* II MIT MAXIMUM NUMBER OF ITERATIONS. 5435* IU NFV ACTUAL NUMBER OF COMPUTED FUNCTION VALUES. 5436* II MFV MAXIMUM NUMBER OF COMPUTED FUNCTION VALUES. 5437* IU NFG ACTUAL NUMBER OF COMPUTED GRADIENT VALUES. 5438* II MFG MAXIMUM NUMBER OF COMPUTED GRADIENT VALUES. 5439* IU NTESX ACTUAL NUMBER OF TESTS ON STEPLENGTH. 5440* II MTESX MAXIMUM NUMBER OF TESTS ON STEPLENGTH. 5441* IU NTESF ACTUAL NUMBER OF TESTS ON FUNCTION DECREASE. 5442* II MTESF MAXIMUM NUMBER OF TESTS ON FUNCTION DECREASE. 5443* II IRES1 RESTART SPECIFICATION. RESTART IS PERFORMED AFTER 5444* IRES1*N+IRES2 ITERATIONS. 5445* II IRES2 RESTART SPECIFICATION. RESTART IS PERFORMED AFTER 5446* IRES1*N+IRES2 ITERATIONS. 5447* IU IREST RESTART INDICATOR. RESTART IS PERFORMED IF IREST>0. 5448* II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. 5449* ITERS=0 FOR ZERO STEP. 5450* IO ITERM TERMINATION INDICATOR. ITERM=1-TERMINATION AFTER MTESX 5451* UNSUFFICIENT STEPLENGTHS. ITERM=2-TERMINATION AFTER MTESF 5452* UNSUFFICIENT FUNCTION DECREASES. ITERM=3-TERMINATION ON LOWER 5453* BOUND FOR FUNCTION VALUE. ITERM=4-TERMINATION ON LOWER BOUND 5454* FOR GRADIENT. ITERM=11-TERMINATION AFTER MAXIMUM NUMBER OF 5455* ITERATIONS. ITERM=12-TERMINATION AFTER MAXIMUM NUMBER OF 5456* COMPUTED FUNCTION VALUES. 5457* 5458 SUBROUTINE PYFUT1(N,F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLB,TOLG,KD, 5459 & NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF,ITES,IRES1, 5460 & IRES2,IREST,ITERS,ITERM) 5461 INTEGER N,KD,NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF, 5462 & ITES,IRES1,IRES2,IREST,ITERS,ITERM 5463 DOUBLE PRECISION F,FO,UMAX,GMAX,DMAX,TOLX,TOLF,TOLG,TOLB 5464 DOUBLE PRECISION TEMP 5465 IF (ITERM.LT.0) RETURN 5466 IF (ITES .LE.0) GO TO 1 5467 IF (ITERS.EQ.0) GO TO 1 5468 IF (NIT.LE.0) FO=F+MIN(SQRT(ABS(F)),ABS(F)/1.0D 1) 5469 IF (F.LE.TOLB) THEN 5470 ITERM = 3 5471 RETURN 5472 END IF 5473 IF (KD.GT.0) THEN 5474 IF (GMAX.LE.TOLG.AND.UMAX.LE.TOLG) THEN 5475 ITERM = 4 5476 RETURN 5477 END IF 5478 END IF 5479 IF (NIT.LE.0) THEN 5480 NTESX = 0 5481 NTESF = 0 5482 END IF 5483 IF (DMAX.LE.TOLX) THEN 5484 ITERM = 1 5485 NTESX = NTESX+1 5486 IF (NTESX.GE.MTESX) RETURN 5487 ELSE 5488 NTESX = 0 5489 END IF 5490 TEMP=ABS(FO-F)/MAX(ABS(F),1.0D 0) 5491 IF (TEMP.LE.TOLF) THEN 5492 ITERM = 2 5493 NTESF = NTESF+1 5494 IF (NTESF.GE.MTESF) RETURN 5495 ELSE 5496 NTESF = 0 5497 END IF 5498 1 IF (NIT.GE.MIT) THEN 5499 ITERM = 11 5500 RETURN 5501 END IF 5502 IF (NFV.GE.MFV) THEN 5503 ITERM = 12 5504 RETURN 5505 END IF 5506 IF (NFG.GE.MFG) THEN 5507 ITERM = 13 5508 RETURN 5509 END IF 5510 ITERM = 0 5511 IF (N.GT.0.AND.NIT-KIT.GE.IRES1*N+IRES2) THEN 5512 IREST=MAX(IREST,1) 5513 END IF 5514 NIT = NIT + 1 5515 RETURN 5516 END 5517* SUBROUTINE PYFUT8 ALL SYSTEMS 98/12/01 5518* PURPOSE : 5519* TERMINATION CRITERIA AND TEST ON RESTART. 5520* 5521* PARAMETERS : 5522* II N ACTUAL NUMBER OF VARIABLES. 5523* RI F NEW VALUE OF THE OBJECTIVE FUNCTION. 5524* RI FO OLD VALUE OF THE OBJECTIVE FUNCTION. 5525* RO GMAX NORM OF THE TRANSFORMED GRADIENT. 5526* RI DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES. 5527* RI RPF3 VALUE OF THE BARRIER PARAMETER. 5528* RI TOLX LOWER BOUND FOR STEPLENGTH. 5529* RI TOLF LOWER BOUND FOR FUNCTION DECREASE. 5530* RI TOLB LOWER BOUND FOR FUNCTION VALUE. 5531* RI TOLG LOWER BOUND FOR GRADIENT. 5532* RI TOLP LOWER BOUND FOR BARRIER PARAMETER. 5533* II KD DEGREE OF REQUIRED DERIVATIVES. 5534* IU NIT ACTUAL NUMBER OF ITERATIONS. 5535* II KIT NUMBER OF THE ITERATION AFTER RESTART. 5536* II MIT MAXIMUM NUMBER OF ITERATIONS. 5537* IU NFV ACTUAL NUMBER OF COMPUTED FUNCTION VALUES. 5538* II MFV MAXIMUM NUMBER OF COMPUTED FUNCTION VALUES. 5539* IU NFG ACTUAL NUMBER OF COMPUTED GRADIENT VALUES. 5540* II MFG MAXIMUM NUMBER OF COMPUTED GRADIENT VALUES. 5541* IU NTESX ACTUAL NUMBER OF TESTS ON STEPLENGTH. 5542* II MTESX MAXIMUM NUMBER OF TESTS ON STEPLENGTH. 5543* IU NTESF ACTUAL NUMBER OF TESTS ON FUNCTION DECREASE. 5544* II MTESF MAXIMUM NUMBER OF TESTS ON FUNCTION DECREASE. 5545* II IRES1 RESTART SPECIFICATION. RESTART IS PERFORMED AFTER 5546* IRES1*N+IRES2 ITERATIONS. 5547* II IRES2 RESTART SPECIFICATION. RESTART IS PERFORMED AFTER 5548* IRES1*N+IRES2 ITERATIONS. 5549* IU IREST RESTART INDICATOR. RESTART IS PERFORMED IF IREST>0. 5550* II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. 5551* ITERS=0 FOR ZERO STEP. 5552* IO ITERM TERMINATION INDICATOR. ITERM=1-TERMINATION AFTER MTESX 5553* UNSUFFICIENT STEPLENGTHS. ITERM=2-TERMINATION AFTER MTESF 5554* UNSUFFICIENT FUNCTION DECREASES. ITERM=3-TERMINATION ON LOWER 5555* BOUND FOR FUNCTION VALUE. ITERM=4-TERMINATION ON LOWER BOUND 5556* FOR GRADIENT. ITERM=11-TERMINATION AFTER MAXIMUM NUMBER OF 5557* ITERATIONS. ITERM=12-TERMINATION AFTER MAXIMUM NUMBER OF 5558* COMPUTED FUNCTION VALUES. 5559* 5560 SUBROUTINE PYFUT8(N,F,FO,GMAX,DMAX,RPF3,TOLX,TOLF,TOLB,TOLG,TOLP, 5561 & KD,NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF,IRES1, 5562 & IRES2,IREST,ITERS,ITERM) 5563 INTEGER N,KD,NIT,KIT,MIT,NFV,MFV,NFG,MFG,NTESX,MTESX,NTESF,MTESF, 5564 & IRES1,IRES2,IREST,ITERS,ITERM 5565 DOUBLE PRECISION F,FO,RPF3,GMAX,DMAX,TOLX,TOLF,TOLG,TOLB,TOLP 5566 DOUBLE PRECISION TEMP 5567 IF (ITERM.LT.0) RETURN 5568 IF (ITERS.EQ.0) GO TO 1 5569 IF (NIT.LE.0) FO=F+MIN(SQRT(ABS(F)),ABS(F)/1.0D 1) 5570 IF (F.LE.TOLB) THEN 5571 ITERM = 3 5572 RETURN 5573 END IF 5574 IF (RPF3.GT.TOLP) GO TO 1 5575 IF (KD.GT.0) THEN 5576 IF (GMAX.LE.TOLG) THEN 5577 ITERM = 4 5578 RETURN 5579 END IF 5580 END IF 5581 IF (NIT.LE.0) THEN 5582 NTESX = 0 5583 NTESF = 0 5584 END IF 5585 IF (DMAX.LE.TOLX) THEN 5586 ITERM = 1 5587 NTESX = NTESX+1 5588 IF (NTESX.GE.MTESX) RETURN 5589 ELSE 5590 NTESX = 0 5591 END IF 5592 TEMP=ABS(FO-F)/MAX(ABS(F),1.0D 0) 5593 IF (TEMP.LE.TOLF) THEN 5594 ITERM = 2 5595 NTESF = NTESF+1 5596 IF (NTESF.GE.MTESF) RETURN 5597 ELSE 5598 NTESF = 0 5599 END IF 5600 1 IF (NIT.GE.MIT) THEN 5601 ITERM = 11 5602 RETURN 5603 END IF 5604 IF (NFV.GE.MFV) THEN 5605 ITERM = 12 5606 RETURN 5607 END IF 5608 IF (NFG.GE.MFG) THEN 5609 ITERM = 13 5610 RETURN 5611 END IF 5612 ITERM = 0 5613 IF (N.GT.0.AND.NIT-KIT.GE.IRES1*N+IRES2) THEN 5614 IREST=MAX(IREST,1) 5615 END IF 5616 NIT = NIT + 1 5617 RETURN 5618 END 5619* SUBROUTINE PYPTSH ALL SYSTEMS 98/12/01 5620* PURPOSE : 5621* POWELL-TOINT GRAPH COLORING ALGORITHM FOR GROUPING COLUMNS OF THE 5622* HESSIAN MATRIX BEFORE NUMERICAL DIFFERENTIATION. 5623* 5624* PARAMETERS : 5625* II NF DECLARED NUMBER OF VARIABLES. 5626* II MMAX MAXIMUM NUMBER OF NONZERO ELEMENTS. 5627* II IH(NF+1) POINTER VECTOR OF SPARSE HESSIAN MATRIX. 5628* II JH(MMAX) INDEX VECTOR OF THE HESSIAN MATRIX. 5629* IO COL(NF) VECTOR DISCERNING GROUPS OF THE HESSIAN COLUMN OF THE 5630* SAME COLOUR. 5631* RA DEG(NF) DEGREES OF THE ADJACENCY GRAPH. 5632* RA ORD(NF) AUXILIARY ARRAY. 5633* RA RADIX(NF+1) AUXILIARY ARRAY. 5634* IA WN11(NF) AUXILIARY VECTOR USED FOR INDICES OF THE COLUMNS 5635* THAT HAVE NOT BEEN COLOURED YET. 5636* IA WN12(NF) AUXILIARY VECTOR. 5637* RA XS(NF) AUXILIARY VECTOR. 5638* IO ITERM TERMINATION INDICATOR. 5639* 5640* SUBPROGRAMS USED : 5641* S PYCSER GROUPING COLUMNS OF THE SPARSE SYMMETRIC MATRIX. 5642* S MXSTG1 WIDTHEN THE STRUCTURE. 5643* S MXSTL1 SHRINK THE STRUCTURE. 5644* S MXVSR2 SORT. 5645* 5646 SUBROUTINE PYPTSH(NF,MMAX,IH,JH,COL,DEG,ORD,RADIX,WN11,WN12,XS, 5647 & ITERM) 5648 INTEGER NF,MMAX,IH(*),JH(*),COL(*) 5649 INTEGER WN11(*),WN12(*),ITERM 5650 DOUBLE PRECISION RADIX(*),ORD(*) 5651 DOUBLE PRECISION XS(*),DEG(*) 5652 INTEGER NCOL,CNM,I,ML,MM,J,K1,L 5653* 5654* SAVE SYMBOLIC STRUCTURE OF FACTOR 5655* 5656 MM=IH(NF+1)-1 5657 IF (2*MM-NF+2.GE.MMAX) THEN 5658 ITERM=-45 5659 RETURN 5660 END IF 5661* 5662* WIDTHEN THE STRUCTURE 5663* 5664 CALL MXSTG1(NF,ML,IH,JH,WN12,WN11) 5665 DO 100 I=1,NF 5666 COL(I)=NF 5667 WN12(I)=0 5668 WN11(I)=I 5669100 CONTINUE 5670* 5671* NUMBER OF THE FREE COLUMNS 5672* 5673 CNM=NF 5674* 5675* NUMBER OF USED COLOURS 5676* 5677 NCOL=1 5678* 5679* DEGREE RECOUNT 5680* 5681 K1=1 5682 DO 110 I=1,NF 5683 L=IH(I+1) 5684 DEG(I)=L-K1 5685 K1=L 5686110 CONTINUE 5687* 5688* COLUMN RESORT 5689* 5690200 CALL MXVSR2(NF,DEG,ORD,RADIX,WN11,CNM) 5691* 5692* ORD REWRITE INTO THE ARRAY INVP 5693* 5694 DO 250 I=1,CNM 5695 WN11(I)=ORD(I) 5696250 CONTINUE 5697* 5698* COLUMNS OF THE NEW COLOUR NCOL 5699* 5700 CALL PYCSER(JH,IH,WN12,XS,DEG,WN11,COL,NCOL,CNM) 5701* 5702* STOP TEST 5703* 5704 IF (CNM.GE.1) THEN 5705 NCOL=NCOL+1 5706 GO TO 200 5707 END IF 5708* 5709* SHRINK THE STRUCTURE 5710* 5711 CALL MXSTL1(NF,ML,IH,JH,WN12) 5712* 5713* INTO COL GIVE INDICES OF THE INDIVIDUAL GROUPS ONE AFTER ANOTHER, 5714* END OF THE GROUP IS MARKED BY THE NEGATIVE INDEX VALUE. 5715* 5716* 5717* READ COL 5718* 5719 DO 300 I=1,NF 5720 WN11(I)=0 5721 300 CONTINUE 5722 DO 400 I=1,NF 5723 J=COL(I) 5724 WN11(J)=WN11(J)+1 5725 400 CONTINUE 5726 WN12(1)=1 5727 L=1 5728 DO 500 I=2,NF 5729 L=L+WN11(I-1) 5730 WN12(I)=L 5731 IF (WN11(I).EQ.0) GO TO 550 5732 500 CONTINUE 5733 550 CONTINUE 5734* 5735* CHANGE COL 5736* 5737 DO 600 I=1,NF 5738 J=COL(I) 5739 WN11(I)=J 5740 600 CONTINUE 5741 DO 700 I=1,NF 5742 J=WN11(I) 5743 COL(WN12(J))=I 5744 WN12(J)=WN12(J)+1 5745 700 CONTINUE 5746 DO 800 I=1,NCOL 5747 L=WN12(I)-1 5748 IF (L.GT.NF) GO TO 900 5749 COL(L)=-COL(L) 5750 800 CONTINUE 5751 900 CONTINUE 5752 RETURN 5753 END 5754* SUBROUTINE PYRMC0 ALL SYSTEMS 98/12/01 5755* PURPOSE : 5756* OLD SIMPLE BOUND IS REMOVED FROM THE ACTIVE SET. TRANSFORMED 5757* GRADIENT OF THE OBJECTIVE FUNCTION IS UPDATED. 5758* 5759* PARAMETERS : 5760* II NF DECLARED NUMBER OF VARIABLES. 5761* II N REDUCED NUMBER OF VARIABLES. 5762* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. 5763* RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 5764* RI EPS8 TOLERANCE FOR CONSTRAINT TO BE REMOVED. 5765* RI UMAX MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER. 5766* RI GMAX NORM OF THE TRANSFORMED GRADIENT. 5767* RO RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER. 5768* II IOLD NUMBER OF REMOVED CONSTRAINTS. 5769* IU IREST RESTART INDICATOR. 5770* 5771 SUBROUTINE PYRMC0(NF,N,IX,G,EPS8,UMAX,GMAX,RMAX,IOLD,IREST) 5772 INTEGER NF,N,IX(*),IOLD,IREST 5773 DOUBLE PRECISION G(*),EPS8,UMAX,GMAX,RMAX 5774 INTEGER I,IXI 5775 IF (N.EQ.0.OR.RMAX.GT.0.0D 0) THEN 5776 IF (UMAX.GT.EPS8*GMAX) THEN 5777 IOLD=0 5778 DO 1 I=1,NF 5779 IXI=IX(I) 5780 IF (IXI.GE.0) THEN 5781 ELSE IF (IXI.LE.-5) THEN 5782 ELSE IF ((IXI.EQ.-1.OR.IXI.EQ.-3).AND.-G(I).LE.0.0D 0) THEN 5783 ELSE IF ((IXI.EQ.-2.OR.IXI.EQ.-4).AND. G(I).LE.0.0D 0) THEN 5784 ELSE 5785 IOLD=IOLD+1 5786 IX(I)=MIN(ABS(IX(I)),3) 5787 IF (RMAX.EQ.0) GO TO 2 5788 END IF 5789 1 CONTINUE 5790 2 IF (IOLD.GT.1) IREST=MAX(IREST,1) 5791 END IF 5792 END IF 5793 RETURN 5794 END 5795* SUBROUTINE PYTCAB ALL SYSTEMS 06/12/01 5796* PURPOSE : 5797* VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED 5798* AND SCALED. TEST VALUE DMAX IS DETERMINED. 5799* 5800* PARAMETERS : 5801* II NC NUMBER OF APPROXIMATED FUNCTIONS. 5802* II MC NUMBER OF NONZERO ELEMENTS IN THE FIELD CG. 5803* RI CG(MC) JACOBIAN MATRIX OF THE APPROXIMATED FUNCTIONS. 5804* RO CGO(MC) SAVED JACOBIAN MATRIX OF THE APPROXIMATED FUNCTIONS. 5805* RI ICG(NC+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD CG. 5806* RI CZ(NC) VECTOR CONTAINING LAGRANGE MULTIPLIERS FOR CONSTRAINTS. 5807* II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. 5808* ITERS=0 FOR ZERO STEP. 5809* II JOB SUBJECTS OF UPDATES. JOB=0-CONSTRAINT FUNCTIONS. 5810* JOB=1-CONSTRAINT FUNCTIONS MULTIPLIED BY SIGNS OF THE 5811* LAGRANGIAN MULTIPLIERS. JOB-2-TERMS OF THE LAGRANGIAN 5812* FUNCTION. 5813* 5814* SUBPROGRAMS USED : 5815* S MXVDIF DIFFERENCE OF TWO VECTORS. 5816* S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE 5817* SUBSTRACTED ONE. 5818* 5819 SUBROUTINE PYTCAB(NC,MC,CG,CGO,ICG,CZ,ITERS,JOB) 5820 INTEGER NC,MC,ICG(*),ITERS,JOB 5821 DOUBLE PRECISION CG(*),CGO(*),CZ(*) 5822 INTEGER J,K,KC,L,M 5823 DOUBLE PRECISION TEMP 5824 IF (ITERS.GT.0) THEN 5825 CALL MXVDIF(MC,CG,CGO,CGO) 5826 ELSE 5827 CALL MXVSAV(MC,CG,CGO) 5828 END IF 5829 DO 4 KC=1,NC 5830 M=ICG(KC) 5831 L=ICG(KC+1)-M 5832 IF (JOB.GT.0) THEN 5833 TEMP=CZ(KC) 5834 IF (JOB.EQ.1) TEMP=SIGN(1.0D 0,TEMP) 5835 K=M 5836 DO 2 J=1,L 5837 CGO(K)=CGO(K)*TEMP 5838 K=K+1 5839 2 CONTINUE 5840 END IF 5841 4 CONTINUE 5842 RETURN 5843 END 5844* SUBROUTINE PYTCUB ALL SYSTEMS 06/12/01 5845* PURPOSE : 5846* VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED 5847* AND SCALED. TEST VALUE DMAX IS DETERMINED. 5848* 5849* PARAMETERS : 5850* II NC NUMBER OF APPROXIMATED FUNCTIONS. 5851* II MC NUMBER OF NONZERO ELEMENTS IN THE FIELD CG. 5852* RI CG(MC) JACOBIAN MATRIX OF THE APPROXIMATED FUNCTIONS. 5853* RO CGO(MC) SAVED JACOBIAN MATRIX OF THE APPROXIMATED FUNCTIONS. 5854* RI ICG(NC+1) POSITION OF THE FIRST ROWS ELEMENTS IN THE FIELD CG. 5855* II IC(NC) VECTOR CONTAINING TYPES OF CONSTRAINTS. 5856* RI CZL(NC) VECTOR CONTAINING LOWER MULTIPLIERS FOR CONSTRAINTS. 5857* RI CZU(NC) VECTOR CONTAINING UPPER MULTIPLIERS FOR CONSTRAINTS. 5858* II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. 5859* ITERS=0 FOR ZERO STEP. 5860* II JOB SUBJECTS OF UPDATES. JOB=0-CONSTRAINT FUNCTIONS. 5861* JOB=1-CONSTRAINT FUNCTIONS MULTIPLIED BY SIGNS OF THE 5862* LAGRANGIAN MULTIPLIERS. JOB-2-TERMS OF THE LAGRANGIAN 5863* FUNCTION. 5864* 5865* SUBPROGRAMS USED : 5866* S MXVDIF DIFFERENCE OF TWO VECTORS. 5867* S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE 5868* SUBSTRACTED ONE. 5869* 5870 SUBROUTINE PYTCUB(NC,MC,CG,CGO,ICG,IC,CZL,CZU,ITERS,JOB) 5871 INTEGER NC,MC,ICG(NC+1),IC(NC),ITERS,JOB 5872 DOUBLE PRECISION CG(*),CGO(*),CZL(*),CZU(*) 5873 INTEGER J,K,KC,KK,L,M 5874 DOUBLE PRECISION TEMP 5875 IF (ITERS.GT.0) THEN 5876 CALL MXVDIF(MC,CG,CGO,CGO) 5877 ELSE 5878 CALL MXVSAV(MC,CG,CGO) 5879 END IF 5880 DO 4 KC=1,NC 5881 M=ICG(KC) 5882 L=ICG(KC+1)-M 5883 IF (JOB.GT.0) THEN 5884 KK=ABS(IC(KC)) 5885 IF (KK.EQ.3.OR.KK.EQ.4) THEN 5886 TEMP= CZU(KC)-CZL(KC) 5887 ELSE IF (KK.EQ.1) THEN 5888 TEMP=-CZL(KC) 5889 ELSE IF (KK.EQ.2) THEN 5890 TEMP= CZU(KC) 5891 ELSE IF (KK.EQ.5) THEN 5892 TEMP= CZL(KC) 5893 END IF 5894 IF (JOB.EQ.1) TEMP=SIGN(1.0D 0,TEMP) 5895 K=M 5896 DO 2 J=1,L 5897 CGO(K)=CGO(K)*TEMP 5898 K=K+1 5899 2 CONTINUE 5900 END IF 5901 4 CONTINUE 5902 RETURN 5903 END 5904* SUBROUTINE PYTRCD ALL SYSTEMS 98/12/01 5905* PURPOSE : 5906* VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED 5907* AND SCALED AND REDUCED. TEST VALUE DMAX IS DETERMINED. 5908* 5909* PARAMETERS : 5910* II NF DECLARED NUMBER OF VARIABLES. 5911* RI X(NF) VECTOR OF VARIABLES. 5912* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. 5913* RU XO(NF) VECTORS OF VARIABLES DIFFERENCE. 5914* RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 5915* RU GO(NF) GRADIENTS DIFFERENCE. 5916* RO R VALUE OF THE STEPSIZE PARAMETER. 5917* RO F NEW VALUE OF THE OBJECTIVE FUNCTION. 5918* RI FO OLD VALUE OF THE OBJECTIVE FUNCTION. 5919* RO P NEW VALUE OF THE DIRECTIONAL DERIVATIVE. 5920* RI PO OLD VALUE OF THE DIRECTIONAL DERIVATIVE. 5921* RO DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES. 5922* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. 5923* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. 5924* IO KD DEGREE OF REQUIRED DERIVATIVES. 5925* IO LD DEGREE OF COMPUTED DERIVATIVES. 5926* II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION. 5927* ITERS=0 FOR ZERO STEP. 5928* 5929* SUBPROGRAMS USED : 5930* S MXVDIF DIFFERENCE OF TWO VECTORS. 5931* S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE 5932* SUBSTRACTED ONE. 5933* 5934 SUBROUTINE PYTRCD(NF,X,IX,XO,G,GO,R,F,FO,P,PO,DMAX,KBF,KD,LD, 5935 & ITERS) 5936 INTEGER NF,IX(*),KBF,KD,LD,ITERS 5937 DOUBLE PRECISION X(*),XO(*),G(*),GO(*),R,F,FO,P,PO,DMAX 5938 INTEGER I 5939 IF (ITERS.GT.0) THEN 5940 CALL MXVDIF(NF,X,XO,XO) 5941 CALL MXVDIF(NF,G,GO,GO) 5942 PO=R*PO 5943 P=R*P 5944 ELSE 5945 F = FO 5946 P = PO 5947 CALL MXVSAV(NF,X,XO) 5948 CALL MXVSAV(NF,G,GO) 5949 LD=KD 5950 END IF 5951 DMAX = 0.0D 0 5952 DO 1 I=1,NF 5953 IF (KBF.GT.0) THEN 5954 IF (IX(I).LT.0) THEN 5955 XO(I)=0.0D 0 5956 GO(I)=0.0D 0 5957 GO TO 1 5958 END IF 5959 END IF 5960 DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0)) 5961 1 CONTINUE 5962 RETURN 5963 END 5964* SUBROUTINE PYTRCG ALL SYSTEMS 99/12/01 5965* PURPOSE : 5966* GRADIENT OF THE OBJECTIVE FUNCTION IS SCALED AND REDUCED. TEST VALUES 5967* GMAX AND UMAX ARE COMPUTED. 5968* 5969* PARAMETERS : 5970* II NF DECLARED NUMBER OF VARIABLES. 5971* II N ACTUAL NUMBER OF VARIABLES. 5972* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. 5973* RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 5974* RI UMAX MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER. 5975* RI GMAX NORM OF THE TRANSFORMED GRADIENT. 5976* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. 5977* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. 5978* II IOLD INDEX OF THE REMOVED CONSTRAINT. 5979* 5980* SUBPROGRAMS USED : 5981* RF MXVMAX L-INFINITY NORM OF A VECTOR. 5982* 5983 SUBROUTINE PYTRCG(NF,N,IX,G,UMAX,GMAX,KBF,IOLD) 5984 INTEGER NF,N,IX(*),KBF,IOLD 5985 DOUBLE PRECISION G(*),UMAX,GMAX 5986 DOUBLE PRECISION TEMP,MXVMAX 5987 INTEGER I 5988 IF (KBF.GT.0) THEN 5989 GMAX = 0.0D 0 5990 UMAX = 0.0D 0 5991 IOLD=0 5992 DO 1 I=1,NF 5993 TEMP=G(I) 5994 IF ( IX(I) .GE. 0) THEN 5995 GMAX=MAX(GMAX,ABS(TEMP)) 5996 ELSE IF (IX(I).LE.-5) THEN 5997 ELSE IF (( IX(I) .EQ. -1 .OR. IX(I) .EQ. -3) 5998 & .AND. UMAX+TEMP .GE. 0.0D 0) THEN 5999 ELSE IF (( IX(I) .EQ. -2 .OR. IX(I) .EQ. -4) 6000 & .AND. UMAX-TEMP .GE. 0.0D 0) THEN 6001 ELSE 6002 IOLD=I 6003 UMAX=ABS(TEMP) 6004 END IF 6005 1 CONTINUE 6006 ELSE 6007 UMAX=0.0D 0 6008 GMAX=MXVMAX(NF,G) 6009 END IF 6010 N=NF 6011 RETURN 6012 END 6013* SUBROUTINE PYTRCS ALL SYSTEMS 98/12/01 6014* PURPOSE : 6015* SCALED AND REDUCED DIRECTION VECTOR IS BACK TRANSFORMED. VECTORS 6016* X,G AND VALUES F,P ARE SAVED. 6017* 6018* PARAMETERS : 6019* II NF DECLARED NUMBER OF VARIABLES. 6020* RI X(NF) VECTOR OF VARIABLES. 6021* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. 6022* RO XO(NF) SAVED VECTOR OF VARIABLES. 6023* RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES. 6024* RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES. 6025* RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION. 6026* RO GO(NF) SAVED GRADIENT OF THE OBJECTIVE FUNCTION. 6027* RO S(NF) DIRECTION VECTOR. 6028* RO RO SAVED VALUE OF THE STEPSIZE PARAMETER. 6029* RO FP PREVIOUS VALUE OF THE OBJECTIVE FUNCTION. 6030* RU FO SAVED VALUE OF THE OBJECTIVE FUNCTION. 6031* RI F VALUE OF THE OBJECTIVE FUNCTION. 6032* RO PO SAVED VALUE OF THE DIRECTIONAL DERIVATIVE. 6033* RI P VALUE OF THE DIRECTIONAL DERIVATIVE. 6034* RO RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER. 6035* RI ETA9 MAXIMUM FOR REAL NUMBERS. 6036* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. 6037* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. 6038* 6039* SUBPROGRAMS USED : 6040* S MXVCOP COPYING OF A VECTOR. 6041* 6042 SUBROUTINE PYTRCS(NF,X,IX,XO,XL,XU,G,GO,S,RO,FP,FO,F,PO,P,RMAX, 6043 & ETA9,KBF) 6044 INTEGER NF,IX(*),KBF 6045 DOUBLE PRECISION X(*),XO(*),XL(*),XU(*),G(*),GO(*),S(*),RO,FP,FO, 6046 & F,PO,P,RMAX,ETA9 6047 INTEGER I 6048 FP = FO 6049 RO = 0.0D 0 6050 FO = F 6051 PO = P 6052 CALL MXVCOP(NF,X,XO) 6053 CALL MXVCOP(NF,G,GO) 6054 IF (KBF.GT.0) THEN 6055 DO 1 I=1,NF 6056 IF (IX(I).LT.0) THEN 6057 S(I)=0.0D 0 6058 ELSE 6059 IF (IX(I).EQ.1.OR.IX(I).GE.3) THEN 6060 IF (S(I).LT.-1.0D 0/ETA9) RMAX=MIN(RMAX,(XL(I)-X(I))/S(I)) 6061 END IF 6062 IF (IX(I).EQ.2.OR.IX(I).GE.3) THEN 6063 IF (S(I).GT. 1.0D 0/ETA9) RMAX=MIN(RMAX,(XU(I)-X(I))/S(I)) 6064 END IF 6065 END IF 6066 1 CONTINUE 6067 END IF 6068 RETURN 6069 END 6070* SUBROUTINE PYTSCH ALL SYSTEMS 99/12/01 6071* PURPOSE : 6072* HESSIAN MATRIX OF THE OBJECTIVE FUNCTION OR ITS APPROXIMATION 6073* IS SCALED. 6074* 6075* PARAMETERS : 6076* II NF DECLARED NUMBER OF VARIABLES. 6077* II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS. 6078* RU H(M) HESSIAN MATRIX OR ITS APPROXIMATION. 6079* II IH(N+1) POINTERS OF THE DIAGONAL ELEMENTS OF H. 6080* II JH(M) INDICES OF THE NONZERO ELEMENTS OF H. 6081* II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS. 6082* KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS. 6083* 6084 SUBROUTINE PYTSCH(NF,IX,H,IH,JH,KBF) 6085 INTEGER NF,IX(*),IH(*),JH(*),KBF 6086 DOUBLE PRECISION H(*) 6087 INTEGER I,J,K,JSTRT,JSTOP 6088 IF (KBF.GT.0) THEN 6089 JSTOP=0 6090 DO 3 I=1,NF 6091 JSTRT=JSTOP+1 6092 JSTOP=IH(I+1)-1 6093 IF (IX(I).GE.0) THEN 6094 DO 1 J=JSTRT,JSTOP 6095 K=JH(J) 6096 IF (K.LT.0) THEN 6097 H(J)=0.0D 0 6098 END IF 6099 1 CONTINUE 6100 ELSE 6101 H(JSTRT)=1.0D 0 6102 DO 2 J=JSTRT+1,JSTOP 6103 H(J)=0.0D 0 6104 2 CONTINUE 6105 END IF 6106 3 CONTINUE 6107 END IF 6108 RETURN 6109 END 6110