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