/* mbuild.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_ /* Table of constant values */ static integer c__1 = 1; static doublereal c_b58 = -1.; static doublereal c_b60 = 1.; static doublereal c_b65 = .5; /* Copyright (C) 2000, International Business Machines */ /* Corporation and others. All Rights Reserved. */ /* Subroutine */ int mdbld_(kappa, g, h__, n, nind, pntint, valint, poly, base, varnt, lpoly, lptint, neqcon, wrk, lwrk, iwrk, liwrk) doublereal *kappa, *g, *h__; integer *n, *nind; doublereal *pntint, *valint, *poly; integer *base, *varnt, *lpoly, *lptint, *neqcon; doublereal *wrk; integer *lwrk, *iwrk, *liwrk; { /* Format strings */ static char fmt_1100[] = "(\002 MDBLD: *** ERROR: LWRK TOO SMALL!\002\ /\002 IT SHOULD BE AT LEAST \002,i5)"; static char fmt_1200[] = "(\002 MDBLD: *** ERROR: LIWRK TOO SMALL!\002\ /\002 IT SHOULD BE AT LEAST \002,i5)"; /* System generated locals */ integer h_dim1, h_offset, i__1; /* Builtin functions */ integer s_wsfe(), do_fio(), e_wsfe(); /* Subroutine */ int s_stop(); /* Local variables */ static logical fail; static integer lenw, imatr, icurw, dd, ilafla, ix; extern /* Subroutine */ int mbldmf_(); static logical froben; extern /* Subroutine */ int mbldnp_(); static integer ivx; /* Fortran I/O blocks */ static cilist io___8 = { 0, 0, 0, fmt_1100, 0 }; static cilist io___9 = { 0, 0, 0, fmt_1200, 0 }; static cilist io___12 = { 0, 0, 0, fmt_1100, 0 }; /* *********************************************************************** */ /* THIS FUNCTION BUILDS A QUADRATIC INTERPOLATION */ /* FUNCTION FOR THE SUPPLIED 'NIND' DATA POINTS THAT ARE */ /* GIVEN IN ARRAY 'PNTINT' WITH VALUES IN ARRAY 'VALINT'. */ /* T T */ /* M(X)= KAPPA + G X + 0.5*X H X */ /* IF 'VARNT' EQUALS 1 AND THE MODEL IS INCOMPLETE QUADRATIC */ /* THEN MINIMUM FROBENIUS MODEL IS BUILT, OTHERWISE THE MODEL */ /* IS BUILT FROM A LINEAR COMBINATION OF CONSTRUCTED NEWTON */ /* FUNDAMENTAL POLYNOMIALS. */ /* PARAMETERS: */ /* N (INPUT) DIMENSION OF THE PROBLEM */ /* NIND (INPUT) NUMBER OF INTERPOLATION POINTS */ /* PNTINT (INPUT) LIST OF 'NIND' DATA POINTS. THE I-TH POINT OCCUPIES */ /* POSITIONS ( I - 1 ) * N + 1 TO I * N. */ /* VALINT (INPUT) VALUES OF THE OBJECTIVE FUNCTION AT THE 'NIND' */ /* DATA POINTS CONTAINED IN PNTINT. */ /* POLY (INPUT) THE ARRAY CONTAINING COEFFICIENTS OF NEWTON FUNDAMENTAL */ /* POLYNOMIALS (AS COMPUTED BY 'NBUILD') */ /* BASE (INPUT) THE INDEX (IN PNTINT) OF THE CURRENT BASE (BEST) POINT. */ /* VARNT (INPUT) VARNT=1 -- FROBENIUS MODEL IS BUILT, VARNT=2 -- SUBQUADRATIC */ /* MODEL IS BUILT */ /* KAPPA (OUTPUT) THE CONSTANT TERM OF THE INTERPOLATION MODEL. */ /* G (OUTPUT) VECTOR OF THE LINEAR TERMS OF THE INTERPOLATION MODEL. */ /* H (OUTPUT) MATRIX OF QUADRATIC TERMS OF THE INTERPOLATION MODEL. */ /* WRK (INPUT) REAL WORKSPACE OF LENGTH LWRK. */ /* IWRK (INPUT) INTEGER WORKSPACE OF LENGTH LIWRK. */ /* ************************************************************************** */ /* SUBROUTINES AND FUNCTIONS CALLED: */ /* APPLICATION: MBLDMF, MBLDNP */ /* LOCAL VARIABLES */ /* Parameter adjustments */ h_dim1 = *n; h_offset = 1 + h_dim1 * 1; h__ -= h_offset; --g; --valint; --poly; --pntint; --wrk; --iwrk; /* Function Body */ dd = (*n + 1) * (*n + 2) / 2; /* SEE IF FROBENIUS MODEL SHOULD BE COMPUTED */ froben = FALSE_; if (*varnt == 1 && *nind > *n + 1 && *nind < dd) { froben = TRUE_; } /* PARTITION REAL WORKSPACE */ if (froben) { imatr = 1; /* Computing 2nd power */ i__1 = *nind; ix = imatr + i__1 * i__1; ivx = ix + (*nind - 1) * *n; icurw = ivx + *nind - 1; lenw = *lwrk - icurw + 1; } /* CHECK IF THE WORKING SPACE IS SUFFICIENT */ if (lenw < 1 && froben) { if (dfocm_1.iprint >= 0) { io___8.ciunit = dfocm_1.iout; s_wsfe(&io___8); i__1 = -lenw + 1; do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer)); e_wsfe(); } s_stop("", (ftnlen)0); } if (*liwrk < 1 && froben) { if (dfocm_1.iprint >= 0) { io___9.ciunit = dfocm_1.iout; s_wsfe(&io___9); i__1 = -(*liwrk) + 1; do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer)); e_wsfe(); } s_stop("", (ftnlen)0); } fail = TRUE_; /* IF APPLICABLE, BUILD MINIMUM FROBENIUS MODEL */ if (froben) { mbldmf_(kappa, &g[1], &h__[h_offset], &pntint[1], &valint[1], base, & wrk[ix], &wrk[ivx], &wrk[imatr], n, nind, neqcon, lptint, & fail, &wrk[icurw], &lenw, &iwrk[1], liwrk); } /* IF MINIMUM FROBENIUS MODEL WAS NOT BUILT, THEN BUILD MODEL BASED ON */ /* NEWTON FUNDAMENTAL POLYNOMIALS */ if (fail) { ilafla = 1; icurw = ilafla + *nind; lenw = *lwrk - icurw + 1; /* CHECK IF THE WORKING SPACE IS SUFFICIENT */ if (lenw < 1) { if (dfocm_1.iprint >= 0) { io___12.ciunit = dfocm_1.iout; s_wsfe(&io___12); i__1 = -lenw + 1; do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer)); e_wsfe(); } s_stop("", (ftnlen)0); } mbldnp_(kappa, &g[1], &h__[h_offset], &pntint[1], &valint[1], &poly[1] , &wrk[ilafla], n, nind, neqcon, lptint, lpoly); } return 0; /* NON-EXECUTABLE STATEMENTS */ } /* mdbld_ */ /* ******************************************************************************* */ /* NEXT SUBROUTINE */ /* ******************************************************************************* */ /* Subroutine */ int mbldnp_(kappa, g, h__, pntint, valint, poly, lafla, n, nind, neqcon, lptint, lpoly) doublereal *kappa, *g, *h__, *pntint, *valint, *poly, *lafla; integer *n, *nind, *neqcon, *lptint, *lpoly; { /* System generated locals */ integer h_dim1, h_offset; /* Local variables */ extern /* Subroutine */ int fd_(), mterms_(); /* *********************************************************************** */ /* THIS FUNCTION BUILDS A QUADRATIC INTERPOLATION */ /* FUNCTION FOR THE SUPPLIED 'NIND' DATA POINTS THAT ARE */ /* GIVEN IN ARRAY 'PNTINT' WITH VALUES IN ARRAY 'VALINT'. */ /* T T */ /* M(X)= KAPPA + G X + 0.5*X H X */ /* THE MODEL IS BUILT AS A LINEAR COMBINATION OF 'NIND' NEWTON */ /* FUNDAMENTAL POLYNOMIALS STORED IN ARRAY 'POLY'. */ /* PARAMETERS: */ /* N (INPUT) DIMENSION OF THE PROBLEM */ /* NIND (INPUT) NUMBER OF INTERPOLATION POINTS */ /* PNTINT (INPUT) LIST OF 'NIND' DATA POINTS. THE I-TH POINT OCCUPIES */ /* POSITIONS ( I - 1 ) * N + 1 TO I * N. */ /* VALINT (INPUT) VALUES OF THE OBJECTIVE FUNCTION AT THE 'NIND' */ /* DATA POINTS CONTAINED IN PNTINT. */ /* POLY (INPUT) THE ARRAY CONTAINING COEFFICIENTS OF NEWTON FUNDAMENTAL */ /* POLYNOMIALS (AS COMPUTED BY 'NBUILD') */ /* LAFLA AN AUXILIARY ARRAY OF LENGTH 'N' FOR FINITE DIFFERENCES */ /* NEQCON (INPUT) NUMBER OF LINEARLY INDEP. EQUALITY CONSTRAINTS */ /* KAPPA (OUTPUT) THE CONSTANT TERM OF THE INTERPOLATION MODEL. */ /* G (OUTPUT) VECTOR OF THE LINEAR TERMS OF THE INTERPOLATION MODEL. */ /* H (OUTPUT) MATRIX OF QUADRATIC TERMS OF THE INTERPOLATION MODEL. */ /* ************************************************************************** */ /* SUBROUTINES AND FUNCTIONS CALLED: */ /* APPLICATION: FD, GTERMS, */ /* BLAS : DCOPY */ /* COMPUTE FINITE DIFFERENCES FROM NEWTON POLYNOMIALS, POINTS AND VALUES */ /* Parameter adjustments */ h_dim1 = *n; h_offset = 1 + h_dim1 * 1; h__ -= h_offset; --g; --lafla; --valint; --pntint; --poly; /* Function Body */ fd_(&lafla[1], &pntint[1], &valint[1], &poly[1], n, nind, neqcon, lptint, lpoly); /* COMPUTE TERMS OF THE MODEL FROM FINITE DIFFERENCES, PNTINT AND NEWTON */ /* POLYNOMIALS */ mterms_(kappa, &g[1], &h__[h_offset], &lafla[1], &poly[1], &pntint[1], nind, n, neqcon, lpoly, lptint); return 0; } /* mbldnp_ */ /* ******************************************************************************* */ /* NEXT SUBROUTINE */ /* ******************************************************************************* */ /* Subroutine */ int mbldmf_(kappa, g, h__, pntint, valint, base, y, vy, matr, n, q, neqcon, lpnt, fail, wrk, lwrk, iwrk, liwrk) doublereal *kappa, *g, *h__, *pntint, *valint; integer *base; doublereal *y, *vy, *matr; integer *n, *q, *neqcon, *lpnt; logical *fail; doublereal *wrk; integer *lwrk, *iwrk, *liwrk; { /* Format strings */ static char fmt_1000[] = "(\002 MBLDMF: *** ERROR: LWRK TOO SMALL !\002\ /\002 IT SHOULD BE AT LEAST \002,i5)"; static char fmt_2000[] = "(\002 MBLDMF: *** ERROR: LIWRK TOO SMALL \ !\002/\002 IT SHOULD BE AT LEAST \002,i5)"; static char fmt_8000[] = "(\002 MBLDMF: Minimum frobenius norm model com\ puted \002,/)"; static char fmt_4000[] = "(\002 MINIMUM FROBENIUS NORM MODEL FAILED\ !\002/\002 FACTORIZATION CANNOT BE COMPUTED\002)"; static char fmt_5000[] = "(\002 MINIMUM FROBENIUS NORM MODEL FAILED\ !\002/\002 CONDITION NUMBER CANNOT BE COMPUTED\002)"; static char fmt_8010[] = "(\002 MBLDMF: Condition number: \002,d14.7,/)"; static char fmt_3000[] = "(\002 MINIMUM FROBENIUS NORM MODEL FAILED\ !\002/\002 TOO ILL-CONDITIONED\002,d14.7,/)"; /* System generated locals */ integer h_dim1, h_offset, matr_dim1, matr_offset, i__1, i__2; doublereal d__1; /* Builtin functions */ integer s_wsfe(), do_fio(), e_wsfe(); /* Subroutine */ int s_stop(); /* Local variables */ extern /* Subroutine */ int dger_(); extern doublereal ddot_(); static integer info, i__, j, jbase; extern /* Subroutine */ int dgemv_(); static doublereal rcond, anorm; extern /* Subroutine */ int dcopy_(), daxpy_(), dtrsv_(); static integer kk, ny; extern /* Subroutine */ int dpocon_(), dpotrf_(), rzrvec_(), rzrmat_(); static doublereal sum; /* Fortran I/O blocks */ static cilist io___13 = { 0, 0, 0, fmt_1000, 0 }; static cilist io___14 = { 0, 0, 0, fmt_2000, 0 }; static cilist io___15 = { 0, 0, 0, fmt_8000, 0 }; static cilist io___21 = { 0, 0, 0, fmt_4000, 0 }; static cilist io___25 = { 0, 0, 0, fmt_5000, 0 }; static cilist io___26 = { 0, 0, 0, fmt_8010, 0 }; static cilist io___28 = { 0, 0, 0, fmt_3000, 0 }; /* *********************************************************************** */ /* THIS FUNCTION BUILDS A QUADRATIC INTERPOLATION */ /* FUNCTION FOR THE SUPPLIED 'Q' DATA POINTS THAT ARE */ /* GIVEN IN ARRAY 'PNTINT' WITH VALUES IN ARRAY 'VALINT'. */ /* T T */ /* M(X)= KAPPA + G X + 0.5*X H X */ /* THE MODEL IS INCOMPLETE QUADRATIC ( 'Q' IS LESS THAN (N+1)(N+2)/2 ) */ /* AND 'H' HAS THE SMALLEST FROBENIUS NORM AMOUNG ALL HESSIANS OF ALL */ /* POSSIBLE QUADRATIC MODELS SATISFYING THE INTERPOLATION CONDITION */ /* FOR 'PNTINT' AND 'VALINT'. */ /* PARAMETERS: */ /* N (INPUT) DIMENSION OF THE PROBLEM */ /* Q (INPUT) NUMBER OF INTERPOLATION POINTS */ /* PNTINT (INPUT) LIST OF 'NIND' DATA POINTS. THE I-TH POINT OCCUPIES */ /* POSITIONS ( I - 1 ) * N + 1 TO I * N. */ /* VALINT (INPUT) VALUES OF THE OBJECTIVE FUNCTION AT THE 'NIND' */ /* DATA POINTS CONTAINED IN PNTINT. */ /* BASE (INPUT) THE INDEX (IN PNTINT) OF THE CURRENT BASE (BEST) POINT. */ /* KAPPA (OUTPUT) THE CONSTANT TERM OF THE INTERPOLATION MODEL. */ /* G (OUTPUT) VECTOR OF THE LINEAR TERMS OF THE INTERPOLATION MODEL. */ /* H (OUTPUT) MATRIX OF QUADRATIC TERMS OF THE INTERPOLATION MODEL. */ /* FAIL (OUTPUT) INDICATES IF THE THE SUBROUTINE HAS FAILED. */ /* Y AUXILIARY ARRAY OF SHIFTED POINTS */ /* VY AUXILIARY ARRAY OF SHIFTED VALUES */ /* MATR AUXILIARY MATRIX QxQ */ /* WRK (INPUT) REAL WORKSPACE OF LENGTH LWRK. */ /* IWRK (INPUT) INTEGER WORKSPACE OF LENGTH LIWRK. */ /* ************************************************************************** */ /* LOCAL VARIABLES */ /* SUBROUTINES AND FUNCTIONS CALLED: */ /* APPLICATION: RZRMAT, RZRVEC */ /* NPSOL : FO6QGT */ /* LAPACK : DPOTRF, DPOCON, DTRSV, */ /* BLAS : DCOPY , DAXPY , DGER , DGEMV */ /* Parameter adjustments */ h_dim1 = *n; h_offset = 1 + h_dim1 * 1; h__ -= h_offset; --g; matr_dim1 = *q; matr_offset = 1 + matr_dim1 * 1; matr -= matr_offset; --vy; --valint; --y; --pntint; --iwrk; --wrk; /* Function Body */ *fail = FALSE_; /* CHECK SUFFICIENCY OF THE WORKSPACE */ if (*lwrk < *q * 3) { if (dfocm_1.iprint >= 0) { io___13.ciunit = dfocm_1.iout; s_wsfe(&io___13); i__1 = *q * 3; do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer)); e_wsfe(); } s_stop("", (ftnlen)0); } if (*liwrk < *q) { if (dfocm_1.iprint >= 0) { io___14.ciunit = dfocm_1.iout; s_wsfe(&io___14); do_fio(&c__1, (char *)&(*q), (ftnlen)sizeof(integer)); e_wsfe(); } s_stop("", (ftnlen)0); } if (dfocm_1.iprint >= 3) { io___15.ciunit = dfocm_1.iout; s_wsfe(&io___15); e_wsfe(); } /* SHIFT ALL THE POINT SO THAT THE BASE IS AT THE ORIGIN AND PUT THE */ /* POINTS, EXCEPT THE BASE, IN Y, AND THEIR VALUES IN VY. IGNORE THE */ /* DUMMY POINTS IF THERE ARE ANY */ jbase = (*base - 1) * *n; ny = 0; i__1 = *q; for (i__ = 1; i__ <= i__1; ++i__) { if (i__ != *base && (i__ <= *n + 1 - *neqcon || i__ > *n + 1)) { i__2 = *n; for (j = 1; j <= i__2; ++j) { y[ny * *n + j] = pntint[(i__ - 1) * *n + j] - pntint[jbase + j]; /* L4: */ } ++ny; } /* L5: */ } /* FORM THE MATRIX 'MATR' BY FORMULAS: IF P IS THE MATRIX WITH */ /* THE INTERPOLATION POINT AS THE COLUMNS THEN */ /* ******************************* */ /* * T 2 * */ /* * M = P P, M=0.5*M +M * */ /* * * */ /* ******************************* */ i__1 = ny; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = ny; for (j = 1; j <= i__2; ++j) { matr[i__ + j * matr_dim1] = ddot_(n, &y[(i__ - 1) * *n + 1], & c__1, &y[(j - 1) * *n + 1], &c__1); /* L10: */ } /* L20: */ } i__1 = ny; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = ny; for (j = 1; j <= i__2; ++j) { /* Computing 2nd power */ d__1 = matr[i__ + j * matr_dim1]; matr[i__ + j * matr_dim1] = d__1 * d__1 * .5 + matr[i__ + j * matr_dim1]; /* L30: */ } /* L40: */ } /* USING LAPACK ROUTINES FIND AN ESTIMATE OF THE CONDITION NUMBER OF 'MATR' */ /* 'MATR' IS SYMMETRIC, SO, FIRST, COMPUTE THE LL^T FACTORIZATION IN 'DPOTRF' */ /* IF FACTORIZATION FAILS, THEN PRINT A MESSAGE AND QUIT */ dpotrf_("U", &ny, &matr[matr_offset], q, &info, (ftnlen)1); if (info != 0) { if (dfocm_1.iprint > 2) { io___21.ciunit = dfocm_1.iout; s_wsfe(&io___21); e_wsfe(); } *fail = TRUE_; return 0; } /* ANORM=F06QGF('1', 'G' , NY , NY, MATR , Q) */ anorm = 0.; i__1 = ny; for (j = 1; j <= i__1; ++j) { sum = 0.; i__2 = ny; for (i__ = 1; i__ <= i__2; ++i__) { sum += (d__1 = matr[i__ + j * matr_dim1], abs(d__1)); /* L70: */ } anorm = max(anorm,sum); /* L80: */ } /* FIND THE ESTIMATE OF THE COND. NUMBER, IF FAIL, THEN QUIT */ dpocon_("U", &ny, &matr[matr_offset], q, &anorm, &rcond, &wrk[1], &iwrk[1] , &info, (ftnlen)1); if (info != 0) { if (dfocm_1.iprint > 2) { io___25.ciunit = dfocm_1.iout; s_wsfe(&io___25); e_wsfe(); } *fail = TRUE_; return 0; } if (dfocm_1.iprint >= 3) { io___26.ciunit = dfocm_1.iout; s_wsfe(&io___26); do_fio(&c__1, (char *)&rcond, (ftnlen)sizeof(doublereal)); e_wsfe(); } /* CHECK IF THE CONDITION NUMBER IS NOT TOO BIG */ /* IF NOT, THEN FORM THE RIGHT HAND SIDE FOR THE LINEAR SYSTEM */ /* AND SOLVE THE SYSTEM USING LAPACK ROUTINES AND FACTORIZATION */ /* COMPUTED BY 'DPOTRF'. STORE THE SOLUTION IN 'VY' */ if (rcond >= dfocm_1.mcheps / 1e4) { kk = 0; i__1 = *q; for (i__ = 1; i__ <= i__1; ++i__) { if (i__ != *base && (i__ <= *n + 1 - *neqcon || i__ > *n + 1)) { ++kk; vy[kk] = valint[i__] - valint[*base]; } /* L50: */ } dtrsv_("U", "T", "N", &ny, &matr[matr_offset], q, &vy[1], &c__1, ( ftnlen)1, (ftnlen)1, (ftnlen)1); dtrsv_("U", "N", "N", &ny, &matr[matr_offset], q, &vy[1], &c__1, ( ftnlen)1, (ftnlen)1, (ftnlen)1); /* FIND THE COEFFICIENT OF THE SHIFTED MODEL BY FORMULAS: */ /* ************************************ */ /* * G <- 0, H <- 0 * */ /* * FOR ALL POINTS P_i IN 'Y' * */ /* * T * */ /* * H <- H + VY(I)*[P_i(P_i) ] * */ /* * * */ /* * G <- G + VY(I)*P_i * */ /* * * */ /* ************************************ */ rzrvec_(&g[1], n); rzrmat_(&h__[h_offset], n, n); i__1 = ny; for (i__ = 1; i__ <= i__1; ++i__) { daxpy_(n, &vy[i__], &y[(i__ - 1) * *n + 1], &c__1, &g[1], &c__1); dger_(n, n, &vy[i__], &y[(i__ - 1) * *n + 1], &c__1, &y[(i__ - 1) * *n + 1], &c__1, &h__[h_offset], n); /* L60: */ } /* SHIFT THE MODEL BACK TO INTERPOLATE THE ORIGINAL POINTS AND FIND KAPPA */ /* ************************************ */ /* * G <- G - H*P_base * */ /* * T * */ /* * KAPPA <- VALINT(BASE)-(P_base) G * */ /* ************************************ */ dgemv_("N", n, n, &c_b58, &h__[h_offset], n, &pntint[jbase + 1], & c__1, &c_b60, &g[1], &c__1, (ftnlen)1); dcopy_(n, &g[1], &c__1, &wrk[1], &c__1); dgemv_("N", n, n, &c_b65, &h__[h_offset], n, &pntint[jbase + 1], & c__1, &c_b60, &wrk[1], &c__1, (ftnlen)1); *kappa = valint[*base] - ddot_(n, &wrk[1], &c__1, &pntint[jbase + 1], &c__1); } else { /* IF THE CONDITION NUMBER IS TOO LARGE, THEN STOP, SINCE IT SHOULD NOT */ /* BE HAPPENING */ if (dfocm_1.iprint > 2) { io___28.ciunit = dfocm_1.iout; s_wsfe(&io___28); do_fio(&c__1, (char *)&rcond, (ftnlen)sizeof(doublereal)); e_wsfe(); } *fail = TRUE_; } return 0; } /* mbldmf_ */