1 * SUBROUTINE PA0GS3 ALL SYSTEMS 91/12/01
3 * NUMERICAL COMPUTATION OF THE GRADIENT OF THE APPROXIMATED
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.
18 * SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION.
20 SUBROUTINE PA0GS3(N,KA,X,FA,GA,IAG,JAG,ETA1,NAV)
21 DOUBLE PRECISION ETA1,FA
23 DOUBLE PRECISION GA(*),X(*)
25 DOUBLE PRECISION ETA,FTEMP,XSTEP,XTEMP
29 DO 10 KVAR = IAG(KA),IAG(KA+1) - 1
34 XSTEP = ETA*MAX(ABS(X(IVAR)),1.0D0)*SIGN(1.0D0,X(IVAR))
36 X(IVAR) = X(IVAR) + XSTEP
37 XSTEP = X(IVAR) - XTEMP
41 * NUMERICAL DIFFERENTIATION
43 GA(IVAR) = (FA-FTEMP)/XSTEP
49 * SUBROUTINE PA0HS3 ALL SYSTEMS 99/12/01
51 * NUMERICAL COMPUTATION OF THE HESSIAN MATRIX OF THE APPROXIMATED
52 * FUNCTION USING ITS VALUES.
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.
71 * SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION.
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
78 INTEGER IVAR,JVAR,KVAR,LVAR,MVAR
79 ETA=ETA1**(1.0D 0/3.0D 0)
82 DO 4 KVAR=MVAR+1,IAG(KA+1)-1
85 IF (IX(IVAR).LE.-5) GO TO 4
91 IF (XTEMPI.GE.0.0D 0) THEN
92 GO(IVAR)= ETA*MAX(ABS(XTEMPI),1.0D 0)
94 GO(IVAR)=-ETA*MAX(ABS(XTEMPI),1.0D 0)
96 X(IVAR)=X(IVAR)+GO(IVAR)
97 GO(IVAR)=X(IVAR)-XTEMPI
104 * NUMERICAL DIFFERENTIATION
106 DO 10 KVAR=MVAR+1,IAG(KA+1)-1
109 IF (IX(IVAR).LE.-5) GO TO 10
112 X(IVAR)=XTEMPI+GO(IVAR)
113 DO 9 LVAR=KVAR,IAG(KA+1)-1
116 IF (IX(JVAR).LE.-5) GO TO 9
119 X(JVAR)=X(JVAR)+GO(JVAR)
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))
133 * SUBROUTINE PA0SQ3 ALL SYSTEMS 92/12/01
135 * COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE FUNCTION
136 * WHICH IS DEFINED AS A SUM OF SQUARES.
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).
157 * SE FUN COMPUTATION OF THE VALUE OF THE APPROXIMATED FUNCTION.
158 * S PA0GS3 NUMERICAL DIFFERENTIATION.
159 * S MXVSET INITIATION OF A VECTOR.
161 SUBROUTINE PA0SQ3(N,X,F,AF,GA,AG,IAG,JAG,G,ETA1,KD,LD,NFV,NFG,
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(*)
168 INTEGER J,JP,K,KA,L,NAV
170 IF (KD.GE.0 .AND. LD.LT.0) THEN
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
180 IF (KD.LT.0) GO TO 30
187 IF (LD.GE.0) GO TO 10
189 10 IF (KD.LT.1) GO TO 30
191 CALL PA0GS3(N,KA,X,FA,GA,IAG,JAG,ETA1,NAV)
199 G(JP) = G(JP) + FA*GA(JP)
204 IF (KD.GE.0 .AND. LD.LT.0) F = 0.5D0*F
205 IF (IDER.EQ.0) NFV=NFV+NAV/N
209 * SUBROUTINE PA1HS3 ALL SYSTEMS 99/12/01
211 * NUMERICAL COMPUTATION OF THE HESSIAN MATRIX OF THE APPROXIMATED
212 * FUNCTION USING ITS GRADIENTS.
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.
231 * SE DFUN COMPUTATION OF THE GRADIENT OF THE APPROXIMATED FUNCTION.
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
238 INTEGER IVAR,JVAR,KVAR,LVAR,MVAR
242 DO 5 KVAR=MVAR+1,IAG(KA+1)-1
245 IF (IX(IVAR).LE.-5) GO TO 5
251 IF (XTEMP.GE.0.0D 0) THEN
252 XSTEP= ETA*MAX(ABS(XTEMP),1.0D 0)
254 XSTEP=-ETA*MAX(ABS(XTEMP),1.0D 0)
258 CALL DFUN(NF,KA,X,GA)
261 * NUMERICAL DIFFERENTIATION
263 DO 4 LVAR=MVAR+1,IAG(KA+1)-1
266 IF (IX(JVAR).LE.-5) GO TO 4
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
274 HA(IJ)=0.5D 0*(HA(IJ)+(GA(JVAR)-GO(JVAR))/XSTEP)
282 * SUBROUTINE PA1SF3 ALL SYSTEMS 97/12/01
284 * COMPUTATION OF THE VALUE AND THE GRADIENT OF THE OBJECTIVE FUNCTION
285 * WHICH IS DEFINED AS A SUM OF SQUARES.
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
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
304 * IU NFV NUMBER OF OBJECTIVE FUNCTION VALUES COMPUTED.
305 * IU NFG NUMBER OF OBJECTIVE FUNCTION GRADIENTS COMPUTED.
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.
312 SUBROUTINE PA1SF3(NF,NA,X,GA,G,AG,IAG,JAG,F,AF,KD,LD,ISNA,
314 INTEGER NF,NA,IAG(*),JAG(*),KD,LD,ISNA,NFV,NFG
315 DOUBLE PRECISION X(*),GA(*),G(*),AG(*),F,AF(*)
319 IF (KD.GE.0.AND.LD.LT.0) THEN
323 IF (KD.GE.1.AND.LD.LT.1) THEN
324 CALL MXVSET(NF,0.0D 0,G)
338 CALL DFUN(NF,KA,X,GA)
344 IF (ISNA.GT.1) AG(K)=GA(JP)
352 * SUBROUTINE PA2SF4 ALL SYSTEMS 97/12/01
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.
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
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.
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
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,
395 DOUBLE PRECISION X(*),GA(*),G(*),GO(*),HA(*),H(*),AF(*),F,ETA1
397 INTEGER J,JP,K,KA,L,NAG
399 IF (KD.GE.0.AND.LD.LT.0) THEN
403 IF (KD.GE.1.AND.LD.LT.1) THEN
404 CALL MXVSET(NF,0.0D 0,G)
407 IF (KD.GE.2.AND.LD.LT.2) CALL MXVSET(IH(NF+1)-1,0.0D 0,H)
419 CALL DFUN(NF,KA,X,GA)
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)
438 * SUBROUTINE PA2SQ4 ALL SYSTEMS 97/12/01
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.
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
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
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.
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
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,
480 DOUBLE PRECISION X(*),GA(*),AG(*),G(*),H(*),AF(*),F,ETA1
481 INTEGER J,JP,K,KA,L,NAV
484 IF (KD.GE.0.AND.LD.LT.0) THEN
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
492 IF (KD.GE.2.AND.LD.LT.2) CALL MXVSET(IH(NF+1)-1,0.0D 0,H)
505 CALL PA0GS3(NF,KA,X,FA,GA,IAG,JAG,ETA1,NAV)
507 CALL DFUN(NF,KA,X,GA)
514 G(JP)=G(JP)+FA*GA(JP)
515 IF (ISNA.GT.1) AG(K)=GA(JP)
518 2 IF (KD.LT.2) GO TO 3
520 CALL PASSH1(H,IH,JH,IAG,JAG,GA,KA,1.0D 0)
522 IF (KD.GE.0.AND.LD.LT.0) F=5.0D-1*F
523 IF (IDER.EQ.0) NFV=NFV+NAV/NA
527 * SUBROUTINE PA2SQ8 ALL SYSTEMS 97/12/01
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.
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
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.
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
570 * S PASSH2 ADDITION OF THE PARTIAL HESSIAN MATRIX TO THE SPARSE
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,
577 DOUBLE PRECISION X(*),GA(*),G(*),GO(*),GS(*),HA(*),H(*),AF(*),F,
579 INTEGER J,JP,K,KA,L,NAV,NAG
582 IF (KD.GE.0.AND.LD.LT.0) THEN
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
590 IF (KD.GE.2.AND.LD.LT.2) CALL MXVSET(IH(NF+1)-1,0.0D 0,H)
604 CALL PA0GS3(NF,KA,X,FA,GA,IAG,JAG,ETA1,NAV)
606 CALL DFUN(NF,KA,X,GA)
613 G(JP)=G(JP)+FA*GA(JP)
621 CALL PA0HS3(NF,KA,X,IX,HA,GO,GS,IAG,JAG,FA,ETA1,KBF,NAV)
623 CALL PA1HS3(NF,KA,X,IX,HA,GO,GA,IAG,JAG,FA,ETA1,KBF,NAG)
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)
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
635 * SUBROUTINE PALNG3 ALL SYSTEMS 97/12/01
637 * COMPUTATION OF THE GRADIENT OF THE LINEAR APPROXIMATED FUNCTION.
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.
646 SUBROUTINE PALNG3(AG,IAG,JAG,GA,KA)
647 DOUBLE PRECISION AG(*),GA(*)
648 INTEGER IAG(*),JAG(*),KA
659 * SUBROUTINE PASED3 ALL SYSTEMS 07/12/01
661 * COMPRESSED SPARSE STRUCTURE OF THE JACOBIAN MATRIX IS COMPUTED FROM
662 * THE COORDINATE FORM.
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.
673 SUBROUTINE PASED3(NA,MA,IAG,JAG,IER)
674 INTEGER NA,MA,IAG(*),JAG(*),IER
677 CALL MXVSR7(MA,IAG,JAG)
678 IF (IAG(1).LT.1.OR.IAG(MA).GT.NA) THEN
682 CALL MXVINS(NA,0,IAG(MA+1))
684 IAG(IAG(J)+MA)=IAG(IAG(J)+MA)+1
688 IAG(KA+1)=IAG(KA)+IAG(KA+MA)
695 CALL MXVSRT(L,JAG(K))
696 IF (JAG(K).LT.1.OR.JAG(K+L-1).GT.NA) THEN
703 IF (J.GT.1.AND.JAG(K).EQ.JAG(K-1)) THEN
711 IAG(NA+1)=IAG(NA+1)-I
715 * SUBROUTINE PASSH1 ALL SYSTEMS 98/12/01
717 * COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE SPARSE JACOBIAN
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
728 * RI GA(NF) GRADIENT OF THE SELECTED FUNCTION.
729 * II KA INDEX OF THE SELECTED FUNCTION (ROW OF THE SPARSE JACOBIAN
731 * RI FACTOR SCALING FACTOR.
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
745 2 IF (ABS(JH(JF)).LT.J) THEN
749 H(JF)=H(JF)+TEMP*GA(J)
754 * SUBROUTINE PASSH2 ALL SYSTEMS 98/12/01
756 * COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE SPARSE JACOBIAN
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
768 * II KA INDEX OF THE SELECTED FUNCTION (ROW OF THE SPARSE JACOBIAN
770 * RI FACTOR SCALING FACTOR.
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
787 2 IF (ABS(JH(JF)).LT.J) THEN
791 H(JF)=H(JF)+FACTOR*HA(K)
799 * SUBROUTINE PASSH3 ALL SYSTEMS 98/12/01
801 * COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE SPARSE JACOBIAN
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
812 * RI GA(NF) GRADIENT OF THE SELECTED FUNCTION.
813 * II KA INDEX OF THE SELECTED FUNCTION (ROW OF THE SPARSE JACOBIAN
815 * RI FACTOR SCALING FACTOR.
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
831 2 IF (ABS(JH(JF)).LT.J) THEN
835 H(JF)=H(JF)+TEMP*GA(J)
840 * SUBROUTINE PCBS04 ALL SYSTEMS 98/12/01
842 * INITIATION OF THE VECTOR CONTAINING TYPES OF CONSTRAINTS.
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.
854 SUBROUTINE PCBS04(NF,X,IX,XL,XU,EPS9,KBF)
856 DOUBLE PRECISION X(*),XL(*),XU(*),EPS9
857 DOUBLE PRECISION TEMP
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)
871 * SUBROUTINE PDSGM1 ALL SYSTEMS 01/09/22
873 * COMPUTATION OF A TRUST-REGION STEP BY THE DOG-LEG METHOD WITH DIRECT
874 * MATRIX DECOMPOSITIONS.
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
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.
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.
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.
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
973 DOUBLE PRECISION B1,B2,B3,D3,S1,S2
974 DOUBLE PRECISION MXSSMQ,MXSPCQ,MXUDOT
977 * DIRECTION DETERMINATION
979 IF (IDEC.LT.0) IDEC=0
981 ELSE IF (IDEC.EQ.1) THEN
987 B2=MXUDOT(NF,G,G,IX,KBF)
990 IF (ALF2*GNORM.LE.XDEL) THEN
993 CALL MXSPCT(NF,MM,MH,MMAX,H,JH,PSL,ITERM)
994 IF (ITERM.NE.0) GO TO 13130
996 * SPARSE GILL-MURRAY DECOMPOSITION
999 CALL MXSPCF(NF,H(MM+1),PSL,JH(MM+1),WN11,WN12,XS,INF,S1,S2)
1004 CALL MXSPCD(NF,H(MM+1),PSL,JH(MM+1),S,INF)
1005 CALL MXVSBP(NF,PERM,S,XS)
1007 * DIRECTION OF NEGATIVE CURVATURE
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)
1016 ELSE IF (GNORM.LE.0.0D 0) THEN
1021 CALL MXVSET(NF,0.0D 0,S)
1026 B1=MXSSMQ(NF,H,IH,JH,G,G)
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)
1033 IF (XDEL.LE.0.0D 0) THEN
1035 * INITIAL TRUST REGION BOUND
1037 IF (B1.LE.0.0D 0) THEN
1042 IF (IEST.EQ.1) XDEL=MIN(XDEL,4.0D 0*(F-FMIN)/GNORM)
1045 IF (B1.LE.0.0D 0.OR.B2*GNORM.GE.B1*XDEL) THEN
1047 * SCALED STEEPEST DESCENT DIRECTION IS ACCEPTED
1049 CALL MXVSCL(NF,-XDEL/GNORM,G,S)
1055 CALL MXSPCT(NF,MM,MH,MMAX,H,JH,PSL,ITERM)
1056 IF (ITERM.NE.0) THEN
1060 * SPARSE GILL-MURRAY DECOMPOSITION
1063 CALL MXSPCF(NF,H(MM+1),PSL,JH(MM+1),WN11,WN12,XS,INF,S1,S2)
1068 * COMPUTATION OF THE NEWTON DIRECTION
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))
1076 * COMPUTATION OF THE STEEPEST DESCENT DIRECTION
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
1087 * NEWTON AND THE STEEPEST DESCENT DIRECTION ARE
1088 * APPROXIMATELY EQUAL
1090 CALL MXUCOP(NF,GO,S,IX,KBF)
1093 ELSE IF (B1.LE.0.0D 0) THEN
1095 * BOUNDARY STEP WITH NEGATIVE INCREMENT
1097 CALL PNSTEP(XDEL,SNORM,-B1,B2,B3)
1098 CALL MXUDIR(NF,-B3,XO,S,S,IX,KBF)
1101 ELSE IF (D3.LE.XDEL) THEN
1103 * NEWTON DIRECTION IS ACCEPTED
1105 CALL MXUCOP(NF,GO,S,IX,KBF)
1110 * DOUBLE DOGLEG STRATEGY
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)
1125 PP=MXSSMQ(NF,H,IH,JH,S,S)*0.5D 0
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
1134 IF (KD.GT.0) P=MXUDOT(NF,G,S,IX,KBF)
1137 * SUBROUTINE PDSGM4 ALL SYSTEMS 01/09/22
1139 * COMPUTATION OF A TRUST-REGION STEP BY THE SHIFTED STEIHAUG-TOINT
1140 * METHOD WITH CONJUGATE GRADIENT ITERATIONS.
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
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
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.
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
1220 * S MXSSMD MATRIX-VECTOR PRODUCT FOLLOWED BY THE ADDITION OF A
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.
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.
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
1257 * DIRECTION DETERMINATION
1263 IF (IDEC.LT.0) IDEC=0
1264 IF (IDEC.NE.0.AND.IDEC.NE.1) THEN
1268 GNORM=SQRT(MXUDOT(NF,G,G,IX,KBF))
1269 IF (GNORM.GE.1.0D 3*GNORMO) EPS=1.0D-6
1271 RHO1=MXUDOT(NF,G,G,IX,KBF)
1272 IF (XDEL.LE.0.0D 0) THEN
1274 * INITIAL TRUST REGION BOUND
1276 RHO2=MXSSMQ(NF,H,IH,JH,G,G)
1277 IF (RHO2.LE.0.0D 0) THEN
1280 XDEL=(GNORM*GNORM/RHO2)*GNORM
1282 IF (IEST.EQ.1) XDEL=MIN(XDEL,4.0D 0*(F-FMIN)/GNORM)
1285 PAR=MIN(EPS,SQRT(GNORM))
1286 IF (PAR.GT.1.0D-2) THEN
1287 PAR=MIN(PAR,1.0D 0/DBLE(NIT))
1295 * INCOMPLETE LANCZOS TRIDIAGONALIZATION
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
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
1321 CALL MXVCOR(NOS1,ETA0,XO)
1326 IF (T.GE.1.0D 5) GO TO 13118
1328 * SOLUTION OF THE TRIDIAGONAL SYSTEM
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)
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)
1342 * MARQUARDT PARAMETER T IS COMPUTED USING THE ONE-DIMENSIONAL
1345 T=T+(RHO/MXVVDP(NOS1,XS,S,S))*((SQRT(RHO)-RHO2)/RHO2)
1349 CALL MXVNEG(NF,G,XO)
1353 * INCOMPLETE GILL-MURRAY DECOMPOSITION
1357 IF (2*M.GE.MMAX) THEN
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
1369 IF (INF.NE.0) NOS2=0
1370 ELSE IF (MOS3.EQ.1) THEN
1371 IF (INF.LT.0) NOS2=0
1377 * PRELIMINARY INEXACT SOLUTION
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
1386 CALL MXVSCL(NF,XDEL/SNORM,S,S)
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
1400 CALL MXVSET(NF,0.0D 0,S)
1401 CALL MXVNEG(NF,G,XS)
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)
1407 RHO=MXUDOT(NF,XS,XO,IX,KBF)
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)
1413 CALL MXSSMM(NF,H,IH,JH,XO,GO)
1415 ALF=MXUDOT(NF,XO,GO,IX,KBF)
1416 IF (ALF.LE.0.0D 0) GO TO 13160
1418 RHO2=SQRT(MXUDEL(NF,ALF,XO,S,IX,KBF))
1419 IF (RHO2.GE.XDEL) GO TO 13160
1423 CALL MXUDIR(NF, ALF,XO,S,S,IX,KBF)
1424 CALL MXUDIR(NF,-ALF,GO,XS,XS,IX,KBF)
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
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)
1435 CALL MXUDIR(NF,ALF,XO,GO,XO,IX,KBF)
1438 CALL MXUDIR(NF,ALF,XO,XS,XO,IX,KBF)
1443 * AN INEXACT SOLUTION IS OBTAINED
1446 SNORM=SQRT(MXUDOT(NF,S,S,IX,KBF))
1450 * BOUNDARY STEP IS COMPUTED
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)
1461 PP=MXSSMQ(NF,H,IH,JH,S,S)*0.5D 0
1462 IF (KD.GT.0) P=MXUDOT(NF,G,S,IX,KBF)
1465 * SUBROUTINE PDSGM7 ALL SYSTEMS 01/09/22
1467 * COMPUTATION OF A TRUST-REGION STEP BY THE MORE-SORENSEN METHOD WITH
1468 * DIRECT MATRIX DECOMPOSITIONS.
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
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.
1532 * SUBPROGRAMS USED :
1533 * S PNSTEP COMPUTATION OF THE BOUNDARY STEP.
1534 * S MXSPCA ADDITION OF THE LEVENBERG-MARQUARDT TERM TO THE SPARSE
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
1550 * S MXSSMG GERSHGORIN BOUNDS FOR EIGENVALUES OF A SPARSE SYMMETRIC
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.
1561 * J.J.MORE, D.C.SORENSEN: COMPUTING A TRUST REGION STEP. REPORT NO.
1562 * ANL-81-83, ARGONNE NATIONAL LAB. 1981.
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
1576 * DIRECTION DETERMINATION
1578 IF (IDEC.LT.0) IDEC=0
1584 GNORM=SQRT(MXUDOT(NF,G,G,IX,KBF))
1585 IF (XDEL.LE.0.0D 0) THEN
1587 * INITIAL TRUST REGION BOUND
1589 RHO1=MXSSMQ(NF,H,IH,JH,G,G)
1591 IF (RHO1.LE.0.0D 0) THEN
1594 XDEL=(RHO2/RHO1)*GNORM
1596 IF (IEST.EQ.1) XDEL=MIN(XDEL,4.0D 0*(F-FMIN)/GNORM)
1600 * INITIAL BOUNDS FOR THE PARAMETER T
1605 E=-MXSSDL(NF,H,IH,JH,INF)
1606 CALL MXSSMG(NF,H,IH,JH,EL,EU,S)
1609 ELSE IF (IDIR.EQ.1) THEN
1611 TL=MAX(TL,GNORM/XDEL-EU)
1613 ELSE IF (IDIR.EQ.2) THEN
1616 TU=MIN(TU,GNORM/XDEL-EL)
1624 IF (T.LE.E.AND.NRED.NE.0) THEN
1626 * THE PARAMETER T IS SHIFTED
1629 T=MAX(T,TL+0.1D 0*(TU-TL))
1630 T=MIN(T,TL+0.9D 0*(TU-TL))
1633 CALL MXSPCT(NF,MM,MH,MMAX,H,JH,PSL,ITERM)
1634 IF (ITERM.NE.0) THEN
1638 * SPARSE GILL-MURRAY DECOMPOSITION
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)
1645 * NEW ESTIMATION E IS COMPUTED (THE MATRIX IS NOT POSITIVE DEFINITE)
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))
1660 * STEP S IS COMPUTED
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))
1669 IF (TU-TL.LE.1.0D-8) THEN
1671 * INTERVAL IS TOO SMALL
1673 IF (T.NE.0.0D 0) THEN
1679 ELSE IF (NRED.GE.20) THEN
1681 * MAXIMUM NUMBER OF OLC REDUCTIONS
1685 ELSE IF (SNORM.GT.DEL2*XDEL) THEN
1691 ELSE IF (SNORM.LT.DEL1*XDEL) THEN
1692 IF (T.NE.0.0D 0) THEN
1699 * STEP IS ACCEPTABLE
1709 * TRYING TO USE BOUNDARY STEP
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)
1727 IF (GNORM.LE.0.0D 0) THEN
1731 * NEW T IS COMPUTED USING ONE STEP OF THE NEWTON METHOD FOR
1732 * NONLINEAR EQUATION
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)
1743 PP=MXSSMQ(NF,H,IH,JH,S,S)*0.5D 0
1745 IF (KD.GT.0) P=MXUDOT(NF,G,S,IX,KBF)
1748 * SUBROUTINE PDSLM1 ALL SYSTEMS 01/09/22
1750 * DIRECTION DETERMINATION FOR LINE SEARCH USING DIRECT MATRIX
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
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.
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.
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
1818 DOUBLE PRECISION ALF,BET
1819 DOUBLE PRECISION MXUDOT
1821 * DIRECTION DETERMINATION
1823 IF (IDEC.LT.0) IDEC=0
1826 CALL MXSPCT(NF,MM,MH,MMAX,H,JH,PSL,ITERM)
1827 IF (ITERM.NE.0) RETURN
1829 * SPARSE GILL-MURRAY DECOMPOSITION
1832 CALL MXSPCF(NF,H(MM+1),PSL,JH(MM+1),WN11,WN12,XO,INF,ALF,BET)
1835 ELSE IF (IDEC.EQ.1) THEN
1840 GNORM=SQRT(MXUDOT(NF,G,G,IX,KBF))
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)
1849 SNORM=SQRT(MXUDOT(NF,S,S,IX,KBF))
1852 * SUBROUTINE PDSLM3 ALL SYSTEMS 01/09/22
1854 * DIRECTION DETERMINATION FOR LINE SEARCH USING CONJUGATE GRADIENT
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
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.
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
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.
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,
1930 INTEGER NOS2,NRED,MMX,INF
1931 DOUBLE PRECISION PAR,ALF,EPS,RHO,RHO1,RHO2,SIG
1932 DOUBLE PRECISION MXUDOT
1935 * DIRECTION DETERMINATION
1941 IF (IDEC.LT.0) IDEC=0
1942 IF (IDEC.NE.0.AND.IDEC.NE.1) THEN
1945 ELSE IF (IDEC.EQ.0) THEN
1948 * INCOMPLETE GILL-MURRAY DECOMPOSITION
1951 IF (2*M.GE.MMAX) THEN
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
1961 IF (INF.NE.0) NOS2=0
1966 RHO1=MXUDOT(NF,G,G,IX,KBF)
1968 PAR=MIN(EPS,SQRT(GNORM))
1969 IF (PAR.GT.1.0D-2) THEN
1970 PAR=MIN(PAR,1.0D 0/DBLE(NIT))
1975 * PRELIMINARY INEXACT SOLUTION
1977 CALL MXVNEG(NF,G,XO)
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))
1994 CALL MXVSET(NF,0.0D 0,S)
1995 CALL MXVNEG(NF,G,XS)
1997 CALL MXVNEG(NF,G,XO)
1998 ELSE IF (MOS2.GT.2) THEN
1999 RHO=MXUDOT(NF,XS,XO,IX,KBF)
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)
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
2011 C IF (ALF.LE.1.0D-8*SIG) THEN
2013 * CG FAILS (THE MATRIX IS NOT POSITIVE DEFINITE)
2028 CALL MXUDIR(NF, ALF,XO,S,S,IX,KBF)
2029 CALL MXUDIR(NF,-ALF,GO,XS,XS,IX,KBF)
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
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)
2040 CALL MXUDIR(NF,ALF,XO,GO,XO,IX,KBF)
2043 CALL MXUDIR(NF,ALF,XO,XS,XO,IX,KBF)
2046 C SIG=RHO2+ALF*ALF*SIG
2050 * SUBROUTINE PF1HS2 ALL SYSTEMS 99/12/01
2052 * NUMERICAL COMPUTATION OF THE HESSIAN MATRIX OF THE MODEL FUNCTION
2053 * USING ITS GRADIENTS - SPARSE VERSION USING DIRECT COLOURING METHOD.
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
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.
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.
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(*),
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
2098 IF (3*MM-NF+ML.GE.M) THEN
2105 CALL MXVCOP(NF,X,XO)
2107 * WIDTHEN THE STRUCTURE
2113 CALL MXSTG1(NF,MX,IH,JH,WN12,WN11)
2114 CALL MXVSET(K,0.0D 0,HF)
2117 IF (IVAR.GT.NF) GO TO 870
2119 IF (COL(J).GE.1) THEN
2130 IF (IX(L).LE.-7) GO TO 400
2135 XS(L)=ETA*MAX(ABS(X(L)),1.0D 0)*SIGN(1.0D 0,X(L))
2144 * NUMERICAL DIFFERENTIATION
2147 * SET AUXILIARY VECTOR DISCERNING THE SINGLETONS IN A GROUP TO ZERO
2153 * DISCERN SINGLETONS OF THE GROUP OF THE SAME COLOR.
2157 DO 500 K=IH(L),IH(L+1)-1
2159 IF (WN11(K1).EQ.0) THEN
2167 * NUMERICAL VALUES COMPUTATION
2171 DO 700 K=IH(L),IH(L+1)-1
2173 IF (WN11(K1).GT.0) THEN
2174 HF(K)=(GF(K1)-GO(K1))/XS(L)
2179 * SET THE ORIGINAL VALUE OF X FOR THE COMPONENTS OF THE ACTUAL COLOR.
2189 * MOVE THE ELEMENTS OF THE HESSIAN APPROXIMATION INTO THE UPPER
2201 IF (HF(L).EQ.0) THEN
2203 ELSE IF (HF(L).NE.0.AND.HF(J).NE.0) THEN
2204 HF(L)=0.5D 0*(HF(J)+HF(L))
2211 * SHRINK THE STRUCTURE
2213 CALL MXSTL2(NF,MX,HF,IH,JH,WN12)
2221 CALL MXVCOP(NF,XO,X)
2226 * SUBROUTINE PFSEB4 ALL SYSTEMS 98/12/01
2228 * COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE PARTITIONED HESSIAN
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.
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
2256 IF (LL.EQ.3.OR.LL.EQ.4) THEN
2257 TEMP= CZU(KC)-CZL(KC)
2258 ELSE IF (LL.EQ.1) THEN
2260 ELSE IF (LL.EQ.2) THEN
2262 ELSE IF (LL.EQ.5) THEN
2265 IF (JOB.EQ.1) TEMP=ABS(TEMP)
2266 ELSE IF (JOB.EQ.2) THEN
2267 IF (ICA(KC).GE.0) GO TO 7
2284 2 IF (JH(JF).LT.J) THEN
2288 B(JF)=B(JF)+TEMP*CH(K)
2297 * SUBROUTINE PFSEB5 ALL SYSTEMS 06/12/01
2299 * COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE PARTITIONED HESSIAN
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.
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
2325 ELSE IF (JOB.EQ.1) THEN
2342 2 IF (JH(JF).LT.J) THEN
2346 B(JF)=B(JF)+TEMP*CH(K)
2355 * SUBROUTINE PFSED3 ALL SYSTEMS 07/12/01
2357 * COMPRESSED SPARSE STRUCTURE OF THE HESSIAN MATRIX IS COMPUTED FROM
2358 * THE COORDINATE FORM.
2361 * II NF NUMBER OF VARIABLES.
2362 * II M NUMBER OF NONZERO ELEMENTS IN THE UPPER PART OF THE SPARSE
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.
2370 SUBROUTINE PFSED3(NF,M,IH,JH,IER)
2371 INTEGER NF,M,IH(*),JH(*),IER
2375 IF (IH(J).GT.JH(J)) THEN
2385 CALL MXVSR7(M+NF,IH,JH)
2386 IF (IH(1).LT.1.OR.IH(M+NF).GT.NF) THEN
2392 IF (IH(J).EQ.K) THEN
2403 CALL MXVSRT(L,JH(K))
2404 IF (JH(K).LT.1.OR.JH(K+L-1).GT.NF) THEN
2411 IF (J.GT.1.AND.JH(K).EQ.JH(K-1)) THEN
2419 IH(NF+1)=IH(NF+1)-LL
2423 * SUBROUTINE PFSET2 ALL SYSTEMS 97/12/01
2425 * COMPUTATION OF THE NUMBER OF NONZERO ELEMENTS OF THE SPARSE
2426 * HESSIAN MATRIX STORED IN THE BLOCKED FORM.
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
2435 SUBROUTINE PFSET2(NA,MB,MC,IAG)
2436 INTEGER NA,MB,MC,IAG(*)
2444 MC=MAX(MC,L*(L+1)/2)
2448 * SUBROUTINE PFSET3 ALL SYSTEMS 97/12/01
2450 * COMPUTATION OF THE SPARSE STRUCTURE OF THE HESSIAN MATRIX FROM THE
2451 * SPARSE STRUCTURE OF THE JACOBIAN MATRIX.
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.
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
2480 2 IF (JH(JF).LT.J.AND.JF.LE.LF) THEN
2482 IF (JF.LE.LF) GO TO 2
2484 IF (JH(JF).GT.J .OR.JF.GT.LF) THEN
2505 * SUBROUTINE PFSET4 ALL SYSTEMS 98/12/01
2507 * COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE PARTITIONED HESSIAN
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.
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
2537 2 IF (JH(JF).LT.J) THEN
2550 * FUNCTION PNFUZ1 ALL SYSTEMS 01/09/22
2552 * COMPUTATION OF LOWER AND UPPER LAGRANGE MULTIPLIERS.
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
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
2568 FUNCTION PNFUZ1(Z,NA,RPF3,AF,AZL,AZU,IEXT)
2570 DOUBLE PRECISION Z,RPF3,AF(*),AZL(*),AZU(*),PNFUZ1
2575 AZU(KA)=RPF3/(Z-AF(KA))
2576 PNFUZ1=PNFUZ1-AZU(KA)
2579 AZL(KA)=RPF3/(Z+AF(KA))
2580 PNFUZ1=PNFUZ1-AZL(KA)
2585 * SUBROUTINE PNINT1 ALL SYSTEMS 91/12/01
2587 * EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH WITH DIRECTIONAL
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
2604 * IO MERR ERROR INDICATOR. MERR=0 FOR NORMAL RETURN.
2607 * EXTRAPOLATION OR INTERPOLATION WITH STANDARD MODEL FUNCTIONS.
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,
2617 IF (MODE.LE.0) RETURN
2618 IF (PL.GE.0.0D 0) THEN
2621 ELSE IF (RU.LE.RL) THEN
2637 ELSE IF (NTYP.EQ.MTYP) THEN
2638 A = (FU-FL)/(PL*(RU-RL))
2643 * QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH ONE DIRECTIONAL
2646 DEN = 2.0D 0*(1.0D 0-A)
2647 ELSE IF (NTYP.EQ.3) THEN
2649 * QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH TWO DIRECTIONAL
2653 ELSE IF (NTYP.EQ.4) THEN
2655 * CUBIC EXTRAPOLATION OR INTERPOLATION
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
2662 ELSE IF (NTYP.EQ.5) THEN
2664 * CONIC EXTRAPOLATION OR INTERPOLATION
2667 IF (DIS.LT.0.0D 0) GO TO 1
2669 IF (DEN.LE.0.0D 0) GO TO 1
2670 DEN = 1.0D 0 - B*(1.0D 0/DEN)**3
2672 IF (MODE.EQ.1.AND.DEN.GT.0.0D 0.AND.DEN.LT.1.0D 0) THEN
2674 * EXTRAPOLATION ACCEPTED
2676 R = RL + (RU-RL)/DEN
2680 ELSE IF (MODE.EQ.2.AND.DEN.GT.1.0D 0) THEN
2682 * INTERPOLATION ACCEPTED
2684 R = RL + (RU-RL)/DEN
2685 IF (RL.EQ.0.0D 0) THEN
2686 R = MAX(R,RL+C2L*(RU-RL))
2688 R = MAX(R,RL+C3L*(RU-RL))
2690 R = MIN(R,RL+C2U*(RU-RL))
2695 * SUBROUTINE PNINT3 ALL SYSTEMS 91/12/01
2697 * EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH WITHOUT DIRECTIONAL
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
2715 * IO MERR ERROR INDICATOR. MERR=0 FOR NORMAL RETURN.
2718 * EXTRAPOLATION OR INTERPOLATION WITH STANDARD MODEL FUNCTIONS.
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
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)
2730 IF (MODE .LE. 0) RETURN
2731 IF (PO .GE. ZERO) THEN
2734 ELSE IF (RU .LE. RL) THEN
2740 DO 1 NTYP = MTYP, 1, -1
2741 IF (NTYP .EQ. 1) THEN
2745 IF (MODE .EQ. 1) THEN
2748 ELSE IF (RI-RL.LE.RU-RI) THEN
2755 ELSE IF (NTYP.EQ.MTYP.AND.L1) THEN
2756 IF (.NOT.L2) AI=(FI-FO)/(RI*PO)
2759 IF (L1.AND.(NTYP.EQ.2.OR.L2)) THEN
2761 * TWO POINT QUADRATIC EXTRAPOLATION OR INTERPOLATION
2763 IF (AU.GE.ONE) GO TO 1
2765 ELSE IF (.NOT.L1.OR..NOT.L2.AND.NTYP.EQ.3) THEN
2767 * THREE POINT QUADRATIC EXTRAPOLATION OR INTERPOLATION
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
2776 * THREE POINT CUBIC EXTRAPOLATION OR INTERPOLATION
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
2784 IF (DEN.EQ.ZERO) GO TO 1
2789 IF (MODE .EQ. 1 .AND. R .GT. RU) THEN
2791 * EXTRAPOLATION ACCEPTED
2796 ELSE IF (MODE .EQ. 2 .AND. R .GT. RL .AND. R .LT. RU) THEN
2798 * INTERPOLATION ACCEPTED
2800 IF (RI.EQ.ZERO.AND.NTYP.NE.4) THEN
2801 R = MAX( R, RL + C2L*(RU-RL))
2803 R = MAX( R, RL + C3L*(RU-RL))
2805 R = MIN( R, RL + C2U*(RU-RL))
2806 IF (R.EQ.RI) GO TO 1
2811 * SUBROUTINE PNNEQ1 ALL SYSTEMS 92/12/01
2813 * SOLUTION OF A SINGLE NONLINEAR EQUATION.
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.
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.
2831 SUBROUTINE PNNEQ1(AA,BB,X,F,EPSX,EPSF,IC,IE,ISYS)
2832 DOUBLE PRECISION AA,BB,X,F,EPSX,EPSF
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
2842 GO TO (1,2,3,4,6) ISYS+1
2845 IF (ITMAX.LE.0) ITMAX=100
2851 IF (ABS(F).LE.EPSF) GO TO 7
2858 IF (ABS(F).LE.EPSF) GO TO 7
2860 IF (FA*FB.GT.0.0D 0) THEN
2875 IF (F1*F2.GT.0.0D 0) THEN
2892 IF (ABS(F1).LT.ABS(F2)) THEN
2899 DEL=EPSX*(ABS(X)+ONE)
2900 IF (ABS(F).LE.EPSF.OR.D.LE.TWO*DEL) GO TO 7
2903 IF (THREE*D.LE.TWO*D1) THEN
2921 IF (ABS(A).LE.1.0D-10) THEN
2922 R=(F2*X1-F1*X2)/(F2-F1)
2925 IF (R.LT.0.0D 0) THEN
2926 R=(F2*X1-F1*X2)/(F2-F1)
2931 IF (ABS(RA-Z).LE.ABS(RB-Z)) THEN
2936 IF (R.LE.MIN(X1,X2).OR.R.GE.MAX(X1,X2)) THEN
2937 R=(F2*X1-F1*X2)/(F2-F1)
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
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
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
2964 IF (F*SIGN(ONE,F1).GE.0.0D 0) THEN
2965 IF (D.LE.ABS(X3-X)) THEN
2983 IF (ITER.LE.ITMAX) GO TO 5
2988 * SUBROUTINE PNSTEP ALL SYSTEMS 89/12/01
2990 * DETERMINATION OF A SCALING FACTOR FOR THE BOUNDARY STEP.
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.
3000 SUBROUTINE PNSTEP(DEL,A,B,C,ALF)
3001 DOUBLE PRECISION DEL, A, B, C, ALF
3002 DOUBLE PRECISION DEN, DIS
3004 DEN = (DEL+A) * (DEL-A)
3005 IF (DEN .LE. 0.0D 0) RETURN
3007 IF (B .GE. 0.0D 0) THEN
3008 ALF = DEN / (SQRT(DIS) + B)
3010 ALF = (SQRT(DIS) - B) / C
3014 * SUBROUTINE PNSTP4 ALL SYSTEMS 99/12/01
3016 * STEPSIZE SELECTION USING POLYHEDRAL APPROXIMATION
3017 * FOR DESCENT STEP IN NONCONVEX VARIABLE METRIC METHOD.
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.
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)
3044 * INITIAL CHOICE OF POSSIBLY ACTIVE LINES
3056 DX = X(I) - AY(JN+I)
3062 IF (MOS3.NE.2) R = R** (DBLE(MOS3)*0.5D0)
3063 ALF = MAX(ABS(ALFL),ETA5*R)
3065 IF (R*R+ (ALF+ALF)/DF.GT.1.0D-6) THEN
3075 IF (BET.GT.0.0D0) TB = MIN(TB,ALF/ (BET-DF))
3076 BETR = MAX(BETR,BET-ALF)
3080 IF (BETR.LE.DF*0.5D0) RETURN
3084 IF (BETR.LE.0.0D0) THEN
3085 IF (T.LT.1.0D0 .OR. BETR.EQ.0.0D0) RETURN
3092 30 IF (LQ.GE.1) THEN
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))
3098 IF (ABS(BETR-BETL)+ABS(ALFR-ALFL).LT.-1.0D-4*DF) RETURN
3099 R = (ALFR-ALFL)/ (BETR-BETL)
3101 IF (ABS(T-R).LT.1.0D-4) RETURN
3108 IF (ALF.LT.0.0D0) GO TO 40
3118 IF (BET.EQ.0.0D0) RETURN
3120 * NEW INTERVAL SELECTION
3123 IF (BET.LT.0.0D0) THEN
3143 * SUBROUTINE PNSTP5 ALL SYSTEMS 99/12/01
3145 * STEPSIZE SELECTION USING POLYHEDRAL APPROXIMATION
3146 * FOR NULL STEP IN NONCONVEX VARIABLE METRIC METHOD.
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.
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
3173 * INITIAL CHOICE OF POSSIBLY ACTIVE PARABOLAS
3185 DX = X(I) - AY(JN+I)
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
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
3224 IF (ABS(T-Q).LT.1.0D-3) RETURN
3230 IF (ALF.LT.0.0D0) GO TO 40
3241 IF (Q.EQ.0.0D0) RETURN
3243 * NEW INTERVAL SELECTION
3246 IF (Q.LT.0.0D0) THEN
3255 * SUBROUTINE PP0BA1 ALL SYSTEMS 05/12/01
3257 * EVALUATION OF THE BARRIER FUNCTION FOR THE SUM OF ABSOLUTE VALUES.
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.
3265 SUBROUTINE PP0BA1(NA,AS,RPF3,F)
3267 DOUBLE PRECISION AS(*),RPF3,F
3269 F=-DBLE(NA)*RPF3*LOG(2.0D 0*RPF3)
3271 F=F+AS(KA)-RPF3*LOG(AS(KA))
3275 * SUBROUTINE PP0BX1 ALL SYSTEMS 05/12/01
3277 * EVALUATION OF THE BARRIER FUNCTION FOR THE MINIMAX OPTIMIZATION.
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
3290 * II IEXT KIND OF THE MINIMAX APPROXIMATION. IEXT=0-CHEBYSHEV
3291 * APPROXIMATION. IEXT=-1-MINIMAX. IEXT=+1-MAXIMIN.
3293 SUBROUTINE PP0BX1(NA,Z,AF,F,FF,PAR,RPF3,MEP,IEXT)
3295 DOUBLE PRECISION Z,AF(*),PAR,RPF3,F,FF
3312 ELSE IF (MEP.EQ.2) THEN
3316 IF (Z-FA.LE.PAR) THEN
3319 F=F+(2.0D 0-0.5D 0*PAR/(Z-FA))*RPF3*PAR/(Z-FA)
3323 IF (Z+FA.LE.PAR) THEN
3326 F=F+(2.0D 0-0.5D 0*PAR/(Z+FA))*RPF3*PAR/(Z+FA)
3330 ELSE IF (MEP.EQ.3) THEN
3334 F=F+RPF3*LOG(1.0D 0/(Z-FA)+1.0D 0)
3337 F=F+RPF3*LOG(1.0D 0/(Z+FA)+1.0D 0)
3340 ELSE IF (MEP.EQ.4) THEN
3344 F=F+RPF3*RPF3/(Z-FA)
3347 F=F+RPF3*RPF3/(Z+FA)
3354 * SUBROUTINE PP1MX3 ALL SYSTEMS 05/12/01
3356 * COMPUTATION OF THE VALUE AND THE GRADIENT OF THE LAGRANGIAN FUNCTION
3357 * FOR THE MINIMAX OPTIMIZATION.
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.
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.
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
3395 IF (KD.LE.LD) RETURN
3396 IF (KD.GE.0.AND.LD.LT.0) THEN
3399 IF (KD.GE.1.AND.LD.LT.1) THEN
3400 CALL MXVSET(NF,0.0D 0,G)
3404 IF (LD.GE.0) GO TO 1
3405 CALL FUN(NF,KA,X,FA)
3406 IF (ISNA.GE.1) AF(KA)=FA
3408 IF (KA.EQ.1) F=ABS(FA)
3410 ELSE IF (IEXT.LT.0) THEN
3413 ELSE IF (IEXT.GT.0) THEN
3417 1 IF (KD.LT.1) GO TO 3
3418 IF (LD.GE.1) GO TO 3
3419 CALL DFUN(NF,KA,X,GA)
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)
3431 IF (ISNA.GE.2) AG(K)=GA(JP)
3437 * SUBROUTINE PP1SA3 ALL SYSTEMS 05/12/01
3439 * COMPUTATION OF THE VALUE AND THE GRADIENT OF THE LAGRANGIAN FUNCTION
3440 * FOR THE SUM OF ABSOLUTE VALUES.
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.
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.
3470 SUBROUTINE PP1SA3(NF,NA,X,GA,AG,IAG,JAG,G,AZ,FA,AF,F,KD,LD,NFV,
3472 INTEGER NF,NA,IAG(*),JAG(*),KD,LD,NFV,NFG,ISNA
3473 DOUBLE PRECISION X(*),GA(*),AG(*),G(*),AZ(*),FA,AF(*),F
3475 IF (KD.LE.LD) RETURN
3476 IF (KD.GE.0.AND.LD.LT.0) THEN
3480 IF (KD.GE.1.AND.LD.LT.1) THEN
3481 CALL MXVSET(NF,0.0D 0,G)
3485 IF (LD.GE.0) GO TO 1
3486 CALL FUN(NF,KA,X,FA)
3487 IF (ISNA.GE.1) AF(KA)=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)
3496 G(JP)=G(JP)+AZ(KA)*GA(JP)
3497 IF (ISNA.GE.2) AG(K)=GA(JP)
3503 * SUBROUTINE PPLAG1 ALL SYSTEMS 05/12/01
3505 * COMPUTATION OF THE LAGRANGE MULTIPLIERS FOR THE SUM OF ABSOLUTE
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.
3515 SUBROUTINE PPLAG1(NA,AF,AS,AZ,RPF3)
3517 DOUBLE PRECISION AF(*),AS(*),AZ(*),RPF3
3522 AS(KA)=RPF3+SQRT(RPF3**2+FA**2)
3527 * SUBROUTINE PS0G01 ALL SYSTEMS 97/12/01
3529 * SIMPLE SEARCH WITH TRUST REGION UPDATE.
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
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.
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.
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,
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
3597 IF (ISYS.EQ.1) GO TO 2
3605 * COMPUTATION OF THE NEW FUNCTION VALUE
3613 IF (KTERS.LT.0.OR.KTERS.GT.5) THEN
3618 IF (DF.LT.EPS4*DFPRED) THEN
3620 * STEP IS TOO LARGE, IT HAS TO BE REDUCED
3624 ELSE IF (MES1.EQ.2) THEN
3625 XDEL=BET2*MIN(0.5D 0*XDEL,SNORM)
3627 XDEL=0.5D 0*PO*SNORM/(PO+DF)
3628 XDEL=MAX(XDEL,BET1*SNORM)
3629 XDEL=MIN(XDEL,BET2*SNORM)
3635 IF (ITERD.GT.2) NRED2=NRED2+1
3637 ELSE IF (DF.LE.EPS5*DFPRED) THEN
3644 * STEP IS TOO SMALL, IT HAS TO BE ENLARGED
3647 XDEL=MAX(XDEL,GAM1*SNORM)
3648 ELSE IF (ITERD.GT.2) THEN
3653 XDEL=MIN(XDEL,XMAX,GAM2*SNORM)
3655 IF (NRED1.GE.MRED) THEN
3665 IF (XDEL.GE.XMAX) MAXST=1
3674 * SUBROUTINE PS0L02 ALL SYSTEMS 97/12/01
3676 * EXTENDED LINE SEARCH WITHOUT DIRECTIONAL DERIVATIVES.
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
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
3724 * IU ISYS CONTROL PARAMETER.
3727 * S PNINT3 EXTRAPOLATION OR INTERPOLATION WITHOUT DIRECTIONAL
3731 * SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH EXTENDED TERMINATION
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,
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
3749 IF (PO.GE.0.0D 0) THEN
3754 IF (RMAX.LE.0.0D 0) THEN
3759 * INITIAL STEPSIZE SELECTION
3761 IF (INITS.GT.0) THEN
3763 ELSE IF (IEST.EQ.0) THEN
3766 RTEMP=MAX(F-FP,FMIN-F)
3772 IF (INIT1.EQ.0) THEN
3773 ELSE IF (INIT1.EQ.1.OR.INITS.GE.1.AND.IEST.EQ.0) THEN
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
3793 * NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION)
3795 2 CALL PNINT3(RO,RL,RU,RI,FO,FL,FU,FI,PO,R,MODE,MTYP,MERR)
3799 ELSE IF (MODE.EQ.1) THEN
3802 ELSE IF (MODE.EQ.2) THEN
3806 * COMPUTATION OF THE NEW FUNCTION VALUE
3813 IF (ITERS.NE.0) GO TO 4
3818 L1=R.LE.RMIN.AND.NIT.NE.KIT
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
3828 * TEST ON TERMINATION
3830 IF (L1.AND..NOT.L3) THEN
3833 ELSE IF (L2.AND..NOT.F.GE.FU) THEN
3839 ELSE IF (L3.AND.L7.AND.KTERS.EQ.5) THEN
3842 ELSE IF (L3.AND.L4.AND.L7.AND.(KTERS.EQ.2.OR.KTERS.EQ.3.OR.
3846 ELSE IF (KTERS.LT.0.OR.KTERS.EQ.6.AND.L7) THEN
3849 ELSE IF (ABS(NRED).GE.MRED) THEN
3857 IF (F.GE.FMAX) MTYP=1
3861 * INTERVAL CHANGE AFTER EXTRAPOLATION
3872 ELSE IF ( MES1 .EQ. 1) THEN
3876 * INTERVAL CHANGE AFTER INTERPOLATION
3878 ELSE IF (R.LE.RI) THEN
3903 * SUBROUTINE PS1L01 ALL SYSTEMS 97/12/01
3905 * STANDARD LINE SEARCH WITH DIRECTIONAL DERIVATIVES.
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
3926 * RO PAR2 PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC
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
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
3959 * IU ISYS CONTROL PARAMETER.
3962 * S PNINT1 EXTRAPOLATION OR INTERPOLATION WITH DIRECTIONAL
3966 * SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH STANDARD TERMINATION
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,
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
3988 IF (PO.GE.0.0D 0) THEN
3993 IF (RMAX.LE.0.0D 0) THEN
3998 * INITIAL STEPSIZE SELECTION
4000 IF (INITS.GT.0) THEN
4002 ELSE IF (IEST.EQ.0) THEN
4005 RTEMP=MAX(F-FP,FMIN-F)
4011 IF (INIT1.EQ.0) THEN
4012 ELSE IF (INIT1.EQ.1.OR.INITS.GE.1.AND.IEST.EQ.0) THEN
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
4028 * NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION)
4030 2 CALL PNINT1(RL,RU,FL,FU,PL,PU,R,MODE,MTYP,MERR)
4034 ELSE IF (MODE.EQ.1) THEN
4037 ELSE IF (MODE.EQ.2) THEN
4041 * COMPUTATION OF THE NEW FUNCTION VALUE AND THE NEW DIRECTIONAL
4053 IF (ITERS.NE.0) GO TO 4
4058 L1=R.LE.RMIN.AND.NIT.NE.KIT
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
4067 M1=ABS(P).LE.CON*ABS(PO).AND.FO-F.GE.(CON1/CON)*ABS(FO)
4071 M2=ABS(P).LE.0.5D 0*ABS(PO).AND.ABS(FO-F).LE.2.0D 0*CON1*ABS(FO)
4078 * TEST ON TERMINATION
4080 IF (L1.AND..NOT.L3) THEN
4083 ELSE IF (L2.AND.L3.AND..NOT.L5) THEN
4086 ELSE IF (M3.AND.MES1.EQ.3) THEN
4089 ELSE IF (L3.AND.L5.AND.L7) THEN
4092 ELSE IF (KTERS.LT.0.OR.KTERS.EQ.6.AND.L7) THEN
4095 ELSE IF (ABS(NRED).GE.MRED) THEN
4104 IF (F.GE.FMAX) MTYP=1
4108 * INTERVAL CHANGE AFTER EXTRAPOLATION
4119 ELSE IF ( MES1 .EQ. 1) THEN
4124 * INTERVAL CHANGE AFTER INTERPOLATION
4140 * SUBROUTINE PS1L18 ALL SYSTEMS 99/12/01
4142 * SPECIAL LINE SEARCH FOR NONSMOOTH NONCONVEX VARIABLE METRIC METHOD.
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
4177 * IU ISYS CONTROL PARAMETER.
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.
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.
4197 * SPECIAL METHOD OF STEP LENGTH DETERMINATION.
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,
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
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)
4216 CALL PNSTP5(N,MA,MAL,X,AF,AG,AY,S,F,PO,T,TB,ETA5,ETA9,MOS3)
4218 10 T = MIN(MAX(T,TMIN),TMAX)
4224 * FUNCTION AND GRADIENT EVALUATION AT A NEW POINT
4226 20 CALL MXVDIR(N,T,S,U,X)
4234 * NULL/DESCENT STEP TEST (ITERS=0/1)
4237 IF (F.LE.FO-T* (EPS1+EPS1)*WK) THEN
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)
4261 * SUBROUTINE PUBBM1 ALL SYSTEMS 97/12/01
4263 * PARTITIONED VARIABLE METRIC UPDATE.
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.
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.
4294 SUBROUTINE PUBBM1(NA,AH,IAG,JAG,S,XO,AGO,ETA0,ETA9,ICOR,NIT,KIT,
4296 INTEGER NA,IAG(*),JAG(*),ICOR,NIT,KIT,ITERH,MET,
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
4302 IF (MET.LE.0) GO TO 22
4303 L1=MET1.GE.3.OR.MET1.EQ.2.AND.NIT.EQ.KIT
4311 * DETERMINATION OF THE PARAMETERS B, C
4313 B=MXWDOT(L,JAG(K),AGO(K),XO,2)
4315 IF (B.LE.1.0D 0/ETA9) GO TO 20
4317 IF (ABS(B).LE.1.0D 0/ETA9) GO TO 20
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
4325 * DETERMINATION OF THE PARAMETER GAM (SELF SCALING)
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)
4342 IF (B.LT.0.0D 0) INEG=INEG+1
4348 IF (ABS(DEN).GT.ETA0*ABS(C)) THEN
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)
4355 ELSE IF (B.LT.0.0D 0) THEN
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)
4367 IF (GAM.NE.1.0D 0) THEN
4368 CALL MXDSMS(L,AH(NB+1),1.0D 0/GAM)
4373 IF (INEG.GE.NA/2) ICOR=1
4377 * SUBROUTINE PUBBM2 ALL SYSTEMS 97/12/01
4379 * PARTITIONED VARIABLE METRIC UPDATE.
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
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.
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.
4411 SUBROUTINE PUBBM2(NA,AH,IAG,JAG,S,XO,AGO,ETA0,ETA9,NIT,KIT,ITERH,
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
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
4427 * DETERMINATION OF THE PARAMETERS B, C
4429 B=MXWDOT(L,JAG(K),AGO(K),XO,2)
4431 IF (B.LE.1.0D 0/ETA9) GO TO 20
4433 IF (ABS(B).LE.1.0D 0/ETA9) GO TO 20
4436 CALL MXBSBM(L,AH(NB+1),JAG(K),XO,S,1)
4437 C=MXWDOT(L,JAG(K),XO,S,1)
4439 IF (ABS(C).LE.1.0D 0/ETA9) GO TO 20
4441 IF (C.LE.1.0D 0/ETA9) GO TO 20
4444 IF (B.LE.0.0D 0) THEN
4446 * POWELL'S CORRECTION
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)
4456 * DETERMINATION OF THE PARAMETER GAM (SELF SCALING)
4460 IF (NIT.NE.KIT) THEN
4461 L3=GAM.LT.CON1.OR.GAM.GT.CON2
4463 ELSE IF (MET1.EQ.4) THEN
4474 ELSE IF (MET.EQ.2) THEN
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)
4485 ELSE IF (MET.EQ.3) THEN
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)
4496 ELSE IF (MET.EQ.4) THEN
4502 IF (ABS(DEN).LE.ETA0*ABS(C)) GO TO 18
4504 IF (DEN.LE.ETA0*C) GO TO 18
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)
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)
4520 IF (GAM.NE.1.0D 0) THEN
4521 CALL MXDSMS(L,AH(NB+1),1.0D 0/GAM)
4528 * SUBROUTINE PUBVI2 ALL SYSTEMS 04/12/01
4530 * NONSMOOTH VARIABLE METRIC UPDATE OF THE INVERSE HESSIAN MATRIX.
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
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.
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.
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
4571 CALL MXVDIF(L,AG(K),AGO(K),U)
4573 * DETERMINATION OF THE PARAMETERS B, C
4575 B=MXWDOT(L,JAG(K),U,XO,2)
4576 IF (ABS(B).LE.1.0D 0/ETA9) GO TO 20
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
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,
4588 IF (B.LT.0.0D 0) INEG=INEG+1
4590 LR=NNK.NE.0.AND.C.LT.GAM*B
4592 IF (B.LT.0.0D 0) GO TO 20
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)
4600 IF (GAM.NE.1.0D 0) THEN
4601 CALL MXDSMS(L,AH(NB+1),1.0D 0/GAM)
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)
4614 * SUBROUTINE PULCI3 ALL SYSTEMS 96/12/01
4616 * LIMITED STORAGE INVERSE COLUMN UPDATE METHODS.
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.
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.
4646 * LIMITED STORAGE VARIABLE METRIC METHODS.
4648 SUBROUTINE PULCI3(N,A,IA,JA,IP,ID,XM,GM,IM,XO,AFO,S,MF,NIT,KIT,
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
4655 DOUBLE PRECISION MXVMX1
4657 MM = MIN(NIT-KIT,MF)
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))
4666 TEMP = MXVMX1(N,AFO,II)
4667 IF (TEMP.LE.0.0D0) THEN
4677 * SUBROUTINE PULSP3 ALL SYSTEMS 02/12/01
4679 * LIMITED STORAGE VARIABLE METRIC UPDATE.
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).
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.
4707 * SHIFTED BFGS METHOD IN THE PRODUCT FORM.
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
4715 IF (B.LE.0.0D 0) THEN
4719 CALL MXDRMM(N,M,XM,GO,GR)
4725 * DETERMINATION OF THE PARAMETER SIG (SHIFT)
4729 IF (A.GT.0.0D 0) 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
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
4738 SIG=MAX(SIG,2.0D-1*POM)
4739 SIG=MIN(SIG,8.0D-1*POM)
4744 * COMPUTATION OF SHIFTED XO AND SHIFTED B
4747 CALL MXVDIR(N,-SIG,GO,XO,XO)
4749 * BFGS-BASED SHIFTED BFGS UPDATE
4752 CALL MXDCMU(N,M,XM,-1.0D 0/BB,XO,GR)
4753 CALL MXVSCL(N,SQRT(PAR/BB),XO,XM(N*M+1))
4759 * SUBROUTINE PULVP3 ALL SYSTEMS 03/12/01
4761 * RANK-TWO LIMITED-STORAGE VARIABLE-METRIC METHODS IN THE PRODUCT FORM.
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
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
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.
4800 * RANK-ONE LIMITED-STORAGE VARIABLE-METRIC METHOD IN THE PRODUCT FORM.
4802 SUBROUTINE PULVP3(N,M,XM,XR,GR,S,SO,XO,GO,R,PO,SIG,ITERH,MET2,
4804 INTEGER N,M,ITERH,MET2,MET3,MET5
4805 DOUBLE PRECISION XM(*),XR(*),GR(*),S(*),SO(*),XO(*),GO(*),
4807 DOUBLE PRECISION MXVDOT
4808 DOUBLE PRECISION DEN,POM,A,B,C,AA,BB,CC,AH,PAR,ZET
4814 IF (B.LE.0.0D 0) THEN
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.
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)
4834 * DETERMINATION OF THE PARAMETER SIG (SHIFT)
4837 IF (A.GT.0.0D 0) 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
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
4846 SIG=MAX(SIG,2.0D-1*POM)
4847 SIG=MIN(SIG,8.0D-1*POM)
4852 * COMPUTATION OF SHIFTED XO AND SHIFTED B
4855 CALL MXVDIR(N,-SIG,GO,XO,XO)
4857 * COMPUTATION OF THE PARAMETER RHO (CORRECTION)
4861 ELSE IF (MET2.EQ.2) THEN
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))
4871 * COMPUTATION OF THE PARAMETER THETA (BFGS)
4873 POM=SIGN(SQRT(PAR*B/CC),BB)
4875 * COMPUTATION OF Q AND P
4879 * RANK ONE UPDATE OF XM
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)
4886 * RANK TWO UPDATE OF XM
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)
4895 * SUBROUTINE PUSMM1 ALL SYSTEMS 97/12/01
4897 * VARIABLE METRIC UPDATE OF A SPARSE SYMMETRIC POSITIVE DEFINITE MATRIX
4898 * USING THE MARWIL PROJECTION.
4901 * II NF NUMBER OF VARIABLES.
4902 * RU H(M) POSITIVE DEFINITE APPROXIMATION OF THE SPARSE HESSIAN
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
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.
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.
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
4942 DOUBLE PRECISION MXUDOT
4943 DOUBLE PRECISION A,B,C,GAM
4947 * DETERMINATION OF THE PARAMETER C AND THE VECTOR S
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)
4955 CALL MXUDIF(NF,GO,G,S,IX,KBF)
4956 CALL MXVSCL(NF,R,S,S)
4964 B=MXUDOT(NF,XO,GO,IX,KBF)
4965 IF (B.GT.0.0D 0.AND.C.GT.0.0D 0) THEN
4967 CALL MXVSCL(MM,1.0D 0/GAM,H,H)
4968 CALL MXVSCL(NF,1.0D 0/GAM,S,S)
4971 CALL MXUDIR(NF,-1.0D 0,S,GO,S,IX,KBF)
4973 * RANK-ONE UPDATE PROJECTED USING MXSSMY
4975 CALL MXSSMY(NF,H,IH,JH,XS,S,XO)
4979 * SUBROUTINE PUSSD5 ALL SYSTEMS 97/12/01
4981 * INITIATION OF A DENSE SYMMETRIC POSITIVE DEFINITE MATRIX
4984 * II NA NUMBER OF APPROXIMATED FUNCTIONS.
4985 * RI AF(NA) VECTOR CONTAINING VALUES OF THE APPROXIMATED
4987 * RU AH(MB) POSITIVE DEFINITE APPROXIMATION OF THE PARTITIONED
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
4993 * II IH(NF+1) POINTERS OF THE DIAGONAL ELEMENTS OF THE SPARSE
4995 * II JH(M) INDICES OF THE NONZERO ELEMENTS OF THE SPARSE HESSIAN
4996 * MATRIX IN THE PACKED ROW FORM.
4998 * SUBPROGRAMS USED :
4999 * S PASSH2 COMPUTATION OF THE SPARSE HESSIAN MATRIX FROM THE
5000 * PARTITIONED HESSIAN MATRIX.
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
5011 CALL PASSH2(H,IH,JH,AH(NB+1),IAG,JAG,KA,AF(KA))
5016 * SUBROUTINE PYABU1 ALL SYSTEMS 04/12/01
5018 * SUBGRADIENT AGGREGATION FOR NONSMOOTH VARIABLE METRIC METHOD.
5021 * II NF NUMBER OF VARIABLES.
5022 * RI H(M) POSITIVE DEFINITE APPROXIMATION OF THE SPARSE HESSIAN
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.
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.
5045 SUBROUTINE PYABU1(NF,H,JH,PSL,PERM,G,GO,GV,S,U,V,ALF,ALFV,RHO,
5047 INTEGER NF,JH(*),PSL(*),PERM(*),JC
5048 DOUBLE PRECISION H(*),G(*),GO(*),GV(*),S(*),U(*),V(*),ALF,ALFV,
5050 DOUBLE PRECISION A,B,ALFM,LAM1,LAM2,PQ,PR,PRQR,QQP,QR,RR,RRP,RRQ,
5053 DOUBLE PRECISION ZERO,ONE,MXVDOT
5054 PARAMETER (ZERO=0.0D 0,ONE=1.0D 0)
5057 * General routine - here always input parameter ALFM=0
5077 PQ=PQ+RHO*(S(I)-U(I))**2
5080 PRQR=PRQR+RHO*U(I)*S(I)
5081 QQP=QQP+RHO+G(I)*(S(I)-U(I))
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)
5090 W1=ONE/H(PSL(I)+I-1)
5091 PQ=PQ+W1*(S(I)-U(I))**2
5094 PRQR=PRQR+W1*U(I)*S(I)
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
5104 IF (W.EQ.ZERO) GO TO 10
5107 IF (LAM1*(LAM1-ONE).LT.ZERO.AND.LAM2*(LAM1+LAM2-ONE).LT.ZERO)
5110 * MINIMUM ON THE BOUNDARY
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
5118 IF (ALFM.LE.ALFV) A=ONE
5119 IF (PR.GT.ZERO) A=MIN(ONE,MAX(ZERO,RRP/PR))
5126 IF (QQP*(QQP-PQ).GE.ZERO) GO TO 40
5127 IF (QR-RRQ-RRQ-QQP*QQP/PQ.GE.W) GO TO 40
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)
5133 ALFV=LAM1*ALFM+LAM2*ALF+A*ALFV
5135 GV(I)=LAM1*GO(I)+LAM2*G(I)+A*GV(I)
5139 * SUBROUTINE PYABU2 ALL SYSTEMS 04/12/01
5141 * SIMPLIFIED AGGREGATION FOR NONSMOOTH VARIABLE METRIC METHOD.
5144 * II NF NUMBER OF VARIABLES.
5145 * RI H(M) POSITIVE DEFINITE APPROXIMATION OF THE SPARSE HESSIAN
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.
5159 * SUBPROGRAMS USED :
5160 * S MXSPCB BACK SUBSTITUTION USING THE SPARSE DECOMPOSITION
5161 * OBTAINED BY MXSPCF.
5162 * S MXVSFP PERMUTATION OF A VECTOR.
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
5169 DOUBLE PRECISION ZERO,ONE
5170 PARAMETER (ZERO=0.0D 0,ONE=1.0D 0)
5183 CALL MXVSFP(NF,PERM,S,V)
5184 CALL MXSPCB(NF,H,PSL,JH,S,1)
5189 LAM=0.5D 0+SIGN(0.5D 0,P)
5190 IF (Q.GT.ZERO) LAM=MIN(ONE,MAX(ZERO,P/Q))
5194 GV(I)=LAM*G(I)+P*GV(I)
5198 * SUBROUTINE PYADC0 ALL SYSTEMS 98/12/01
5200 * NEW SIMPLE BOUNDS ARE ADDED TO THE ACTIVE SET.
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.
5211 SUBROUTINE PYADC0(NF,N,X,IX,XL,XU,INEW)
5212 INTEGER NF,N,IX(NF),INEW
5213 DOUBLE PRECISION X(*),XL(*),XU(*)
5222 ELSE IF ((IXI.EQ.1.OR.IXI.EQ.3.OR.IXI.EQ.4).AND.X(I).LE.XL(I))
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))
5241 IF (II.GT.0) INEW=INEW+1
5246 * SUBROUTINE PYBUN1 ALL SYSTEMS 97/12/01
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
5263 * SUBPROGRAMS USED :
5264 * S MXVCOP COPYING OF A VECTOR.
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
5280 IF (G(I).NE.AG(KN+I)) GO TO 2
5288 AF(K+MB*3)=AF(K+1+MB*3)
5290 CALL MXVCOP(N,AG(KN),AG(KN-N))
5291 CALL MXVCOP(N,AY(KN),AY(KN-N))
5298 IF (L.GT.0.AND.KN.EQ.0) THEN
5300 AF(3*MB+NB+1)=AF(3*MB+NB)
5302 CALL MXVCOP(N,AG(KN-N),AG(KN))
5303 CALL MXVCOP(N,AY(KN-N),AY(KN))
5310 CALL MXVCOP(N,G,AG(K))
5311 CALL MXVCOP(N,X,AY(K))
5314 * SUBROUTINE PYCSER ALL SYSTEMS 98/12/01
5316 * GROUP OF THE SAME COLOUR FOR THE POWELL-TOINT ALGORITHM FOR SPARSE
5317 * HESSIANS APPROXIMATIONS IS CREATED.
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
5329 * IU NCOL NUMBER OF COLOURS USED SO FAR.
5330 * IU CNM NUMBER OF COLUMNS THAT HAVE NOT BEEN COLOURED SO FAR.
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
5338 * DEFINITION OF THE INCIDENCE ARRAY A
5342 * ELEMENT WAS MARKED THAT IT IS INSERTED
5344 DO 100 I=IH(L),IH(L+1)-1
5347 * COLUMN OF THIS NUMBER HAS APPEARED IN ONE OF THE PREVIOUS GROUPS
5349 IF (COL(K).LT.NCOL) GO TO 100
5354 * COLUMN IS INSERTED
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
5363 IF (CNM.EQ.1) GO TO 250
5366 * TRANSFORMATION OF THE EXAMINED COLUMN I IS
5370 DO 300 J=IH(L),IH(L+1)-1
5372 IF (COL(K).LT.NCOL) GO TO 300
5373 IF (WN02(K).GE.NCOL) GO TO 200
5379 * COPY OF THE WN03 ARRAY INTO WN02 FOR THE COLUMN WAS ACCEPTED
5382 WN02(INT(WN03(K)))=NCOL
5383 DEG(INT(WN03(K)))=DEG(INT(WN03(K)))-1
5387 * INSERT THE COLUMN INTO THE PROCESSED GROUP
5391 * END OF THE MAIN CYCLE
5404 IF (COL(L).EQ.NCOL) THEN
5416 * SUBROUTINE PYFUT1 ALL SYSTEMS 98/12/01
5418 * TERMINATION CRITERIA AND TEST ON RESTART.
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.
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)
5474 IF (GMAX.LE.TOLG.AND.UMAX.LE.TOLG) THEN
5483 IF (DMAX.LE.TOLX) THEN
5486 IF (NTESX.GE.MTESX) RETURN
5490 TEMP=ABS(FO-F)/MAX(ABS(F),1.0D 0)
5491 IF (TEMP.LE.TOLF) THEN
5494 IF (NTESF.GE.MTESF) RETURN
5498 1 IF (NIT.GE.MIT) THEN
5502 IF (NFV.GE.MFV) THEN
5506 IF (NFG.GE.MFG) THEN
5511 IF (N.GT.0.AND.NIT-KIT.GE.IRES1*N+IRES2) THEN
5517 * SUBROUTINE PYFUT8 ALL SYSTEMS 98/12/01
5519 * TERMINATION CRITERIA AND TEST ON RESTART.
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.
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)
5574 IF (RPF3.GT.TOLP) GO TO 1
5576 IF (GMAX.LE.TOLG) THEN
5585 IF (DMAX.LE.TOLX) THEN
5588 IF (NTESX.GE.MTESX) RETURN
5592 TEMP=ABS(FO-F)/MAX(ABS(F),1.0D 0)
5593 IF (TEMP.LE.TOLF) THEN
5596 IF (NTESF.GE.MTESF) RETURN
5600 1 IF (NIT.GE.MIT) THEN
5604 IF (NFV.GE.MFV) THEN
5608 IF (NFG.GE.MFG) THEN
5613 IF (N.GT.0.AND.NIT-KIT.GE.IRES1*N+IRES2) THEN
5619 * SUBROUTINE PYPTSH ALL SYSTEMS 98/12/01
5621 * POWELL-TOINT GRAPH COLORING ALGORITHM FOR GROUPING COLUMNS OF THE
5622 * HESSIAN MATRIX BEFORE NUMERICAL DIFFERENTIATION.
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
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.
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.
5646 SUBROUTINE PYPTSH(NF,MMAX,IH,JH,COL,DEG,ORD,RADIX,WN11,WN12,XS,
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
5654 * SAVE SYMBOLIC STRUCTURE OF FACTOR
5657 IF (2*MM-NF+2.GE.MMAX) THEN
5662 * WIDTHEN THE STRUCTURE
5664 CALL MXSTG1(NF,ML,IH,JH,WN12,WN11)
5671 * NUMBER OF THE FREE COLUMNS
5675 * NUMBER OF USED COLOURS
5690 200 CALL MXVSR2(NF,DEG,ORD,RADIX,WN11,CNM)
5692 * ORD REWRITE INTO THE ARRAY INVP
5698 * COLUMNS OF THE NEW COLOUR NCOL
5700 CALL PYCSER(JH,IH,WN12,XS,DEG,WN11,COL,NCOL,CNM)
5709 * SHRINK THE STRUCTURE
5711 CALL MXSTL1(NF,ML,IH,JH,WN12)
5713 * INTO COL GIVE INDICES OF THE INDIVIDUAL GROUPS ONE AFTER ANOTHER,
5714 * END OF THE GROUP IS MARKED BY THE NEGATIVE INDEX VALUE.
5731 IF (WN11(I).EQ.0) GO TO 550
5748 IF (L.GT.NF) GO TO 900
5754 * SUBROUTINE PYRMC0 ALL SYSTEMS 98/12/01
5756 * OLD SIMPLE BOUND IS REMOVED FROM THE ACTIVE SET. TRANSFORMED
5757 * GRADIENT OF THE OBJECTIVE FUNCTION IS UPDATED.
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.
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
5775 IF (N.EQ.0.OR.RMAX.GT.0.0D 0) THEN
5776 IF (UMAX.GT.EPS8*GMAX) 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
5786 IX(I)=MIN(ABS(IX(I)),3)
5787 IF (RMAX.EQ.0) GO TO 2
5790 2 IF (IOLD.GT.1) IREST=MAX(IREST,1)
5795 * SUBROUTINE PYTCAB ALL SYSTEMS 06/12/01
5797 * VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED
5798 * AND SCALED. TEST VALUE DMAX IS DETERMINED.
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
5814 * SUBPROGRAMS USED :
5815 * S MXVDIF DIFFERENCE OF TWO VECTORS.
5816 * S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
5819 SUBROUTINE PYTCAB(NC,MC,CG,CGO,ICG,CZ,ITERS,JOB)
5820 INTEGER NC,MC,ICG(*),ITERS,JOB
5821 DOUBLE PRECISION CG(*),CGO(*),CZ(*)
5823 DOUBLE PRECISION TEMP
5824 IF (ITERS.GT.0) THEN
5825 CALL MXVDIF(MC,CG,CGO,CGO)
5827 CALL MXVSAV(MC,CG,CGO)
5834 IF (JOB.EQ.1) TEMP=SIGN(1.0D 0,TEMP)
5844 * SUBROUTINE PYTCUB ALL SYSTEMS 06/12/01
5846 * VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED
5847 * AND SCALED. TEST VALUE DMAX IS DETERMINED.
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
5865 * SUBPROGRAMS USED :
5866 * S MXVDIF DIFFERENCE OF TWO VECTORS.
5867 * S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
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)
5878 CALL MXVSAV(MC,CG,CGO)
5885 IF (KK.EQ.3.OR.KK.EQ.4) THEN
5886 TEMP= CZU(KC)-CZL(KC)
5887 ELSE IF (KK.EQ.1) THEN
5889 ELSE IF (KK.EQ.2) THEN
5891 ELSE IF (KK.EQ.5) THEN
5894 IF (JOB.EQ.1) TEMP=SIGN(1.0D 0,TEMP)
5904 * SUBROUTINE PYTRCD ALL SYSTEMS 98/12/01
5906 * VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED
5907 * AND SCALED AND REDUCED. TEST VALUE DMAX IS DETERMINED.
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.
5929 * SUBPROGRAMS USED :
5930 * S MXVDIF DIFFERENCE OF TWO VECTORS.
5931 * S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
5934 SUBROUTINE PYTRCD(NF,X,IX,XO,G,GO,R,F,FO,P,PO,DMAX,KBF,KD,LD,
5936 INTEGER NF,IX(*),KBF,KD,LD,ITERS
5937 DOUBLE PRECISION X(*),XO(*),G(*),GO(*),R,F,FO,P,PO,DMAX
5939 IF (ITERS.GT.0) THEN
5940 CALL MXVDIF(NF,X,XO,XO)
5941 CALL MXVDIF(NF,G,GO,GO)
5947 CALL MXVSAV(NF,X,XO)
5948 CALL MXVSAV(NF,G,GO)
5954 IF (IX(I).LT.0) THEN
5960 DMAX=MAX(DMAX,ABS(XO(I))/MAX(ABS(X(I)),1.0D 0))
5964 * SUBROUTINE PYTRCG ALL SYSTEMS 99/12/01
5966 * GRADIENT OF THE OBJECTIVE FUNCTION IS SCALED AND REDUCED. TEST VALUES
5967 * GMAX AND UMAX ARE COMPUTED.
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.
5980 * SUBPROGRAMS USED :
5981 * RF MXVMAX L-INFINITY NORM OF A VECTOR.
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
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
6013 * SUBROUTINE PYTRCS ALL SYSTEMS 98/12/01
6015 * SCALED AND REDUCED DIRECTION VECTOR IS BACK TRANSFORMED. VECTORS
6016 * X,G AND VALUES F,P ARE SAVED.
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.
6039 * SUBPROGRAMS USED :
6040 * S MXVCOP COPYING OF A VECTOR.
6042 SUBROUTINE PYTRCS(NF,X,IX,XO,XL,XU,G,GO,S,RO,FP,FO,F,PO,P,RMAX,
6044 INTEGER NF,IX(*),KBF
6045 DOUBLE PRECISION X(*),XO(*),XL(*),XU(*),G(*),GO(*),S(*),RO,FP,FO,
6052 CALL MXVCOP(NF,X,XO)
6053 CALL MXVCOP(NF,G,GO)
6056 IF (IX(I).LT.0) THEN
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))
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))
6070 * SUBROUTINE PYTSCH ALL SYSTEMS 99/12/01
6072 * HESSIAN MATRIX OF THE OBJECTIVE FUNCTION OR ITS APPROXIMATION
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.
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
6093 IF (IX(I).GE.0) THEN
6102 DO 2 J=JSTRT+1,JSTOP