1 #ifndef lint 2 static char vcid[] = "$Id: newtr1.c,v 1.8 1995/02/22 00:52:41 curfman Exp $"; 3 #endif 4 5 #include <math.h> 6 #include "nonlin/nlall.h" 7 #include "nonlin/snes/nlepriv.h" 8 9 /*D 10 NLE_NTR1 - Implements Newton's Method with a trust region 11 approach for solving systems of nonlinear equations. 12 13 Input parameters: 14 . nlP - nonlinear context obtained from NLCreate() 15 16 Returns: 17 Number of global iterations until termination. The precise type of 18 termination can be examined by calling NLGetTerminationType() after 19 NLSolve(). 20 21 Calling sequence: 22 $ nlP = NLCreate(NLE_NTR1,0) 23 $ NLCreateDVectors() 24 $ NLSet***() 25 $ NLSetUp() 26 $ NLSolve() 27 $ NLDestroy() 28 29 Notes: 30 See NLCreate() and NLSetUp() for information on the definition and 31 initialization of the nonlinear solver context. 32 33 The basic algorithm is taken from "The Minpack Project", by More', 34 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 35 of Mathematical Software", Wayne Cowell, editor. See the examples 36 in nonlin/examples. 37 D*/ 38 /* 39 This is intended as a model implementation, since it does not 40 necessarily have many of the bells and whistles of other 41 implementations. 42 43 The code is DATA-STRUCTURE NEUTRAL and can be called RECURSIVELY. 44 The following context variable is used: 45 NLCtx *nlP - The nonlinear solver context, which is created by 46 calling NLCreate(NLE_NTR1). 47 48 The step_compute routine must return two values: 49 1) ynorm - the norm of the step 50 2) gpnorm - the predicted value for the residual norm at the new 51 point, assuming a local linearization. The value is 0 if the 52 step lies within the trust region and is > 0 otherwise. 53 */ 54 int NLENewtonTR1Solve( nlP ) 55 NLCtx *nlP; 56 { 57 NLENewtonTR1Ctx *neP = (NLENewtonTR1Ctx *) nlP->MethodPrivate; 58 void *X, *F, *Y, *G, *TMP; 59 int maxits, i, iters, history_len, nlconv; 60 double rho, fnorm, gnorm, gpnorm, xnorm, delta; 61 double *history, ynorm; 62 FILE *fp = nlP->fp; 63 NLMonCore *mc = &nlP->mon.core; 64 VECntx *vc = NLGetVectorCtx(nlP); 65 66 CHKCOOKIEN(nlP,NL_COOKIE); 67 nlconv = 0; /* convergence monitor */ 68 history = nlP->conv_hist; /* convergence history */ 69 history_len = nlP->conv_hist_len; /* convergence history length */ 70 maxits = nlP->max_its; /* maximum number of iterations */ 71 X = nlP->vec_sol; /* solution vector */ 72 F = nlP->vec_rg; /* residual vector */ 73 Y = nlP->work[0]; /* work vectors */ 74 G = nlP->work[1]; 75 76 INITIAL_GUESS( X ); /* X <- X_0 */ 77 VNORM( vc, X, &xnorm ); /* xnorm = || X || */ 78 79 RESIDUAL( X, F ); /* Evaluate (+/-) F(X) */ 80 VNORM( vc, F, &fnorm ); /* fnorm <- || F || */ 81 nlP->norm = fnorm; 82 if (history && history_len > 0) history[0] = fnorm; 83 delta = neP->delta0*fnorm; 84 neP->delta = delta; 85 mc->nvectors += 4; mc->nscalars += 3; 86 MONITOR( X, F, &fnorm ); /* Monitor progress */ 87 88 for ( i=0; i<maxits; i++ ) { 89 nlP->iter = i+1; 90 91 STEP_SETUP( X ); /* Step set-up phase */ 92 while(1) { 93 iters = STEP_COMPUTE( X, F, Y, &fnorm, &delta, 94 &(nlP->trunctol), &gpnorm, &ynorm, (void *)0 ); 95 CHKERRV(1,-(NL)); /* Step compute phase */ 96 VAXPY( vc, 1.0, X, Y ); /* Y <- X + Y */ 97 RESIDUAL( Y, G ); /* Evaluate (+/-) G(Y) */ 98 VNORM( vc, G, &gnorm ); /* gnorm <- || g || */ 99 if (fnorm == gpnorm) rho = 0.0; 100 else rho = (fnorm*fnorm - gnorm*gnorm)/ 101 (fnorm*fnorm - gpnorm*gpnorm); 102 103 /* Update size of trust region */ 104 if (rho < neP->mu) delta *= neP->delta1; 105 else if (rho < neP->eta) delta *= neP->delta2; 106 else delta *= neP->delta3; 107 108 if (fp) fprintf(fp,"%d: f=%g, g=%g, ynorm=%g\n", 109 i, fnorm, gnorm, ynorm ); 110 if (fp) fprintf(fp," gpred=%g, rho=%g, delta=%g, iters=%d\n", 111 gpnorm, rho, delta, iters); 112 113 neP->delta = delta; 114 mc->nvectors += 4; mc->nscalars += 8; 115 if (rho > neP->sigma) break; 116 neP->itflag = 0; 117 if (CONVERGED( &xnorm, &ynorm, &fnorm )) { 118 /* We're not progressing, so return with the current iterate */ 119 if (X != nlP->vec_sol) VCOPY( vc, X, nlP->vec_sol ); 120 return i; 121 } 122 nlP->mon.nunsuc++; 123 } 124 STEP_DESTROY(); /* Step destroy phase */ 125 fnorm = gnorm; 126 nlP->norm = fnorm; 127 if (history && history_len > i+1) history[i+1] = fnorm; 128 TMP = F; F = G; G = TMP; 129 TMP = X; X = Y; Y = TMP; 130 VNORM( vc, X, &xnorm ); /* xnorm = || X || */ 131 mc->nvectors += 2; 132 mc->nscalars++; 133 MONITOR( X, F, &fnorm ); /* Monitor progress */ 134 135 /* Test for convergence */ 136 neP->itflag = 1; 137 if (CONVERGED( &xnorm, &ynorm, &fnorm )) { 138 /* Verify solution is in corect location */ 139 if (X != nlP->vec_sol) VCOPY( vc, X, nlP->vec_sol ); 140 break; 141 } 142 } 143 if (i == maxits) i--; 144 return i+1; 145 } 146 /* -------------------------------------------------------------*/ 147 void NLENewtonTR1Create( nlP ) 148 NLCtx *nlP; 149 { 150 NLENewtonTR1Ctx *neP; 151 152 CHKCOOKIE(nlP,NL_COOKIE); 153 nlP->method = NLE_NTR1; 154 nlP->method_type = NLE; 155 nlP->setup = NLENewtonTR1SetUp; 156 nlP->solver = NLENewtonTR1Solve; 157 nlP->destroy = NLENewtonTR1Destroy; 158 nlP->set_param = NLENewtonTR1SetParameter; 159 nlP->get_param = NLENewtonTR1GetParameter; 160 nlP->usr_monitor = NLENewtonDefaultMonitor; 161 nlP->converged = NLENewtonTR1DefaultConverged; 162 nlP->term_type = NLENewtonTR1DefaultConvergedType; 163 164 neP = NEW(NLENewtonTR1Ctx); CHKPTR(neP); 165 nlP->MethodPrivate = (void *) neP; 166 neP->mu = 0.25; 167 neP->eta = 0.75; 168 neP->delta = 0.0; 169 neP->delta0 = 0.2; 170 neP->delta1 = 0.3; 171 neP->delta2 = 0.75; 172 neP->delta3 = 2.0; 173 neP->sigma = 0.0001; 174 neP->itflag = 0; 175 } 176 /*------------------------------------------------------------*/ 177 /*ARGSUSED*/ 178 void NLENewtonTR1SetUp( nlP ) 179 NLCtx *nlP; 180 { 181 CHKCOOKIE(nlP,NL_COOKIE); 182 nlP->nwork = 2; 183 nlP->work = VGETVECS( nlP->vc, nlP->nwork ); CHKPTR(nlP->work); 184 NLiBasicSetUp( nlP, "NLENewtonTR1SetUp" ); CHKERR(1); 185 } 186 /*------------------------------------------------------------*/ 187 /*ARGSUSED*/ 188 void NLENewtonTR1Destroy( nlP ) 189 NLCtx *nlP; 190 { 191 CHKCOOKIE(nlP,NL_COOKIE); 192 VFREEVECS( nlP->vc, nlP->work, nlP->nwork ); 193 NLiBasicDestroy( nlP ); CHKERR(1); 194 } 195 /*------------------------------------------------------------*/ 196 /*ARGSUSED*/ 197 /*@ 198 NLENewtonTR1DefaultConverged - Default test for monitoring the 199 convergence of the method NLENewtonTR1Solve. 200 201 Input Parameters: 202 . nlP - nonlinear context obtained from NLCreate() 203 . xnorm - 2-norm of current iterate 204 . pnorm - 2-norm of current step 205 . fnorm - 2-norm of residual 206 207 Returns: 208 $ 1 if ( delta < xnorm*deltatol ), 209 $ 2 if ( fnorm < atol ), 210 $ 3 if ( pnorm < xtol*xnorm ), 211 $ -2 if ( nres > max_res ), 212 $ -1 if ( delta < xnorm*epsmch ), 213 $ 0 otherwise, 214 215 where 216 $ atol - absolute residual norm tolerance, 217 $ set with NLSetAbsConvergenceTol() 218 $ delta - trust region paramenter 219 $ deltatol - trust region size tolerance, 220 $ set with NLSetTrustRegionTol() 221 $ epsmch - machine epsilon 222 $ max_res - maximum number of residual evaluations, 223 $ set with NLSetMaxResidualEvaluations() 224 $ nres - number of residual evaluations 225 $ xtol - relative residual norm tolerance, 226 $ set with NLSetRelConvergenceTol() 227 228 Note: 229 Call NLGetConvergenceType() after calling NLSolve() to obtain 230 information about the type of termination which occurred for the 231 nonlinear solver. 232 @*/ 233 int NLENewtonTR1DefaultConverged( nlP, xnorm, pnorm, fnorm ) 234 NLCtx *nlP; 235 double *xnorm; 236 double *pnorm; 237 double *fnorm; 238 { 239 NLENewtonTR1Ctx *neP = (NLENewtonTR1Ctx *)nlP->MethodPrivate; 240 double epsmch = 1.0e-14; /* This must be fixed */ 241 242 CHKCOOKIEN(nlP,NL_COOKIE); 243 if (nlP->method_type != NLE) { 244 SETERRC(1,"Compatible with NLE component only"); 245 return 0; 246 } 247 nlP->conv_info = 0; 248 if (neP->delta < *xnorm * nlP->deltatol) nlP->conv_info = 1; 249 if (neP->itflag) { 250 if (nlP->conv_info) return nlP->conv_info; 251 nlP->conv_info = 252 NLENewtonDefaultConverged( nlP, xnorm, pnorm, fnorm ); 253 } 254 if (neP->delta < *xnorm * epsmch) nlP->conv_info = -1; 255 return nlP->conv_info; 256 } 257 /*------------------------------------------------------------*/ 258 /* 259 NLENewtonTR1DefaultConvergedType - Returns information regarding 260 the type of termination which occurred within the 261 NLENewtonTR1DefaultConverged() test. 262 263 Input Parameter: 264 . nlP - nonlinear context obtained from NLCreate() 265 266 Returns: 267 Character string - message detailing the type of termination which 268 occurred. 269 */ 270 char *NLENewtonTR1DefaultConvergedType( nlP ) 271 NLCtx *nlP; 272 { 273 char *mesg; 274 275 CHKCOOKIEN(nlP,NL_COOKIE); 276 if ((int)nlP->converged != (int)NLENewtonTR1DefaultConverged) { 277 mesg = "Compatible only with NLENewtonTR1DefaultConverged.\n"; 278 SETERRC(1,"Compatible only with NLENewtonTR1DefaultConverged."); 279 } else { 280 switch (nlP->conv_info) { 281 case 1: 282 mesg = "Trust region parameter satisfies the trust region tolerance.\n"; 283 break; 284 case -1: 285 mesg = "Machine epsilon tolerance exceeds the trust region parameter.\n"; 286 break; 287 default: 288 mesg = NLENewtonDefaultConvergedType( nlP ); 289 } } 290 return mesg; 291 } 292 /*------------------------------------------------------------*/ 293 /* 294 NLENewtonTR1SetParameter - Sets a chosen parameter used by the 295 NLE_NTR1 method to the desired value. 296 297 Note: 298 Possible parameters for the NLE_NTR1 method are 299 $ param = "mu" - used to compute trust region parameter 300 $ param = "eta" - used to compute trust region parameter 301 $ param = "sigma" - used to determine termination 302 $ param = "delta0" - used to initialize trust region parameter 303 $ param = "delta1" - used to compute trust region parameter 304 $ param = "delta2" - used to compute trust region parameter 305 $ param = "delta3" - used to compute trust region parameter 306 */ 307 void NLENewtonTR1SetParameter( nlP, param, value ) 308 NLCtx *nlP; 309 char *param; 310 double *value; 311 { 312 NLENewtonTR1Ctx *ctx = (NLENewtonTR1Ctx *)nlP->MethodPrivate; 313 314 CHKCOOKIE(nlP,NL_COOKIE); 315 if (nlP->method != NLE_NTR1) { 316 SETERRC(1,"Compatible only with NLE_NTR1 method"); 317 return; 318 } 319 if (!strcmp(param,"mu")) ctx->mu = *value; 320 else if (!strcmp(param,"eta")) ctx->eta = *value; 321 else if (!strcmp(param,"sigma")) ctx->sigma = *value; 322 else if (!strcmp(param,"delta0")) ctx->delta0 = *value; 323 else if (!strcmp(param,"delta1")) ctx->delta1 = *value; 324 else if (!strcmp(param,"delta2")) ctx->delta2 = *value; 325 else if (!strcmp(param,"delta3")) ctx->delta3 = *value; 326 else SETERRC(1,"Invalid parameter name for NLE_NTR1"); 327 } 328 /*------------------------------------------------------------*/ 329 /* 330 NLENewtonTR1GetParameter - Returns the value of a chosen parameter 331 used by the NLE_NTR1 method. 332 333 Note: 334 Possible parameters for the NLE_NTR1 method are 335 $ param = "mu" - used to compute trust region parameter 336 $ param = "eta" - used to compute trust region parameter 337 $ param = "sigma" - used to determine termination 338 $ param = "delta0" - used to initialize trust region parameter 339 $ param = "delta1" - used to compute trust region parameter 340 $ param = "delta2" - used to compute trust region parameter 341 $ param = "delta3" - used to compute trust region parameter 342 */ 343 double NLENewtonTR1GetParameter( nlP, param ) 344 NLCtx *nlP; 345 char *param; 346 { 347 NLENewtonTR1Ctx *ctx = (NLENewtonTR1Ctx *)nlP->MethodPrivate; 348 double value = 0.0; 349 350 CHKCOOKIEN(nlP,NL_COOKIE); 351 if (nlP->method != NLE_NTR1) { 352 SETERRC(1,"Compatible only with NLE_NTR1 method"); 353 return value; 354 } 355 if (!strcmp(param,"mu")) value = ctx->mu; 356 else if (!strcmp(param,"eta")) value = ctx->eta; 357 else if (!strcmp(param,"sigma")) value = ctx->sigma; 358 else if (!strcmp(param,"delta0")) value = ctx->delta0; 359 else if (!strcmp(param,"delta1")) value = ctx->delta1; 360 else if (!strcmp(param,"delta2")) value = ctx->delta2; 361 else if (!strcmp(param,"delta3")) value = ctx->delta3; 362 else SETERRC(1,"Invalid parameter name for NLE_NTR1"); 363 return value; 364 } 365