1 #ifndef lint 2 static char vcid[] = "$Id: tr.c,v 1.1 1995/04/12 20:36:40 bsmith Exp bsmith $"; 3 #endif 4 5 #include <math.h> 6 #include "tr.h" 7 8 /* 9 NLE_NTR1 - Implements Newton's Method with a trust region 10 approach for solving systems of nonlinear equations. 11 12 Input parameters: 13 . nlP - nonlinear context obtained from NLCreate() 14 15 Returns: 16 Number of global iterations until termination. The precise type of 17 termination can be examined by calling NLGetTerminationType() after 18 NLSolve(). 19 20 Calling sequence: 21 $ nlP = NLCreate(NLE_NTR1,0) 22 $ NLCreateDVectors() 23 $ NLSet***() 24 $ NLSetUp() 25 $ NLSolve() 26 $ NLDestroy() 27 28 Notes: 29 See NLCreate() and NLSetUp() for information on the definition and 30 initialization of the nonlinear solver context. 31 32 The basic algorithm is taken from "The Minpack Project", by More', 33 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 34 of Mathematical Software", Wayne Cowell, editor. See the examples 35 in nonlin/examples. 36 */ 37 /* 38 This is intended as a model implementation, since it does not 39 necessarily have many of the bells and whistles of other 40 implementations. 41 42 The code is DATA-STRUCTURE NEUTRAL and can be called RECURSIVELY. 43 The following context variable is used: 44 NLCtx *nlP - The nonlinear solver context, which is created by 45 calling NLCreate(NLE_NTR1). 46 47 The step_compute routine must return two values: 48 1) ynorm - the norm of the step 49 2) gpnorm - the predicted value for the residual norm at the new 50 point, assuming a local linearization. The value is 0 if the 51 step lies within the trust region and is > 0 otherwise. 52 */ 53 static int SNESSolve_TR(SNES snes, int *its ) 54 { 55 SNES_TR *neP = (SNES_TR *) snes->data; 56 Vec X, F, Y, G, TMP; 57 int maxits, i, iters, history_len, nlconv,ierr,lits; 58 double rho, fnorm, gnorm, gpnorm, xnorm, delta,norm; 59 double *history, ynorm; 60 Scalar one = 1.0; 61 62 nlconv = 0; /* convergence monitor */ 63 history = snes->conv_hist; /* convergence history */ 64 history_len = snes->conv_hist_len; /* convergence history length */ 65 maxits = snes->max_its; /* maximum number of iterations */ 66 X = snes->vec_sol; /* solution vector */ 67 F = snes->vec_res; /* residual vector */ 68 Y = snes->work[0]; /* work vectors */ 69 G = snes->work[1]; 70 71 ierr = SNESComputeInitialGuess(snes,X); CHKERR(ierr); /* X <- X_0 */ 72 VecNorm(X, &xnorm ); /* xnorm = || X || */ 73 74 ierr = SNESComputeResidual(snes,X,F); CHKERR(ierr); /* (+/-) F(X) */ 75 VecNorm(F, &fnorm ); /* fnorm <- || F || */ 76 snes->norm = fnorm; 77 if (history && history_len > 0) history[0] = fnorm; 78 delta = neP->delta0*fnorm; 79 neP->delta = delta; 80 if (snes->Monitor)(*snes->Monitor)(snes,0,X,F,fnorm,snes->monP); 81 82 for ( i=0; i<maxits; i++ ) { 83 snes->iter = i+1; 84 85 (*snes->ComputeJacobian)(X,&snes->jacobian,snes->jacP); 86 while(1) { 87 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian,0); 88 ierr = SLESSolve(snes->sles,F,Y,&lits); CHKERR(ierr); 89 /* Scale Y if need be and predict new value of F norm */ 90 VecNorm( Y, &norm ); 91 if (norm > delta) { 92 norm = delta/norm; 93 gpnorm = (1.0 - norm)*fnorm; 94 VecScale( &norm, Y ); 95 PLogInfo((PetscObject)snes, "Scaling direction by %g \n",norm ); 96 ynorm = delta; 97 } else { 98 gpnorm = 0.0; 99 PLogInfo((PetscObject)snes,"Direction is in Trust Region \n" ); 100 ynorm = norm; 101 } 102 VecAXPY(&one, X, Y ); /* Y <- X + Y */ 103 ierr = SNESComputeResidual(snes,Y,G); CHKERR(ierr); /* (+/-) F(X) */ 104 VecNorm( G, &gnorm ); /* gnorm <- || g || */ 105 if (fnorm == gpnorm) rho = 0.0; 106 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 107 108 /* Update size of trust region */ 109 if (rho < neP->mu) delta *= neP->delta1; 110 else if (rho < neP->eta) delta *= neP->delta2; 111 else delta *= neP->delta3; 112 113 PLogInfo((PetscObject)snes,"%d: f_norm=%g, g_norm=%g, ynorm=%g\n", 114 i, fnorm, gnorm, ynorm ); 115 PLogInfo((PetscObject)snes,"gpred=%g, rho=%g, delta=%g,iters=%d\n", 116 gpnorm, rho, delta, lits); 117 118 neP->delta = delta; 119 if (rho > neP->sigma) break; 120 neP->itflag = 0; 121 if ((*snes->Converged)( snes, xnorm, ynorm, fnorm,snes->cnvP)) { 122 /* We're not progressing, so return with the current iterate */ 123 if (X != snes->vec_sol) VecCopy( X, snes->vec_sol ); 124 return i; 125 } 126 } 127 fnorm = gnorm; 128 snes->norm = fnorm; 129 if (history && history_len > i+1) history[i+1] = fnorm; 130 TMP = F; F = G; G = TMP; 131 TMP = X; X = Y; Y = TMP; 132 VecNorm(X, &xnorm ); /* xnorm = || X || */ 133 if (snes->Monitor) (*snes->Monitor)(snes,0,X,F,fnorm,snes->monP); 134 135 /* Test for convergence */ 136 neP->itflag = 1; 137 if ((*snes->Converged)( snes, xnorm, ynorm, fnorm,snes->cnvP )) { 138 /* Verify solution is in corect location */ 139 if (X != snes->vec_sol) VecCopy(X, snes->vec_sol ); 140 break; 141 } 142 } 143 if (i == maxits) *its = i-1; else *its = i; 144 return 0; 145 } 146 /* -------------------------------------------------------------*/ 147 148 /*------------------------------------------------------------*/ 149 static int SNESSetUp_TR( SNES snes ) 150 { 151 int ierr; 152 snes->nwork = 2; 153 ierr = VecGetVecs(snes->vec_sol,snes->nwork,&snes->work ); CHKERR(ierr); 154 return 0; 155 } 156 /*------------------------------------------------------------*/ 157 static int SNESDestroy_TR(PetscObject obj ) 158 { 159 SNES snes = (SNES) obj; 160 VecFreeVecs(snes->work, snes->nwork ); 161 return 0; 162 } 163 /*------------------------------------------------------------*/ 164 /*@ 165 SNESTRDefaultConverged - Default test for monitoring the 166 convergence of the method NLENewtonTR1Solve. 167 168 Input Parameters: 169 . snes - nonlinear context obtained from NLCreate() 170 . xnorm - 2-norm of current iterate 171 . pnorm - 2-norm of current step 172 . fnorm - 2-norm of residual 173 174 Returns: 175 $ 1 if ( delta < xnorm*deltatol ), 176 $ 2 if ( fnorm < atol ), 177 $ 3 if ( pnorm < xtol*xnorm ), 178 $ -2 if ( nres > max_res ), 179 $ -1 if ( delta < xnorm*epsmch ), 180 $ 0 otherwise, 181 182 where 183 $ atol - absolute residual norm tolerance, 184 $ set with NLSetAbsConvergenceTol() 185 $ delta - trust region paramenter 186 $ deltatol - trust region size tolerance, 187 $ set with NLSetTrustRegionTol() 188 $ epsmch - machine epsilon 189 $ max_res - maximum number of residual evaluations, 190 $ set with NLSetMaxResidualEvaluations() 191 $ nres - number of residual evaluations 192 $ xtol - relative residual norm tolerance, 193 $ set with NLSetRelConvergenceTol() 194 195 Note: 196 Call NLGetConvergenceType() after calling NLSolve() to obtain 197 information about the type of termination which occurred for the 198 nonlinear solver. 199 @*/ 200 int SNESTRDefaultConverged(SNES snes, double xnorm, double pnorm, 201 double fnorm, void *ctx ) 202 { 203 SNES_TR *neP = (SNES_TR *)snes->data; 204 double epsmch = 1.0e-14; /* This must be fixed */ 205 206 if (neP->delta < xnorm * neP->deltatol) return 1; 207 if (neP->itflag) { 208 return SNESDefaultConverged( snes, xnorm, pnorm, fnorm,ctx ); 209 } 210 if (neP->delta < xnorm * epsmch) return -1; 211 return 0; 212 } 213 214 static int SNESSetFromOptions_TR(SNES snes) 215 { 216 SNES_TR *ctx = (SNES_TR *)snes->data; 217 double tmp; 218 219 if (OptionsGetDouble(0,snes->prefix,"-mu",&tmp)) {ctx->mu = tmp;} 220 if (OptionsGetDouble(0,snes->prefix,"-eta",&tmp)) {ctx->eta = tmp;} 221 if (OptionsGetDouble(0,snes->prefix,"-sigma",&tmp)) {ctx->sigma = tmp;} 222 if (OptionsGetDouble(0,snes->prefix,"-delta0",&tmp)) {ctx->delta0 = tmp;} 223 if (OptionsGetDouble(0,snes->prefix,"-delta1",&tmp)) {ctx->delta1 = tmp;} 224 if (OptionsGetDouble(0,snes->prefix,"-delta2",&tmp)) {ctx->delta2 = tmp;} 225 if (OptionsGetDouble(0,snes->prefix,"-delta3",&tmp)) {ctx->delta3 = tmp;} 226 return 0; 227 } 228 229 static int SNESPrintHelp_TR(SNES snes) 230 { 231 SNES_TR *ctx = (SNES_TR *)snes->data; 232 char *prefix = "-"; 233 if (snes->prefix) prefix = snes->prefix; 234 fprintf(stderr,"%smu mu (default %g)\n",prefix,ctx->mu); 235 fprintf(stderr,"%seta eta (default %g)\n",prefix,ctx->eta); 236 fprintf(stderr,"%ssigma sigma (default %g)\n",prefix,ctx->sigma); 237 fprintf(stderr,"%sdelta0 delta0 (default %g)\n",prefix,ctx->delta0); 238 fprintf(stderr,"%sdelta1 delta1 (default %g)\n",prefix,ctx->delta1); 239 fprintf(stderr,"%sdelta2 delta2 (default %g)\n",prefix,ctx->delta2); 240 fprintf(stderr,"%sdelta3 delta3 (default %g)\n",prefix,ctx->delta3); 241 return 0; 242 } 243 244 int SNESCreate_TR(SNES snes ) 245 { 246 SNES_TR *neP; 247 248 snes->type = SNES_NTR; 249 snes->Setup = SNESSetUp_TR; 250 snes->Solver = SNESSolve_TR; 251 snes->destroy = SNESDestroy_TR; 252 snes->Monitor = 0; 253 snes->Converged = SNESTRDefaultConverged; 254 snes->PrintHelp = SNESPrintHelp_TR; 255 snes->SetFromOptions = SNESSetFromOptions_TR; 256 257 neP = NEW(SNES_TR); CHKPTR(neP); 258 snes->data = (void *) neP; 259 neP->mu = 0.25; 260 neP->eta = 0.75; 261 neP->delta = 0.0; 262 neP->delta0 = 0.2; 263 neP->delta1 = 0.3; 264 neP->delta2 = 0.75; 265 neP->delta3 = 2.0; 266 neP->sigma = 0.0001; 267 neP->itflag = 0; 268 } 269