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