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