5 #define max(a,b) ((a) > (b) ? (a) : (b))
6 #define min(a,b) ((a) < (b) ? (a) : (b))
7 #define iabs(a) ((a) < 0 ? -(a) : (a))
9 /* subroutines extracted from pssubs.for */
10 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
11 /* SUBROUTINE PCBS04 ALL SYSTEMS 98/12/01
13 * INITIATION OF THE VECTOR CONTAINING TYPES OF CONSTRAINTS.
16 * II NF NUMBER OF VARIABLES.
17 * RI X(NF) VECTOR OF VARIABLES.
18 * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS.
19 * RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
20 * RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
21 * RI EPS9 TOLERANCE FOR ACTIVE CONSTRAINTS.
22 * II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
23 * KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
25 void luksan_pcbs04__(int *nf, double *x, int *ix,
26 double *xl, double *xu, double *eps9, int *kbf)
28 /* System generated locals */
36 /* Parameter adjustments */
45 for (i__ = 1; i__ <= i__1; ++i__) {
47 ixi = (i__2 = ix[i__], iabs(i__2));
49 d__2 = (d__1 = xl[i__], fabs(d__1));
50 if ((ixi == 1 || ixi == 3 || ixi == 4) && x[i__] <= xl[i__] + *
51 eps9 * max(d__2,temp)) {
55 d__2 = (d__1 = xu[i__], fabs(d__1));
56 if ((ixi == 2 || ixi == 3 || ixi == 4) && x[i__] >= xu[i__] - *
57 eps9 * max(d__2,temp)) {
64 } /* luksan_pcbs04__ */
66 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
67 /* SUBROUTINE PNINT1 ALL SYSTEMS 91/12/01 */
69 /* EXTRAPOLATION OR INTERPOLATION FOR LINE SEARCH WITH DIRECTIONAL */
73 /* RI RL LOWER VALUE OF THE STEPSIZE PARAMETER. */
74 /* RI RU UPPER VALUE OF THE STEPSIZE PARAMETER. */
75 /* RI FL VALUE OF THE OBJECTIVE FUNCTION FOR R=RL. */
76 /* RI FU VALUE OF THE OBJECTIVE FUNCTION FOR R=RU. */
77 /* RI PL DIRECTIONAL DERIVATIVE FOR R=RL. */
78 /* RI PU DIRECTIONAL DERIVATIVE FOR R=RU. */
79 /* RO R VALUE OF THE STEPSIZE PARAMETER OBTAINED. */
80 /* II MODE MODE OF LINE SEARCH. */
81 /* II MTYP METHOD SELECTION. MTYP=1-BISECTION. MTYP=2-QUADRATIC */
82 /* INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE). */
83 /* MTYP=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL */
84 /* DERIVATIVES). MTYP=4-CUBIC INTERPOLATION. MTYP=5-CONIC */
86 /* IO MERR ERROR INDICATOR. MERR=0 FOR NORMAL RETURN. */
89 /* EXTRAPOLATION OR INTERPOLATION WITH STANDARD MODEL FUNCTIONS. */
91 void luksan_pnint1__(double *rl, double *ru, double *fl,
92 double *fu, double *pl, double *pu, double *r__,
93 int *mode, int *mtyp, int *merr)
95 /* System generated locals */
99 double a, b, c__, d__, den, dis;
109 } else if (*ru <= *rl) {
113 for (ntyp = *mtyp; ntyp >= 1; --ntyp) {
122 *r__ = (*rl + *ru) * .5;
125 } else if (ntyp == *mtyp) {
126 a = (*fu - *fl) / (*pl * (*ru - *rl));
131 /* QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH ONE DIRECTIONAL */
135 } else if (ntyp == 3) {
137 /* QUADRATIC EXTRAPOLATION OR INTERPOLATION WITH TWO DIRECTIONAL */
141 } else if (ntyp == 4) {
143 /* CUBIC EXTRAPOLATION OR INTERPOLATION */
145 c__ = b - a * 2. + 1.;
146 d__ = b - a * 3. + 2.;
147 dis = d__ * d__ - c__ * 3.;
151 den = d__ + sqrt(dis);
152 } else if (ntyp == 5) {
154 /* CONIC EXTRAPOLATION OR INTERPOLATION */
164 /* Computing 3rd power */
166 den = 1. - b * (d__1 * (d__1 * d__1));
168 if (*mode == 1 && den > 0. && den < 1.) {
170 /* EXTRAPOLATION ACCEPTED */
172 *r__ = *rl + (*ru - *rl) / den;
174 d__1 = *r__, d__2 = *ru * 1.1;
175 *r__ = max(d__1,d__2);
177 d__1 = *r__, d__2 = *ru * 1e3;
178 *r__ = min(d__1,d__2);
180 } else if (*mode == 2 && den > 1.) {
182 /* INTERPOLATION ACCEPTED */
184 *r__ = *rl + (*ru - *rl) / den;
187 d__1 = *r__, d__2 = *rl + (*ru - *rl) * .01;
188 *r__ = max(d__1,d__2);
191 d__1 = *r__, d__2 = *rl + (*ru - *rl) * .1;
192 *r__ = max(d__1,d__2);
195 d__1 = *r__, d__2 = *rl + (*ru - *rl) * .9;
196 *r__ = min(d__1,d__2);
203 } /* luksan_pnint1__ */
205 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
206 /* SUBROUTINE PS1L01 ALL SYSTEMS 97/12/01
208 * STANDARD LINE SEARCH WITH DIRECTIONAL DERIVATIVES.
211 * RO R VALUE OF THE STEPSIZE PARAMETER.
212 * RO RP PREVIOUS VALUE OF THE STEPSIZE PARAMETER.
213 * RO F VALUE OF THE OBJECTIVE FUNCTION.
214 * RI FO INITIAL VALUE OF THE OBJECTIVE FUNCTION.
215 * RO FP PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
216 * RO P VALUE OF THE DIRECTIONAL DERIVATIVE.
217 * RI PO INITIAL VALUE OF THE DIRECTIONAL DERIVATIVE.
218 * RO PP PREVIOUS VALUE OF THE DIRECTIONAL DERIVATIVE.
219 * RI FMIN LOWER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
220 * RI MAXF UPPER BOUND FOR VALUE OF THE OBJECTIVE FUNCTION.
221 * RI RMIN MINIMUM VALUE OF THE STEPSIZE PARAMETER
222 * RI RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER
223 * RI TOLS TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE
224 * CHANGE OF THE FUNCTION VALUE).
225 * RI TOLP TERMINATION TOLERANCE FOR LINE SEARCH (IN TEST ON THE
226 * CHANGE OF THE DIRECTIONAL DERIVATIVE).
227 * RO PAR1 PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC
229 * RO PAR2 PARAMETER FOR CONTROLLED SCALING OF VARIABLE METRIC
231 * II KD DEGREE OF REQUIRED DERIVATIVES.
232 * IO LD DEGREE OF PREVIOUSLY COMPUTED DERIVATIVES OF OBJECTIVE
233 * II NIT ACTUAL NUMBER OF ITERATIONS.
234 * II KIT NUMBER OF THE ITERATION AFTER LAST RESTART.
235 * IO NRED ACTUAL NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
236 * II MRED MAXIMUM NUMBER OF EXTRAPOLATIONS OR INTERPOLATIONS.
237 * IO MAXST MAXIMUM STEPSIZE INDICATOR. MAXST=0 OR MAXST=1 IF MAXIMUM
238 * STEPSIZE WAS NOT OR WAS REACHED.
239 * II IEST LOWER BOUND SPECIFICATION. IEST=0 OR IEST=1 IF LOWER BOUND
240 * IS NOT OR IS GIVEN.
241 * II INITS CHOICE OF THE INITIAL STEPSIZE. INITS=0-INITIAL STEPSIZE
242 * IS SPECIFIED IN THE CALLING PROGRAM. INITS=1-UNIT INITIAL
243 * STEPSIZE. INITS=2-COMBINED UNIT AND QUADRATICALLY ESTIMATED
244 * INITIAL STEPSIZE. INITS=3-QUADRATICALLY ESTIMATED INITIAL
246 * IO ITERS TERMINATION INDICATOR. ITERS=0-ZERO STEP. ITERS=1-PERFECT
247 * LINE SEARCH. ITERS=2 GOLDSTEIN STEPSIZE. ITERS=3-CURRY
248 * STEPSIZE. ITERS=4-EXTENDED CURRY STEPSIZE.
249 * ITERS=5-ARMIJO STEPSIZE. ITERS=6-FIRST STEPSIZE.
250 * ITERS=7-MAXIMUM STEPSIZE. ITERS=8-UNBOUNDED FUNCTION.
251 * ITERS=-1-MRED REACHED. ITERS=-2-POSITIVE DIRECTIONAL
252 * DERIVATIVE. ITERS=-3-ERROR IN INTERPOLATION.
253 * II KTERS TERMINATION SELECTION. KTERS=1-PERFECT LINE SEARCH.
254 * KTERS=2-GOLDSTEIN STEPSIZE. KTERS=3-CURRY STEPSIZE.
255 * KTERS=4-EXTENDED CURRY STEPSIZE. KTERS=5-ARMIJO STEPSIZE.
256 * KTERS=6-FIRST STEPSIZE.
257 * II MES METHOD SELECTION. MES=1-BISECTION. MES=2-QUADRATIC
258 * INTERPOLATION (WITH ONE DIRECTIONAL DERIVATIVE).
259 * MES=3-QUADRATIC INTERPOLATION (WITH TWO DIRECTIONAL
260 * DERIVATIVES). MES=4-CUBIC INTERPOLATION. MES=5-CONIC
262 * IU ISYS CONTROL PARAMETER.
265 * S PNINT1 EXTRAPOLATION OR INTERPOLATION WITH DIRECTIONAL
269 * SAFEGUARDED EXTRAPOLATION AND INTERPOLATION WITH STANDARD TERMINATION
272 void luksan_ps1l01__(double *r__, double *rp,
273 double *f, double *fo, double *fp, double *p,
274 double *po, double *pp, double *minf, double *maxf,
275 double *rmin, double *rmax, double *tols, double *
276 tolp, double *par1, double *par2, int *kd, int *ld,
277 int *nit, int *kit, int *nred, int *mred, int *
278 maxst, int *iest, int *inits, int *iters, int *kters,
281 /* System generated locals */
284 /* Local variables */
285 unsigned l1, l2, l3, m1, l5, m2, l7, m3;
286 static double fl, fu, pl, rl, pu, ru;
287 static int mes1, mes2, mes3, mode;
310 /* INITIAL STEPSIZE SELECTION */
314 } else if (*iest == 0) {
318 d__1 = *f - *fp, d__2 = *minf - *f;
319 rtemp = max(d__1,d__2);
321 init1 = iabs(*inits);
326 } else if (init1 == 1 || (*inits >= 1 && *iest == 0)) {
328 } else if (init1 == 2) {
330 d__1 = 1., d__2 = rtemp * 4. / *po;
331 *r__ = min(d__1,d__2);
332 } else if (init1 == 3) {
334 d__1 = 1., d__2 = rtemp * 2. / *po;
335 *r__ = min(d__1,d__2);
336 } else if (init1 == 4) {
337 *r__ = rtemp * 2. / *po;
339 *r__ = max(*r__,*rmin);
340 *r__ = min(*r__,*rmax);
346 /* NEW STEPSIZE SELECTION (EXTRAPOLATION OR INTERPOLATION) */
349 luksan_pnint1__(&rl, &ru, &fl, &fu, &pl, &pu, r__, &mode, &mtyp, &merr);
353 } else if (mode == 1) {
355 *r__ = min(*r__,*rmax);
356 } else if (mode == 2) {
360 /* COMPUTATION OF THE NEW FUNCTION VALUE AND THE NEW DIRECTIONAL */
379 l1 = *r__ <= *rmin && *nit != *kit;
381 l3 = *f - *fo <= *tols * *r__ * *po;
382 l5 = *p >= *tolp * *po || (mes2 == 2 && mode == 2);
383 l7 = mes2 <= 2 || mode != 0;
388 m1 = fabs(*p) <= fabs(*po) * .01 && *fo - *f >= fabs(*fo) *
389 9.9999999999999994e-12;
393 m2 = fabs(*p) <= fabs(*po) * .5 && (d__1 = *fo - *f, fabs(d__1)) <=
394 fabs(*fo) * 2.0000000000000001e-13;
403 /* TEST ON TERMINATION */
408 } else if (l2 && l3 && ! l5) {
411 } else if (m3 && mes1 == 3) {
414 } else if (l3 && l5 && l7) {
417 } else if (*kters < 0 || (*kters == 6 && l7)) {
420 } else if (iabs(*nred) >= *mred) {
435 /* INTERVAL CHANGE AFTER EXTRAPOLATION */
446 } else if (mes1 == 1) {
451 /* INTERVAL CHANGE AFTER INTERPOLATION */
467 } /* luksan_ps1l01__ */
469 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
470 /* SUBROUTINE PULSP3 ALL SYSTEMS 02/12/01
472 * LIMITED STORAGE VARIABLE METRIC UPDATE.
475 * II N NUMBER OF VARIABLES (NUMBER OF ROWS OF XM).
476 * II M NUMBER OF COLUMNS OF XM.
477 * II MF MAXIMUM NUMBER OF COLUMNS OF XM.
478 * RI XM(N*M) RECTANGULAR MATRIX IN THE PRODUCT FORM SHIFTED BROYDEN
479 * METHOD (STORED COLUMNWISE): H-SIGMA*I=XM*TRANS(XM)
480 * RO GR(M) MATRIX TRANS(XM)*GO.
481 * RU XO(N) VECTORS OF VARIABLES DIFFERENCE XO AND VECTOR XO-TILDE.
482 * RU GO(N) GRADIENT DIFFERENCE GO AND VECTOR XM*TRANS(XM)*GO.
483 * RI R STEPSIZE PARAMETER.
484 * RI PO OLD DIRECTIONAL DERIVATIVE (MULTIPLIED BY R)
485 * RU SIG SCALING PARAMETER (ZETA AND SIGMA).
486 * IO ITERH TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION.
487 * ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS.
488 * II MET3 CHOICE OF SIGMA (1-CONSTANT, 2-QUADRATIC EQUATION).
491 * S MXDRMM MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR
492 * MATRIX BY A VECTOR.
493 * S MXDCMU UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX.
494 * WITH CONTROLLING OF POSITIVE DEFINITENESS.
495 * S MXVDIR VECTOR AUGMENTED BY A SCALED VECTOR.
496 * RF MXVDOT DOT PRODUCT OF VECTORS.
497 * S MXVSCL SCALING OF A VECTOR.
500 * SHIFTED BFGS METHOD IN THE PRODUCT FORM.
502 void luksan_pulsp3__(int *n, int *m, int *mf,
503 double *xm, double *gr, double *xo, double *go,
504 double *r__, double *po, double *sig, int *iterh,
507 /* System generated locals */
508 double d__1, d__2, d__3, d__4;
510 /* Builtin functions */
512 /* Local variables */
513 double a, b, c__, aa, bb, ah, den, par, pom;
515 /* Parameter adjustments */
525 b = luksan_mxvdot__(n, &xo[1], &go[1]);
530 luksan_mxdrmm__(n, m, &xm[1], &go[1], &gr[1]);
531 ah = luksan_mxvdot__(n, &go[1], &go[1]);
532 aa = luksan_mxvdot__(m, &gr[1], &gr[1]);
536 /* DETERMINATION OF THE PARAMETER SIG (SHIFT) */
541 den = luksan_mxvdot__(n, &xo[1], &xo[1]);
544 d__1 = 0., d__2 = 1. - aa / a;
546 d__3 = 0., d__4 = 1. - b * b / (den * ah);
547 *sig = sqrt((max(d__1,d__2))) / (sqrt((max(d__3,d__4))) + 1.) *
551 d__1 = 0., d__2 = *sig * ah / a;
553 d__3 = 0., d__4 = 1. - b * b / (den * ah);
554 *sig = sqrt((max(d__1,d__2))) / (sqrt((max(d__3,d__4))) + 1.) *
558 d__1 = *sig, d__2 = pom * .2;
559 *sig = max(d__1,d__2);
561 d__1 = *sig, d__2 = pom * .8;
562 *sig = min(d__1,d__2);
567 /* COMPUTATION OF SHIFTED XO AND SHIFTED B */
571 luksan_mxvdir__(n, &d__1, &go[1], &xo[1], &xo[1]);
573 /* BFGS-BASED SHIFTED BFGS UPDATE */
577 luksan_mxdcmu__(n, m, &xm[1], &d__1, &xo[1], &gr[1]);
578 d__1 = sqrt(par / bb);
579 luksan_mxvscl__(n, &d__1, &xo[1], &xm[*n * *m + 1]);
584 } /* luksan_pulsp3__ */
586 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
587 /* SUBROUTINE PULVP3 ALL SYSTEMS 03/12/01
589 * RANK-TWO LIMITED-STORAGE VARIABLE-METRIC METHODS IN THE PRODUCT FORM.
592 * II N NUMBER OF VARIABLES (NUMBER OF ROWS OF XM).
593 * II M NUMBER OF COLUMNS OF XM.
594 * RI XM(N*M) RECTANGULAR MATRIX IN THE PRODUCT FORM SHIFTED BROYDEN
595 * METHOD (STORED COLUMNWISE): H-SIGMA*I=XM*TRANS(XM)
596 * RO XR(M) VECTOR TRANS(XM)*H**(-1)*XO.
597 * RO GR(M) MATRIX TRANS(XM)*GO.
598 * RA S(N) AUXILIARY VECTORS (H**(-1)*XO AND U).
599 * RA SO(N) AUXILIARY VECTORS ((H-SIGMA*I)*H**(-1)*XO AND V).
600 * RU XO(N) VECTORS OF VARIABLES DIFFERENCE XO AND VECTOR XO-TILDE.
601 * RU GO(N) GRADIENT DIFFERENCE GO AND VECTOR XM*TRANS(XM)*GO.
602 * RI R STEPSIZE PARAMETER.
603 * RI PO OLD DIRECTIONAL DERIVATIVE (MULTIPLIED BY R)
604 * RU SIG SCALING PARAMETER (ZETA AND SIGMA).
605 * IO ITERH TERMINATION INDICATOR. ITERH<0-BAD DECOMPOSITION.
606 * ITERH=0-SUCCESSFUL UPDATE. ITERH>0-NONPOSITIVE PARAMETERS.
607 * II MET2 CHOICE OF THE CORRECTION PARAMETER (1-THE UNIT VALUE,
608 * 2-THE BALANCING VALUE, 3-THE SQUARE ROOT, 4-THE GEOMETRIC
610 * II MET3 CHOICE OF THE SHIFT PARAMETER (4-THE FIRST FORMULA,
611 * 5-THE SECOND FORMULA).
612 * II MET5 CHOICE OF THE METHOD (1-RANK-ONE METHOD, 2-RANK-TWO
616 * S MXDRMM MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR
617 * MATRIX BY A VECTOR.
618 * S MXDCMU UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX.
619 * WITH CONTROLLING OF POSITIVE DEFINITENESS. RANK-ONE FORMULA.
620 * S MXDCMV UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX.
621 * WITH CONTROLLING OF POSITIVE DEFINITENESS. RANK-TWO FORMULA.
622 * S MXVDIR VECTOR AUGMENTED BY A SCALED VECTOR.
623 * RF MXVDOT DOT PRODUCT OF VECTORS.
624 * S MXVLIN LINEAR COMBINATION OF TWO VECTORS.
625 * S MXVSCL SCALING OF A VECTOR.
628 * RANK-ONE LIMITED-STORAGE VARIABLE-METRIC METHOD IN THE PRODUCT FORM.
630 void luksan_pulvp3__(int *n, int *m, double *xm,
631 double *xr, double *gr, double *s, double *so,
632 double *xo, double *go, double *r__, double *po,
633 double *sig, int *iterh, int *met2, int *met3,
636 /* System generated locals */
637 double d__1, d__2, d__3, d__4;
639 /* Builtin functions */
641 /* Local variables */
642 double a, b, c__, aa, bb, cc, ah, den, par, pom, zet;
644 /* Parameter adjustments */
656 /* COMPUTATION OF B */
658 b = luksan_mxvdot__(n, &xo[1], &go[1]);
664 /* COMPUTATION OF GR=TRANS(XM)*GO, XR=TRANS(XM)*H**(-1)*XO */
665 /* AND S=H**(-1)*XO, SO=(H-SIGMA*I)*H**(-1)*XO. COMPUTATION */
666 /* OF AA=GR*GR, BB=GR*XR, CC=XR*XR. COMPUTATION OF A AND C. */
668 luksan_mxdrmm__(n, m, &xm[1], &go[1], &gr[1]);
669 luksan_mxvscl__(n, r__, &s[1], &s[1]);
670 luksan_mxdrmm__(n, m, &xm[1], &s[1], &xr[1]);
672 luksan_mxvdir__(n, &d__1, &s[1], &xo[1], &so[1]);
673 ah = luksan_mxvdot__(n, &go[1], &go[1]);
674 aa = luksan_mxvdot__(m, &gr[1], &gr[1]);
675 bb = luksan_mxvdot__(m, &gr[1], &xr[1]);
676 cc = luksan_mxvdot__(m, &xr[1], &xr[1]);
680 /* DETERMINATION OF THE PARAMETER SIG (SHIFT) */
684 den = luksan_mxvdot__(n, &xo[1], &xo[1]);
687 d__1 = 0., d__2 = 1. - aa / a;
689 d__3 = 0., d__4 = 1. - b * b / (den * ah);
690 *sig = sqrt((max(d__1,d__2))) / (sqrt((max(d__3,d__4))) + 1.) *
694 d__1 = 0., d__2 = *sig * ah / a;
696 d__3 = 0., d__4 = 1. - b * b / (den * ah);
697 *sig = sqrt((max(d__1,d__2))) / (sqrt((max(d__3,d__4))) + 1.) *
701 d__1 = *sig, d__2 = pom * .2;
702 *sig = max(d__1,d__2);
704 d__1 = *sig, d__2 = pom * .8;
705 *sig = min(d__1,d__2);
710 /* COMPUTATION OF SHIFTED XO AND SHIFTED B */
714 luksan_mxvdir__(n, &d__1, &go[1], &xo[1], &xo[1]);
716 /* COMPUTATION OF THE PARAMETER RHO (CORRECTION) */
720 } else if (*met2 == 2) {
722 } else if (*met2 == 3) {
723 par = sqrt(1. - aa / a);
724 } else if (*met2 == 4) {
725 par = sqrt(sqrt(1. - aa / a) * (*sig * ah / b));
727 par = zet / (zet + *sig);
730 /* COMPUTATION OF THE PARAMETER THETA (BFGS) */
732 d__1 = sqrt(par * b / cc);
733 pom = copysign(d__1, bb);
735 /* COMPUTATION OF Q AND P */
739 /* RANK ONE UPDATE OF XM */
741 luksan_mxvdir__(m, &pom, &xr[1], &gr[1], &xr[1]);
742 luksan_mxvlin__(n, &par, &xo[1], &pom, &so[1], &s[1]);
743 d__1 = -1. / (par * b + pom * bb);
744 luksan_mxdcmu__(n, m, &xm[1], &d__1, &s[1], &xr[1]);
747 /* RANK TWO UPDATE OF XM */
749 d__1 = par / pom - bb / b;
750 luksan_mxvdir__(n, &d__1, &xo[1], &so[1], &s[1]);
753 luksan_mxdcmv__(n, m, &xm[1], &d__1, &xo[1], &gr[1], &d__2, &s[1], &xr[1]);
758 } /* luksan_pulvp3__ */
760 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
761 /* SUBROUTINE PYADC0 ALL SYSTEMS 98/12/01
763 * NEW SIMPLE BOUNDS ARE ADDED TO THE ACTIVE SET.
766 * II NF DECLARED NUMBER OF VARIABLES.
767 * II N REDUCED NUMBER OF VARIABLES.
768 * RI X(NF) VECTOR OF VARIABLES.
769 * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS.
770 * RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
771 * RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
772 * IO INEW NUMBER OF ACTIVE CONSTRAINTS.
774 void luksan_pyadc0__(int *nf, int *n, double *x,
775 int *ix, double *xl, double *xu, int *inew)
777 /* System generated locals */
780 /* Local variables */
783 /* Parameter adjustments */
793 for (i__ = 1; i__ <= i__1; ++i__) {
798 } else if ((ixi == 1 || ixi == 3 || ixi == 4) && x[i__] <= xl[i__]) {
809 } else if ((ixi == 2 || ixi == 3 || ixi == 4) && x[i__] >= xu[i__]) {
824 } /* luksan_pyadc0__ */
826 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
827 /* SUBROUTINE PYFUT1 ALL SYSTEMS 98/12/01
829 * TERMINATION CRITERIA AND TEST ON RESTART.
832 * II N ACTUAL NUMBER OF VARIABLES.
833 * RI F NEW VALUE OF THE OBJECTIVE FUNCTION.
834 * RI FO OLD VALUE OF THE OBJECTIVE FUNCTION.
835 * RI UMAX MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
836 * RO GMAX NORM OF THE TRANSFORMED GRADIENT.
837 * RI DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
838 * RI TOLX LOWER BOUND FOR STEPLENGTH.
839 * RI TOLF LOWER BOUND FOR FUNCTION DECREASE.
840 * RI TOLB LOWER BOUND FOR FUNCTION VALUE.
841 * RI TOLG LOWER BOUND FOR GRADIENT.
842 * II KD DEGREE OF REQUIRED DERIVATIVES.
843 * IU NIT ACTUAL NUMBER OF ITERATIONS.
844 * II KIT NUMBER OF THE ITERATION AFTER RESTART.
845 * II MIT MAXIMUM NUMBER OF ITERATIONS.
846 * IU NFV ACTUAL NUMBER OF COMPUTED FUNCTION VALUES.
847 * II MFV MAXIMUM NUMBER OF COMPUTED FUNCTION VALUES.
848 * IU NFG ACTUAL NUMBER OF COMPUTED GRADIENT VALUES.
849 * II MFG MAXIMUM NUMBER OF COMPUTED GRADIENT VALUES.
850 * IU NTESX ACTUAL NUMBER OF TESTS ON STEPLENGTH.
851 * II MTESX MAXIMUM NUMBER OF TESTS ON STEPLENGTH.
852 * IU NTESF ACTUAL NUMBER OF TESTS ON FUNCTION DECREASE.
853 * II MTESF MAXIMUM NUMBER OF TESTS ON FUNCTION DECREASE.
854 * II IRES1 RESTART SPECIFICATION. RESTART IS PERFORMED AFTER
855 * IRES1*N+IRES2 ITERATIONS.
856 * II IRES2 RESTART SPECIFICATION. RESTART IS PERFORMED AFTER
857 * IRES1*N+IRES2 ITERATIONS.
858 * IU IREST RESTART INDICATOR. RESTART IS PERFORMED IF IREST>0.
859 * II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
860 * ITERS=0 FOR ZERO STEP.
861 * IO ITERM TERMINATION INDICATOR. ITERM=1-TERMINATION AFTER MTESX
862 * UNSUFFICIENT STEPLENGTHS. ITERM=2-TERMINATION AFTER MTESF
863 * UNSUFFICIENT FUNCTION DECREASES. ITERM=3-TERMINATION ON LOWER
864 * BOUND FOR FUNCTION VALUE. ITERM=4-TERMINATION ON LOWER BOUND
865 * FOR GRADIENT. ITERM=11-TERMINATION AFTER MAXIMUM NUMBER OF
866 * ITERATIONS. ITERM=12-TERMINATION AFTER MAXIMUM NUMBER OF
867 * COMPUTED FUNCTION VALUES.
869 void luksan_pyfut1__(int *n, double *f, double *
870 fo, double *umax, double *gmax, double *dmax__,
871 double *tolx, double *tolf, double *tolb, double *
872 tolg, int *kd, int *nit, int *kit, int *mit, int *
873 nfv, int *mfv, int *nfg, int *mfg, int *ntesx,
874 int *mtesx, int *ntesf, int *mtesf, int *ites,
875 int *ires1, int *ires2, int *irest, int *iters,
878 /* System generated locals */
881 /* Builtin functions */
883 /* Local variables */
897 d__1 = sqrt((fabs(*f))), d__2 = fabs(*f) / 10.;
898 *fo = *f + min(d__1,d__2);
905 if (*gmax <= *tolg && *umax <= *tolg) {
914 if (*dmax__ <= *tolx) {
917 if (*ntesx >= *mtesx) {
925 temp = (d__1 = *fo - *f, fabs(d__1)) / max(d__2,1.);
929 if (*ntesf >= *mtesf) {
949 if (*n > 0 && *nit - *kit >= *ires1 * *n + *ires2) {
950 *irest = max(*irest,1);
954 } /* luksan_pyfut1__ */
956 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
957 /* SUBROUTINE PYRMC0 ALL SYSTEMS 98/12/01
959 * OLD SIMPLE BOUND IS REMOVED FROM THE ACTIVE SET. TRANSFORMED
960 * GRADIENT OF THE OBJECTIVE FUNCTION IS UPDATED.
963 * II NF DECLARED NUMBER OF VARIABLES.
964 * II N REDUCED NUMBER OF VARIABLES.
965 * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS.
966 * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION.
967 * RI EPS8 TOLERANCE FOR CONSTRAINT TO BE REMOVED.
968 * RI UMAX MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
969 * RI GMAX NORM OF THE TRANSFORMED GRADIENT.
970 * RO RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER.
971 * II IOLD NUMBER OF REMOVED CONSTRAINTS.
972 * IU IREST RESTART INDICATOR.
974 void luksan_pyrmc0__(int *nf, int *n, int *ix,
975 double *g, double *eps8, double *umax, double *gmax,
976 double *rmax, int *iold, int *irest)
978 /* System generated locals */
979 int i__1, i__2, i__3;
981 /* Local variables */
984 /* Parameter adjustments */
989 if (*n == 0 || *rmax > 0.) {
990 if (*umax > *eps8 * *gmax) {
993 for (i__ = 1; i__ <= i__1; ++i__) {
996 } else if (ixi <= -5) {
997 } else if ((ixi == -1 || ixi == -3) && -g[i__] <= 0.) {
998 } else if ((ixi == -2 || ixi == -4) && g[i__] <= 0.) {
1002 i__3 = (i__2 = ix[i__], iabs(i__2));
1003 ix[i__] = min(i__3,3);
1012 *irest = max(*irest,1);
1017 } /* luksan_pyrmc0__ */
1019 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
1020 /* SUBROUTINE PYTRCD ALL SYSTEMS 98/12/01
1022 * VECTORS OF VARIABLES DIFFERENCE AND GRADIENTS DIFFERENCE ARE COMPUTED
1023 * AND SCALED AND REDUCED. TEST VALUE DMAX IS DETERMINED.
1026 * II NF DECLARED NUMBER OF VARIABLES.
1027 * RI X(NF) VECTOR OF VARIABLES.
1028 * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS.
1029 * RU XO(NF) VECTORS OF VARIABLES DIFFERENCE.
1030 * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION.
1031 * RU GO(NF) GRADIENTS DIFFERENCE.
1032 * RO R VALUE OF THE STEPSIZE PARAMETER.
1033 * RO F NEW VALUE OF THE OBJECTIVE FUNCTION.
1034 * RI FO OLD VALUE OF THE OBJECTIVE FUNCTION.
1035 * RO P NEW VALUE OF THE DIRECTIONAL DERIVATIVE.
1036 * RI PO OLD VALUE OF THE DIRECTIONAL DERIVATIVE.
1037 * RO DMAX MAXIMUM RELATIVE DIFFERENCE OF VARIABLES.
1038 * II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
1039 * KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
1040 * IO KD DEGREE OF REQUIRED DERIVATIVES.
1041 * IO LD DEGREE OF COMPUTED DERIVATIVES.
1042 * II ITERS TERMINATION INDICATOR FOR STEPLENGTH DETERMINATION.
1043 * ITERS=0 FOR ZERO STEP.
1045 * SUBPROGRAMS USED :
1046 * S MXVDIF DIFFERENCE OF TWO VECTORS.
1047 * S MXVSAV DIFFERENCE OF TWO VECTORS WITH COPYING AND SAVING THE
1050 void luksan_pytrcd__(int *nf, double *x, int *ix,
1051 double *xo, double *g, double *go, double *r__,
1052 double *f, double *fo, double *p, double *po,
1053 double *dmax__, int *kbf, int *kd, int *ld, int *
1056 /* System generated locals */
1058 double d__1, d__2, d__3, d__4, d__5;
1060 /* Local variables */
1063 /* Parameter adjustments */
1072 luksan_mxvdif__(nf, &x[1], &xo[1], &xo[1]);
1073 luksan_mxvdif__(nf, &g[1], &go[1], &go[1]);
1079 luksan_mxvsav__(nf, &x[1], &xo[1]);
1080 luksan_mxvsav__(nf, &g[1], &go[1]);
1085 for (i__ = 1; i__ <= i__1; ++i__) {
1095 d__5 = (d__2 = x[i__], fabs(d__2));
1096 d__3 = *dmax__, d__4 = (d__1 = xo[i__], fabs(d__1)) / max(d__5,1.);
1097 *dmax__ = max(d__3,d__4);
1102 } /* luksan_pytrcd__ */
1104 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
1105 /* SUBROUTINE PYTRCG ALL SYSTEMS 99/12/01
1107 * GRADIENT OF THE OBJECTIVE FUNCTION IS SCALED AND REDUCED. TEST VALUES
1108 * GMAX AND UMAX ARE COMPUTED.
1111 * II NF DECLARED NUMBER OF VARIABLES.
1112 * II N ACTUAL NUMBER OF VARIABLES.
1113 * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS.
1114 * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION.
1115 * RI UMAX MAXIMUM ABSOLUTE VALUE OF THE NEGATIVE LAGRANGE MULTIPLIER.
1116 * RI GMAX NORM OF THE TRANSFORMED GRADIENT.
1117 * II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
1118 * KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
1119 * II IOLD INDEX OF THE REMOVED CONSTRAINT.
1121 * SUBPROGRAMS USED :
1122 * RF MXVMAX L-INFINITY NORM OF A VECTOR.
1124 void luksan_pytrcg__(int *nf, int *n, int *ix,
1125 double *g, double *umax, double *gmax, int *kbf,
1128 /* System generated locals */
1132 /* Local variables */
1136 /* Parameter adjustments */
1146 for (i__ = 1; i__ <= i__1; ++i__) {
1150 d__1 = *gmax, d__2 = fabs(temp);
1151 *gmax = max(d__1,d__2);
1152 } else if (ix[i__] <= -5) {
1153 } else if ((ix[i__] == -1 || ix[i__] == -3) && *umax + temp >= 0.)
1155 } else if ((ix[i__] == -2 || ix[i__] == -4) && *umax - temp >= 0.)
1165 *gmax = luksan_mxvmax__(nf, &g[1]);
1169 } /* luksan_pytrcg__ */
1171 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
1172 /* SUBROUTINE PYTRCS ALL SYSTEMS 98/12/01
1174 * SCALED AND REDUCED DIRECTION VECTOR IS BACK TRANSFORMED. VECTORS
1175 * X,G AND VALUES F,P ARE SAVED.
1178 * II NF DECLARED NUMBER OF VARIABLES.
1179 * RI X(NF) VECTOR OF VARIABLES.
1180 * II IX(NF) VECTOR CONTAINING TYPES OF BOUNDS.
1181 * RO XO(NF) SAVED VECTOR OF VARIABLES.
1182 * RI XL(NF) VECTOR CONTAINING LOWER BOUNDS FOR VARIABLES.
1183 * RI XU(NF) VECTOR CONTAINING UPPER BOUNDS FOR VARIABLES.
1184 * RI G(NF) GRADIENT OF THE OBJECTIVE FUNCTION.
1185 * RO GO(NF) SAVED GRADIENT OF THE OBJECTIVE FUNCTION.
1186 * RO S(NF) DIRECTION VECTOR.
1187 * RO RO SAVED VALUE OF THE STEPSIZE PARAMETER.
1188 * RO FP PREVIOUS VALUE OF THE OBJECTIVE FUNCTION.
1189 * RU FO SAVED VALUE OF THE OBJECTIVE FUNCTION.
1190 * RI F VALUE OF THE OBJECTIVE FUNCTION.
1191 * RO PO SAVED VALUE OF THE DIRECTIONAL DERIVATIVE.
1192 * RI P VALUE OF THE DIRECTIONAL DERIVATIVE.
1193 * RO RMAX MAXIMUM VALUE OF THE STEPSIZE PARAMETER.
1194 * RI ETA9 MAXIMUM FOR REAL NUMBERS.
1195 * II KBF SPECIFICATION OF SIMPLE BOUNDS. KBF=0-NO SIMPLE BOUNDS.
1196 * KBF=1-ONE SIDED SIMPLE BOUNDS. KBF=2=TWO SIDED SIMPLE BOUNDS.
1198 * SUBPROGRAMS USED :
1199 * S MXVCOP COPYING OF A VECTOR.
1201 void luksan_pytrcs__(int *nf, double *x, int *ix,
1202 double *xo, double *xl, double *xu, double *g,
1203 double *go, double *s, double *ro, double *fp,
1204 double *fo, double *f, double *po, double *p,
1205 double *rmax, double *eta9, int *kbf)
1207 /* System generated locals */
1211 /* Local variables */
1214 /* Parameter adjustments */
1229 luksan_mxvcop__(nf, &x[1], &xo[1]);
1230 luksan_mxvcop__(nf, &g[1], &go[1]);
1233 for (i__ = 1; i__ <= i__1; ++i__) {
1237 if (ix[i__] == 1 || ix[i__] >= 3) {
1238 if (s[i__] < -1. / *eta9) {
1240 d__1 = *rmax, d__2 = (xl[i__] - x[i__]) / s[i__];
1241 *rmax = min(d__1,d__2);
1244 if (ix[i__] == 2 || ix[i__] >= 3) {
1245 if (s[i__] > 1. / *eta9) {
1247 d__1 = *rmax, d__2 = (xu[i__] - x[i__]) / s[i__];
1248 *rmax = min(d__1,d__2);
1256 } /* luksan_pytrcs__ */