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