/* xgnew.f -- translated by f2c (version 20000817). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ #include "f2c.h" #include "blaswrap.h" /* Common Block Declarations */ struct { integer iout, iprint; doublereal mcheps, cnstol; } dfocm_; #define dfocm_1 dfocm_ struct { doublereal gmod[200], hmod[40000] /* was [200][200] */; } mdlpar_; #define mdlpar_1 mdlpar_ /* Table of constant values */ static integer c__1 = 1; /* Copyright (C) 2000, International Business Machines */ /* Corporation and others. All Rights Reserved. */ /* Subroutine */ int ptnew_(xgnew, ipoly, vnew, pntint, lptint, n, poly, lpoly, nind, base, delta, lb, ub, a, lda, ldcj, nclin, ncnln, pivval, pivthr, xchthr, wrk, lwrk, iwrk, liwrk, fail, info, minmax, x0) doublereal *xgnew; integer *ipoly; doublereal *vnew, *pntint; integer *lptint, *n; doublereal *poly; integer *lpoly, *nind, *base; doublereal *delta, *lb, *ub, *a; integer *lda, *ldcj, *nclin, *ncnln; doublereal *pivval, *pivthr, *xchthr, *wrk; integer *lwrk, *iwrk, *liwrk; logical *fail; integer *info, *minmax; doublereal *x0; { /* Format strings */ static char fmt_1000[] = "(\002 PTNEW: *** ERROR: LWRK TOO SMALL!\002\ /\002 IT SHOULD BE AT LEAST \002,i12,/)"; static char fmt_8000[] = "(\002 PTNEW: A new point is found by maximizin\ g the pivot,\002,/\002 polynomial: \002,i4,\002 pivot value: \002,d14.\ 7,/)"; static char fmt_8020[] = "(\002 PTNEW: A new point replaced the \002,i4\ ,\002-th point\002,/,\002 old pivot: \002,d14.7,\002 new pivot: \002,d\ 14.7,/)"; static char fmt_8030[] = "(\002 PTNEW: A new point DID NOT replace the\ \002,i4,\002-th point\002,/\002 old pivot: \002,d14.7,\002 new pivot: \ \002,d14.7,/)"; static char fmt_8010[] = "(\002 PTNEW: No point which can be replaced w\ as found\002,/)"; /* System generated locals */ integer i__1, i__2; doublereal d__1; /* Builtin functions */ integer s_wsfe(), do_fio(), e_wsfe(); /* Subroutine */ int s_stop(); /* Local variables */ static doublereal kmod, mval; static integer lenw; static doublereal vmax; static integer i__, j; static doublereal kappa; extern /* Subroutine */ int getnp_(), dcopy_(), shift_(), mintr_(); static integer icurw, dd, ig, ih, ixbase; static doublereal pivmin; static integer minpiv; extern /* Subroutine */ int unshft_(), nextnp_(); static integer np1; /* Fortran I/O blocks */ static cilist io___6 = { 0, 0, 0, fmt_1000, 0 }; static cilist io___15 = { 0, 0, 0, fmt_8000, 0 }; static cilist io___18 = { 0, 0, 0, fmt_8020, 0 }; static cilist io___19 = { 0, 0, 0, fmt_8030, 0 }; static cilist io___20 = { 0, 0, 0, fmt_8010, 0 }; /* ***************************************************************************** */ /* THIS SUBROUTINE TRIES TO FIND A NEW POINT THAT CAN BE INCLUDED IN */ /* INTERPOLATION TO IMPROVE GEOMETRY. IF THE INTERPOLATION SET IS NOT */ /* COMPLETE, THEN WE LOOK FOR A 'GOOD' POINT TO ADD TO THE SET; IF */ /* THE INTERPOLATION SET IS COMPLETE, THEN WE CHOOSE A 'BAD' POINT THAT */ /* WE WOULD LIKE TO REPLACE AND TRY TO FIND A 'GOOD' POINT TO REPLACE IT. */ /* IF WE FIND SUCH A POINT THEN FAIL=FASLE AND THE NEW POINT IS XGNEW. */ /* IF WE DO NOT FIND IT THEN FAIL=TRUE. */ /* PARAMETERS */ /* POLY (INPUT) THE ARRAY OF NEWTON FUNDAMENTAL POLYNOMIALS */ /* PNTINT (INPUT) THE ARRAY OF INTERPOLATION POINTS */ /* NIND (INPUT) CARDINALITY OF THE INTERPOLATION SET */ /* LB (INPUT) LOWER BOUNDS OF THE PROBLEM */ /* LU (INPUT) UPPER BOUNDS OF THE PROBLEM */ /* DELTA (INPUT) TRUST REGION RADIUS */ /* A (INPUT) MATRIX OF LINEAR CONSTRAINTS OF THE PROBLEM */ /* PIVVAL (INPUT) ARRAY OF PIVOT VALUES ASSOCIATED WITH THE */ /* INTERPOLATION POINTS */ /* PIVTHR (INPUT) THRESHOLD FOR ACCEPTABLE PIVOT VALUE */ /* XCHTHR (INPUT) THRESHOLD FOR THE MINIMUM ACCEPTABLE IMPROVEMENT */ /* IN THE PIVOT VALUE DUE TO REPLACING A POINT */ /* X0 (INPUT) CURRENT SHIFT OF THE POINTS FROM THE ORIGINAL POSITION */ /* XGNEW (OUTPUT) THE NEW POINT FOUND IN ORDER TO IMPROVE GEOMETRY */ /* VNEW (OUTPUT) THE PIVOT VALUE CORRESPONDING TO XGNEW, IN CASE IT */ /* CAN BE ADDED */ /* IPOLY (OUTPUT) THE INDEX WHICH SHOULD BE ASSIGNED TO XGNEW, WHEN */ /* IT IS INCLUDED IN THE INTERPOLATION SET. */ /* (IF XGNEW IS ADDED, THEN IPOLY=NIND+1, OTHERWISE */ /* IPOLY IS THE INDEX OF THE POINT WHICH IS TO BE */ /* REPLACED BY XGNEW */ /* FAIL (OUTPUT) LOGICAL VARIABLE INDICATING IF THE SEARCH SUCCEEDED */ /* FAIL=TRUE INDICATES THAT WE COULD NEITHER FIND ANY */ /* POINT TO BE ADDED TO THE SET, NOR WE COULD FIND A */ /* POINT WHICH WOULD REPLACE ANOTHER POINT IN THE SET */ /* WITH SIGNIFICANT IMPROVEMENT IN GEOMETRY. */ /* (NOTICE: WE DID NOT FIND SUCH POINTS, DOES NOT MEAN */ /* THAT THEY DO NOT EXISTS) */ /* MINMAX (INPUT) INDICATOR THAT ENABLES USER TO MINIMIZE THE NEXT */ /* POLYNOMIAL, IF ON THE PREVIOUS CALL IT WAS MAXIMIZED */ /* 0 IF WE SHOULD MAXIMIZE OF ABSOLUTE VALUE OF PIVOT */ /* -1 IF WE SHOULD MINIMIZE THE REAL VALUE OF PIVOT */ /* 1 IF WE SHOULD MAXIMIZE THE REAL VALUE OF PIVOT */ /* (OUTPUT) 1 IF THE INPUT VALUE WAS 0 AND ANSWER ACHIEVED */ /* BY MAXIMIZATION */ /* -1 IF THE INPUT VALUE WAS 0 AND ANSWER ACHIEVED */ /* BY MINIMIZATION */ /* 0 IF THE INPUT VALUE WAS 1 OR -1 */ /* ******************************************************************** */ /* COMMON VARIABLES */ /* PRINTOUT PARAMETERS */ /* MODEL PARAMETERS */ /* EXTERNAL SUBROUTINES */ /* APPLICATIONS: MINTR, GETNP, NEXTNP */ /* LOCAL VARIABLES */ /* CONSTANTS */ /* PARTITION THE REAL SPACE */ /* Parameter adjustments */ --pntint; --x0; --xgnew; --poly; --pivval; --a; --ub; --lb; --wrk; --iwrk; /* Function Body */ ig = 1; ih = ig + *n; ixbase = ih + *n * *n; icurw = ixbase + *n; lenw = *lwrk - icurw + 1; /* CHECK IF REAL SPACE IS SUFFICIENT */ if (lenw < 1) { if (dfocm_1.iprint >= 0) { io___6.ciunit = dfocm_1.iout; s_wsfe(&io___6); i__1 = -lenw + 1; do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer)); e_wsfe(); } s_stop("", (ftnlen)0); } np1 = *n + 1; dd = np1 * (*n + 2) / 2; /* NO POINT WAS FOUND YET, THEREFORE FAIL IS TRUE */ *fail = TRUE_; vmax = 0.; /* IF THE INTERPOLATION SET IS INCOMPLETE THEN WE LOOK FOR */ /* A POINT TO BE ADDED. WE DO IT BY MAXIMIZING THE NEXT PIVOT. */ /* THE NEXT PIVOT IS THE VALUE AT A CERTAIN POINT OF THE 'NEXT' */ /* (NIND+1ST) POLYNOMIAL, UPDATED SO, THAT IS IS ZERO AT ALL */ /* POINT IN THE SET. */ if (*nind < dd) { /* UPDATE THE 'NEXT POLYNOMIAL, SO THAT IT IS ZERO AT ALL POINTS */ /* OF THE INTERPOLATION SET */ i__1 = *nind + 1; nextnp_(&i__1, &poly[1], &pntint[1], nind, n, lpoly, lptint); /* PUT THE COEFFICIENTS OF THE POLYNOMIAL IN THE FORM OF QUADRATIC FORM */ i__1 = *nind + 1; getnp_(&i__1, &poly[1], lpoly, n, &kappa, &wrk[ig], &wrk[ih]); if (*minmax != 1) { /* SET PARAMETERS FOR TRUST-REGION MINIMIZATION */ dcopy_(n, &pntint[(*base - 1) * *n + 1], &c__1, &wrk[ixbase], & c__1); kmod = kappa; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { mdlpar_1.gmod[i__ - 1] = wrk[ig + i__ - 1]; kmod -= mdlpar_1.gmod[i__ - 1] * x0[i__]; i__2 = *n; for (j = 1; j <= i__2; ++j) { mdlpar_1.hmod[i__ + j * 200 - 201] = wrk[ih + (i__ - 1) * *n + j - 1]; mdlpar_1.gmod[i__ - 1] -= mdlpar_1.hmod[i__ + j * 200 - 201] * x0[j]; kmod += x0[i__] * .5 * mdlpar_1.hmod[i__ + j * 200 - 201] * x0[j]; /* L10: */ } /* L20: */ } /* MINIMIZE THE 'NEXT' POLYNOMIAL OVER THE TRUST-REGION */ unshft_(n, &x0[1], &wrk[ixbase]); mintr_(n, &wrk[ixbase], &mval, delta, &lb[1], &ub[1], &a[1], lda, ldcj, nclin, ncnln, &wrk[icurw], &lenw, &iwrk[1], liwrk, info); shift_(n, &x0[1], &wrk[ixbase]); /* STORE THE VALUE AND THE POINT */ if (*info == 0) { *vnew = mval + kmod; vmax = abs(*vnew); dcopy_(n, &wrk[ixbase], &c__1, &xgnew[1], &c__1); } } /* SET PARAMETERS FOR MAXIMIZATION OVER THE TRUST-REGION (MINIMIZATION */ /* WITH NEGATIVE SIGN) */ if (*minmax != -1) { kmod = -kappa; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { mdlpar_1.gmod[i__ - 1] = -wrk[ig + i__ - 1]; kmod -= mdlpar_1.gmod[i__ - 1] * x0[i__]; i__2 = *n; for (j = 1; j <= i__2; ++j) { mdlpar_1.hmod[i__ + j * 200 - 201] = -wrk[ih + (i__ - 1) * *n + j - 1]; mdlpar_1.gmod[i__ - 1] -= mdlpar_1.hmod[i__ + j * 200 - 201] * x0[j]; kmod += x0[i__] * .5 * mdlpar_1.hmod[i__ + j * 200 - 201] * x0[j]; /* L30: */ } /* L40: */ } dcopy_(n, &pntint[(*base - 1) * *n + 1], &c__1, &wrk[ixbase], & c__1); /* MAXIMIZE THE 'NEXT' POLYNOMIAL OVER THE TRUST-REGION */ unshft_(n, &x0[1], &wrk[ixbase]); mintr_(n, &wrk[ixbase], &mval, delta, &lb[1], &ub[1], &a[1], lda, ldcj, nclin, ncnln, &wrk[icurw], &lenw, &iwrk[1], liwrk, info); shift_(n, &x0[1], &wrk[ixbase]); } /* CHOOSE THE LARGER ABSOLUTE VALUE BETWEEN THE MAXIMUM AND THE MINIMUM */ /* AND PICK APPROPRIATE POINT */ if (*minmax == 0) { if ((d__1 = mval + kmod, abs(d__1)) > vmax && *info == 0) { *vnew = -mval - kmod; *minmax = 1; vmax = abs(*vnew); dcopy_(n, &wrk[ixbase], &c__1, &xgnew[1], &c__1); } else { *minmax = -1; } } else { *minmax = 0; if ((d__1 = mval + kmod, abs(d__1)) > vmax && *info == 0) { *vnew = -mval - kmod; vmax = abs(*vnew); dcopy_(n, &wrk[ixbase], &c__1, &xgnew[1], &c__1); } } /* IF THE PIVOT VALUE IS ACCEPTABLE, THEN WE ARE DONE */ if (vmax > *pivthr) { *fail = FALSE_; *ipoly = *nind + 1; *info = 0; if (dfocm_1.iprint >= 3) { io___15.ciunit = dfocm_1.iout; s_wsfe(&io___15); do_fio(&c__1, (char *)&(*ipoly), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&vmax, (ftnlen)sizeof(doublereal)); e_wsfe(); } } } /* IF WE DID NOT MANAGE TO FIND A POINT TO ADD (BECAUSE THE */ /* INTERPOLATION SET WAS FULL OR PIVOT VALUE TOO SMALL), THEN */ /* WE TRY TO FIND A POINT TO INCLUDE BY REPLACING SOME OTHER POINT */ if (*fail) { /* FIRST CHOOSE A POINT WHICH WE WANT TO REPLACE. IT WILL BE THE POINT */ /* WITH THE SMALLEST ASSOCIATED PIVOT VALUE */ pivmin = 1. - dfocm_1.cnstol * *delta; minpiv = 0; i__1 = *nind; for (i__ = 1; i__ <= i__1; ++i__) { if ((d__1 = pivval[i__], abs(d__1)) < pivmin) { pivmin = (d__1 = pivval[i__], abs(d__1)); minpiv = i__; } /* L45: */ } /* IF THE THERE SMALLEST PIVOT IS REASONABLY SMALL AND THE CHOSEN POINT */ /* IS NOT THE BASE, THEN WE SET IPOLY EQUAL TO MINPIV - THE INDEX OF */ /* THE CANDIDATE FOR REPLACEMENT. */ if (minpiv > 0 && minpiv != *base && (minpiv > np1 || *nind <= np1)) { *ipoly = minpiv; /* MAXIMIZE THE ABSOLUTE VALUE OF THE NEWTON POLYNOMIAL WITH INDEX IPOLY */ getnp_(ipoly, &poly[1], lpoly, n, &kappa, &wrk[ig], &wrk[ih]); dcopy_(n, &pntint[(*base - 1) * *n + 1], &c__1, &wrk[ixbase], & c__1); kmod = kappa; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { mdlpar_1.gmod[i__ - 1] = wrk[ig + i__ - 1]; kmod -= mdlpar_1.gmod[i__ - 1] * x0[i__]; i__2 = *n; for (j = 1; j <= i__2; ++j) { mdlpar_1.hmod[i__ + j * 200 - 201] = wrk[ih + (i__ - 1) * *n + j - 1]; mdlpar_1.gmod[i__ - 1] -= mdlpar_1.hmod[i__ + j * 200 - 201] * x0[j]; kmod += x0[i__] * .5 * mdlpar_1.hmod[i__ + j * 200 - 201] * x0[j]; /* L50: */ } /* L60: */ } /* DO TRUST-REGION MINIMIZATION */ unshft_(n, &x0[1], &wrk[ixbase]); mintr_(n, &wrk[ixbase], &mval, delta, &lb[1], &ub[1], &a[1], lda, ldcj, nclin, ncnln, &wrk[icurw], &lenw, &iwrk[1], liwrk, info); shift_(n, &x0[1], &wrk[ixbase]); if (*info == 0) { *vnew = mval + kmod; vmax = abs(*vnew); dcopy_(n, &wrk[ixbase], &c__1, &xgnew[1], &c__1); } kmod = -kappa; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { mdlpar_1.gmod[i__ - 1] = -wrk[ig + i__ - 1]; kmod -= mdlpar_1.gmod[i__ - 1] * x0[i__]; i__2 = *n; for (j = 1; j <= i__2; ++j) { mdlpar_1.hmod[i__ + j * 200 - 201] = -wrk[ih + (i__ - 1) * *n + j - 1]; mdlpar_1.gmod[i__ - 1] -= mdlpar_1.hmod[i__ + j * 200 - 201] * x0[j]; kmod += x0[i__] * .5 * mdlpar_1.hmod[i__ + j * 200 - 201] * x0[j]; /* L70: */ } /* L80: */ } dcopy_(n, &pntint[(*base - 1) * *n + 1], &c__1, &wrk[ixbase], & c__1); /* DO TRUST-REGION MAXIMIZATION */ unshft_(n, &x0[1], &wrk[ixbase]); mintr_(n, &wrk[ixbase], &mval, delta, &lb[1], &ub[1], &a[1], lda, ldcj, nclin, ncnln, &wrk[icurw], &lenw, &iwrk[1], liwrk, info); shift_(n, &x0[1], &wrk[ixbase]); /* CHOOSE THE BETTER POINT BETWEEN MAXIMIZER AND MINIMIZER */ if ((d__1 = mval + kmod, abs(d__1)) > vmax && *info == 0) { *vnew = -mval - kmod; vmax = abs(*vnew); dcopy_(n, &wrk[ixbase], &c__1, &xgnew[1], &c__1); } /* CHECK IF THE NEW PIVOT GIVES AT LEAST 'XCHTHR' TIMES IMPROVEMENT */ /* OVER THE OLD PIVOT VALUE. IF IT DOES, WE ACCEPT THE POINT BY SETTING */ /* FAIL=FALSE */ if (vmax / *delta > *xchthr / pivmin) { *fail = FALSE_; *info = 0; if (dfocm_1.iprint >= 3) { io___18.ciunit = dfocm_1.iout; s_wsfe(&io___18); do_fio(&c__1, (char *)&(*ipoly), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&pivmin, (ftnlen)sizeof(doublereal)) ; do_fio(&c__1, (char *)&vmax, (ftnlen)sizeof(doublereal)); e_wsfe(); } } else { if (dfocm_1.iprint >= 3) { io___19.ciunit = dfocm_1.iout; s_wsfe(&io___19); do_fio(&c__1, (char *)&(*ipoly), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&pivmin, (ftnlen)sizeof(doublereal)) ; do_fio(&c__1, (char *)&vmax, (ftnlen)sizeof(doublereal)); e_wsfe(); } } } else { if (dfocm_1.iprint >= 3) { io___20.ciunit = dfocm_1.iout; s_wsfe(&io___20); e_wsfe(); } } } return 0; } /* ptnew_ */