1 #ifndef lint 2 static char vcid[] = "$Id: newls1.c,v 1.9 1995/02/22 00:52:02 curfman Exp $"; 3 #endif 4 5 #include <math.h> 6 #include "nonlin/nlall.h" 7 #include "nonlin/snes/nlepriv.h" 8 #include "system/system.h" 9 #include "system/time/usec.h" 10 11 /*D 12 NLE_NLS1 - Implements a line search variant of Newton's Method 13 for solving systems of nonlinear equations. 14 15 Input parameters: 16 . nlP - nonlinear context obtained from NLCreate() 17 18 Returns: 19 Number of global iterations until termination. The precise type of 20 termination can be examined by calling NLGetTerminationType() after 21 NLSolve(). 22 23 Calling sequence: 24 $ nlP = NLCreate(NE_NLS1,0); 25 $ NLCreateDVectors() 26 $ NLSetXXX() 27 $ NLSetUp() 28 $ NLSolve() 29 $ NLDestroy() 30 31 Notes: 32 See NLCreate() and NLSetUp() for information on the definition and 33 initialization of the nonlinear solver context. 34 35 This implements essentially a truncated Newton method with a 36 line search. By default a cubic backtracking line search 37 is employed, as described in the text "Numerical Methods for 38 Unconstrained Optimization and Nonlinear Equations" by Dennis 39 and Schnabel. See the examples in nonlin/snes/examples. 40 D*/ 41 /* 42 This is intended as a model implementation, since it does not 43 necessarily have many of the bells and whistles of other 44 implementations. 45 46 The code is DATA-STRUCTURE NEUTRAL and can be called RECURSIVELY. 47 The following context variable is used: 48 NLCtx *nlP - The nonlinear solver context, which is created by 49 calling NLCreate(NLE_NLS1). 50 */ 51 52 int NLENewtonLS1Solve( nlP ) 53 NLCtx *nlP; 54 { 55 NLENewtonLS1Ctx *neP = (NLENewtonLS1Ctx *) nlP->MethodPrivate; 56 int maxits, i, iters, line, nlconv, history_len; 57 double fnorm, gnorm, gpnorm, xnorm, ynorm, *history; 58 double t_elp, t_cpu; 59 void *Y, *X, *F, *G, *W, *TMP; 60 FILE *fp = nlP->fp; 61 NLMonCore *mc = &nlP->mon.core; 62 VECntx *vc = NLGetVectorCtx(nlP); 63 64 CHKCOOKIEN(nlP,NL_COOKIE); 65 nlconv = 0; /* convergence monitor */ 66 history = nlP->conv_hist; /* convergence history */ 67 history_len = nlP->conv_hist_len; /* convergence history length */ 68 maxits = nlP->max_its; /* maximum number of iterations */ 69 X = nlP->vec_sol; /* solution vector */ 70 F = nlP->vec_rg; /* residual vector */ 71 Y = nlP->work[0]; /* work vectors */ 72 G = nlP->work[1]; 73 W = nlP->work[2]; 74 75 INITIAL_GUESS( X ); /* X <- X_0 */ 76 VNORM( vc, X, &xnorm ); /* xnorm = || X || */ 77 RESIDUAL( X, F ); /* Evaluate (+/-) F(X) */ 78 VNORM( vc, F, &fnorm ); /* fnorm <- || F || */ 79 nlP->norm = fnorm; 80 if (history && history_len > 0) history[0] = fnorm; 81 mc->nvectors += 4; mc->nscalars += 2; 82 MONITOR( X, F, &fnorm ); /* Monitor progress */ 83 84 for ( i=0; i<maxits; i++ ) { 85 nlP->iter = i+1; 86 87 /* Solve J Y = -F, where J is Jacobian matrix */ 88 STEP_SETUP( X ); /* Step set-up phase */ 89 iters = STEP_COMPUTE( X, F, Y, &fnorm, &(neP->maxstep), 90 &(nlP->trunctol), &gpnorm, &ynorm, (void *)0 ); 91 CHKERRV(1,-(NL)); /* Step compute phase, 92 excluding line search */ 93 94 /* Line search should really be part of step compute phase */ 95 line = (*neP->line_search)( nlP, X, F, G, Y, W, fnorm, 96 &ynorm, &gnorm ); CHKERRV(1,-(NL)); 97 98 if (!line) nlP->mon.nunsuc++; 99 if (fp) fprintf(fp,"%d: fnorm=%g, gnorm=%g, ynorm=%g, iters=%d\n", 100 nlP->iter, fnorm, gnorm, ynorm, iters ); 101 TMP = F; F = G; G = TMP; 102 TMP = X; X = Y; Y = TMP; 103 fnorm = gnorm; 104 105 STEP_DESTROY(); /* Step destroy phase */ 106 nlP->norm = fnorm; 107 if (history && history_len > i+1) history[i+1] = fnorm; 108 VNORM( vc, X, &xnorm ); /* xnorm = || X || */ 109 mc->nvectors += 2; 110 mc->nscalars++; 111 MONITOR( X, F, &fnorm ); /* Monitor progress */ 112 113 /* Test for convergence */ 114 if (CONVERGED( &xnorm, &ynorm, &fnorm )) { 115 /* Verify that solution is in corect location */ 116 if (X != nlP->vec_sol) VCOPY( vc, X, nlP->vec_sol ); 117 break; 118 } 119 } 120 if (i == maxits) i--; 121 return i+1; 122 } 123 /* ------------------------------------------------------------ */ 124 /*ARGSUSED*/ 125 void NLENewtonLS1SetUp( nlP ) 126 NLCtx *nlP; 127 { 128 NLENewtonLS1Ctx *ctx = (NLENewtonLS1Ctx *)nlP->MethodPrivate; 129 130 CHKCOOKIE(nlP,NL_COOKIE); 131 nlP->nwork = 3; 132 nlP->work = VGETVECS( nlP->vc, nlP->nwork ); CHKPTR(nlP->work); 133 if (!ctx->line_search) { 134 SETERRC(1,"NLENewtonLS1SetUp needs line search routine!\n"); 135 return; 136 } 137 NLiBasicSetUp( nlP, "NLENewtonLS1SetUp" ); CHKERR(1); 138 } 139 /* ------------------------------------------------------------ */ 140 void NLENewtonLS1Create( nlP ) 141 NLCtx *nlP; 142 { 143 NLENewtonLS1Ctx *neP; 144 145 CHKCOOKIE(nlP,NL_COOKIE); 146 nlP->method = NLE_NLS1; 147 nlP->method_type = NLE; 148 nlP->setup = NLENewtonLS1SetUp; 149 nlP->solver = NLENewtonLS1Solve; 150 nlP->destroy = NLENewtonLS1Destroy; 151 nlP->set_param = NLENewtonLS1SetParameter; 152 nlP->get_param = NLENewtonLS1GetParameter; 153 nlP->usr_monitor = NLENewtonDefaultMonitor; 154 nlP->converged = NLENewtonDefaultConverged; 155 nlP->term_type = NLENewtonDefaultConvergedType; 156 157 neP = NEW(NLENewtonLS1Ctx); CHKPTR(neP); 158 nlP->MethodPrivate = (void *) neP; 159 neP->line_search = NLStepDefaultLineSearch; 160 neP->alpha = 1.e-4; 161 neP->maxstep = 1.e8; 162 neP->steptol = 1.e-12; 163 } 164 /* ------------------------------------------------------------ */ 165 /*ARGSUSED*/ 166 void NLENewtonLS1Destroy( nlP ) 167 NLCtx *nlP; 168 { 169 VFREEVECS( nlP->vc, nlP->work, nlP->nwork ); 170 NLiBasicDestroy( nlP ); CHKERR(1); 171 } 172 /* ------------------------------------------------------------ */ 173 /*ARGSUSED*/ 174 /*@ 175 NLStepSimpleLineSearch - This routine is not a line search at all; 176 it simply uses the full Newton step. Thus, this routine is intended 177 to serve as a template and is not recommended for general use. 178 179 Input Parameters: 180 . nlP - nonlinear context 181 . x - current iterate 182 . f - residual evaluated at x 183 . y - search direction (contains new iterate on output) 184 . w - work vector 185 . fnorm - 2-norm of f 186 187 Output parameters: 188 . g - residual evaluated at new iterate y 189 . y - new iterate (contains search direction on input) 190 . gnorm - 2-norm of g 191 . ynorm - 2-norm of search length 192 193 Returns: 194 1, indicating success of the step. 195 196 Notes: 197 Use either NLSetStepLineSearchRoutines() or NLSetLineSearchRoutine() 198 to set this routine within the NLE_NLS1 method. 199 @*/ 200 int NLStepSimpleLineSearch( nlP, x, f, g, y, w, fnorm, 201 ynorm, gnorm ) 202 NLCtx *nlP; 203 void *x; 204 void *f; 205 void *g; 206 void *y; 207 void *w; 208 double fnorm; 209 double *ynorm; 210 double *gnorm; 211 { 212 NLMonCore *mc = &nlP->mon.core; 213 VECntx *vc = NLGetVectorCtx(nlP); 214 215 CHKCOOKIEN(nlP,NL_COOKIE); 216 VNORM( vc, y, ynorm ); /* ynorm = || y || */ 217 VAXPY( vc, 1.0, x, y ); /* y <- x + y */ 218 RESIDUAL( y, g ); /* Evaluate (+/-) g(y) */ 219 VNORM( vc, g, gnorm ); /* gnorm = || g || */ 220 mc->nvectors += 6; 221 mc->nscalars += 2; 222 return 1; 223 } 224 /* ------------------------------------------------------------------ */ 225 /*@ 226 NLStepDefaultLineSearch - This routine performs a cubic line search. 227 228 Input Parameters: 229 . nlP - nonlinear context 230 . x - current iterate 231 . f - residual evaluated at x 232 . y - search direction (contains new iterate on output) 233 . w - work vector 234 . fnorm - 2-norm of f 235 236 Output parameters: 237 . g - residual evaluated at new iterate y 238 . y - new iterate (contains search direction on input) 239 . gnorm - 2-norm of g 240 . ynorm - 2-norm of search length 241 242 Returns: 243 1 if the line search succeeds; 0 if the line search fails. 244 245 Notes: 246 Use either NLSetStepLineSearchRoutines() or NLSetLineSearchRoutine() 247 to set this routine within the NLE_NLS1 method. 248 249 This line search is taken from "Numerical Methods for Unconstrained 250 Optimization and Nonlinear Equations" by Dennis and Schnabel, page 325. 251 252 @*/ 253 int NLStepDefaultLineSearch( nlP, x, f, g, y, w, fnorm, ynorm, gnorm ) 254 NLCtx *nlP; 255 void *x; 256 void *f; 257 void *g; 258 void *y; 259 void *w; 260 double fnorm; 261 double *ynorm; 262 double *gnorm; 263 { 264 double alpha, maxstep, steptol, initslope; 265 double minlambda, lambda, lambdaprev, gnormprev, lambdatemp; 266 double a, b, d, t1, t2; 267 int count; 268 FILE *fp = nlP->fp; 269 NLENewtonLS1Ctx *neP = (NLENewtonLS1Ctx *) nlP->MethodPrivate; 270 NLMonCore *mc = &nlP->mon.core; 271 VECntx *vc = NLGetVectorCtx(nlP); 272 273 CHKCOOKIEN(nlP,NL_COOKIE); 274 alpha = neP->alpha; 275 maxstep = neP->maxstep; 276 steptol = neP->steptol; 277 278 VNORM( vc, y, ynorm ); 279 if (*ynorm > maxstep) { /* Step too big, so scale back */ 280 VSCALE( vc, maxstep/(*ynorm), y ); 281 *ynorm = maxstep; 282 mc->nvectors++; 283 } 284 minlambda = steptol/(*ynorm); 285 VDOT( vc, f, y, &initslope ); 286 if (initslope > 0.0) initslope = -initslope; 287 if (initslope == 0.0) initslope = -1.0; 288 289 VCOPY( vc, y, w ); 290 VAXPY( vc, 1.0, x, w ); 291 RESIDUAL( w, g ); /* Evaluate (+/-) g(w) */ 292 VNORM( vc, g, gnorm ); 293 if (*gnorm <= fnorm + alpha*initslope) { /* Sufficient reduction */ 294 if (fp) fprintf(fp,"Taking full Newton step\n"); 295 VCOPY( vc, w, y ); 296 mc->nvectors += 8; mc->nscalars += 3; 297 return 1; 298 } 299 300 /* Fit points with quadratic */ 301 lambda = 1.0; count = 0; 302 lambdatemp = -initslope/(2.0*(*gnorm - fnorm - initslope)); 303 lambdaprev = lambda; 304 gnormprev = *gnorm; 305 if (lambdatemp <= .1*lambda) { 306 lambda = .1*lambda; 307 mc->nscalars++; 308 } else lambda = lambdatemp; 309 VCOPY( vc, x, w ); 310 VAXPY( vc, lambda, y, w ); 311 RESIDUAL( w, g ); /* Evaluate (+/-) g(w) */ 312 VNORM( vc, g, gnorm ); 313 if (*gnorm <= fnorm + alpha*initslope) { /* sufficient reduction */ 314 if (fp) fprintf(fp,"Taking Newton step from quadratic \n"); 315 VCOPY( vc, w, y ); 316 mc->nvectors += 12; mc->nscalars += 8; 317 return 1; 318 } 319 320 /* Fit points with cubic */ 321 count = 1; 322 while (1) { 323 if (lambda <= minlambda) { /* bad luck; use full step */ 324 fprintf(stderr,"Unable to find good step length! %d \n",count); 325 fprintf(stderr, "f %g fnew %g ynorm %g lambda %g \n", 326 fnorm,*gnorm, *ynorm,lambda); 327 VCOPY( vc, w, y ); 328 mc->nvectors += 12 + 4*(count-1); 329 mc->nscalars += 8 + 28*(count-1); 330 return 0; 331 } 332 t1 = *gnorm - fnorm - lambda*initslope; 333 t2 = gnormprev - fnorm - lambdaprev*initslope; 334 a = (t1/(lambda*lambda) - 335 t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 336 b = (-lambdaprev*t1/(lambda*lambda) + 337 lambda*t2/(lambdaprev*lambdaprev))/(lambda-lambdaprev); 338 d = b*b - 3*a*initslope; 339 if (d < 0.0) d = 0.0; 340 if (a == 0.0) { 341 lambdatemp = -initslope/(2.0*b); 342 mc->nscalars += 2; 343 } else { 344 lambdatemp = (-b + sqrt(d))/(3.0*a); 345 mc->nscalars += 4; 346 } 347 if (lambdatemp > .5*lambda) { 348 lambdatemp = .5*lambda; 349 mc->nscalars++; 350 } 351 lambdaprev = lambda; 352 gnormprev = *gnorm; 353 if (lambdatemp <= .1*lambda) { 354 lambda = .1*lambda; 355 mc->nscalars++; 356 } 357 else lambda = lambdatemp; 358 VCOPY( vc, x, w ); 359 VAXPY( vc, lambda, y, w ); 360 RESIDUAL( w, g ); /* Evaluate (+/-) g(w) */ 361 VNORM( vc, g, gnorm ); 362 if (*gnorm <= fnorm + alpha*initslope) { /* is reduction enough */ 363 if (fp) fprintf(fp,"Taking Newton step from cubic %d\n",count); 364 VCOPY( vc, w, y ); 365 mc->nvectors += 12 + 4*count; 366 mc->nscalars += 8 + 28*count; 367 return 1; 368 } 369 count++; 370 } 371 } 372 /* ------------------------------------------------------------ */ 373 /* 374 NLENewtonLS1SetParameter - Sets a chosen parameter used by the 375 NLE_NLS1 method to the specified value. 376 377 Note: 378 Possible parameters for the NLE_NLS1 method are 379 $ param = "alpha" - used to determine sufficient reduction 380 $ param = "maxstep" - maximum step size 381 $ param = "steptol" - step comvergence tolerance 382 */ 383 void NLENewtonLS1SetParameter( nlP, param, value ) 384 NLCtx *nlP; 385 char *param; 386 double *value; 387 { 388 NLENewtonLS1Ctx *ctx = (NLENewtonLS1Ctx *)nlP->MethodPrivate; 389 390 CHKCOOKIE(nlP,NL_COOKIE); 391 if (nlP->method != NLE_NLS1) { 392 SETERRC(1,"Compatible only with NLE_NLS1 method"); 393 return; 394 } 395 if (!strcmp(param,"alpha")) ctx->alpha = *value; 396 else if (!strcmp(param,"maxstep")) ctx->maxstep = *value; 397 else if (!strcmp(param,"steptol")) ctx->steptol = *value; 398 else SETERRC(1,"Invalid parameter name for NLE_NLS1 method"); 399 } 400 /* ------------------------------------------------------------ */ 401 /* 402 NLENewtonLS1GetParameter - Returns the value of a chosen parameter 403 used by the NLE_NLS1 method. 404 405 Note: 406 Possible parameters for the NLE_NLS1 method are 407 $ param = "alpha" - used to determine sufficient reduction 408 $ param = "maxstep" - maximum step size 409 $ param = "steptol" - step comvergence tolerance 410 */ 411 double NLENewtonLS1GetParameter( nlP, param ) 412 NLCtx *nlP; 413 char *param; 414 { 415 NLENewtonLS1Ctx *ctx = (NLENewtonLS1Ctx *)nlP->MethodPrivate; 416 double value = 0.0; 417 418 CHKCOOKIEN(nlP,NL_COOKIE); 419 if (nlP->method != NLE_NLS1) { 420 SETERRC(1,"Compatible only with NLE_NLS1 method"); 421 return value; 422 } 423 if (!strcmp(param,"alpha")) value = ctx->alpha; 424 else if (!strcmp(param,"maxstep")) value = ctx->maxstep; 425 else if (!strcmp(param,"steptol")) value = ctx->steptol; 426 else SETERRC(1,"Invalid parameter name for NLE_NLS1 method"); 427 return value; 428 } 429 /* ------------------------------------------------------------ */ 430 /*@C 431 NLSetLineSearchRoutine - Sets the line search routine to be used 432 by the method NLE_NLS1. 433 434 Input Parameters: 435 . nlP - nonlinear context obtained from NLCreate() 436 . func - pointer to int function 437 438 Possible routines: 439 NLStepDefaultLineSearch() - default line search 440 NLStepSimpleLineSearch() - the full Newton step (actually not a 441 line search) 442 443 Calling sequence of func: 444 . func (NLCtx *nlP, void *x, void *f, void *g, void *y, 445 void *w, double fnorm, double *ynorm, double *gnorm ) 446 447 Input parameters for func: 448 . nlP - nonlinear context 449 . x - current iterate 450 . f - residual evaluated at x 451 . y - search direction (contains new iterate on output) 452 . w - work vector 453 . fnorm - 2-norm of f 454 455 Output parameters for func: 456 . g - residual evaluated at new iterate y 457 . y - new iterate (contains search direction on input) 458 . gnorm - 2-norm of g 459 . ynorm - 2-norm of search length 460 461 Returned by func: 462 1 if the line search succeeds; 0 if the line search fails. 463 @*/ 464 void NLSetLineSearchRoutine( nlP, func ) 465 NLCtx *nlP; 466 int (*func)(); 467 { 468 CHKCOOKIE(nlP,NL_COOKIE); 469 if ((nlP)->method == NLE_NLS1) 470 ((NLENewtonLS1Ctx *)(nlP->MethodPrivate))->line_search = func; 471 } 472