/* ptinit.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 { doublereal gmod[200], hmod[40000] /* was [200][200] */; } mdlpar_; #define mdlpar_1 mdlpar_ struct { integer lpoly, lpnts, lvalue, lptint, lvlint; } rpart_; #define rpart_1 rpart_ struct { integer iout, iprint; doublereal mcheps, cnstol; } dfocm_; #define dfocm_1 dfocm_ struct { integer npmin, layer, effort; } opti_; #define opti_1 opti_ /* Table of constant values */ static integer c__1 = 1; static integer c__2 = 2; static integer c__0 = 0; /* Copyright (C) 2000, International Business Machines */ /* Corporation and others. All Rights Reserved. */ /* Subroutine */ int ptinit_(n, x, ldx, fx, nx, np0, nf, poly, points, values, pntint, valint, sp2in, in2sp, nind, base, dist, delta, delmax, pivthr, pivval, lb, ub, a, lda, ldcj, nclin, ncnln, scale, scal, wrk, lwrk, iwrk, liwrk, inform__, neqcon) integer *n; doublereal *x; integer *ldx; doublereal *fx; integer *nx, *np0, *nf; doublereal *poly, *points, *values, *pntint, *valint; integer *sp2in, *in2sp, *nind, *base; doublereal *dist, *delta, *delmax, *pivthr, *pivval, *lb, *ub, *a; integer *lda, *ldcj, *nclin, *ncnln, *scale; doublereal *scal, *wrk; integer *lwrk, *iwrk, *liwrk, *inform__, *neqcon; { /* Format strings */ static char fmt_1000[] = "(\002 WARNING: THE FIRST INITIAL POINT IS OUT \ OF BOUNDS\002,/10x,\002IT WILL BE PROJECTED ON THESE BOUNDS\002,/)"; static char fmt_2000[] = "(\002 SOME PARAMETER IN THE PROBLEM FORMULATIO\ N \002,/\002 HAS ILLEGAL VALUE OR DERIVATIVES OF THE CONSTRAINTS\002/\002 AR\ E WRONG. THE PROGRAM WILL STOP\002,/)"; static char fmt_2010[] = "(\002 FEASIBLE SET SEEMS TO BE EMPTY \002,/\ \002 THE PROGRAM WILL STOP\002,/)"; static char fmt_1010[] = "(\002 WARNING: THE FIRST INITIAL POINT IS NOT \ FEASIBLE\002,/10x,\002IT WILL BE PROJECTED ON THE FEASIBLE SET\002,/)"; static char fmt_2020[] = "(\002 FUNCTION VALUE WAS NOT FOUND FOR INITIAL\ POINT\002,/\002 OR ITS PROJECTION. THE PROGRAM WILL STOP\002,/)"; static char fmt_1020[] = "(\002 WARNING: THE \002,i4,\002-TH INITIAL POI\ NT IS OUT OF BOUNDS\002,/10x,\002IT WILL BE IGNORED\002,/)"; static char fmt_1030[] = "(\002 WARNING: THE \002,i4,\002-TH INITIAL POI\ NT DOES NOT SATISFY\002,/10x,\002LINEAR CONSTRAINTS, IT WILL BE IGNORED\002,\ /)"; static char fmt_1040[] = "(\002 WARNING: THE \002,i4,\002-TH INITIAL POI\ NT IS NOT FEASIBLE\002,/10x,\002IT WILL BE IGNORED\002,/)"; static char fmt_8000[] = "(\002 PTINIT: Getting an auxiliary point \002/)" ; static char fmt_1050[] = "(\002 WARNING: AT LEAST ONE INITIAL POINT IS \ TOO FAR FROM\002,/\002 THE BASE, IT WILL NOT BE IN THE INTERPOLATION SET \ \002,/\002 TO INCLUDE IT, INCREASE PARAMETER DELTA OR LAYER\002,/)"; /* System generated locals */ integer i__1, i__2, i__3; doublereal d__1, d__2; /* Builtin functions */ integer s_wsfe(), e_wsfe(), do_fio(); /* Local variables */ extern doublereal ddot_(); static integer i__, j; static doublereal distb; static logical iferr; extern /* Subroutine */ int shift_(), dcopy_(), unscl_(), mintr_(); static integer np, ix; static logical infeas; extern /* Subroutine */ int nbuild_(); static logical bdvltd; extern /* Subroutine */ int getdis_(), intdim_(), funcon_(), ranlux_(); static doublereal del; extern /* Subroutine */ int scl_(); static doublereal val; static integer inp; extern /* Subroutine */ int fun_(); static integer iix; /* Fortran I/O blocks */ static cilist io___6 = { 0, 0, 0, fmt_1000, 0 }; static cilist io___9 = { 0, 0, 0, fmt_2000, 0 }; static cilist io___10 = { 0, 0, 0, fmt_2010, 0 }; static cilist io___11 = { 0, 0, 0, fmt_1010, 0 }; static cilist io___12 = { 0, 0, 0, fmt_2020, 0 }; static cilist io___17 = { 0, 0, 0, fmt_1020, 0 }; static cilist io___18 = { 0, 0, 0, fmt_1030, 0 }; static cilist io___19 = { 0, 0, 0, fmt_2000, 0 }; static cilist io___20 = { 0, 0, 0, fmt_1040, 0 }; static cilist io___21 = { 0, 0, 0, fmt_1040, 0 }; static cilist io___22 = { 0, 0, 0, fmt_8000, 0 }; static cilist io___24 = { 0, 0, 0, fmt_2000, 0 }; static cilist io___25 = { 0, 0, 0, fmt_2010, 0 }; static cilist io___26 = { 0, 0, 0, fmt_2000, 0 }; static cilist io___27 = { 0, 0, 0, fmt_2010, 0 }; static cilist io___28 = { 0, 0, 0, fmt_1050, 0 }; /* Added by Sergej V. Aksenov; integer flag to see if fun_ had internal error and it is necessary to abort dfo */ int funflag = 0; /* ****************************************************************** */ /* THIS SUBROUTINE BUILDS THE INITIAL INTERPOLATION SET. */ /* GIVEN ONE STARTING POINT AND A TOTAL NUMBER OF POINT */ /* IN INITIAL MODEL (NP0 - SPECIFIED BY THE USER) ADDITIONAL */ /* POINTS ARE CONSTRUCTED IN SOME NEIGHBORHOOD OF THE INITIAL POINT. */ /* AT LEAST 2 POINTS ARE NEEDED TO BUILD INITIAL MODEL. */ /* AFTER DETERMINING THE POINTS AND THEIR FUNCTION VALUES WE */ /* BUILD THE BASIS OF NEWTON FUNDAMENTAL POLYNOMIALS FOR THE */ /* INITIAL INTERPOLATION SET */ /* PARAMETERS */ /* N (INPUT) PROBLEM DIMENSION */ /* NP0 (INPUT) NUMBER OF POINTS REQUIRED FOR INITIAL MODEL BY THE */ /* USER */ /* (OUTPUT) NUMBER OF POINTS FOR THE INITIAL MODEL SUPPLIED BY */ /* THE PROGRAM */ /* X (INPUT) ARRAY OF LENGTH LDX*NX CONTAINING THE 'NX 'STARTING */ /* POINTS, PROVIDED BY THE USER. */ /* (OUTPUT) CURRENT BEST POINT, POSSIBLY AFTER PROJECTION ONTO */ /* FEASIBLE SET, IN FIRST N ENTREES */ /* LDX (INPUT) LEADING DIMENSION OF ARRAY X */ /* FX (INPUT) ARRAY OF 'NX' FUNCTION VALUES AT POINTS IN 'X' */ /* (OUTPUT) THE BEST CURRENT VALUE, IN THE FIRST ENTREE OF 'FX' */ /* NX (INPUT) NUMBER OF INITIAL POINTS PROVIDED */ /* DELTA (INPUT) TRUST REGION RADIUS */ /* DELMAX (INPUT) THE MAXIMUM TRUST REGION RADIUS ALLOWED */ /* PIVTHR (INPUT) PIVOT THRESHOLD VALUE */ /* LB (INPUT) ARRAY OF LENGTH N+NCLIN+NCNLN OF LOWER BOUNDS */ /* UB (INPUT) '' '' '' UPPER '' */ /* NCLIN (INPUT) NUMBER OF LINEAR ANALYTIC CONSTRAINTS */ /* A (INPUT) (LDA X N) MATRIX OF LINEAR ANALYTIC CONSTRAINTS */ /* LDA (INPUT) LEADING DIMENSION OF MATRIX A */ /* LDCJ (INPUT) LEADING DIMENSION OF THE JACOBIAN OF NONLINEAR */ /* CONSTRAINTS, AS COMPUTED BY THE USER ROUTINE 'FUNCON' */ /* NCNLN (INPUT) NUMBER OF NONLINEAR INEQUALITIES */ /* POINTS (OUTPUT) ARRAY (N,NP0) WITH THE POOL OF POTENTIAL INTERPOLATION */ /* (SAMPLE) POINTS */ /* VALUE (OUTPUT) ARRAY (NP0) OF VALUES AT THE POINTS */ /* NF (OUTPUT) NUMBER OF FUNCTION CALLS DONE SO FAR */ /* BASE (OUTPUT) INDEX OF THE BASE POINT */ /* DISTP (OUTPUT) ARRAY (LVALUE) OF DISTANCES OF POINTS TO THE BASE POINTS */ /* POLY (OUTPUT) LONG ARRAY CONTAINING ALL NEWTON FUNDAMENTAL POLYNOMIALS */ /* NIND (OUTPUT) NUMBER OF POINTS INCLUDED IN THE INTERPOLATION */ /* PNTINT (OUTPUT) ARRAY (N*NIND) OF POINTS INCLUDED IN THE INTERPOLATION */ /* IN THE ORDER OF INCLUSION */ /* VALINT (OUTPUT) ARRAY (NIND) OF VALUES OF POINTS IN 'PNTINT' */ /* PIVVAL (OUTPUT) ARRAY (NIND) OF VALUES OF THE PIVOTS PRODUCED IN NBUILD */ /* SP2IN (OUTPUT) ARRAY OF 0/1 WHICH INDICATES WHICH POINTS IN "POINTS" */ /* ARE INCLUDED IN THE INTERPOLATION */ /* IN2SP (OUTPUT) ARRAY (NIND) OF INDICES WHICH FOR EVERY POINT IN "PNTINT" */ /* INDICATES ITS POSITION IN "POINTS" */ /* WRK REAL SPACE WORKING ARRAY */ /* IWRK INTEGER SPACE WORKING ARRAY */ /* INFORM (OUTPUT) INFORMATION ON EXIT */ /* 0 SUCCESSFUL MINIMIZATION */ /* 1 THE DERIVATIVES OF THE CONSTRAINT OR SOME PARAMETER */ /* SET BY THE USER IS INCORRECT */ /* 2 PROBLEM IS PROBABLY INFEASIBLE */ /* -1 CANNOT COMPUTE FUNCTION VALUE AT ONE OF THE */ /* GENERATED POINTS OR ITS PROJECTION */ /* ********************************************************************** */ /* COMMON VARIABLES: */ /* MODEL PARAMETERS */ /* LENGTH OF ARRAYS */ /* PRINTOUT PARAMETERS */ /* INTERPOLATION CONTROL PARAMETERS */ /* EXTERNAL ROUTINES */ /* LOCAL VARIABLES */ /* SUBROUTINES AND FUNCTIONS CALLED: */ /* APPLICATION: MINTR , FUN , GETDIS, SCL , UNSCL, */ /* SHIFT , RANLUX, FUNCON, NBUILD */ /* FORTRAN SUPPLIED: MIN , ABS */ /* BLAS: DCOPY , DDOT */ /* ************************************************************** */ /* PROCESS THE FIRST POINT, PROVIDED BY THE USER, TO MAKE SURE IT */ /* IS INCLUDED IN THE SAMPLE SET */ /* ************************************************************** */ /* Parameter adjustments */ --scal; --fx; --x; --dist; --sp2in; --poly; --points; --values; --pntint; --valint; --in2sp; --pivval; --a; --ub; --lb; --wrk; /* Function Body */ iferr = FALSE_; bdvltd = FALSE_; infeas = FALSE_; /* COPY THE POINT IN THE SET OF CURRENT SAMPLE POINTS */ dcopy_(n, &x[1], &c__1, &points[1], &c__1); values[1] = fx[1]; val = 0.; /* **************************************************** */ /* CHECK FEASIBILITY OF PROVIDED POINT */ /* **************************************************** */ /* CHECK IF THE INITIAL POINT IS FEASIBLE FOR SIMPLE BOUNDS */ /* IF IT IS NOT, THEN PROJECT IT ON THE BOUNDS */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { if (points[i__] < lb[i__]) { points[i__] = lb[i__]; bdvltd = TRUE_; } else if (points[i__] > ub[i__]) { points[i__] = ub[i__]; bdvltd = TRUE_; } /* L10: */ } /* IF PROJECTED, PRINT WARNING MESSAGE */ if (bdvltd) { if (dfocm_1.iprint >= 0) { io___6.ciunit = dfocm_1.iout; s_wsfe(&io___6); e_wsfe(); } dcopy_(n, &points[1], &c__1, &x[1], &c__1); } /* CHECK FEASIBILITY OF FIRST POINT WRT LINEAR CONSTRAINTS */ i__1 = *nclin; for (i__ = 1; i__ <= i__1; ++i__) { val = ddot_(n, &a[i__], lda, &points[1], &c__1); if (val > ub[*n + i__] || val < lb[*n + i__]) { infeas = TRUE_; } /* L20: */ } val = 0.; /* ****************************************************************** */ /* IF THERE ARE NONLINEAR CONSTRAINTS OR IF THE POINT VIOLATES LINEAR */ /* CONSTRAINTS, THEN PROJECT IT ONTO FEASIBLE SET */ /* (IF IT IS ALREADY FEASIBLE, IS REMAINS THE SAME). */ /* IF F_del DENOTES INTERSECTION OF THE FEASIBLE REGION */ /* AND TRUST REGION WITH CENTER AT X_1 AND RADIUS DEL, THEN WE SOLVE */ /* 2 */ /* MIN ||X-X_1|| S.T. { X IN F_del } */ /* ********************************************************** */ del = *delta; if (*ncnln > 0 || infeas) { i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { mdlpar_1.gmod[i__ - 1] = x[i__] * -2.; mdlpar_1.hmod[i__ + i__ * 200 - 201] = 2.; i__2 = *n; for (j = i__ + 1; j <= i__2; ++j) { mdlpar_1.hmod[i__ + j * 200 - 201] = 0.; mdlpar_1.hmod[j + i__ * 200 - 201] = 0.; /* L30: */ } /* L40: */ } /* THIS MINIMIZATION PRODUCES PROJECTION OF THE POINT ONTO */ /* FEASIBLE REGION, INTERSECTED WITH TRUST REGION WITH RADIUS DEL */ /* FIND THE PROJECTION */ L50: mintr_(n, &points[1], &val, &del, &lb[1], &ub[1], &a[1], lda, ldcj, nclin, ncnln, &wrk[1], lwrk, iwrk, liwrk, inform__); if (*inform__ == 1) { if (dfocm_1.iprint > 0) { io___9.ciunit = dfocm_1.iout; s_wsfe(&io___9); e_wsfe(); } return 0; } else if (*inform__ == 2) { /* IF FEASIBLE SOLUTION WAS NOT FOUND, TRY TO INCREASE */ /* TRUST REGION RADIUS */ if (del < *delmax) { del *= 2; goto L50; } else { /* IF THE TRUST REGION RADIUS IS AT IT'S MAXIMUM VALUE */ /* THEN ASSUME THE PROBLEM IS INFEASIBLE AND QUIT */ if (dfocm_1.iprint > 0) { io___10.ciunit = dfocm_1.iout; s_wsfe(&io___10); e_wsfe(); } return 0; } } /* IF THE POINT HAS CHANGED, THEN THE STARTING POINT IS NOT */ /* FEASIBLE AND WE FOUND ITS PROJECTION. PRINT A WARNING. */ val = 0.; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { val += (d__1 = points[i__] - x[i__], abs(d__1)); /* L60: */ } if (dfocm_1.iprint >= 0) { if (val > *n * 100 * dfocm_1.mcheps) { io___11.ciunit = dfocm_1.iout; s_wsfe(&io___11); e_wsfe(); } } } if (val > *n * 100 * dfocm_1.mcheps || bdvltd) { /* Changed by Sergej V. Aksenov to read the output flag of fun_ into int funflag */ /* COMPUTE FUNCTION VALUE FOR THE FIRST POINT IF IT */ /* HAD TO BE PROJECTED. IF SUCH COMPUTATION */ /* FAILS, THEN QUIT THE PROGRAM */ ++(*nf); if (*scale != 0) { unscl_(n, &points[1], &scal[1]); } funflag = fun_(n, &points[1], &values[1], &iferr); if (*scale != 0) { scl_(n, &points[1], &scal[1]); } if (iferr) { if (dfocm_1.iprint > 0) { io___12.ciunit = dfocm_1.iout; s_wsfe(&io___12); e_wsfe(); } *inform__ = -1; return 0; } } /* Added by Sergej V. Aksenov to get out of dfo if funflag=-1 the output code of ptinit is set to 99 */ if (funflag) { *inform__ = 99; return 0; } /* INITIALIZE POINTERS TO CURRENT SAMPLE POINT (IN 'POINTS') */ np = 1; inp = *n + 1; /* CYCLE OVER ALL OTHER POINTS PROVIDED BY USER */ i__1 = *nx; for (ix = 2; ix <= i__1; ++ix) { /* POINTER TO CURRENT POINT PROVIDED BY USER (IN 'X') */ iix = (ix - 1) * *ldx + 1; /* COPY THE POINT IN THE SET OF CURRENT SAMPLE POINTS */ dcopy_(n, &x[iix], &c__1, &points[inp], &c__1); /* **************************************************** */ /* CHECK FEASIBILITY OF PROVIDED POINT */ /* **************************************************** */ /* CHECK IF THE INITIAL POINT IS FEASIBLE FOR SIMPLE BOUNDS */ /* IF IT IS NOT, SKIP THIS POINT, AND MOVE TO THE NEXT ONE */ /* PRINT A WARNING */ i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { if (points[inp + i__ - 1] < lb[i__] - dfocm_1.cnstol || points[ inp + i__ - 1] > ub[i__] + dfocm_1.cnstol) { if (dfocm_1.iprint >= 0) { io___17.ciunit = dfocm_1.iout; s_wsfe(&io___17); do_fio(&c__1, (char *)&ix, (ftnlen)sizeof(integer)); e_wsfe(); } goto L120; } /* L70: */ } /* CHECK IF THE CURRENT POINT IS FEASIBLE WRT LINEAR CONSTRAINTS */ /* IF IT IS NOT, SKIP THIS POINT, AND MOVE TO THE NEXT ONE */ /* PRINT A WARNING */ i__2 = *nclin; for (i__ = 1; i__ <= i__2; ++i__) { val = ddot_(n, &a[i__], lda, &points[inp], &c__1); if (val > ub[*n + i__] + dfocm_1.cnstol || val < lb[*n + i__] - dfocm_1.cnstol) { if (dfocm_1.iprint >= 0) { io___18.ciunit = dfocm_1.iout; s_wsfe(&io___18); do_fio(&c__1, (char *)&ix, (ftnlen)sizeof(integer)); e_wsfe(); } goto L120; } /* L80: */ } /* CHECK IF THE POINT IS FEASIBLE WRT NONLINEAR CONSTRAINTS */ /* BY PROJECTING IT ONTO THE FEASIBLE REGION (SEE HOW IT IS */ /* DONE FOR THE FIRST POINT) AND CHECKING IF THE PROJECTION */ /* IS DIFFERENT FROM THE POINT */ val = 0.; if (*ncnln > 0) { i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { mdlpar_1.gmod[i__ - 1] = x[iix + i__ - 1] * -2.; mdlpar_1.hmod[i__ + i__ * 200 - 201] = 2.; i__3 = *n; for (j = i__ + 1; j <= i__3; ++j) { mdlpar_1.hmod[i__ + j * 200 - 201] = 0.; mdlpar_1.hmod[j + i__ * 200 - 201] = 0.; /* L90: */ } /* L100: */ } /* THIS MINIMIZATION PRODUCES PROJECTION OF THE POINT ONTO */ /* FEASIBLE REGION, INTERSECTED WITH TRUST REGION WITH RADIUS DEL */ /* FIND THE PROJECTION */ mintr_(n, &points[inp], &val, &del, &lb[1], &ub[1], &a[1], lda, ldcj, nclin, ncnln, &wrk[1], lwrk, iwrk, liwrk, inform__); if (*inform__ == 1) { if (dfocm_1.iprint > 0) { io___19.ciunit = dfocm_1.iout; s_wsfe(&io___19); e_wsfe(); } return 0; } else if (*inform__ == 2) { /* IF NO FEASIBLE SOLUTION WAS FOUND THEN SKIP THIS POINTS AND */ /* MOVE TO THE NEXT ONE */ if (dfocm_1.iprint > 0) { io___20.ciunit = dfocm_1.iout; s_wsfe(&io___20); do_fio(&c__1, (char *)&ix, (ftnlen)sizeof(integer)); e_wsfe(); } goto L120; } /* IF THE POINT HAS CHANGED, THEN THE STARTING POINT IS NOT */ /* FEASIBLE, WE SKIP IT AND MOVE TO THE NEXT POINT */ val = 0.; i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { val += (d__1 = points[inp + i__ - 1] - x[iix + i__ - 1], abs( d__1)); /* L110: */ } if (val > *n * 100 * dfocm_1.mcheps) { if (dfocm_1.iprint > 0) { io___21.ciunit = dfocm_1.iout; s_wsfe(&io___21); do_fio(&c__1, (char *)&ix, (ftnlen)sizeof(integer)); e_wsfe(); } goto L120; } } /* IF THE POINT PASSED ALL FEASIBILITY TESTS, THEN ACCEPT IT AS */ /* A SAMPLE POINT AND RECORD ITS FUNCTION VALUE */ ++np; inp += *n; values[np] = fx[ix]; L120: ; } *np0 = np; /* -------------------------------------------------------------- */ /* IF THERE IS ONLY ONE POINT IN THE SAMPLE SET, TRY TO FIND ANOTHER */ /* -------------------------------------------------------------- */ /* GENERATE A POINT RANDOMLY WITHIN DELTA DISTANCE FROM X */ /* AND SATISFYING THE SIMPLE BOUNDS */ /* WRITE THE POINT IN ARRAY 'POINTS' */ if (*np0 == 1) { if (dfocm_1.iprint >= 2) { io___22.ciunit = dfocm_1.iout; s_wsfe(&io___22); e_wsfe(); } dcopy_(n, &points[1], &c__1, &x[1], &c__1); /* GENERATE N RANDOM NUMBERS BETWEEN 0 AND 1 AND RECORD THEM TO */ /* POINTS STARTING FROM N+1-ST ENTREE */ ranlux_(&points[*n + 1], n); i__1 = *n; for (j = 1; j <= i__1; ++j) { distb = 0.; /* Computing MIN */ d__1 = *delta, d__2 = ub[j] - x[j]; distb = min(d__1,d__2); if (distb > dfocm_1.mcheps) { points[*n + j] = distb * points[*n + j] + x[j]; } if (distb <= dfocm_1.mcheps) { /* Computing MIN */ d__1 = *delta, d__2 = x[j] - lb[j]; distb = min(d__1,d__2); points[*n + j] = -distb * points[*n + j] + x[j]; } /* L130: */ } /* CHECK FEASIBILITY OF AUXILIARY POINT WRT LINEAR CONSTRAINTS */ infeas = FALSE_; i__1 = *nclin; for (i__ = 1; i__ <= i__1; ++i__) { val = ddot_(n, &a[i__], lda, &points[*n + 1], &c__1); if (val > ub[*n + i__] || val < lb[*n + i__]) { infeas = TRUE_; } /* L140: */ } val = 0.; /* ------------------------------------------------------- */ /* FIND THE SECOND POINT FOR INTERPOLATION */ /* ------------------------------------------------------- */ if (*ncnln > 0 || infeas) { i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { mdlpar_1.gmod[i__ - 1] = points[*n + i__] * -2.; mdlpar_1.hmod[i__ + i__ * 200 - 201] = 2.; i__2 = *n; for (j = i__ + 1; j <= i__2; ++j) { mdlpar_1.hmod[i__ + j * 200 - 201] = 0.; mdlpar_1.hmod[j + i__ * 200 - 201] = 0.; /* L150: */ } /* L160: */ } mintr_(n, &x[1], &val, delta, &lb[1], &ub[1], &a[1], lda, ldcj, nclin, ncnln, &wrk[1], lwrk, iwrk, liwrk, inform__); if (*inform__ == 1) { if (dfocm_1.iprint > 0) { io___24.ciunit = dfocm_1.iout; s_wsfe(&io___24); e_wsfe(); } return 0; } else if (*inform__ == 2) { if (dfocm_1.iprint > 0) { io___25.ciunit = dfocm_1.iout; s_wsfe(&io___25); e_wsfe(); } return 0; } /* ------------------------------------------------------- */ /* IF FIRST AND SECOND POINTS COINCIDE, FIND A DIFFERENT SECOND POINT */ /* ------------------------------------------------------- */ val = 0.; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { val += (d__1 = x[i__] - points[i__], abs(d__1)); /* L170: */ } if (val < *n * 10 * *pivthr) { i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { mdlpar_1.gmod[i__ - 1] = points[*n + i__] * 2.; mdlpar_1.hmod[i__ + i__ * 200 - 201] = -2.; i__2 = *n; for (j = i__ + 1; j <= i__2; ++j) { mdlpar_1.hmod[i__ + j * 200 - 201] = 0.; mdlpar_1.hmod[j + i__ * 200 - 201] = 0.; /* L180: */ } /* L190: */ } mintr_(n, &x[1], &val, &del, &lb[1], &ub[1], &a[1], lda, ldcj, nclin, ncnln, &wrk[1], lwrk, iwrk, liwrk, inform__); if (*inform__ == 1) { if (dfocm_1.iprint > 3) { io___26.ciunit = dfocm_1.iout; s_wsfe(&io___26); e_wsfe(); } return 0; } else if (*inform__ == 2) { if (dfocm_1.iprint > 3) { io___27.ciunit = dfocm_1.iout; s_wsfe(&io___27); e_wsfe(); } return 0; } } dcopy_(n, &x[1], &c__1, &points[*n + 1], &c__1); } /* -------------------------------------------------------------- */ /* INTERPOLATION POINTS ARE COMPUTED */ /* -------------------------------------------------------------- */ /* COMPUTE FUNCTION VALUE FOR THE AUXILIARY POINT. */ /* IF THERE ARE NONLINEAR CONSTRAINS, AND IF FUNCTION EVALUATION */ /* FAILS FOR THE AUXILIARY POINT, WE QUIT. */ /* IF THERE ARE NO NONLINEAR CONSTRAINTS, THEN IF FUNCTION EVALUATION */ /* FAILS FOR AUXILIARY POINT, WE SUBSTITUTE IT BY ITS CONVEX */ /* COMBINATION WITH THE FIRST POINT: */ /* X_k <-- 1/2(X_k+X_1) */ /* Changed by Sergej V. Aksenov to read the output flag of fun_ into int funflag */ L200: if (*scale != 0) { unscl_(n, &points[*n + 1], &scal[1]); } funflag = fun_(n, &points[*n + 1], &values[2], &iferr); if (*scale != 0) { scl_(n, &points[*n + 1], &scal[1]); } ++(*nf); if (iferr) { i__1 = *n; for (j = 1; j <= i__1; ++j) { points[*n + j] = (points[j] + points[*n + j]) * .5; /* L210: */ } /* Added by Sergej V. Aksenov to get out of dfo if funflag=-1 the output code of ptinit is set to 99 */ if (funflag) { *inform__ = 99; return 0; } /* IF THE SECOND POINT GETS TOO CLOSE TO FIRST POINT, QUIT */ getdis_(n, &c__2, &c__2, &points[1], &c__1, &dist[1], &wrk[1], lwrk); if (dist[2] < *pivthr) { *inform__ = -1; return 0; } /* IF THE SECOND POINTS BECOMES NON-FEASIBLE, QUIT */ if (*ncnln > 0) { funcon_(&c__1, ncnln, n, ldcj, iwrk, &points[*n + 1], &wrk[1], &wrk[*ncnln + 1], &c__1); i__1 = *ncnln; for (j = 1; j <= i__1; ++j) { if (wrk[j] < lb[*n + *nclin + j] - dfocm_1.cnstol || wrk[ j] > ub[*n + *nclin + j] + dfocm_1.cnstol) { *inform__ = -1; return 0; } /* L220: */ } } goto L200; } /* SET THE NUMBER OF SAMPLE POINTS TO EQUAL 2 */ *np0 = 2; } /* CHOOSE THE BASE POINT */ *base = 1; i__1 = *np0; for (i__ = 1; i__ <= i__1; ++i__) { if (values[i__] < values[*base] + dfocm_1.mcheps * 100.) { *base = i__; } /* L230: */ } /* COPY THE BASE POINT INTO X */ dcopy_(n, &points[(*base - 1) * *n + 1], &c__1, &x[1], &c__1); /* THE MAXIMUM POSSIBLE NUMBER OF INDEPENDENT INTERPOLATION POINTS */ /* DEPENDS ON THE ACTUAL DIMENSION OF THE FEASIBLE SET, RATHER THEN */ /* ON THE DIMENSION OF THE WHOLE SPACE. SUBROUTINE 'INTDIM' COMPUTES */ /* THE DIMENSION OF THE SPACE SPANNED BY EQUALITY CONSTRAINTS AT */ /* THE BASE. WE CONSIDER A CONSTRAINT AS EQUALITY IF THE DIFFERENCE */ /* BETWEEN ITS UPPER AND LOWER BOUNDS IS LESS THAN PIVOT THRESHOLD */ /* NEQCON - NUMBER OF LINEARLY INDEPENDENT EQUALITY CONSTRAINTS */ intdim_(&x[1], n, nclin, ncnln, neqcon, &a[1], lda, ldcj, &lb[1], &ub[1], pivthr, &wrk[1], lwrk, iwrk, liwrk); /* THE BASE IS CHOSEN, SHIFT THE POINTS SO THAT THE BASE IS AT THE ORIGIN */ i__1 = *np0; for (i__ = 1; i__ <= i__1; ++i__) { shift_(n, &x[1], &points[(i__ - 1) * *n + 1]); /* L240: */ } /* GET THE DISTANCES OF ALL POINTS IN 'POINTS' TO THE BASE */ getdis_(n, np0, &c__0, &points[1], base, &dist[1], &wrk[1], lwrk); /* CHECK IF ANY OF THE POINTS ARE TOO FAR, AND WOULD NOT BE INCLUDED */ /* IN INTERPOLATION SET */ bdvltd = FALSE_; i__1 = *np0; for (i__ = 1; i__ <= i__1; ++i__) { if (dist[i__] > opti_1.layer * *delta) { bdvltd = TRUE_; } /* L250: */ } if (bdvltd) { if (dfocm_1.iprint >= 0) { io___28.ciunit = dfocm_1.iout; s_wsfe(&io___28); e_wsfe(); } } /* --------------------------------------------------------------------- */ /* COMPUTE THE BASIS OF NEWTON FUNDAMENTAL POLYNOMIALS */ /* --------------------------------------------------------------------- */ if (*base == 1) { dcopy_(n, &points[1], &c__1, &pntint[1], &c__1); valint[1] = values[1]; } in2sp[*base] = *base; nbuild_(&poly[1], &points[1], &values[1], &pntint[1], &valint[1], &sp2in[ 1], &in2sp[1], nind, n, base, &dist[1], delta, pivthr, &pivval[1], neqcon); *inform__ = 0; return 0; /* L2030: */ } /* ptinit_ */