/***************************************************************************** * RESOLUTION OF A SET OF NON-LINEAR EQUATIONS * * -------------------------------------------------------------------------- * * COPYRIGHT 1991, BY R.S. BAIN * * Originally wrapped by bain@luther.che.wisc.edu on Tue Mar 30 16:27:48 1993 * * * * As of this writing, we do not know how to contact Rod Bain, the * * author of NNES. * * * * -- David M. Gay (dmg@bell-labs.com) * * Bell Labs, Murray Hill * * 8 February 1999 * * * * Visual C++ Release By J-P Moreau, Paris. * * (PART 3/3) * * (www.jpmoreau.fr) * *****************************************************************************/ #include #include //Headers of functions used below void ascalf (int, REAL, REAL *, REAL **, REAL *); void ascalx (int, REAL, REAL **, REAL *); void ataov (bool *, int, int, int, int, REAL **, REAL **, REAL *); void atvov (bool *, int, int, int, int, REAL **, REAL *, REAL *); void avmul (int, int, int, int, REAL **, REAL *, REAL *); void broyfa (bool, bool *, bool, bool, int, int, int, int, REAL, REAL **, REAL *, REAL *, REAL *, REAL **, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *); void broyun (bool *, int, int, int, int, REAL, REAL *, REAL *, REAL **, REAL *, REAL *, REAL *); void cholde (int, REAL *, REAL *, REAL, REAL **, REAL **); void chsolv (bool, bool *, int, int, int, int, REAL **, REAL *, REAL *); void condno (bool, bool *, int, int, int, int, REAL *, REAL **, REAL *); void delcau (bool, bool, bool *, int, int maxexp, int, int, int, REAL *, REAL *, REAL *, REAL, REAL *, REAL *, REAL *, REAL **, REAL *, REAL *); void dogleg (bool *, bool *, bool, bool *, int *, int, int *, int, int, REAL *, REAL, REAL *, REAL, REAL, REAL, REAL *, REAL *, REAL *, REAL *n, REAL *, REAL *); REAL Dot_Product(int, REAL *, REAL *); void fcn1(bool *, int, REAL *, REAL *); void fcnevl(bool *, int, int, int, int, REAL, REAL *, REAL *, REAL *); void gradf (bool, bool *, bool, bool, bool, int, int, int, int, int, int, REAL *, REAL *, REAL **, REAL *, REAL *); int IMAX(int, int); int IMIN(int, int); void initch (bool *, bool *, bool *, bool *, bool *, bool *, int, int *, int, int, int, int, int, int, int, int, int, REAL *, REAL *, REAL, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *); void innerp (bool, bool *, int, int, int, int, int, int, REAL *, REAL *, REAL *); void jacobi (bool, bool *, bool *, int, int, int, int, REAL, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL **, REAL **, REAL *, REAL *); void line (bool *, bool *, bool, bool *, bool, bool, bool *, bool *, bool *, bool *, bool *, bool *, int *, int, int *, int, int *, int, int, int *, int, int, int, int, int *, int, int, int, int *, int *, REAL *, REAL *, REAL, REAL *, REAL *, REAL, REAL, REAL, REAL, REAL *, REAL *, REAL **, REAL *, REAL *, REAL *, REAL *, REAL *, REAL **, REAL *, REAL **, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *); void Line0(FILE *); void Line1(FILE *); void llfa (bool, bool *, bool *, bool *, int, int, int, int, int, REAL, REAL *, REAL **, REAL *, REAL *, REAL *, REAL **, REAL **, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *); void llun (bool, bool *, int, int, int, int, int, REAL, REAL, REAL *, REAL *, REAL **, REAL **, REAL *, REAL *, REAL *, REAL *); void matcop (int, int, int, int, int, int, REAL **, REAL **); void matprt (int, int, REAL **); REAL MAX(REAL, REAL); void maxst (bool *, int, int, int, int, REAL, REAL *, REAL, REAL *, REAL *); REAL MIN(REAL, REAL); void mmulv(int, REAL **, REAL *, REAL *); void Msg(FILE *, char *); void nersl (bool, bool, bool, bool, int, int, int, int, int, REAL, REAL *, REAL *); void nestop (bool *, bool, bool *, bool *, bool *, int, int *, int, int *, int *, int *, int *, int *, int, int, int, int *, REAL *, REAL, REAL, REAL *, REAL, REAL *, REAL *, REAL *, REAL *, REAL *); void qrsolv(bool, bool *, int, int, int, int, REAL **, REAL *, REAL *, REAL *); void rsolv (bool, bool *, int, int, int, int, REAL **, REAL *, REAL *); REAL Sign(REAL, REAL); REAL Sqr(REAL); void twonrm (bool *, int, int, REAL, REAL *, REAL *); void update (int, int *, int, int *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *); void uvmul (int, int, int, int, REAL **, REAL *, REAL *); void qform (int n, REAL **a, REAL *hhpi, REAL **jac) { /*---------------------------------------------------- ! FEB. 14, 1991 ! ! FORM Q^ FROM THE HOUSEHOLDER MATRICES STORED IN ! MATRICES A AND HHPI AND STORE IT IN JAC. !---------------------------------------------------*/ REAL tau; int i,j,k; REAL V1[10], V2[10]; for(i=1; i<=n; i++) for(j=1; j<=n; j++) jac[i][j] = ZERO; for(j=1; j<=n; j++) jac[j][j] = ONE; for (k=1; k 0) for (kk=1; kk<=ncc32; kk++) { k=k+32; sum=sum + amat[k-31][i]*amat[k-31][j]+amat[k-30][i]*amat[k-30][ j]+amat[k-29][i]*amat[k-29][j]+amat[k-28][i]*amat[k-28][j]+ amat[k-27][i]*amat[k-27][j]+amat[k-26][i]*amat[k-26][j]+ amat[k-25][i]*amat[k-25][j]+amat[k-24][i]*amat[k-24][j]; sum=sum + amat[k-23][i]*amat[k-23][j]+amat[k-22][i]*amat[k-22][ j]+amat[k-21][i]*amat[k-21][j]+amat[k-20][i]*amat[k-20][j]+ amat[k-19][i]*amat[k-19][j]+amat[k-18][i]*amat[k-18][j]+ amat[k-17][i]*amat[k-17][j]+amat[k-16][i]*amat[k-16][j]; sum=sum + amat[k-15][i]*amat[k-15][j]+amat[k-14][i]*amat[k-14][ j]+amat[k-13][i]*amat[k-13][j]+amat[k-12][i]*amat[k-12][j]+ amat[k-11][i]*amat[k-11][j]+amat[k-10][i]*amat[k-10][j]+ amat[k-9][i]*amat[k-9][j]+amat[k-8][i]*amat[k-8][j]; sum=sum + amat[k-7][i]*amat[k-7][j]+amat[k-6][i]*amat[k-6][j]+ amat[k-5][i]*amat[k-5][j]+amat[k-4][i]*amat[k-4][j]+amat[k-3][ i]*amat[k-3][j]+amat[k-2][i]*amat[k-2][j]+amat[k-1][i]* amat[k-1][j]+amat[k][i]*amat[k][j]; } if (ncc16 > 0) for (kk=1; kk<=ncc16; k++) { k=k+16; sum=sum + amat[k-15][i]*amat[k-15][j]+amat[k-14][i]*amat[k-14][ j]+amat[k-13][i]*amat[k-13][j]+amat[k-12][i]*amat[k-12][j]+ amat[k-11][i]*amat[k-11][j]+amat[k-10][i]*amat[k-10][j]+ amat[k-9][i]*amat[k-9][j]+amat[k-8][i]*amat[k-8][j]; sum=sum + amat[k-7][i]*amat[k-7][j]+amat[k-6][i]*amat[k-6][j]+ amat[k-5][i]*amat[k-5][j]+amat[k-4][i]*amat[k-4][j]+amat[k-3][ i]*amat[k-3][j]+amat[k-2][i]*amat[k-2][j]+amat[k-1][i]* amat[k-1][j]+amat[k][i]*amat[k][j]; } if (ncc8 > 0) for (kk=1; kk<=ncc8; kk++) { k=k+8; sum=sum + amat[k-7][i]*amat[k-7][j]+amat[k-6][i]*amat[k-6][j]+ amat[k-5][i]*amat[k-5][j]+amat[k-4][i]*amat[k-4][j]+amat[k-3][ i]*amat[k-3][j]+amat[k-2][i]*amat[k-2][j]+amat[k-1][i]* amat[k-1][j]+amat[k][i]*amat[k][j]; } if (ncc4 > 0) for (kk=1; kk<=ncc4; kk++) { k=k+4; sum=sum + amat[k-3][i]*amat[k-3][j]+amat[k-2][i]*amat[k-2][j] + amat[k-1][i]*amat[k-1][j]+amat[k][i]*amat[k][j]; } if (ncc4r > 0) for (kk=1; kk<=ncc4r; kk++) { k++; sum += amat[k][i]*amat[k][j]; } bmat[i][j]=sum; if (i != j) bmat[j][i]=bmat[i][j]; } // j loop } // i loop } // utumul() void rtrmul (int n, REAL **a, REAL **h, REAL *rdiag) { /*------------------------------------------------------- ! SEPT. 4, 1991 ! ! FIND R^R FOR QR-DECOMPOSED JACOBIAN. ! ! R IS STORED IN STRICT UPPER TRIANGLE OF A AND RDIAG. !------------------------------------------------------*/ int i; REAL wv1[10]; // TEMPORARILY REPLACE DIAGONAL OF R IN A (A IS RESTORED LATER). for (i=1; i<=n; i++) { wv1[i]=a[i][i]; a[i][i]=rdiag[i]; } utumul (n, n, n, n, n, n, a, h); for (i=1; i<=n; i++) a[i][i]=wv1[i]; } // rtrmul() void onenrm (bool *abort, bool *pertrb, int n, int nunit, int output, REAL epsmch, REAL *h1norm, REAL **h, REAL *scalex) { /*--------------------------------------------------------------------------- ! FEB. 23, 1992 ! ! FIND 1-NORM OF H MATRIX IF PERTURBATION IS DESIRED AND PERTURB DIAGONAL. !--------------------------------------------------------------------------*/ REAL sqrtep, temp; int i, j; char s[80]; extern FILE *fp; sqrtep=SQRT(epsmch); if (output > 4) { Line0(fp); Msg(fp," DIAGONAL OF MATRIX H (= JAC^JAC) BEFORE BEING PERTURBED:"); Line0(fp); for (i=1; i<=n; i++) { sprintf(s," H(%3d,%3d) = %12.3f", i, i, h[i][i]); Msg(fp,s); } } *h1norm=ZERO; for (j=1; j<=n; j++) *h1norm = *h1norm + ABS(h[1][j])/scalex[j]; *h1norm = *h1norm / scalex[1]; for (i=2; i<=n; i++) { temp=ZERO; for (j=1; j<=i; j++) temp += ABS(h[j][i])/scalex[j]; for (j=i+1; j<=n; j++) temp += ABS(h[i][j])/scalex[j]; *h1norm=MAX(*h1norm, temp/scalex[i]); } if (output > 4) { Line0(fp); sprintf(s," 1-NORM OF MATRIX H: %11.3f", *h1norm); Msg(fp,s); } if (*h1norm < epsmch) { if (output > 0) { Line1(fp); Line0(fp); Msg(fp," PROGRAM FAILS AS 1-NORM OF JACOBIAN IS TOO SMALL."); Line0(fp); Line1(fp); } *abort=TRUE; return; } else { // PERTURB DIAGONAL OF MATRIX H - USE THIS TO FIND "SN" *pertrb=TRUE; for (i=1; i<=n; i++) h[i][i] += SQRT(ONE*n)*sqrtep*(*h1norm)*scalex[i]*scalex[i]; if (output > 4) { Line0(fp); Line0(fp); Msg(fp," PERTURBED H MATRIX:"); matprt (n, n, h); } *abort = FALSE; } } // onenrm() void nstpun (bool *abort, bool *linesr, bool overch, bool *overfl, bool *qrsing, bool *sclfch, bool sclxch, int itnum, int maxexp, int n, int nunit, int output, REAL epsmch, REAL **a, REAL *delf, REAL *fvecc, REAL **h, REAL *hhpi, REAL **jac, REAL *rdiag, REAL *scalef, REAL *scalex, REAL *sn) { /*------------------------------------------------------------------------------- ! FEB. 23, 1992 ! ! THIS SUBROUTINE FINDS THE NEWTON STEP. ! ! IF THE JACOBIAN IS DETECTED AS SINGULAR OR IF THE ESTIMATED CONDITION ! NUMBER IS TOO HIGH (GREATER THAN EPSMCH^(-2/3) ) THEN H=J^J IS FORMED ! AND THE DIAGONAL IS PERTURBED BY ADDING SQRT(N*EPSMCH)*H1NORM*SCALEX(I)^2 ! TO THE CORRESPONDING ELEMENT. A CHOLESKY DECOMPOSITION IS PERFORMED ON ! THIS MODIFIED MATRIX PRODUCING A PSEUDO-NEWTON STEP. ! NOTE: THIS PROCEDURE MAY BE BE BYPASSED FOR ILL-CONDITIONED JACOBIANS ! BY SETTING THE LOGICAL VARIABLE BYPASS TO TRUE IN THE DRIVER. ! ! ABORT IF THE 1-NORM OF MATRIX H BECOMES TOO SMALL ! ALTHOUGH BUT NOT AT A SOLUTION ! BYPASS ALLOWS BYPASSING OF THE SPECIAL TREATMENT FOR ! BADLY CONDITIONED JACOBIANS ! QRSING INDICATES SINGULAR JACOBIAN DETECTED. !------------------------------------------------------------------------------*/ REAL connum, h1norm, maxadd, maxffl, scalfi, scalxj, sqrtep; REAL wv1[10]; bool pertrb; int i,j,newstm; char s[80]; extern FILE *fp; extern bool matsup; extern bool bypass; newstm = 0; *abort=FALSE; *overfl=FALSE; pertrb=FALSE; sqrtep=SQRT(epsmch); if (output > 3) { Line0(fp); Line0(fp); Msg(fp," SOLUTION OF LINEAR SYSTEM FOR NEWTON STEP, SN"); } // STORE (POSSIBLY SCALED) JACOBIAN IN MATRIX A if (!(*sclfch)) matcop (n, n, n, n, n, n, jac, a); else for (i=1; i<=n; i++) if (scalef[i] != ONE) { scalfi=scalef[i]; for (j=1; j<=n; j++) a[i][j]=jac[i][j]*scalfi; } else for (j=1; j<=n; j++) a[i][j]=jac[i][j]; // SCALED JACOBIAN IS PRINTED ONLY IF AT LEAST ONE SCALING // FACTOR IS NOT 1.0 if (output > 4 && (*sclfch)) { Line0(fp); Line0(fp); Msg(fp," SCALED JACOBIAN MATRIX:"); matprt (n, n, a); } // QR DECOMPOSITION OF (POSSIBLY SCALED) JACOBIAN qrdcom (qrsing, n, epsmch, a, hhpi, rdiag); // SAVE MATRIX A FOR BACK SUBSTITUTION TO CHECK DEUFLHARDS"S SECOND // STEP ACCEPTANCE CRITERION IN LINE SEARCH OR TRUST REGION METHOD matcop (n, n, n, n, n, n, a, h); if (output > 4 && n > 1 && (!matsup)) { Line0(fp); if (!(*sclfch)) Msg(fp," QR DECOMPOSITION OF JACOBIAN MATRIX:"); else Msg(fp," QR DECOMPOSITION OF SCALED JACOBIAN MATRIX:"); matprt (n, n, a); Line0(fp); Msg(fp," DIAGONAL OF R PI FACTORS FROM QR DECOMPOSITION"); Line0(fp); for (i=1; i 1) { if (sclxch) // SET UP FOR CONDITION NUMBER ESTIMATOR - SCALE R WRT X"S for(j=1; j<=n; j++) { if (scalex[j] != ONE) { scalxj=scalex[j]; rdiag[j] /= scalxj; for (i=1; i ONE/pow(sqrtep, 1.333)) || newstm == 77) { /* FORM H=(DF*JAC]^(DF*JAC) WHERE DF=DIAG(SCALEF). USE PREVIOUSLY COMPUTED QR DECOMPOSITION OF (SCALED) JACOBIAN WHERE R IS STORED IN THE UPPER TRIANGLE OF A AND RDIAG. */ if (overch) { *overfl=FALSE; ataov (overfl, maxexp, n, nunit, output, jac, h, scalef); } else rtrmul (n, a, h, rdiag); if (output > 3) { Line0(fp); if ((*qrsing) && (!(*overfl))) Msg(fp," SINGULAR JACOBIAN DETECTED: JACOBIAN PERTURBED."); else // NOTE: IF OVERFL IS TRUE THEN QRSING MUST BE TRUE if (*overfl) { Msg(fp," POTENTIAL OVERFLOW DETECTED IN CONDITION NUMBER ESTIMATOR."); Msg(fp," MATRIX ""ASSIGNED"" AS SINGULAR AND JACOBIAN PERTURBED."); } else { sprintf(s," CONDITION NUMBER TOO HIGH: %12.3f, JACOBIAN PERTURBED.", connum); Msg(fp,s); } } *overfl=FALSE; if (newstm != 77) { // FIND 1-NORM OF H MATRIX AND PERTURB DIAGONAL onenrm (abort, &pertrb, n, nunit, output, epsmch, &h1norm, h, scalex); if (*abort) return; } // CHOLESKY DECOMPOSITION OF H MATRIX - MAXFFL=0 IMPLIES // THAT H IS KNOWN TO BE POSITIVE DEFINITE. maxffl=ZERO; cholde (n, &maxadd, &maxffl, sqrtep, h, a); if (output > 4) { Line0(fp); Line0(fp); Msg(fp," CHOLESKY DECOMPOSITION OF H MATRIX:"); matprt (n, n, a); } /* FIND NEWTON STEP FROM CHOLESKY DECOMPOSITION. IF THE DIAGONAL HAS BEEN PERTURBED THEN THIS IS NOT THE ACTUAL NEWTON STEP BUT ONLY AN APPORXIMATION THEREOF. */ for (i=1; i<=n; i++) wv1[i]=-delf[i]; chsolv (overch, overfl, maxexp, n, nunit, output, a, wv1, sn); *overfl=FALSE; if (output > 3) if (!sclxch) { Line0(fp); if (!pertrb) Msg(fp," NEWTON STEP FROM CHOLESKY DECOMPOSITION"); else Msg(fp," APPROXIMATE NEWTON STEP FROM PERTURBED JACOBIAN"); Line0(fp); for(i=1; i<=n; i++) { sprintf(s," SN(%3d) = %12.3f", i, sn[i]); Msg(fp,s); } } else { Line0(fp); if (!pertrb) Msg(fp," NEWTON STEP FROM CHOLESKY DECOMPOSITION IN SCALED UNITS:"); else Msg(fp," APPROXIMATE NEWTON STEP FROM PERTURBED JACOBIAN IN SCALED UNITS:"); Line0(fp); for(i=1; i<=n; i++) { sprintf(s," SN(%3d) = %12.3f SN(%3d) = %12.3f", i, sn[i], i, scalex[i]*sn[i]); Msg(fp,s); } } /* SET QRSING TO TRUE SO THAT THE CORRECT MATRIX FACTORIZATION IS USED IN THE BACK-CALCULATION OF SBAR FOR DEUFLHARD RELAXATION FACTOR INITIALIZATION IN LINE SEARCH (ONLY MATTERS WHEN JACOBIAN IS ILL- CONDITIONED BUT NOT SINGULAR). */ *qrsing=TRUE; } else { if (output > 3 && n > 1) { if ((!bypass) && connum <= ONE/pow(sqrtep,1.333)) { Line0(fp); sprintf(s," CONDITION NUMBER ACCEPTABLE: %9.2f, JACOBIAN NOT PERTURBED.", connum); Msg(fp,s); } if ((bypass) && connum > ONE/pow(sqrtep,1.333)) { Line0(fp); sprintf(s," CONDITION NUMBER HIGH: %9.2f, JACOBIAN NOT PERTURBED AS", connum); Msg(fp,s); Msg(fp," BYPASS IS TRUE."); } } // NOTE: HERE SN STORES THE R.H.S. - IT IS OVERWRITTEN for(i=1; i<=n; i++) sn[i] = -fvecc[i]*scalef[i]; if ((!bypass) && (sclxch)) // R WAS SCALED BEFORE THE CONDITION NUMBER ESTIMATOR - // THIS CONVERTS IT BACK TO THE UNSCALED FORM. for(j=1; j<=n; j++) if (scalex[j] != ONE) { scalxj=scalex[j]; rdiag[j] *= scalxj; for (i= 1; i 3) if (!sclxch) { Line0(fp); Msg(fp," NEWTON STEP FROM QR DECOMPOSITION:"); Line0(fp); for(i=1; i<=n; i++) { sprintf(s," SN(%3d) = %12.3f", i, sn[i]); Msg(fp,s); } } else { Line0(fp); Msg(fp," NEWTON STEP FROM QR DECOMPOSITION IN SCALED UNITS:"); Line0(fp); for(i=1; i<=n; i++) { sprintf(s," SN(%3d) = %12.3f SN(%3d = %12.3f", i, sn[i], i, scalex[i]*sn[i]); Msg(fp,s); } } // TRANSFORM MATRICES FOR SUBSEQUENT CALCULATIONS IN TRUST // REGION METHODS (A IS STORED ABOVE IN H). if (!(*linesr)) for (i=1; i<=n; i++) { a[i][i]=rdiag[i]; for (j=1; j 3) { Line0(fp); Line0(fp); Msg(fp," SOLUTION OF LINEAR SYSTEM FOR NEWTON STEP, SN"); } // STORE (POSSIBLY SCALED) JACOBIAN IN MATRIX A if (!sclfch) matcop (n, n, n, n, n, n, jac, a); else for (i=1; i<=n; i++) if (scalef[i] != ONE) { scalfi=scalef[i]; for(j=1; j<=n; j++) a[i][j]=jac[i][j]*scalfi; } else for (j=1; j<=n; j++) a[i][j]=jac[i][j]; // SCALED JACOBIAN IS PRINTED ONLY IF AT LEAST ONE SCALING FACTOR IS NOT 1.0 if (output > 4 && (sclfch)) { Line0(fp); Line0(fp); Msg(fp," SCALED JACOBIAN MATRIX:"); matprt (n, n, a); } // QR DECOMPOSITION OF (POSSIBLY SCALED) JACOBIAN qrdcom (&qrsing, n, epsmch, a, hhpi, rdiag); if (output > 4 && n > 1 && (!matsup)) { Line0(fp); if (!sclfch) Msg(fp," QR DECOMPOSITION OF JACOBIAN MATRIX:"); else Msg(fp," QR DECOMPOSITION OF SCALED JACOBIAN MATRIX:"); matprt (n, n, a); Line0(fp); Msg(fp," DIAGONAL OF R PI FACTORS FROM QR DECOMPOSITION"); Line0(fp); for (i=1; i 1) { if (sclxch) // SET UP FOR CONDITION NUMBER ESTIMATION - SCALE R WRT X'S for (j=1; j<=n; j++) if (scalex[j] != ONE) { scalxj=scalex[j]; rdiag[j] /= scalxj; for (i=1; i ONE/pow(sqrtep,1.333)) || *newstm == 77) { /* FORM H=(DF*JAC)^(DF*JAC) WHERE DF=DIAG(SCALEF). USE PREVIOUSLY COMPUTED QR DECOMPOSITION OF (SCALED) JACOBIAN WHERE R IS STORED IN THE UPPER TRIANGLE OF A AND RDIAG. */ overfl=FALSE; if (overch) ataov (&overfl, maxexp, n, nunit, output, jac, h, scalef); else rtrmul (n, a, h, rdiag); if (output > 3) { Line0(fp); if ((qrsing) && (!overfl)) Msg(fp," SINGULAR JACOBIAN DETECTED: JACOBIAN PERTURBED."); else if (overfl) { Msg(fp," POTENTIAL OVERFLOW DETECTED IN CONDITION NUMBER ESTIMATOR,"); Msg(fp," MATRIX ""ASSIGNED"" AS SINGULAR AND JACOBIAN PERTURBED."); } else { sprintf(s," CONDITION NUMBER TOO HIGH: %12.3f, JACOBIAN PERTURBED.", connum); Msg(fp,s); } } overfl=FALSE; if (*newstm != 77) // FIND 1-NORM OF H MATRIX AND PERTURB DIAGONAL onenrm (abort, &pertrb, n, nunit, output, epsmch, &h1norm, h, scalex); // CHOLESKY DECOMPOSITION OF H MATRIX - MAXFFL=0 INDICATES // THAT H IS KNOWN TO BE POSITIVE DEFINITE. maxffl=ZERO; cholde (n, &maxadd, &maxffl, sqrtep, h, a); if (output > 4) { Line0(fp); Line0(fp); Msg(fp," CHOLESKY DECOMPOSITION OF H MATRIX:"); matprt (n, n, a); } /* FIND NEWTON STEP FROM CHOLESKY DECOMPOSITION. IF THE DIAGONAL HAS BEEN PERTURBED THEN THIS IS NOT THE ACTUAL NEWTON STEP BUT ONLY AN APPROXIMATION THEREOF. */ for (i=1; i<=n; i++) wv1[i] = -delf[i]; chsolv (overch, &overfl, maxexp, n, nunit, output, a, wv1, sn); overfl=FALSE; if (output > 3) if (!sclxch) { if (!pertrb) Msg(fp," NEWTON STEP FROM CHOLESKY DECOMPOSITION"); else Msg(fp," APPROXIMATE NEWTON STEP FROM PERTURBED JACOBIAN"); Line0(fp); for(i=1; i<=n; i++) { sprintf(s," SN(%3d) = %12.3f", i, sn[i]); Msg(fp,s); } } else { Line0(fp); if (! pertrb) Msg(fp," NEWTON STEP FROM CHOLESKY DECOMPOSITION IN SCALED UNITS"); else Msg(fp," APPROXIMATE NEWTON STEP FROM PERTURBED JACOBIAN IN SCALED UNITS"); Line0(fp); for(i=1; i<=n; i++) { sprintf(s," SN(%3d) = %12.3f SN(%3d) = %12.3f", i, sn[i], i, scalex[i]*sn[i]); Msg(fp,s); } } /* SET QRSING TO TRUE SO THAT THE CORRECT MATRIX FACTORIZATION IS USED IN THE BACK-CALCULATION OF SBAR FOR DEUFLHARD RELAXATION FACTOR INITIALIZATION. */ qrsing=TRUE; } else { if (output > 3 && n > 1) { if ((!bypass) && connum <= ONE/pow(sqrtep,1.333)) { Line0(fp); sprintf(s," CONDITION NUMBER ACCEPTABLE: %9.2f, JACOBIAN NOT PERTURBED.",connum); Msg(fp,s); } if ((bypass) && connum > ONE/pow(sqrtep,1.333)) { Line0(fp); sprintf(s," CONDITION NUMBER HIGH: %9.2f, JACOBIAN NOT PERTURBED AS", connum); Msg(fp,s); Msg(fp," BYPASS IS TRUE."); } } for (i=1; i<=n; i++) { sum=ZERO; for (j=1; j<=n; j++) sum -= jac[i][j]*scalef[j]*fvecc[j]; sn[i]=sum; } rsolv (overch, &overfl, maxexp, n, nunit, output, a, rdiag, sn); overfl=FALSE; if (output > 3) { if (!sclxch) { Line0(fp); Msg(fp," NEWTON STEP FROM QR DECOMPOSITION"); Line0(fp); for (i=1; i<=n; i++) { sprintf(s," SN(%3d) = %12.3f", i, sn[i]); Msg(fp,s); } } else { Line0(fp); Msg(fp," NEWTON STEP FROM QR DECOMPOSITION IN SCALED UNITS"); Line0(fp); for (i=1; i<=n; i++) { sprintf(s," SN(%3d) = %12.3f SN(%3d) = %12.3f", i, sn[i], i, scalex[i]*sn[i]); Msg(fp,s); } } } // TRANSFORM MATRICES FOR SUBSEQUENT CALCULATIONS IN TRUST REGION METHOD if (!linesr) for (i=2; i<=n; i++) for (j=1; j 0) { sprintf(s," ANALYTICAL JACOBIAN USED, CHECKED NUMERICALLY, NJACCH: %5d", njacch); Msg(fp,s); } else Msg(fp," ANALYTICAL JACOBIAN USED; NOT CHECKED."); else if (jactyp == 1) Msg(fp," JACOBIAN ESTIMATED USING FORWARD DIFFERENCES."); else if (jactyp == 2) Msg(fp," JACOBIAN ESTIMATED USING BACKWARD DIFFERENCES."); else Msg(fp," JACOBIAN ESTIMATED USING CENTRAL DIFFERENCES."); Line0(fp); Line1(fp); } else { if (linesr) { Line0(fp); if (deuflh) Msg(fp," DEUFLHARD RELAXATION FACTOR INITIALIZATION IN EFFECT."); else Msg(fp," DEUFLHARD RELAXATION FACTOR INITIALIZATION NOT IN EFFECT."); } else { if (etafac == ONE) Msg(fp," METHOD: TRUST REGION USING SINGLE DOGLEG STEPS."); else Msg(fp," METHOD: TRUST REGION USING DOUBLE DOGLEG STEPS."); Line0(fp); if (cauchy) Msg(fp," INITIAL STEP CONSTRAINED BY SCALED CAUCHY STEP."); else Msg(fp," INITIAL STEP CONSTRAINED BY SCALED NEWTON STEP."); } if (geoms) Msg(fp," METHOD: GEOMETRIC SEARCH."); else Msg(fp," METHOD: SEARCH BASED ON SUCCESSIVE MINIMIZATIONS."); if (overch) Msg(fp," OVERLOW CHECKING IN USE."); else Msg(fp," OVERFLOW CHECKING NOT IN USE."); if (jupdm == 0) Msg(fp," NO QUASI-NEWTON UPDATE USED."); if (jupdm == 1) if (qnupdm == 0) Msg(fp," BROYDEN QUASI-NEWTON UPDATE OF UNFACTORED JACOBIAN."); else Msg(fp," BROYDEN QUASI-NEWTON UPDATE OF FACTORED JACOBIAN."); else if (jupdm == 2) if (qnupdm == 0) Msg(fp," LEE AND LEE QUASI-NEWTON UPDATE OF UNFACTORED JACOBIAN."); else Msg(fp," LEE AND LEE QUASI-NEWTON UPDATE OF FACTORED JACOBIAN."); if (jactyp == 0) if (njacch > 0) { sprintf(s," ANALYTICAL JACOBIAN USED, CHECKED NUMERICALLY, NJACCH: %5d", njacch); Msg(fp,s); } else Msg(fp," ANALYTICAL JACOBIAN USED; NOT CHECKED."); else if (jactyp == 1) Msg(fp," JACOBIAN ESTIMATED USING FORWARD DIFFERENCES."); else if (jactyp == 2) Msg(fp," JACOBIAN ESTIMATED USING BACKWARD DIFFERENCES."); else Msg(fp," JACOBIAN ESTIMATED USING CENTRAL DIFFERENCES."); if (!linesr) { Line0(fp); if (trupdm == 0 && jupdm > 0) Msg(fp," TRUST REGION UPDATED USING POWELL STRATEGY."); else Msg(fp," TRUST REGION UPDATED USING DENNIS AND SCHNABEL STRATEGY."); } Line0(fp); Line1(fp); Line0(fp); if (itsclf != 0) { sprintf(s," ADAPTIVE FUNCTION SCALING STARTED AT ITERATION: ..........%6d", itsclf); Msg(fp,s); Line0(fp); } if (itsclx != 0) { sprintf(s," ADAPTIVE VARIABLE SCALING STARTED AT ITERATION: ..........%6d", itsclx); Msg(fp,s); Line0(fp); } if (linesr) if (jupdm == 0) { sprintf(s," MAXIMUM NUMBER OF STEPS IN LINE SEARCH, MAXNS: ...........%6d", maxns); Msg(fp,s); } else { sprintf(s," MAXIMUM NUMBER OF NEWTON LINE SEARCH STEPS, MAXNS: .......%6d", maxns); Msg(fp,s); sprintf(s," MAXIMUM NUMBER OF QUASI-NEWTON LINE SEARCH STEPS, MAXQNS: %6d", maxqns); Msg(fp,s); } else if (jupdm == 0) { sprintf(s," MAXIMUM NUMBER OF TRUST REGION UPDATES, MAXNS: ...........%6d", maxns); Msg(fp,s); } else { sprintf(s," MAXIMUM NO. OF NEWTON TRUST REGION UPDATES, MAXNS: .......%6d", maxns); Msg(fp,s); sprintf(s," MAXIMUM NO. OF QUASI-NEWTON TRUST REGION UPDATES, MAXQNS: %6d", maxqns); Msg(fp,s); } if (narmij < maxit) { Line0(fp); sprintf(s," NUMBER OF OBJECTIVE FUNCTION VALUES COMPARED, MGLL: ......%6d", mgll); Msg(fp,s); } if (jupdm > 0) { if (narmij == maxit) Line0(fp); sprintf(s," MINIMUM NUMBER OF STEPS BETWEEN JACOBIAN UPDATES, MINQNS: %6d", minqns); Msg(fp,s); sprintf(s," NUMBER OF NON-QUASI-NEWTON STEPS AT START, NINITN: .......%6d", ninitn); Msg(fp,s); } sprintf(s," NUMBER OF ARMIJO STEPS AT START, NARMIJ: .................%6d", narmij); Msg(fp,s); } Line0(fp); if (stopcr == 3) { sprintf(s," FUNCTION AND STEP SIZE STOPPING CRITERIA, STOPCR: ........%6d", stopcr); Msg(fp,s); } else if (stopcr == 12) { sprintf(s," FUNCTION OR STEP SIZE STOPPING CRITERIA, STOPCR: .........%6d", stopcr); Msg(fp,s); } else if (stopcr == 1) { sprintf(s," STEP SIZE STOPPING CRITERION, STOPCR: ....................%6d", stopcr); Msg(fp,s); } else { sprintf(s," FUNCTION STOPPING CRITERION, STOPCR: .....................%6d", stopcr); Msg(fp,s); } if (!newton) { Line0(fp); if (acptcr == 12) { sprintf(s," FUNCTION AND STEP SIZE ACCEPTANCE CRITERIA, ACPTCR: ......%6d", acptcr); Msg(fp,s); } else if (acptcr == 2) { sprintf(s," STEP SIZE ACCEPTANCE CRITERION, ACPTCR: ..................%6d", acptcr); Msg(fp,s); } else { sprintf(s," FUNCTION ACCEPTANCE CRITERION, ACPTCR: ...................%6d", acptcr); Msg(fp,s); } if (contyp != 0) { Line0(fp); sprintf(s," CONSTRAINTS IN USE, CONTYP: ..............................%6d", contyp); Msg(fp,s); } } Line0(fp); Line1(fp); Line0(fp); sprintf(s," ESTIMATED MACHINE EPSILON, EPSMCH: ......%e", epsmch); Msg(fp,s); Line0(fp); sprintf(s," FACTOR TO ESTABLISH MAXIMUM STEP SIZE, MSTPF : .......%10.3f", mstpf); Msg(fp,s); sprintf(s," CALCULATED MAXIMUM STEP SIZE, MAXSTP: ................%10.3f", maxstp); Msg(fp,s); if (!linesr) { if (delta < ZERO) { Line0(fp); Msg(fp," INITIAL TRUST REGION NOT PROVIDED."); } else { Line0(fp); sprintf(s," INITIAL TRUST REGION SIZE, DELTA: ....................%10.3f", delta); Msg(fp,s); } if (etafac < ONE) { sprintf(s," FACTOR TO SET DIRECTION OF TRUST REGION STEP, ETAFAC: ... %6.4f", etafac); Msg(fp,s); } sprintf(s," TRUST REGION UPDATING FACTOR, DELFAC: ................%10.3f", delfac); Msg(fp,s); } if (!newton) { Line0(fp); sprintf(s," FACTOR IN OBJECTIVE FUNCTION COMPARISON, ALPHA: ......%10.5f", alpha); Msg(fp,s); if (linesr && (!newton)) { Line0(fp); sprintf(s," REDUCTION FACTOR FOR RELAXATION FACTOR, SIGMA: .......%10.3f", sigma); Msg(fp,s); } if (jupdm != 0) { Line0(fp); sprintf(s," REDUCTION REQUIRED IN OBJ. FUNCTION FOR QN STEP, RATIOF: %6.4f", ratiof); Msg(fp,s); } if (jupdm == 2) { Line0(fp); sprintf(s," FACTOR IN LEE AND LEE UPDATE, OMEGA: .....................%6.4f", omega); Msg(fp,s); } } Line0(fp); if (stopcr != 2) { sprintf(s," STOPPING TOLERANCE FOR STEP SIZE, STPTOL (E-10): .....%10.3f", 1e10*stptol); Msg(fp,s); sprintf(s," STOPPING TOLERANCE FOR NEWTON STEP, NSTTOL (E-10): ...%10.3f", 1e10*nsttol); Msg(fp,s); } if (stopcr != 1) { sprintf(s," STOPPING TOLERANCE FOR OBJECTIVE FUNCTION, FTOL(E-5)..%10.3f", 1e5*ftol); Msg(fp,s); } if ((linesr) && (!newton) && lam0 < ONE) { Line0(fp); sprintf(s," INITIAL LAMBDA IN LINE SEARCH, LAM0: .................%10.3f", lam0); Msg(fp,s); } if (contyp > 0) { Line0(fp); sprintf(s," FACTOR TO ENSURE STEP WITHIN CONSTRAINTS, CONFAC: ........%6.4f", confac); Msg(fp,s); } Line0(fp); Line1(fp); Line0(fp); Msg(fp," SCALING FACTORS"); Line0(fp); Msg(fp," COMPONENT VALUES FUNCTION VALUES"); Line0(fp); for(i=1; i<=n; i++) { sprintf(s," SCALEX(%3d) = %10.3f SCALEF(%3d) = %10.3f", i, scalex[i], i, scalef[i]); Msg(fp,s); } if (contyp > 0) { Line0(fp); Line1(fp); Line0(fp); Msg(fp," LOWER AND UPPER BOUNDS"); Line0(fp); Msg(fp," LOWER BOUNDS UPPER BOUNDS"); Line0(fp); for(i=1; i<=n; i++) { if (ABS(boundl[i])<100000 && ABS(boundu[i])<100000) sprintf(s," BOUNDL(%3d) = %12.3f BOUNDU(%3d) = %12.3f", i, boundl[i], i, boundu[i]); else sprintf(s," BOUNDL(%3d) = %e BOUNDU(%3d) = %e", i, boundl[i], i, boundu[i]); Msg(fp,s); } } Line0(fp); e100:if (output == 2) Line1(fp); Line1(fp); Line0(fp); Msg(fp," INITIAL ESTIMATES INITIAL FUNCTION VALUES"); Line0(fp); for(i=1; i<=n; i++) { if (ABS(fvecc[i]) < 100000) sprintf(s," X(%3d) = %12.3f F(%3d) = %14.3f", i, xc[i], i, fvecc[i]); else sprintf(s," X(%3d) = %12.3f F(%3d) = %e", i, xc[i], i, fvecc[i]); Msg(fp,s); } Line0(fp); if (ABS(fcnold) < 100000) sprintf(s," INITIAL OBJECTIVE FUNCTION VALUE = %10.3f", fcnold); else sprintf(s," INITIAL OBJECTIVE FUNCTION VALUE = %e", fcnold); Msg(fp,s); Line0(fp); Line1(fp); Line1(fp); } //title() void qmin (int nunit, int output, REAL delfts, REAL *delta, REAL deltaf, REAL stplen) { /*---------------------------------------------------------------------- ! FEB. 9, 1991 ! ! SET THE NEW TRUST REGION SIZE, DELTA, BASED ON A QUADRATIC ! MINIMIZATION WHERE DELTA IS THE INDEPENDENT VARIABLE. ! ! DELTAF IS THE DIFFERENCE IN THE SUM-OF-SQUARES OBJECTIVE ! FUNCTION VALUE AND DELFTS IS THE DIRECTIONAL DERIVATIVE IN ! THE DIRECTION OF THE CURRENT STEP, S, WHICH HAS STEP LENGTH STPLEN. !---------------------------------------------------------------------*/ #define point1 0.1 REAL deltmp; char s[80]; extern FILE *fp; if (deltaf-delfts != ZERO) { // CALCULATE DELTA WHERE MINIMUM WOULD OCCUR - DELTMP. // THIS IS PROVISIONAL AS IT MUST BE WITHIN CERTAIN LIMITS TO BE ACCEPTED deltmp=-delfts*stplen/(TWO*(deltaf-delfts)); if (output > 4) { Line0(fp); sprintf(s," TEMPORARY DELTA FROM QUADRATIC MINIMIZATION: %12.3f", deltmp); Msg(fp,s); sprintf(s," VERSUS CURRENT DELTA: %12.3f", *delta); Msg(fp,s); } /* REDUCE DELTA DEPENDING ON THE MAGNITUDE OF DELTMP. IT MUST BE WITHIN (0.1*DELTA,0.5*DELTA) TO BE ACCEPTED - OTHERWISE THE NEAREST ENDPOINT OF THE INTERVAL IS USED */ if (deltmp < point1*(*delta)) { *delta=point1*(*delta); if (output > 4) { Line0(fp); Msg(fp," NEW DELTA SET TO 0.1 CURRENT DELTA"); } } else if (deltmp > *delta/TWO) { *delta=*delta/TWO; if (output > 4) { Line0(fp); Msg(fp," NEW DELTA SET TO 0.5 CURRENT DELTA"); } } else { *delta=deltmp; if (output > 4) { Line0(fp); Msg(fp," NEW DELTA SET TO DELTMP"); } } } else { if (output > 4) { Line0(fp); Msg(fp," TO AVOID OVERFLOW NEW DELTA SET TO 0.5 CURRENT DELTA."); } *delta=*delta/TWO; } } // qmin() void trstup (bool geoms, bool *newtkn, bool overch, bool *overfl, bool qrsing, bool *sclfch, bool *sclxch, int *acpcod, int *acpstr, int acptcr, int contyp, int isejac, int jupdm, int *maxexp, int mgll, int mnew, int n, int narmij, int *nfunc, int notrst, int nunit, int output, int qnupdm, int *retcod, int *trupdm, REAL alpha, REAL confac, REAL delfac, REAL *delstr, REAL *delta, REAL epsmch, REAL *fcnmax, REAL *fcnnew, REAL fcnold, REAL *fcnpre, REAL *maxstp, REAL newlen, REAL *newmax, REAL *powtau, REAL *rellen, REAL stptol, REAL **a, REAL **astore, REAL *boundl, REAL *boundu, REAL *delf, REAL *fplpre, REAL *ftrack, REAL *fvec, REAL *fvecc, REAL *hhpi, REAL **jac, REAL *rdiag, REAL *rhs, REAL *s, REAL *sbar, REAL *scalef, REAL *scalex, REAL *strack, REAL *xc, REAL *xplpre, REAL *xplus) { /*------------------------------------------------------------ ! FEB. 28, 1992 ! ! THIS SUBROUTINE CHECKS FOR ACCEPTANCE OF A TRUST REGION ! STEP GENERATED BY THE DOUBLE DOGLEG METHOD. !-----------------------------------------------------------*/ // Label: e10 #define pt5 0.5 #define threeq 0.75 #define onept1 1.1 bool convio; int i, j, k; REAL delfpr, delfts, deltaf, dlftsm, dmult, powlam, powmu, ratio, ratiom, sbrnrm, sp, ss, stplen, sum; REAL *V1, *wv3; REAL **temp1, **temp2; void *vmblock = NULL; char s1[80]; extern FILE *fp; vmblock = vminit(); V1 = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); wv3 = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); temp1 = (REAL **) vmalloc(vmblock, MATRIX, n+1, n+1); temp2 = (REAL **) vmalloc(vmblock, MATRIX, n+1, n+1); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return; } // NOTE: ACCEPTANCE CODE, ACPCOD, IS 0 ON ENTRANCE TO TRSTUP convio=FALSE; *overfl=FALSE; *retcod = 0; if (output > 3) { Line0(fp); Line0(fp); if ((!(*sclfch)) && (!(*sclxch))) Msg(fp," TRUST REGION UPDATING"); else Msg(fp," TRUST REGION UPDATING (ALL X""S AND F""S IN UNSCALED UNITS)"); Line0(fp); } // CHECK TO MAKE SURE "S" IS A DESCENT DIRECTION - FIND DIRECTIONAL // DERIVATIVE AT CURRENT XC USING S GENERATED BY DOGLEG SUBROUTINE. innerp (overch, overfl, *maxexp, n, n, n, nunit, output, &delfts, delf, s); if (output > 3) { Line0(fp); sprintf(s1," INNER PRODUCT OF DELF AND S, DELFTS: ........%13.4f", delfts); Msg(fp,s1); } if (delfts > ZERO) { if (output > 3) { Line0(fp); Msg(fp," DIRECTIONAL DERIVATIVE POSITIVE, SEARCH DIRECTION REVERSED."); } for (i=1; i<=n; i++) s[i] = -s[i]; } /* FIND MAXIMUM OBJECTIVE FUNCTION VALUE AND MAXIMIUM STEP LENGTH FOR NONMONOTONIC SEARCH. THIS HAS TO BE DONE ONLY ONCE DURING EACH ITERATION (WHERE NOTRST=1). */ if (notrst == 1) { *newmax=newlen; *fcnmax=fcnold; if (isejac > narmij) if (isejac < narmij+mgll) for (j=1; j<=mnew; j++) { *fcnmax=MAX(*fcnmax, ftrack[j-1]); *newmax=MAX(*newmax, strack[j-1]); } else for (j=0; j<=mnew; j++) { *fcnmax=MAX(*fcnmax, ftrack[j]); *newmax=MAX(*newmax, strack[j]); } } // TEST TRIAL POINT - FIND XPLUS AND TEST FOR CONSTRAINT // VIOLATIONS IF CONTYP DOES NOT EQUAL 0. for (i=1; i<=n; i++) { wv3[i]=-ONE; // WV3 IS A MARKER FOR "VIOLATORS" - IT CHANGES TO 1.0 OR 2.0 xplus[i]=xc[i] + s[i]; if (contyp > 0) if (xplus[i] < boundl[i]) { convio=TRUE; wv3[i]=ONE; } else if (xplus[i] > boundu[i]) { convio=TRUE; wv3[i]=TWO; } } // IF CONSTRAINT IS VIOLATED... if (convio) { if (output > 3) { Line0(fp); Msg(fp," CONSTRAINT VIOLATED:"); Line0(fp); Msg(fp," TRIAL ESTIMATES (VIOLATIONS MARKED):"); Line0(fp); for (i=1; i<=n; i++) if (wv3[i] > ZERO) { sprintf(s1," XPLUS(%3d) = %12.3f *", i, xplus[i]); Msg(fp,s1); } else { sprintf(s1," XPLUS(%3d) = %12.3f", i, xplus[i]); Msg(fp,s1); } } /* FIND STEP WITHIN CONSTRAINED REGION. FIND THE RATIO OF THE DISTANCE FROM THE (I)TH COMPONENT TO ITS CONSTRAINT TO THE LENGTH OF THE PROPOSED STEP, XPLUS(I)-XC(I). MULTIPLY THIS BY CONFAC (DEFAULT 0.95) TO ENSURE THE NEW STEP STAYS WITHIN THE ACCEPTABLE REGION UNLESS XC IS CLOSE TO THE BOUNDARY (RATIO <= 1/2). IN SUCH CASES A FACTOR OF 0.5*CONFAC IS USED. NOTE: ONLY THE VIOLATING COMPONENTS ARE REDUCED. */ ratiom=ONE; // RATIOM STORES THE MINIMUM VALUE OF RATIO for (i=1; i<=n; i++) { if (wv3[i] == ONE) ratio=(boundl[i]-xc[i])/s[i]; else if (wv3[i] == TWO) ratio=(boundu[i]-xc[i])/s[i]; if (wv3[i] > ZERO) { // NOTE: RATIO IS STORED IN WV3 FOR OUTPUT ONLY wv3[i]=ratio; ratiom=MIN(ratiom,ratio); if (ratio > pt5) xplus[i]=xc[i] + confac*ratio*s[i]; else // WITHIN BUFFER ZONE xplus[i]=xc[i] + confac*ratio*s[i]/TWO; s[i]=xplus[i] - xc[i]; } } if (output > 3) { Line0(fp); Line0(fp); Msg(fp," NEW S AND XPLUS VECTORS (WITH RATIOS FOR VIOLATIONS)"); Line0(fp); Msg(fp," NOTE: RATIOS ARE RATIO OF LENGTH TO BOUNDARY FROM CURRENT"); Msg(fp," X VECTOR TO MAGNITUDE OF CORRESPONDING PROPOSED STEP."); Line0(fp); for (i=1; i<=n; i++) if (wv3[i] < ZERO ) { sprintf(s1," S(%3d) = %12.3f XPLUS(%3d) = %12.3f", i, s[i], i, xplus[i]); Msg(fp,s1); } else { sprintf(s1," S(%3d) = %12.3f XPLUS(%3d) = %12.3f %9.3f", i, s[i], i, xplus[i], wv3[i]); Msg(fp,s1); } Line0(fp); sprintf(s1," MINIMUM OF RATIOS, RATIOM: %12.3f", ratiom); Msg(fp,s1); } // THE NEW POINT, XPLUS, IS NOT NECESSARILY IN A DESCENT DIRECTION. // CHECK DIRECTIONAL DERIVATIVE FOR MODIFIED STEP, DLFTSM. innerp (overch, overfl, *maxexp, n, n, n, nunit, output, &dlftsm, delf, s); if (output > 3) { Line0(fp); sprintf(s1," INNER PRODUCT OF DELF AND MODIFIED S, DL FTSM: %12.3f", dlftsm); Msg(fp,s1); } // IF DLFTSM IS POSITIVE REDUCE TRUST REGION. IF NOT, TEST NEW POINT if (dlftsm > ZERO) { *delta=confac*ratiom*(*delta); *retcod=8; goto e10; } } // if convio // CONSTRAINTS NOT (OR NO LONGER) VIOLATED - TEST NEW POINT fcn1 (overfl, n, fvec, xplus); *nfunc = *nfunc + 1; // IF OVERFLOW AT NEW POINT REDUCE TRUST REGION AND RETURN if (*overfl) { /* IF THE OVERFLOW COMES AS A RESULT OF INCREASING DELTA WITHIN THE CURRENT ITERATION (IMPLYING DELSTR IS POSITIVE) AND DIVIDING DELTA BY 10 WOULD PRODUCE A DELTA WHICH IS SMALLER THAN THAT AT THE STORED POINT, THEN USE STORED POINT AS THE UPDATED RESULT. */ if (*delstr > *delta/TEN) { for (i=1; i<=n; i++) { temp1[i][1]=xplpre[i]; temp2[i][1]=xplus[i]; } matcop (n, n, 1, 1, n, 1, temp1, temp2); for(i=1; i<=n; i++) xplus[i]=temp2[i][1]; for (i=1; i<=n; i++) { temp1[i][1]=fplpre[i]; temp2[i][1]=fvec[i]; } matcop (n, n, 1, 1, n, 1, temp1, temp2); for(i=1; i<=n; i++) fvec[i]=temp2[i][1]; acpcod=acpstr; *delta=*delstr; fcnnew=fcnpre; *retcod=1; } else { *delta=*delta/TEN; *retcod=9; } goto e10; } // NO OVERFLOW IN RESIDUAL VECTOR if (output > 3) { Line0(fp); Msg(fp," TRIAL ESTIMATES FUNCTION VALUES"); Line0(fp); for(i=1; i<=n; i++) { sprintf(s1," XPLUS(%3d) = %12.3f FVEC(%3d) = %12.3f", i, xplus[i], i, fvec[i]); Msg(fp,s1); } } // IF NO OVERFLOW WITHIN RESIDUAL VECTOR FIND OBJECTIVE FUNCTION. fcnevl (overfl, *maxexp, n, nunit, output, epsmch, fcnnew, fvec, scalef); // IF OVERFLOW IN OBJECTIVE FUNCTION EVALUATION REDUCE TRUST REGION AND RETURN. if (*overfl) { /* IF THE OVERFLOW COMES AS A RESULT OF INCREASING DELTA WITHIN THE CURRENT ITERATION (SO THAT DELSTR IS POSITIVE) AND DIVIDING DELTA BY 10 WOULD PRODUCE A DELTA WHICH IS SMALLER THAN THAT AT THE STORED POINT THEN USE STORED POINT AS THE UPDATED RESULT. */ if (*delstr > *delta/TEN) { for(i=1; i<=n; i++) { temp1[i][1]=xplpre[i]; temp2[i][1]=xplus[i]; } matcop (n, n, 1, 1, n, 1, temp1, temp2); for(i=1; i<=n; i++) xplus[i]=temp2[i][1]; for(i=1; i<=n; i++) { temp1[i][1]=fplpre[i]; temp2[i][1]=fvec[i]; } matcop (n, n, 1, 1, n, 1, temp1, temp2); for(i=1; i<=n; i++) fvec[i]=temp2[i][1]; acpcod=acpstr; *delta=*delstr; fcnnew=fcnpre; *retcod=2; } else { *delta=*delta/TEN; *retcod=10; } goto e10; } else { // NO OVERFLOW AT TRIAL POINT - COMPARE OBJECTIVE FUNCTION TO FCNMAX if (output > 3) { Line0(fp); if (!*sclfch) { sprintf(s1," OBJECTIVE FUNCTION AT XPLUS, FCNNEW: .........%12.4f", *fcnnew); Msg(fp,s1); } else { sprintf(s1," SCALED OBJECTIVE FUNCTION AT XPLUS, FCNNEW: ..%12.4f", *fcnnew); Msg(fp,s1); } sprintf(s1," COMPARE TO FCNMAX+ALPHA*DELFTS: ..............%12.4f", (*fcnmax + alpha*delfts)); Msg(fp,s1); } } /* IF ACPTCR=12 CHECK SECOND DEUFLHARD STEP ACCEPTANCE TEST BY FINDING 2-NORM OF SBAR. THERE ARE FOUR POSSIBILITIES DEPENDING ON WHETHER THE JACOBIAN IS OR IS NOT SINGULAR AND WHETHER QNUPDM IS 0 OR 1 */ if (acptcr == 12) { if (qrsing) { // FORM -J^F AS RIGHT HAND SIDE - METHOD DEPENDS ON // WHETHER QNUPDM EQUALS 0 OR 1 (UNFACTORED OR FACTORED) if (qnupdm == 0) { // UNSCALED JACOBIAN IS IN MATRIX JAC for (k=1; k<=n; k++) wv3[k] = -fvec[k]*scalef[k]*scalef[k]; mmulv(n, jac, wv3, rhs); } else { // R IN UPPER TRIANGLE OF A PLUS RDIAG AND Q^ IN JAC // FROM QR DECOMPOSITION OF SCALED JACOBIAN. for (i=1; i<=n; i++) { sum=ZERO; for (j=1; j<=n; j++) sum -= jac[i][j]*fvec[j]*scalef[j]; wv3[i]=sum; } rhs[1]=rdiag[1]*wv3[1]; for (j=2; j<=n; j++) { for (k=1; k 4) { Line0(fp); if (!(*sclxch)) { Msg(fp," DEUFLHARD SBAR VECTOR"); Line0(fp); for (i=1; i<=n; i++) { sprintf(s1," SBAR(%3d) = %12.3f", i, sbar[i]); Msg(fp,s1); } } else { Msg(fp," DEUFLHARD SBAR VECTOR IN SCALED X UNITS"); Line0(fp); for (i=1; i<=n; i++) { sprintf(s1," SBAR(%3d) =%12.3f SBAR(%3d) = %12.3f", i, sbar[i], i, scalex[i]*sbar[i]); Msg(fp,s1); } } } if (output > 3) { Line0(fp); if (!(*sclxch)) { sprintf(s1," VALUE OF SBRNRM AT XPLUS: .................%12.4f", sbrnrm); Msg(fp,s1); } else { sprintf(s1," VALUE OF SCALED SBRNRM AT XPLUS: ..........%12.4f", sbrnrm); Msg(fp,s1); } sprintf(s1," NEWMAX: ...................................%12.4f", *newmax); Msg(fp,s1); } if (sbrnrm < *newmax) *acpcod=2; // FUNCTION VALUE ACCEPTANCE IS ALSO CHECKED REGARDLESS // OF WHETHER SECOND STEP ACCEPTANCE CRITERION WAS MET } // if acptor=12 // ESTABLISH DELTAF FOR USE IN COMPARISON TO PREDICTED // CHANGE IN OBJECTIVE FUNCTION, DELFPR, LATER deltaf = *fcnnew - fcnold; if (*fcnnew >= *fcnmax + alpha*delfts) { // FAILURE OF FIRST STEP ACCEPTANCE TEST. TEST LENGTH OF // STEP TO ENSURE PROGRESS IS STILL BEING MADE *rellen=ZERO; for (i=1; i<=n; i++) *rellen=MAX(*rellen, ABS(s[i])/MAX((ABS(xplus[i])),ONE/scalex[i])); if (*rellen < stptol) { // NO PROGRESS BEING MADE - RETCOD = 7 STOPS PROGRAM for(i=1; i<=n; i++) { temp1[i][1]=xc[i]; temp2[i][1]=xplus[i]; } matcop (n, n, 1, 1, n, 1, temp1, temp2); for (i=1; i<=n; i++) xplus[i]=temp2[i][1]; *retcod=7; goto e10; } else { /* FAILURE OF STEP BY OBJECTIVE FUNCTION CRITERION. ESTABLISH A NEW DELTA FROM EITHER SIMPLE DIVISION BY DELFAC OR BY FINDING THE MINIMUM OF A QUADRATIC MODEL */ if (geoms) *delta=*delta/delfac; else { // FIRST FIND LENGTH OF TRUST REGION STEP for (i=1; i<=n; i++) wv3[i]=s[i]*scalex[i]; twonrm (overfl, *maxexp, n, epsmch, &stplen, wv3); qmin (nunit, output, delfts, delta, deltaf, stplen); } if (*delta < *delstr) { /* IF DELTA HAS BEEN INCREASED AT THIS ITERATION AND THE DELTA FROM QMIN IS LESS THAN THE DELTA AT THE PREVIOUSLY ACCEPTED (STORED) POINT THEN RETURN TO THAT POINT AND ACCEPT IT AS THE UPDATED ITERATE. */ for(i=1; i<=n; i++) { temp1[i][1]=xplpre[i]; temp2[i][1]=xplus[i]; } matcop (n, n, 1, 1, n, 1, temp1, temp2); for(i=1; i<=n; i++) xplus[i]=temp2[i][1]; for(i=1; i<=n; i++) { temp1[i][1]=fplpre[i]; temp2[i][1]=fvec[i]; } matcop (n, n, 1, 1, n, 1, temp1, temp2); for(i=1; i<=n; i++) fvec[i]=temp2[i][1]; *acpcod=*acpstr; *delta=*delstr; *fcnnew=*fcnpre; *retcod=3; goto e10; } /* IF THE SECOND ACCEPTANCE TEST HAS BEEN PASSED RETURN WITH NEW TRUST REGION AND CONTINUE ON TO NEXT ITERATION; OTHERWISE TRY A NEW STEP WITH REDUCED DELTA. */ if (*acpcod == 2) *retcod=4; else // FAILURE OF FIRST STEP ACCEPTANCE TEST *retcod=11; goto e10; } } else { /* OBJECTIVE FUNCTION MEETS FIRST ACCEPTANCE CRITERION. IN NONMONOTONIC SEARCHES IT MAY BE GREATER THAN THE PREVIOUS OBJECTIVE FUNCTION VALUE - CONSIDER THIS CASE FIRST. */ if (deltaf >= alpha*delfts) { /* AN ACCEPTABLE STEP HAS BEEN FOUND FOR THE NONMONOTONIC SEARCH BUT THE OBJECTIVE FUNCTION VALUE IS NOT A "DECREASE" FROM THE PREVIOUS ITERATION (ACTUALLY IT MIGHT BE BETWEEN ZERO AND ALPHA*DELFTS). ACCEPT STEP BUT REDUCE DELTA. */ *delta=*delta/delfac; *retcod=5; if (*acpcod == 2) *acpcod=12; else *acpcod=1; goto e10; } /* COMPARE DELTAF TO DELTAF PREDICTED, DELFPR, TO DETERMINE NEXT TRUST REGION SIZE. NOTE: DELTAF MUST BE LESS THAN ALPHA*DELFTS (IN ESSENCE NEGATIVE) TO HAVE REACHED THIS POINT IN TRSTUP. R IS IN UPPER TRIANGLE OF MATRIX A SO THE FOLLOWING CODE FINDS: DELFPR = DELF^S + 1/2 S^J^JS = DELF^S + 1/2 S^R^RS */ uvmul (n, n, n, n, a, s, wv3); delfpr = delfts + Dot_Product(n, wv3, wv3) / TWO; if (output > 4) { Line0(fp); sprintf(s1," PREDICTED CHANGE IN OBJECTIVE FUNCTION, DELFPR:%12.3f", delfpr); Msg(fp,s1); sprintf(s1," ACTUAL CHANGE IN OBJECTIVE FUNCTION, DELTAF:%12.3f", deltaf); Msg(fp,s1); } if ((*retcod <= 6 && ABS(delfpr-deltaf) <= ABS(deltaf)/TEN) || (deltaf <= delfts && (!newtkn) && (!convio) && *delstr == ZERO)) { if (MIN(newlen, *maxstp)/(*delta) > onept1) { // PROMISING STEP - INCREASE TRUST REGION. // STORE CURRENT POINT. for(i=1; i<=n; i++) { temp1[i][1]=xplus[i]; temp2[i][1]=xplpre[i]; } matcop (n, n, 1, 1, n, 1, temp1, temp2); for(i=1; i<=n; i++) xplpre[i]=temp2[i][1]; for(i=1; i<=n; i++) { temp1[i][1]=fvec[i]; temp2[i][1]=fplpre[i]; } matcop (n, n, 1, 1, n, 1, temp1, temp2); for(i=1; i<=n; i++) fplpre[i]=temp2[i][1]; *delstr=*delta; fcnpre=fcnnew; // IF NONMONOTONIC STEPS ARE BEING USED EXPAND TRUST // REGION TO NEWLEN, OTHERWISE EXPAND BY DELFAC. if (isejac > narmij) *delta=MIN(newlen, *maxstp); else *delta=MIN(delfac*(*delta), *maxstp); *retcod=12; if (*acpcod == 2) *acpstr=12; else *acpstr=1; *acpcod=0; } else { *retcod=0; if (*acpcod == 2) *acpcod=12; else *acpcod=1; } goto e10; } else { // CHANGE TRUST REGION SIZE DEPENDING ON DELTAF AND DELFPR *retcod=6; if (*acpcod == 2) *acpcod=12; else *acpcod=1; if (deltaf >= delfpr/TEN) { *delta=*delta/delfac; if (output > 3) { Line0(fp); Msg(fp," CHANGE IN F, DELTAF, IS > 0.1 DELFPR - REDUCE DELTA."); } } else if (trupdm == 0 && jupdm > 0) { // POWELL"S UPDATING SCHEME - FIND JAC S FIRST if (qnupdm == 0) { // UNSCALED JACOBIAN IN JAC for (i=1; i<=n; i++) rhs[i] = s[i]*scalef[i]; avmul (n, n, n, n, jac, rhs, wv3); } else { // MULTIPLY BY R FIRST uvmul (n, n, n, n, a, s, rhs); // THEN Q (IN JAC^) mmulv (n, jac, rhs, wv3); } dmult=delfpr/TEN - deltaf; sp=ZERO; ss=ZERO; for (k=1; k<=n; k++) { wv3[k] += fvecc[k]; sp += ABS(fvec[k]*(fvec[k] - wv3[k])); ss += (fvec[k]-wv3[k])*(fvec[k] - wv3[k]); } if (sp + SQRT(sp*sp + dmult*ss) < epsmch) powlam=TEN; else powlam=ONE + dmult/(sp + SQRT(sp*sp + dmult*ss)); powlam=SQRT(powlam); powmu=MIN(MIN(delfac,powlam), *powtau); *powtau=powlam/powmu; if (output > 3) { Line0(fp); Msg(fp," FACTORS IN POWELL UPDATING SCHEME:"); Line0(fp); sprintf(s1," LAMBDA: %12.3f MU: %12.3f TAU: %12.3f", powlam, powmu, *powtau); Msg(fp,s1); Msg(fp," DELTA IS MINIMUM OF MU*DELTA AND MAXSTP."); } *delta=MIN(powmu*(*delta), *maxstp); } else { if (deltaf < threeq*delfpr) { *delta=MIN(delfac*(*delta), *maxstp); if (output > 3) { Line0(fp); Msg(fp," CHANGE IN F, DELTAF, IS LESS THAN 0.75 X PREDICTED."); sprintf(s1," DELTA INCREASED TO: %12.3f", delta); Msg(fp,s1); } } else if (output > 3) { Line0(fp); Msg(fp," DELTAF BETWEEN 0.1 AND 0.75 DELFPR - LEAVE DELTA UNCHANGED"); } } } } // if fcnnew... e10: vmfree(vmblock); //free memory } //trstup() void rcdprt (int nunit, int retcod, REAL delta, REAL rellen, REAL stptol) { /*------------------------------------------------------------------------ ! FEB. 14, 1991 ! ! DESCRIBE MEANING OF RETURN CODES, RETCOD, FROM TRUST REGION UPDATING. !-----------------------------------------------------------------------*/ char s[80]; extern FILE *fp; Line0(fp); sprintf(s," RETCOD, FROM TRUST REGION UPDATING: %5d", retcod); Msg(fp,s); Line0(fp); if (retcod == 1) { Msg(fp," PROMISING STEP FOUND; DELTA HAS BEEN INCREASED TO NEWLEN BUT"); Msg(fp," BECAUSE OF OVERFLOWS IN THE FUNCTION VECTOR(S) IN SUBSEQUENT"); Msg(fp," STEP(S) THE PROJECTED DELTA IS LESS THAN THAT AT THE ALREADY"); Msg(fp," SUCCESSFUL STEP - RETURN TO SUCCESSFUL STEP AND ACCEPT AS"); Msg(fp," NEW POINT."); } else if (retcod == 2) { Msg(fp," PROMISING STEP FOUND; DELTA HAS BEEN INCREASED TO NEWLEN BUT"); Msg(fp," BECAUSE OF OVERFLOWS IN THE OBJECTIVE FUNCTION IN SUBSEQUENT"); Msg(fp," STEP(S) THE PROJECTED DELTA IS LESS THAN THAT AT THE ALREADY"); Msg(fp," SUCCESSFUL STEP - RETURN TO SUCCESSFUL STEP AND ACCEPT AS"); Msg(fp," NEW POINT."); } else if (retcod == 3) { Msg(fp," PROMISING STEP FOUND; DELTA HAS BEEN INCREASED TO NEWLEN BUT"); Msg(fp," BECAUSE OF SUBSEQUENT FAILURES IN THE STEP ACCEPTANCE TEST[S)"); Msg(fp," THE PROJECTED DELTA IS LESS THAN THAT AT THE ALREADY"); Msg(fp," SUCCESSFUL STEP - RETURN TO SUCCESSFUL STEP AND ACCEPT."); } else if (retcod == 4) Msg(fp," STEP ACCEPTED BY STEP SIZE CRITERION ONLY - DELTA REDUCED"); else if (retcod == 5) { Msg(fp," STEP ACCEPTED - NEW FUNCTION VALUE > PREVIOUS =>"); Msg(fp," REDUCE TRUST REGION."); } else if (retcod == 6) Msg(fp," STEP ACCEPTED - DELTA CHANGED AS DETAILED ABOVE."); else if (retcod == 7) { Msg(fp," NO PROGRESS MADE: RELATIVE STEP SIZE IS TOO SMALL"); sprintf(s," REL. STEP SIZE, RELLEN = %12.3f, STPTOL = %12.3f", rellen, stptol); Msg(fp,s); } else if (retcod == 8) { Msg(fp," POINT MODifIED BY CONSTRAINTS NOT A DESCENT DIRECTION"); Msg(fp," DELTA REDUCED TO CONFAC*RATIOM*DELTA."); } else if (retcod == 9) Msg(fp," OVERFLOW DETECTED IN FUNCTION VECTOR - DELTA REDUCED."); else if (retcod == 10) Msg(fp," OVERFLOW IN OBJECTIVE FUNCTION - DELTA REDUCED."); else if (retcod == 11) { Msg(fp," STEP NOT ACCEPTED - REDUCE TRUST REGION SIZE BY MINIMIZATION"); Msg(fp," OF QUADRATIC MODEL IN STEP DIRECTION."); } else Msg(fp," PROMISING STEP - INCREASE DELTA TO NEWLEN AND TRY A NEW STEP."); Line0(fp); sprintf(s," DELTA ON RETURN FROM TRUST REGION UPDATING: %11.3f", delta); Msg(fp,s); } // rcdprt() void nnes (bool absnew, bool cauchy, bool deuflh, bool geoms, bool linesr, bool newton, bool overch, int acptcr, int *itsclf, int *itsclx, int jactyp, int jupdm, int *maxexp, int maxit, int maxns, int maxqns, int mgll, int minqns, int n, int narmij, int niejev, int njacch, int *njetot, int nunit, int *output, int qnupdm, int *stopcr, int supprs, int *trmcod, int trupdm, REAL *alpha, REAL *confac, REAL *delta, REAL *delfac, REAL epsmch, REAL *etafac, REAL *fcnnew, REAL *fdtolj, REAL *ftol, REAL *lam0, REAL *mstpf, REAL *nsttol, REAL *omega, REAL *ratiof, REAL *sigma, REAL *stptol, REAL **a, REAL * boundl, REAL *boundu, REAL *delf, REAL *fsave, REAL *ftrack, REAL *fvec, REAL *fvecc, REAL **h, REAL *hhpi, REAL **jac, REAL **plee, REAL *rdiag, REAL *s, REAL *sbar, REAL *scalef, REAL *scalex, REAL *sn, REAL *ssdhat, REAL *strack, REAL *vhat, REAL *xc, REAL *xplus, REAL *xsave, char *help) { //!---------------------------------------------------------------------------- //! FEB. 28, 1992 //! *** Solve a non-linear system - Main function *** //!---------------------------------------------------------------------------- //Labels: e310, e320, e330, e400 #define pt01 0.01 REAL beta, caulen, delstr, delta0, fcnmax, fcnmin, fcnold, fcnpre, fstore, maxstp, newlen, newmax, powtau, rellen, sbrnrm, sqrtz, stpmax, xnorm; REAL *wv1, *wv2, *wv4; int acpcod, acpstr, contyp, countr, i, isejac, itemp, itstr, j, maxlin, maxtrs, mnew, mold, nac1, nac2, nac12, nfail, nfestr, nfetot, nosubt, notrst, outtmp, retcod; bool abort, checkj, frstdg, instop, jacerr, newtkn, overfl, qnfail, qrsing, restrt, savest, sclfch, sclxch; int newstm; char ss[80]; void *vmblock = NULL; extern FILE *fp; extern itnum, nfunc; extern bool matsup; newstm=0; vmblock = vminit(); wv1 = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); wv2 = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); wv4 = (REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return; } qnfail=FALSE; restrt=TRUE; savest=FALSE; countr=0; isejac=0; mnew=0; nac1=0; nac2=0; nac12=0; nfunc=1; njetot=0; outtmp=*output; *trmcod=0; retcod = 0; nfetot=1; // ESTABLISH INITIAL FUNCTION VALUE AND CHECK FOR STARTING ESTIMATE WHICH // IS A SOLUTION. ALSO, CHECK FOR INCOMPATIBILITIES IN INPUT PARAMETERS. initch (&instop, &linesr, &newton, &overfl, &sclfch, &sclxch, acptcr, &contyp, jactyp, jupdm, *maxexp, n, nunit, *output, qnupdm, *stopcr, trupdm, &epsmch, &fcnold, *ftol, boundl, boundu, fvecc, scalef, scalex, xc); // IF A FATAL ERROR IS DETECTED IN INITCH RETURN TO MAIN PROGRAM. // NOTE: SOME INCOMPATIBILITIES ARE CORRECTED WITHIN INITCH AND EXECUTION // CONTINUES. WARNINGS ARE GENERATED WITHIN INITCH. if (instop) goto e400; // ESTABLISH MAXIMUM STEP LENGTH ALLOWED (USUALLY THIS IS MUCH // LARGER THAN ACTUAL STEP SIZES - IT IS ONLY TO PREVENT // EXCESSIVELY LARGE STEPS). THE FACTOR MSTPF CONTROLS THE // MAGNITUDE OF MAXSTP AND IS SET BY THE USER (DEFAULT=1000). maxst (&overfl, *maxexp, n, nunit, *output, epsmch, &maxstp, *mstpf, scalex, xc); // WRITE TITLE AND RECORD PARAMETERS. title (cauchy, deuflh, geoms, linesr, newton, overch, acptcr, contyp, *itsclf, *itsclx, jactyp, jupdm, maxit, maxns, maxqns, mgll, minqns, n, narmij, niejev, njacch, nunit, *output, qnupdm, *stopcr, trupdm, *alpha, *confac, *delfac, *delta, epsmch, *etafac, fcnold, *ftol, *lam0, maxstp, *mstpf, *nsttol, *omega, *ratiof, *sigma, *stptol, boundl, boundu, fvecc, scalef, scalex, xc); // INITIALIZE FTRACK AND STRACK VECTORS (FTRACK STORES (TRACKS) THE FUNCTION // VALUES FOR THE NONMONOTONIC COMPARISON AND STRACK, SIMILARLY, STORES THE // LENGTH OF THE NEWTON STEPS - WHICH ARE USED IN CONJUNCTION WITH // DEUFLHARD'S SECOND ACCEPTANCE CRITERION). for (j=0; j 2) { Line1(fp); Line0(fp); sprintf(ss," ITERATION NUMBER: %5d", itnum); Msg(fp,ss); } // UPDATE ITERATION COUNTER, ISEJAC, FOR QUASI-NEWTON METHODS // ONLY (I.E. JUPDM > 0). if (jupdm > 0 && (restrt) && (!newton)) isejac=1; else isejac++; if (*output > 4 && jupdm > 0 && (!newton)) { Line0(fp); if (restrt) if (itnum > niejev) Msg(fp," RESTRT IS TRUE, ISEJAC SET TO 1."); else { sprintf(ss," # OF ITERATIONS SINCE EXPLICIT JACOBIAN, ISEJAC, INCREASED TO %4d", isejac); Msg(fp,ss); } } // WHEN AN EXPLICIT JACOBIAN IS BEING USED IN QUASI-NEWTON METHODS THEN // MAXNS STEPS ARE ALLOWED IN THE LINE SEARCH. FOR STEPS BASED ON A // QUASI-NEWTON APPROXIMATION MAXQNS ARE ALLOWED. SIMILARLY MAXNS AND // MAXQNS STEPS ARE ALLOWED IN TRUST REGION METHODS RESPECTIVELY. // THIS IS AN ATTEMPT TO AVOID AN EXCESSIVE NUMBER OF FUNCTION // EVALUATIONS IN A DIRECTION WHICH WOULD NOT LEAD TO A SIGNIFICANT // REDUCTION. if (!newton) { if (linesr) if (jupdm > 0) if (isejac == 1) // JACOBIAN UPDATED EXPLICITLY maxlin=maxns; else // QUASI-NEWTON UPDATE maxlin=maxqns; else maxlin=maxns; else if (jupdm > 0) if (isejac == 1) // JACOBIAN UPDATED EXPLICITLY maxtrs=maxns; else // QUASI-NEWTON UPDATE maxtrs=maxqns; else maxtrs=maxns; if (jupdm > 0 && *output > 4) { Line0(fp); if (linesr) { sprintf(ss," MAXLIN SET TO: %5d", maxlin); Msg(fp,ss); } else { sprintf(ss," MAXTRS SET TO: %5d", maxtrs); Msg(fp,ss); } } } // ESTABLISH WHETHER JACOBIAN IS TO BE CHECKED NUMERICALLY - // NJACCH ESTABLISHES THE NUMBER OF ITERATIONS FOR WHICH JACOBIAN // CHECKING IS DESIRED. IF CHECKJ IS TRUE, A FORWARD DIFFERENCE NUMERICAL // APPROXIMATION OF THE JACOBIAN IS COMPARED TO THE ANALYTICAL VERSION. // STORE THE NUMBER OF FUNCTION EVALUATIONS SO THAT THESE "EXTRA" ARE NOT // INCLUDED IN OVERALL STATISTICS. if (jactyp == 0) if (itnum > njacch) checkj=FALSE; else { checkj=TRUE; nfestr=nfetot; } else checkj=FALSE; // EVALUATE JACOBIAN AT FIRST STEP OR IF NO QUASI-NEWTON // UPDATE IS BEING USED (RESTRT IS FALSE ONLY IN QUASI- // NEWTON METHODS WHEN THE QUASI-NEWTON UPDATE IS BEING USED). if ((restrt) || jupdm == 0) { // IF MORE THAN ONE DAMPED NEWTON STEP HAS BEEN REQUESTED AT THE START, // IDENTIFY THIS AS THE REASON FOR EXPLICIT JACOBIAN EVALUATION. if (jupdm > 0 && itnum <= niejev && itnum > 1 && *output > 4) { Line0(fp); Msg(fp," AS ITNUM <= NIEJEV JACOBIAN EVALUATED EXPLICITLY."); } // OTHERWISE JUPDM IS 0 OR RESTRT IS TRUE AND ITNUM IS > NIEJEV if (jupdm > 0 && itnum > 1 && *output > 4 && itnum > niejev) { Line0(fp); Msg(fp," RESTRT IS TRUE - JACOBIAN EVALUATED EXPLICITLY."); } // NOTE: MATRIX H IS USED HERE TO HOLD THE FINITE DIFFERENCE // ESTIMATION USED IN CHECKING THE ANALYTICAL JACOBIAN IF CHECKJ // IS TRUE. VECTORS WV1 AND, FOR CENTRAL DIFFERENCES, WV2 // TEMPORARILY HOLD THE FINITE DIFFERENCE FUNCTION EVALUATIONS. jacobi (checkj, &jacerr, &overfl, jactyp, n, nunit, *output, epsmch, fdtolj, boundl, boundu, fvecc, wv1, wv2, jac, h, scalex, xc); njetot++; // RETURN IF ANALYTICAL AND NUMERICAL JACOBIANS DON'T AGREE // (APPLICABLE ONLY IF CHECKJ IS TRUE AND A DISCREPANCY IS FOUND // WITHIN SUBROUTINE JACOBI). A WARNING IS GIVEN FROM WITHIN JACOBI. if (jacerr) return; // RESET TOTAL NUMBER OF FUNCTION EVALUATIONS TO NEGLECT // THOSE USED IN CHECKING ANALYTICAL JACOBIAN. if (checkj) nfetot=nfestr; if (jupdm > 0) { // FCNMIN IS THE MINIMUM OF THE OBJECTIVE FUNCTION FOUND SINCE // THE LAST EXPLICIT JACOBIAN EVALUATION. IT IS USED TO ESTABLISH // WHICH STEP THE PROGRAM RETURNS TO WHEN A QUASI-NEWTON STEP FAILS. fcnmin=fcnold; // POWTAU IS THE TAU FROM POWELL"S TRUST REGION UPDATING SCHEME // (USED WHEN JUPDM=1). IT IS RESET TO 1.0 AT EVERY EXPLICIT // JACOBIAN EVALUATION IN QUASI-NEWTON METHODS. powtau=ONE; // AT EVERY EXPLICIT JACOBIAN EVALUATION EXCEPT POSSIBLY THE FIRST, // IN QUASI-NEWTON METHODS, A NEW TRUST REGION IS CALCULATED // INTERNALLY USING EITHER THE NEWTON STEP OR THE CAUCHY STEP AS // SPECIFIED BY THE USER IN LOGICAL VARIABLE "CAUCHY" IN // SUBROUTINE DELCAU. THIS IS FORCED BY SETTING DELTA TO A // NEGATIVE NUMBER. if (itnum > 1) *delta=-ONE; // RESET COUNTER FOR NUMBER OF FAILURES OF RATIO TEST (IS // FCNNEW/FCNOLD < RATIOF?) SINCE LAST EXPLICIT JACOBIAN UPDATE. nfail=0; if (*output > 4 && itnum > niejev) { Line0(fp); Msg(fp," NUMBER OF FAILURES OF RATIO TEST, NFAIL, SET BACK TO 0."); } // ESTABLISH A NEW MAXIMUM STEP LENGTH ALLOWED if (itnum > 1) maxst (&overfl, *maxexp, n, nunit, *output, epsmch, &maxstp, *mstpf, scalex, xc); // SET "P" MATRIX TO IDENTITY FOR LEE AND LEE UPDATE (JUPDM=2). if (jupdm == 2 && itnum >= niejev) for (j=1; j<=n; j++) { for (i=1; i<=n; i++) plee[i][j]=ZERO; plee[j][j]=ONE; } } } if (((restrt) || qnupdm == 0) && *output > 4 && (!matsup)) { // WRITE JACOBIAN MATRIX Line0(fp); Msg(fp," JACOBIAN MATRIX:"); matprt (n, n, jac); } // ESTABLISH SCALING MATRICES IF DESIRED (ITSCLF=0 => NO ADAPTIVE // SCALING, WHILE ITSCLF > 0 MEANS ADAPTIVE SCALING STARTS AFTER THE // [ITSCLF]TH ITERATION). SIMILARLY FOR COMPONENT SCALING... // NOTE: SCALING FACTORS ARE UPDATED ONLY WHEN THE JACOBIAN IS UPDATED // EXPLICITLY IN QUASI-NEWTON METHODS. // FUNCTION SCALING. if ((restrt) && *itsclf > 0 && *itsclf <= itnum) { ascalf (n, epsmch, fvecc, jac, scalef); if (*output > 4) { Line0(fp); Msg(fp," FUNCTION SCALING VECTOR:"); Line0(fp); for(i=1; i<=n; i++) { sprintf(ss," SCALEF(%3d) = %12.3f", i, scalef[i]); Msg(fp,ss); } } // RECALCULATE OBJECTIVE FUNCTION WITH NEW SCALING FACTORS. // THIS AVOIDS PREMATURE FAILURES IF THE CHANGE IS SCALING FACTORS // WOULD MAKE THE PREVIOUS OBJECTIVE FUNCTION VALUE SMALLER. fcnevl (&overfl, *maxexp, n, nunit, *output, epsmch, &fcnold, fvecc, scalef); } // COMPONENT SCALING if ((restrt) && *itsclx > 0 && *itsclx <= itnum) { ascalx (n, epsmch, jac, scalex); if (*output > 4) { Line0(fp); Msg(fp," COMPONENT SCALING VECTOR:"); Line0(fp); for (i=1; i<=n; i++) { sprintf(ss," SCALEX(%3d) = %12.3f", i, scalex[i]); Msg(fp,ss); } } } // FIND GRADIENT OF 1/2 FVECC^FVECC (NOT USED IF DECOMPOSED // MATRIX IS UPDATED IN WHICH CASE THE GRADIENT IS FOUND // WITHIN THAT SUBROUTINE - CALL IS MADE FOR OUTPUT ONLY). gradf (overch, &overfl, restrt, sclfch, sclxch, jupdm, *maxexp, n, nunit, *output, qnupdm, delf, fvecc, jac, scalef, scalex); // FIND NEWTON STEP USING QR DECOMPOSITION. // IF JUPDM = 0 OR QNUPDM = 0 THEN THE UNFACTORED FORM // FOR CALCULATING THE NEWTON STEP IS USED. if (jupdm == 0 || qnupdm == 0) // NEWTON STEP - UNFACTORED FORM BEING UPDATED. nstpun(&abort, &linesr, overch, &overfl, &qrsing, &sclfch, sclxch, itnum, *maxexp, n, nunit, *output, epsmch, a, delf, fvecc, h, hhpi, jac, rdiag, scalef, scalex, sn); else // NEWTON STEP - FACTORED FORM BEING UPDATED. nstpfa (&abort, linesr, overch, overfl, qrsing, restrt, sclfch, sclxch, itnum, *maxexp, n, &newstm, nunit, *output, epsmch, a, delf, fvecc, h, hhpi, jac, rdiag, scalef, scalex, sn); // RUN IS ABORTED IF THE JACOBIAN BECOMES ESSENTIALLY ALL ZEROS. if (abort) return; // CHECK FOR CONVERGENCE ON LENGTH OF NEWTON STEP IF STOPCR = 1, 12 OR 3. if (*stopcr != 2) { stpmax=ZERO; for(i=1; i<=n; i++) { stpmax=MAX(stpmax, ABS(sn[i])/MAX(xc[i], scalex[i])); wv1[i]=scalex[i]*xc[i]; } twonrm (&overfl, *maxexp, n, epsmch, &xnorm, wv1); if (stpmax <= *nsttol*(ONE+xnorm)) { *trmcod=1; // IF STOPCR=3 THEN OBJECTIVE FUNCTION VALUE MUST BE // DETERMINED AS WELL - OTHERWISE A SOLUTION HAS BEEN FOUND. if (*stopcr != 3) goto e330; // NOTE: STATEMENT 202 PRECEDES CONVERGENCE CHECKING SUBROUTINE } } // FIND LENGTH OF (SCALED) NEWTON STEP, NEWLEN for(i=1; i<=n; i++) wv1[i]=scalex[i]*sn[i]; twonrm (&overfl, *maxexp, n, epsmch, &newlen, wv1); // FOR ITERATIONS AFTER THE ARMIJO STEPS HAVE BEEN COMPLETED // AT THE BEGINNING (IN OTHER WORDS THE MONOTONIC STEPS) // STORE THE FUNCTION AND, POSSIBLY, THE NEWTON STEP LENGTHS // IN THE FTRACK AND STRACK VECTORS, RESPECTIVELY. if (isejac >= narmij) { if (isejac == 1) { strack[0]=newlen; // NEWMAX IS USED TO KEEP A BOUND ON THE ENTRIES IN // THE STRACK VECTOR. newmax=newlen; } else strack[countr]=MIN(newmax,newlen); // THE OBJECTIVE FUNCTION VALUE IS STORED EVEN IF IT IS // GREATER THAN ANY PRECEEDING FUNCTION VALUE. ftrack[countr]=fcnold; // WRITE FTRACK AND STRACK VECTORS IF DESIRED. SINCE ONLY THE LAST // MGLL VALUES ARE NEEDED THE COUNTER CIRCULATES THROUGH THE VECTOR // CAUSING ONLY THE MGLL MOST RECENT VALUES TO BE KEPT. // NOTE: THESE VECTORS ARE NOT APPLICABLE IF NEWTON"S METHOD IS BEING // USED. if ((!newton) && *output > 4) { Line0(fp); Line0(fp); // IF ONLY THE FUNCTION VALUE ACCEPTANCE TEST IS BEING USED, // THUS ACPTCR=1, THEN ONLY THE FTRACK VECTOR IS APPLICABLE. if (acptcr == 1) { sprintf(ss," CURRENT FTRACK VECTOR; LATEST CHANGE: ELEMENT %4d", countr); Msg(fp,ss); Line0(fp); for (j=0; j 0) { Line0(fp); Line1(fp); } return; } // NOTE: 201 PRECEDES PRINTING OF ITERATION RESULTS if (newton) goto e320; } else { // TRUST REGION METHOD // ESTABLISH INITIAL TRUST REGION SIZE, DELTA, AND/OR // FIND LENGTH OF SCALED DESCENT STEP, CAULEN. delcau (cauchy, overch, &overfl, isejac, *maxexp, n, nunit, *output, &beta, &caulen, delta, epsmch, &maxstp, &newlen, &sqrtz, a, delf, scalex); frstdg=TRUE; if (*output > 3) { Line0(fp); Line0(fp); Msg(fp," SUMMARY OF TRUST REGION METHOD USING (DOUBLE) DOGLEG STEP."); } // MAIN INTERNAL LOOP FOR TRUST REGION METHOD. // THE TRUST REGION SIZE IS STORED FOR COMPARISON // LATER TO SET THE PARAMETER POWTAU USED IN POWELL'S // TRUST REGION UPDATING SCHEME (QUASI-NEWTON, TRUPDM=1) delta0=*delta; for (notrst=1; notrst<=maxtrs; notrst++) { dogleg (&frstdg, &newtkn, overch, &overfl, maxexp, n, ¬rst, nunit, *output, &beta, caulen, delta, *etafac, newlen, sqrtz, delf, s, scalex, sn, ssdhat, vhat); // NOTE: WV1 AND WV4 HOLD THE COMPONENT AND RESIDUAL VECTOR // RESPECTIVELY FOR A TRIAL POINT WHICH HAS BEEN FOUND // TO BE ACCEPTABLE WHILE THE TRUST REGION IS EXPANDED // AND A NEW TRIAL POINT TESTED. // WV2 AND WV3 ARE WORK VECTORS. // H IS CALLED ASTORE INSIDE TRSTUP. trstup (geoms, &newtkn, overch, &overfl, qrsing, &sclfch, &sclxch, &acpcod, &acpstr, acptcr, contyp, isejac, jupdm, maxexp, mgll, mnew, n, narmij, &nfunc, notrst, nunit, *output, qnupdm, &retcod, &trupdm, *alpha, *confac, *delfac, &delstr, delta, epsmch, &fcnmax, fcnnew, fcnold, &fcnpre, &maxstp, newlen, &newmax, &powtau, &rellen, *stptol, a, h, boundl, boundu, delf, wv1, ftrack, fvec, fvecc, hhpi, jac, rdiag, wv2, s, sbar, scalef, scalex, strack, xc, wv4, xplus); if ((*output > 4 || retcod == 7) && *output > 2) rcdprt (nunit, retcod, *delta, rellen, *stptol); // IF NO PROGRESS WAS BEING MADE (RETCOD=7) IN A QUASI-NEWTON // STEP RETRY WITH AN EXPLICIT JACOBIAN EVALUATION. if (retcod == 7 && (!restrt)) qnfail=TRUE; // RETURN CODE LESS THAN 8 EXITS FROM TRUST REGION LOOP if (retcod < 8) goto e310; } // IF NO SUCCESSFUL STEP FOUND IN A QUASI-NEWTON STEP, // RETRY WITN AN EXPLICIT JACOBIAN EVALUATION. if (!restrt) qnfail=TRUE; } // else e310: if (!linesr) { if (*delta < delta0) powtau=ONE; *delta=MAX(*delta, 1e-10); } // IF RETCOD=7 AND STOPCR=2, RESET STOPPING CRITERION TO // AVOID HANGING IN TRUST REGION METHOD. (RETCOD=7 MEANS // THE RELATIVE STEP LENGTH WAS LESS THAN STPTOL). if ((!linesr) && retcod == 7 && *stopcr == 2 && (!qnfail)) *stopcr=12; // RETAIN NUMBER OF STEPS ACCEPTED BY EACH CRITERION FOR PERFORMANCE EVALUATION. if (!newton) if (acpcod == 1) nac1++; else if (acpcod == 2) nac2++; else if (acpcod == 12) nac12++; // *** PRINT RESULTS OF ITERATION *** e320:if (*output > 2) nersl (newton, restrt, sclfch, sclxch, acpcod, jupdm, n, nunit, *output, *fcnnew, fvec, xplus); // CHECK FOR CONVERGENCE. STATEMENT 202 IS USED IF THE // STEP SIZE OF THE NEWTON STEP IS FOUND TO BE WITHIN // THE SPECIFIED TOLERANCE AND STOPCR IS 1 OR 12 e330:if (!qnfail) // IF QNFAIL IS TRUE THE QUASI-NEWTON SEARCH FAILED TO // FIND A SATISFACTORY STEP - SINCE THE JACOBIAN IS TO // BE RE-EVALUATED AVOID PREMATURE STOPPAGES IN NESTOP. nestop (&absnew, linesr, &newton, &sclfch, &sclxch, acptcr, &itnum, n, &nac1, &nac2, &nac12, &nfunc, njetot, nunit, *output, *stopcr, trmcod, fcnnew, *ftol, *nsttol, &stpmax, *stptol, fvec, scalef, scalex, xc, xplus); // IF THE TERMINATION CODE, TRMCOD, IS GREATER THAN 0 THEN // CONVERGENCE HAS BEEN REACHED. if (*trmcod > 0) return; //** normal exit ** // QUASI-NEWTON UPDATING - JUPDM > 0 if (jupdm > 0) { // QNFAIL MEANS A FAILURE IN THE QUASI-NEWTON SEARCH. // RE-EVALUATE JACOBIAN AND TRY A DAMPED NEWTON STEP. // MAXLIN IS CHANGED FROM MAXQNS TO MAXNS OR MAXTRS IS // CHANGED, SIMILARLY, AT THE START OF THE LOOP. if (qnfail) if (*output > 4) { Line0(fp); Msg(fp," FAILURE IN QUASI-NEWTON SEARCH: QNFAIL IS TRUE."); } else if (!newton) { // QNFAIL IS FALSE - THE STEP HAS BEEN ACCEPTED. if (*output > 4) { Line0(fp); sprintf(ss," FCNNEW= %11.3f FCNOLD= %11.3f RATIO= %11.3f", *fcnnew, fcnold, *fcnnew/fcnold); Msg(fp,ss); } if (*fcnnew/fcnold > *ratiof) { // STEP ACCEPTED BUT NOT A SIGNIFICANT IMPROVEMENT. nfail++; if (*output > 4) { Line0(fp); sprintf(ss," RATIO > RATIOF SO NFAIL INCREASED TO: %5d", nfail); Msg(fp,ss); Line0(fp); } } else { // STEP ACCEPTED WITH A SIGNIFICANT IMPROVEMENT. if (*fcnnew/fcnold > pt01) nosubt=0; else nosubt=1; // ITEMPT IS USED LOCALLY FOR OUTPUT CONTROL. itemp=nfail; nfail=IMAX(nfail-nosubt,0); if (*output > 4) { Line0(fp); if (itemp == nfail) { sprintf(ss," NFAIL STAYS AT: %5d", nfail); Msg(fp,ss); } else { sprintf(ss," NFAIL CHANGED TO: %5d", nfail); Msg(fp,ss); } } } // SAVE THE RESULTS FOR RESTART IF A FAILURE IN THE // QUASI-NEWTON METHOD OCCURS - ESSENTIALLY THIS // FINDS THE BEST POINT SO FAR. if (isejac == 1 || nfail == 1 || (nfail <= minqns && *fcnnew/fcnmin < ONE)) { savest=TRUE; fcnmin=*fcnnew; } if (savest) { savest=FALSE; fstore=*fcnnew; if (*output > 4) { Line0(fp); Msg(fp," STEP IS SAVED"); Line0(fp); Msg(fp," SAVED COMPONENT AND FUNCTION VALUES."); Line0(fp); } for (i=1; i<=n; i++) { xsave[i]=xplus[i]; fsave[i]=fvec[i]; if (*output > 4) { sprintf(ss," XSAVE(%3d) = %12.3f FSAVE(%3d) = %12.3f", i, xsave[i], i, fsave[i]); Msg(fp,ss); } } if (*output > 4) Line0(fp); itstr=itnum; } } // NOTE: IF QNFAIL IS TRUE THEN NFAIL CANNOT HAVE INCREASED IMPLYING // THAT NFAIL CANNOT NOW BE GREATER THAN MINQNS. if ((qnfail) || nfail > minqns) { // RESTART FROM BEST POINT FOUND SO FAR. restrt=TRUE; if (*output > 4) { Line0(fp); Line0(fp); Msg(fp," RESTRT IS TRUE."); } for (j=0; j 4) { Line0(fp); sprintf(ss," RETURN TO ITERATION: %5d WHERE:", itstr); Msg(fp,ss); Line0(fp); } for (i=1; i<=n; i++) { xc[i]=xsave[i]; fvecc[i]=fsave[i]; if (*output > 4) { sprintf(ss," XC(%3d) = %12.3f FVECC(%3d) = %12.3f", i, xc[i], i, fvecc[i]); Msg(fp,ss); } } if (*output > 4) Line0(fp); } else if (itnum >= niejev) restrt=FALSE; } // UPDATE JACOBIAN IF DESIRED: // QNUPDM = 0 ==> ACTUAL JACOBIAN BEING UPDATED // QNUPDM = 1 ==> FACTORED JACOBIAN BEING UPDATED. if (!restrt) { if (qnupdm == 0) { if (jupdm == 1) // USE BROYDEN UPDATE broyun (&overfl, *maxexp, n, nunit, *output, epsmch, fvec, fvecc, jac, scalex, xc, xplus); else if (jupdm == 2) // USE LEE AND LEE UPDATE llun (overch, &overfl, isejac, *maxexp, n, nunit, *output, epsmch, *omega, fvec, fvecc, jac, plee, s, scalex, xc, xplus); // THE FACTORED FORM OF THE JACOBIAN IS UPDATED. if (jupdm == 1) broyfa (overch, &overfl, sclfch, sclxch, *maxexp, n, nunit, *output, epsmch, a, delf, fvec, fvecc, jac, rdiag, s, scalef, scalex, wv1, wv2, xc, xplus); else if (jupdm == 2) llfa (overch, &overfl, &sclfch, &sclxch, isejac, *maxexp, n, nunit, *output, epsmch, omega, a, delf, fvec, fvecc, jac, plee, rdiag, s, scalef, scalex, xc, xplus); } } // UPDATE CURRENT VALUES - RESET TRMCOD TO ZERO. // UPDATE M "VECTOR" (ACTUALLY ONLY THE LATEST VALUE IS NEEDED). mold=mnew; if (isejac < narmij) mnew=0; else mnew=IMIN(mold+1,mgll-1); if (jupdm == 0 || (jupdm > 0 && (!restrt))) update (mnew, &mold, n, trmcod, fcnnew, &fcnold, fvec, fvecc, xc, xplus); } // *** end of itnum main loop *** if (*output > 0) { Line1(fp); Line0(fp); sprintf(ss," NO SOLUTION FOUND AFTER %6d ITERATION(S).", itnum-1); Msg(fp,ss); Line0(fp); Msg(fp," FINAL ESTIMATES FINAL FUNCTION VALUES *"); Line0(fp); for(i=1; i<=n; i++) { sprintf(ss," X(%3d) = %15.6f F(%3d) = %15.6f", i, xplus[i], i, fvec[i]); Msg(fp,ss); } Line0(fp); sprintf(ss," FINAL OBJECTIVE FUNCTION VALUE = %15.6f", fcnnew); Msg(fp,ss); Line0(fp); Line1(fp); } e400: vmfree(vmblock); } // nnes() // end of file unnes2.cpp