1 #ifndef lint 2 static char vcid[] = "$Id: tr.c,v 1.7 1995/05/12 04:18:49 bsmith Exp bsmith $"; 3 #endif 4 5 #include <math.h> 6 #include "tr.h" 7 8 9 /* 10 Implements Newton's Method with a very simple trust region 11 approach for solving systems of nonlinear equations. 12 13 Input parameters: 14 . nlP - nonlinear context obtained from SNESCreate() 15 16 The basic algorithm is taken from "The Minpack Project", by More', 17 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 18 of Mathematical Software", Wayne Cowell, editor. See the examples 19 in nonlin/examples. 20 */ 21 /* 22 This is intended as a model implementation, since it does not 23 necessarily have many of the bells and whistles of other 24 implementations. 25 26 */ 27 static int SNESSolve_TR(SNES snes, int *its ) 28 { 29 SNES_TR *neP = (SNES_TR *) snes->data; 30 Vec X, F, Y, G, TMP, Ytmp; 31 int maxits, i, history_len, nlconv,ierr,lits; 32 MatStructure flg; 33 double rho, fnorm, gnorm, gpnorm, xnorm, delta,norm; 34 double *history, ynorm; 35 Scalar one = 1.0,cnorm; 36 double epsmch = 1.0e-14; /* This must be fixed */ 37 38 nlconv = 0; /* convergence monitor */ 39 history = snes->conv_hist; /* convergence history */ 40 history_len = snes->conv_hist_len; /* convergence history length */ 41 maxits = snes->max_its; /* maximum number of iterations */ 42 X = snes->vec_sol; /* solution vector */ 43 F = snes->vec_func; /* residual vector */ 44 Y = snes->work[0]; /* work vectors */ 45 G = snes->work[1]; 46 Ytmp = snes->work[2]; 47 48 ierr = SNESComputeInitialGuess(snes,X); CHKERR(ierr); /* X <- X_0 */ 49 VecNorm(X, &xnorm ); /* xnorm = || X || */ 50 51 ierr = SNESComputeFunction(snes,X,F); CHKERR(ierr); /* (+/-) F(X) */ 52 VecNorm(F, &fnorm ); /* fnorm <- || F || */ 53 snes->norm = fnorm; 54 if (history && history_len > 0) history[0] = fnorm; 55 delta = neP->delta0*fnorm; 56 neP->delta = delta; 57 if (snes->Monitor)(*snes->Monitor)(snes,0,fnorm,snes->monP); 58 59 for ( i=0; i<maxits; i++ ) { 60 snes->iter = i+1; 61 62 (*snes->ComputeJacobian)(snes,X,&snes->jacobian,&snes->jacobian_pre, 63 &flg,snes->jacP); 64 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg); 65 ierr = SLESSolve(snes->sles,F,Ytmp,&lits); CHKERR(ierr); 66 VecNorm( Ytmp, &norm ); 67 while(1) { 68 VecCopy(Ytmp,Y); 69 /* Scale Y if need be and predict new value of F norm */ 70 71 if (norm >= delta) { 72 norm = delta/norm; 73 gpnorm = (1.0 - norm)*fnorm; 74 cnorm = norm; 75 VecScale( &cnorm, Y ); 76 norm = gpnorm; 77 PLogInfo((PetscObject)snes, "Scaling direction by %g \n",norm ); 78 ynorm = delta; 79 } else { 80 gpnorm = 0.0; 81 PLogInfo((PetscObject)snes,"Direction is in Trust Region \n" ); 82 ynorm = norm; 83 } 84 VecAXPY(&one, X, Y ); /* Y <- X + Y */ 85 ierr = SNESComputeFunction(snes,Y,G); CHKERR(ierr); /* (+/-) F(X) */ 86 VecNorm( G, &gnorm ); /* gnorm <- || g || */ 87 if (fnorm == gpnorm) rho = 0.0; 88 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 89 90 /* Update size of trust region */ 91 if (rho < neP->mu) delta *= neP->delta1; 92 else if (rho < neP->eta) delta *= neP->delta2; 93 else delta *= neP->delta3; 94 95 PLogInfo((PetscObject)snes,"%d: f_norm=%g, g_norm=%g, ynorm=%g\n", 96 i, fnorm, gnorm, ynorm ); 97 PLogInfo((PetscObject)snes,"gpred=%g, rho=%g, delta=%g,iters=%d\n", 98 gpnorm, rho, delta, lits); 99 100 neP->delta = delta; 101 if (rho > neP->sigma) break; 102 PLogInfo((PetscObject)snes,"Trying again in smaller region\n"); 103 /* check to see if progress is hopeless */ 104 if (neP->delta < xnorm * epsmch) return -1; 105 } 106 fnorm = gnorm; 107 snes->norm = fnorm; 108 if (history && history_len > i+1) history[i+1] = fnorm; 109 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 110 TMP = X; X = Y;snes->vec_sol_always = X; Y = TMP; 111 VecNorm(X, &xnorm ); /* xnorm = || X || */ 112 if (snes->Monitor) (*snes->Monitor)(snes,i,fnorm,snes->monP); 113 114 /* Test for convergence */ 115 if ((*snes->Converged)( snes, xnorm, ynorm, fnorm,snes->cnvP )) { 116 /* Verify solution is in corect location */ 117 if (X != snes->vec_sol) { 118 VecCopy(X, snes->vec_sol ); 119 snes->vec_sol_always = snes->vec_sol; 120 snes->vec_func_always = snes->vec_func; 121 } 122 break; 123 } 124 } 125 if (i == maxits) *its = i-1; else *its = i; 126 return 0; 127 } 128 /* -------------------------------------------------------------*/ 129 130 /*------------------------------------------------------------*/ 131 static int SNESSetUp_TR( SNES snes ) 132 { 133 int ierr; 134 snes->nwork = 3; 135 ierr = VecGetVecs(snes->vec_sol,snes->nwork,&snes->work ); CHKERR(ierr); 136 return 0; 137 } 138 /*------------------------------------------------------------*/ 139 static int SNESDestroy_TR(PetscObject obj ) 140 { 141 SNES snes = (SNES) obj; 142 VecFreeVecs(snes->work, snes->nwork ); 143 return 0; 144 } 145 /*------------------------------------------------------------*/ 146 147 #include "options.h" 148 static int SNESSetFromOptions_TR(SNES snes) 149 { 150 SNES_TR *ctx = (SNES_TR *)snes->data; 151 double tmp; 152 153 if (OptionsGetDouble(0,snes->prefix,"-mu",&tmp)) {ctx->mu = tmp;} 154 if (OptionsGetDouble(0,snes->prefix,"-eta",&tmp)) {ctx->eta = tmp;} 155 if (OptionsGetDouble(0,snes->prefix,"-sigma",&tmp)) {ctx->sigma = tmp;} 156 if (OptionsGetDouble(0,snes->prefix,"-delta0",&tmp)) {ctx->delta0 = tmp;} 157 if (OptionsGetDouble(0,snes->prefix,"-delta1",&tmp)) {ctx->delta1 = tmp;} 158 if (OptionsGetDouble(0,snes->prefix,"-delta2",&tmp)) {ctx->delta2 = tmp;} 159 if (OptionsGetDouble(0,snes->prefix,"-delta3",&tmp)) {ctx->delta3 = tmp;} 160 return 0; 161 } 162 163 static int SNESPrintHelp_TR(SNES snes) 164 { 165 SNES_TR *ctx = (SNES_TR *)snes->data; 166 char *prefix = "-"; 167 if (snes->prefix) prefix = snes->prefix; 168 fprintf(stderr,"%smu mu (default %g)\n",prefix,ctx->mu); 169 fprintf(stderr,"%seta eta (default %g)\n",prefix,ctx->eta); 170 fprintf(stderr,"%ssigma sigma (default %g)\n",prefix,ctx->sigma); 171 fprintf(stderr,"%sdelta0 delta0 (default %g)\n",prefix,ctx->delta0); 172 fprintf(stderr,"%sdelta1 delta1 (default %g)\n",prefix,ctx->delta1); 173 fprintf(stderr,"%sdelta2 delta2 (default %g)\n",prefix,ctx->delta2); 174 fprintf(stderr,"%sdelta3 delta3 (default %g)\n",prefix,ctx->delta3); 175 return 0; 176 } 177 178 int SNESCreate_TR(SNES snes ) 179 { 180 SNES_TR *neP; 181 182 snes->type = SNES_NTR; 183 snes->setup = SNESSetUp_TR; 184 snes->solve = SNESSolve_TR; 185 snes->destroy = SNESDestroy_TR; 186 snes->Converged = SNESDefaultConverged; 187 snes->printhelp = SNESPrintHelp_TR; 188 snes->setfromoptions = SNESSetFromOptions_TR; 189 190 neP = NEW(SNES_TR); CHKPTR(neP); 191 snes->data = (void *) neP; 192 neP->mu = 0.25; 193 neP->eta = 0.75; 194 neP->delta = 0.0; 195 neP->delta0 = 0.2; 196 neP->delta1 = 0.3; 197 neP->delta2 = 0.75; 198 neP->delta3 = 2.0; 199 neP->sigma = 0.0001; 200 neP->itflag = 0; 201 return 0; 202 } 203