1 2 #include <private/snesimpl.h> /*I "petscsnes.h" I*/ 3 #include <petscdmshell.h> /*I "petscdmshell.h" I*/ 4 5 PetscBool SNESRegisterAllCalled = PETSC_FALSE; 6 PetscFList SNESList = PETSC_NULL; 7 8 /* Logging support */ 9 PetscClassId SNES_CLASSID; 10 PetscLogEvent SNES_Solve, SNES_LineSearch, SNES_FunctionEval, SNES_JacobianEval, SNES_GSEval; 11 12 #undef __FUNCT__ 13 #define __FUNCT__ "SNESDMComputeJacobian" 14 /* 15 Translates from a SNES call to a DM call in computing a Jacobian 16 */ 17 PetscErrorCode SNESDMComputeJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr) 18 { 19 PetscErrorCode ierr; 20 DM dm; 21 22 PetscFunctionBegin; 23 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 24 ierr = DMComputeJacobian(dm,X,*J,*B,flag);CHKERRQ(ierr); 25 PetscFunctionReturn(0); 26 } 27 28 #undef __FUNCT__ 29 #define __FUNCT__ "SNESSetErrorIfNotConverged" 30 /*@ 31 SNESSetErrorIfNotConverged - Causes SNESSolve() to generate an error if the solver has not converged. 32 33 Logically Collective on SNES 34 35 Input Parameters: 36 + snes - iterative context obtained from SNESCreate() 37 - flg - PETSC_TRUE indicates you want the error generated 38 39 Options database keys: 40 . -snes_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false) 41 42 Level: intermediate 43 44 Notes: 45 Normally PETSc continues if a linear solver fails to converge, you can call SNESGetConvergedReason() after a SNESSolve() 46 to determine if it has converged. 47 48 .keywords: SNES, set, initial guess, nonzero 49 50 .seealso: SNESGetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged() 51 @*/ 52 PetscErrorCode SNESSetErrorIfNotConverged(SNES snes,PetscBool flg) 53 { 54 PetscFunctionBegin; 55 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 56 PetscValidLogicalCollectiveBool(snes,flg,2); 57 snes->errorifnotconverged = flg; 58 59 PetscFunctionReturn(0); 60 } 61 62 #undef __FUNCT__ 63 #define __FUNCT__ "SNESGetErrorIfNotConverged" 64 /*@ 65 SNESGetErrorIfNotConverged - Will SNESSolve() generate an error if the solver does not converge? 66 67 Not Collective 68 69 Input Parameter: 70 . snes - iterative context obtained from SNESCreate() 71 72 Output Parameter: 73 . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE 74 75 Level: intermediate 76 77 .keywords: SNES, set, initial guess, nonzero 78 79 .seealso: SNESSetErrorIfNotConverged(), KSPGetErrorIfNotConverged(), KSPSetErrorIFNotConverged() 80 @*/ 81 PetscErrorCode SNESGetErrorIfNotConverged(SNES snes,PetscBool *flag) 82 { 83 PetscFunctionBegin; 84 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 85 PetscValidPointer(flag,2); 86 *flag = snes->errorifnotconverged; 87 PetscFunctionReturn(0); 88 } 89 90 #undef __FUNCT__ 91 #define __FUNCT__ "SNESSetFunctionDomainError" 92 /*@ 93 SNESSetFunctionDomainError - tells SNES that the input vector to your FormFunction is not 94 in the functions domain. For example, negative pressure. 95 96 Logically Collective on SNES 97 98 Input Parameters: 99 . snes - the SNES context 100 101 Level: advanced 102 103 .keywords: SNES, view 104 105 .seealso: SNESCreate(), SNESSetFunction() 106 @*/ 107 PetscErrorCode SNESSetFunctionDomainError(SNES snes) 108 { 109 PetscFunctionBegin; 110 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 111 snes->domainerror = PETSC_TRUE; 112 PetscFunctionReturn(0); 113 } 114 115 116 #undef __FUNCT__ 117 #define __FUNCT__ "SNESGetFunctionDomainError" 118 /*@ 119 SNESGetFunctionDomainError - Gets the status of the domain error after a call to SNESComputeFunction; 120 121 Logically Collective on SNES 122 123 Input Parameters: 124 . snes - the SNES context 125 126 Output Parameters: 127 . domainerror Set to PETSC_TRUE if there's a domain error; PETSC_FALSE otherwise. 128 129 Level: advanced 130 131 .keywords: SNES, view 132 133 .seealso: SNESSetFunctionDomainError, SNESComputeFunction() 134 @*/ 135 PetscErrorCode SNESGetFunctionDomainError(SNES snes, PetscBool *domainerror) 136 { 137 PetscFunctionBegin; 138 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 139 PetscValidPointer(domainerror, 2); 140 *domainerror = snes->domainerror; 141 PetscFunctionReturn(0); 142 } 143 144 145 #undef __FUNCT__ 146 #define __FUNCT__ "SNESView" 147 /*@C 148 SNESView - Prints the SNES data structure. 149 150 Collective on SNES 151 152 Input Parameters: 153 + SNES - the SNES context 154 - viewer - visualization context 155 156 Options Database Key: 157 . -snes_view - Calls SNESView() at end of SNESSolve() 158 159 Notes: 160 The available visualization contexts include 161 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 162 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 163 output where only the first processor opens 164 the file. All other processors send their 165 data to the first processor to print. 166 167 The user can open an alternative visualization context with 168 PetscViewerASCIIOpen() - output to a specified file. 169 170 Level: beginner 171 172 .keywords: SNES, view 173 174 .seealso: PetscViewerASCIIOpen() 175 @*/ 176 PetscErrorCode SNESView(SNES snes,PetscViewer viewer) 177 { 178 SNESKSPEW *kctx; 179 PetscErrorCode ierr; 180 KSP ksp; 181 PetscBool iascii,isstring; 182 183 PetscFunctionBegin; 184 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 185 if (!viewer) { 186 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&viewer);CHKERRQ(ierr); 187 } 188 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 189 PetscCheckSameComm(snes,1,viewer,2); 190 191 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 192 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 193 if (iascii) { 194 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)snes,viewer,"SNES Object");CHKERRQ(ierr); 195 if (snes->ops->view) { 196 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 197 ierr = (*snes->ops->view)(snes,viewer);CHKERRQ(ierr); 198 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 199 } 200 ierr = PetscViewerASCIIPrintf(viewer," maximum iterations=%D, maximum function evaluations=%D\n",snes->max_its,snes->max_funcs);CHKERRQ(ierr); 201 ierr = PetscViewerASCIIPrintf(viewer," tolerances: relative=%G, absolute=%G, solution=%G\n", 202 snes->rtol,snes->abstol,snes->xtol);CHKERRQ(ierr); 203 ierr = PetscViewerASCIIPrintf(viewer," total number of linear solver iterations=%D\n",snes->linear_its);CHKERRQ(ierr); 204 ierr = PetscViewerASCIIPrintf(viewer," total number of function evaluations=%D\n",snes->nfuncs);CHKERRQ(ierr); 205 if (snes->ksp_ewconv) { 206 kctx = (SNESKSPEW *)snes->kspconvctx; 207 if (kctx) { 208 ierr = PetscViewerASCIIPrintf(viewer," Eisenstat-Walker computation of KSP relative tolerance (version %D)\n",kctx->version);CHKERRQ(ierr); 209 ierr = PetscViewerASCIIPrintf(viewer," rtol_0=%G, rtol_max=%G, threshold=%G\n",kctx->rtol_0,kctx->rtol_max,kctx->threshold);CHKERRQ(ierr); 210 ierr = PetscViewerASCIIPrintf(viewer," gamma=%G, alpha=%G, alpha2=%G\n",kctx->gamma,kctx->alpha,kctx->alpha2);CHKERRQ(ierr); 211 } 212 } 213 if (snes->lagpreconditioner == -1) { 214 ierr = PetscViewerASCIIPrintf(viewer," Preconditioned is never rebuilt\n");CHKERRQ(ierr); 215 } else if (snes->lagpreconditioner > 1) { 216 ierr = PetscViewerASCIIPrintf(viewer," Preconditioned is rebuilt every %D new Jacobians\n",snes->lagpreconditioner);CHKERRQ(ierr); 217 } 218 if (snes->lagjacobian == -1) { 219 ierr = PetscViewerASCIIPrintf(viewer," Jacobian is never rebuilt\n");CHKERRQ(ierr); 220 } else if (snes->lagjacobian > 1) { 221 ierr = PetscViewerASCIIPrintf(viewer," Jacobian is rebuilt every %D SNES iterations\n",snes->lagjacobian);CHKERRQ(ierr); 222 } 223 } else if (isstring) { 224 const char *type; 225 ierr = SNESGetType(snes,&type);CHKERRQ(ierr); 226 ierr = PetscViewerStringSPrintf(viewer," %-3.3s",type);CHKERRQ(ierr); 227 } 228 if (snes->pc && snes->usespc) { 229 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 230 ierr = SNESView(snes->pc, viewer);CHKERRQ(ierr); 231 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 232 } 233 if (snes->usesksp) { 234 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 235 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 236 ierr = KSPView(ksp,viewer);CHKERRQ(ierr); 237 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 238 } 239 ierr = PetscViewerASCIIPrintf(viewer, " line search variant: %s\n",((PetscObject)snes->linesearch)->type_name);CHKERRQ(ierr); 240 PetscFunctionReturn(0); 241 } 242 243 /* 244 We retain a list of functions that also take SNES command 245 line options. These are called at the end SNESSetFromOptions() 246 */ 247 #define MAXSETFROMOPTIONS 5 248 static PetscInt numberofsetfromoptions; 249 static PetscErrorCode (*othersetfromoptions[MAXSETFROMOPTIONS])(SNES); 250 251 #undef __FUNCT__ 252 #define __FUNCT__ "SNESAddOptionsChecker" 253 /*@C 254 SNESAddOptionsChecker - Adds an additional function to check for SNES options. 255 256 Not Collective 257 258 Input Parameter: 259 . snescheck - function that checks for options 260 261 Level: developer 262 263 .seealso: SNESSetFromOptions() 264 @*/ 265 PetscErrorCode SNESAddOptionsChecker(PetscErrorCode (*snescheck)(SNES)) 266 { 267 PetscFunctionBegin; 268 if (numberofsetfromoptions >= MAXSETFROMOPTIONS) { 269 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Too many options checkers, only %D allowed", MAXSETFROMOPTIONS); 270 } 271 othersetfromoptions[numberofsetfromoptions++] = snescheck; 272 PetscFunctionReturn(0); 273 } 274 275 extern PetscErrorCode SNESDefaultMatrixFreeCreate2(SNES,Vec,Mat*); 276 277 #undef __FUNCT__ 278 #define __FUNCT__ "SNESSetUpMatrixFree_Private" 279 static PetscErrorCode SNESSetUpMatrixFree_Private(SNES snes, PetscBool hasOperator, PetscInt version) 280 { 281 Mat J; 282 KSP ksp; 283 PC pc; 284 PetscBool match; 285 PetscErrorCode ierr; 286 287 PetscFunctionBegin; 288 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 289 290 if(!snes->vec_func && (snes->jacobian || snes->jacobian_pre)) { 291 Mat A = snes->jacobian, B = snes->jacobian_pre; 292 ierr = MatGetVecs(A ? A : B, PETSC_NULL,&snes->vec_func);CHKERRQ(ierr); 293 } 294 295 if (version == 1) { 296 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 297 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 298 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 299 } else if (version == 2) { 300 if (!snes->vec_func) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"SNESSetFunction() must be called first"); 301 #if !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 302 ierr = SNESDefaultMatrixFreeCreate2(snes,snes->vec_func,&J);CHKERRQ(ierr); 303 #else 304 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP, "matrix-free operator rutines (version 2)"); 305 #endif 306 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "matrix-free operator rutines, only version 1 and 2"); 307 308 ierr = PetscInfo1(snes,"Setting default matrix-free operator routines (version %D)\n", version);CHKERRQ(ierr); 309 if (hasOperator) { 310 /* This version replaces the user provided Jacobian matrix with a 311 matrix-free version but still employs the user-provided preconditioner matrix. */ 312 ierr = SNESSetJacobian(snes,J,0,0,0);CHKERRQ(ierr); 313 } else { 314 /* This version replaces both the user-provided Jacobian and the user- 315 provided preconditioner matrix with the default matrix free version. */ 316 void *functx; 317 ierr = SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);CHKERRQ(ierr); 318 ierr = SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,functx);CHKERRQ(ierr); 319 /* Force no preconditioner */ 320 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 321 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); 322 ierr = PetscTypeCompare((PetscObject)pc,PCSHELL,&match);CHKERRQ(ierr); 323 if (!match) { 324 ierr = PetscInfo(snes,"Setting default matrix-free preconditioner routines\nThat is no preconditioner is being used\n");CHKERRQ(ierr); 325 ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr); 326 } 327 } 328 ierr = MatDestroy(&J);CHKERRQ(ierr); 329 PetscFunctionReturn(0); 330 } 331 332 #undef __FUNCT__ 333 #define __FUNCT__ "SNESSetUpMatrices" 334 /*@ 335 SNESSetUpMatrices - ensures that matrices are available for SNES, to be called by SNESSetUp_XXX() 336 337 Collective 338 339 Input Arguments: 340 . snes - snes to configure 341 342 Level: developer 343 344 .seealso: SNESSetUp() 345 @*/ 346 PetscErrorCode SNESSetUpMatrices(SNES snes) 347 { 348 PetscErrorCode ierr; 349 DM dm; 350 SNESDM sdm; 351 352 PetscFunctionBegin; 353 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 354 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 355 if (!sdm->computejacobian && snes->dm) { 356 Mat J,B; 357 ierr = DMCreateMatrix(snes->dm,MATAIJ,&B);CHKERRQ(ierr); 358 if (snes->mf_operator) { 359 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 360 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 361 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 362 } else { 363 J = B; 364 ierr = PetscObjectReference((PetscObject)J);CHKERRQ(ierr); 365 } 366 ierr = SNESSetJacobian(snes,J,B,SNESDMComputeJacobian,PETSC_NULL);CHKERRQ(ierr); 367 ierr = MatDestroy(&J);CHKERRQ(ierr); 368 ierr = MatDestroy(&B);CHKERRQ(ierr); 369 } else if (!snes->jacobian && sdm->computejacobian == MatMFFDComputeJacobian) { 370 Mat J; 371 void *functx; 372 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 373 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 374 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 375 ierr = SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);CHKERRQ(ierr); 376 ierr = SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,functx);CHKERRQ(ierr); 377 ierr = MatDestroy(&J);CHKERRQ(ierr); 378 } else if (snes->dm && snes->mf_operator && !snes->jacobian_pre && !snes->jacobian) { 379 Mat J,B; 380 void *functx; 381 ierr = MatCreateSNESMF(snes,&J);CHKERRQ(ierr); 382 ierr = MatMFFDSetOptionsPrefix(J,((PetscObject)snes)->prefix);CHKERRQ(ierr); 383 ierr = MatSetFromOptions(J);CHKERRQ(ierr); 384 ierr = DMCreateMatrix(snes->dm,MATAIJ,&B);CHKERRQ(ierr); 385 ierr = SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);CHKERRQ(ierr); 386 ierr = SNESSetJacobian(snes,J,B,SNESDMComputeJacobian,functx);CHKERRQ(ierr); 387 ierr = MatDestroy(&J);CHKERRQ(ierr); 388 ierr = MatDestroy(&B);CHKERRQ(ierr); 389 } else if (snes->dm && !snes->jacobian_pre) { 390 Mat J,B; 391 J = snes->jacobian; 392 ierr = DMCreateMatrix(snes->dm,MATAIJ,&B);CHKERRQ(ierr); 393 ierr = SNESSetJacobian(snes,J?J:B,B,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 394 ierr = MatDestroy(&B);CHKERRQ(ierr); 395 } 396 PetscFunctionReturn(0); 397 } 398 399 #undef __FUNCT__ 400 #define __FUNCT__ "SNESSetFromOptions" 401 /*@ 402 SNESSetFromOptions - Sets various SNES and KSP parameters from user options. 403 404 Collective on SNES 405 406 Input Parameter: 407 . snes - the SNES context 408 409 Options Database Keys: 410 + -snes_type <type> - ls, tr, ngmres, ncg, richardson, qn, vi, fas 411 . -snes_stol - convergence tolerance in terms of the norm 412 of the change in the solution between steps 413 . -snes_atol <abstol> - absolute tolerance of residual norm 414 . -snes_rtol <rtol> - relative decrease in tolerance norm from initial 415 . -snes_max_it <max_it> - maximum number of iterations 416 . -snes_max_funcs <max_funcs> - maximum number of function evaluations 417 . -snes_max_fail <max_fail> - maximum number of line search failures allowed before stopping, default is none 418 . -snes_max_linear_solve_fail - number of linear solver failures before SNESSolve() stops 419 . -snes_lag_preconditioner <lag> - how often preconditioner is rebuilt (use -1 to never rebuild) 420 . -snes_lag_jacobian <lag> - how often Jacobian is rebuilt (use -1 to never rebuild) 421 . -snes_trtol <trtol> - trust region tolerance 422 . -snes_no_convergence_test - skip convergence test in nonlinear 423 solver; hence iterations will continue until max_it 424 or some other criterion is reached. Saves expense 425 of convergence test 426 . -snes_monitor <optional filename> - prints residual norm at each iteration. if no 427 filename given prints to stdout 428 . -snes_monitor_solution - plots solution at each iteration 429 . -snes_monitor_residual - plots residual (not its norm) at each iteration 430 . -snes_monitor_solution_update - plots update to solution at each iteration 431 . -snes_monitor_draw - plots residual norm at each iteration 432 . -snes_fd - use finite differences to compute Jacobian; very slow, only for testing 433 . -snes_mf_ksp_monitor - if using matrix-free multiply then print h at each KSP iteration 434 - -snes_converged_reason - print the reason for convergence/divergence after each solve 435 436 Options Database for Eisenstat-Walker method: 437 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 438 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 439 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 440 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 441 . -snes_ksp_ew_gamma <gamma> - Sets gamma 442 . -snes_ksp_ew_alpha <alpha> - Sets alpha 443 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 444 - -snes_ksp_ew_threshold <threshold> - Sets threshold 445 446 Notes: 447 To see all options, run your program with the -help option or consult 448 the <A href="../../docs/manual.pdf#nameddest=ch_snes">SNES chapter of the users manual</A>. 449 450 Level: beginner 451 452 .keywords: SNES, nonlinear, set, options, database 453 454 .seealso: SNESSetOptionsPrefix() 455 @*/ 456 PetscErrorCode SNESSetFromOptions(SNES snes) 457 { 458 PetscBool flg,mf,mf_operator,pcset; 459 PetscInt i,indx,lag,mf_version,grids; 460 MatStructure matflag; 461 const char *deft = SNESLS; 462 const char *convtests[] = {"default","skip"}; 463 SNESKSPEW *kctx = NULL; 464 char type[256], monfilename[PETSC_MAX_PATH_LEN]; 465 const char *optionsprefix; 466 PetscViewer monviewer; 467 PetscErrorCode ierr; 468 469 PetscFunctionBegin; 470 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 471 472 if (!SNESRegisterAllCalled) {ierr = SNESRegisterAll(PETSC_NULL);CHKERRQ(ierr);} 473 ierr = PetscObjectOptionsBegin((PetscObject)snes);CHKERRQ(ierr); 474 if (((PetscObject)snes)->type_name) { deft = ((PetscObject)snes)->type_name; } 475 ierr = PetscOptionsList("-snes_type","Nonlinear solver method","SNESSetType",SNESList,deft,type,256,&flg);CHKERRQ(ierr); 476 if (flg) { 477 ierr = SNESSetType(snes,type);CHKERRQ(ierr); 478 } else if (!((PetscObject)snes)->type_name) { 479 ierr = SNESSetType(snes,deft);CHKERRQ(ierr); 480 } 481 /* not used here, but called so will go into help messaage */ 482 ierr = PetscOptionsName("-snes_view","Print detailed information on solver used","SNESView",0);CHKERRQ(ierr); 483 484 ierr = PetscOptionsReal("-snes_stol","Stop if step length less than","SNESSetTolerances",snes->xtol,&snes->xtol,0);CHKERRQ(ierr); 485 ierr = PetscOptionsReal("-snes_atol","Stop if function norm less than","SNESSetTolerances",snes->abstol,&snes->abstol,0);CHKERRQ(ierr); 486 487 ierr = PetscOptionsReal("-snes_rtol","Stop if decrease in function norm less than","SNESSetTolerances",snes->rtol,&snes->rtol,0);CHKERRQ(ierr); 488 ierr = PetscOptionsInt("-snes_max_it","Maximum iterations","SNESSetTolerances",snes->max_its,&snes->max_its,PETSC_NULL);CHKERRQ(ierr); 489 ierr = PetscOptionsInt("-snes_max_funcs","Maximum function evaluations","SNESSetTolerances",snes->max_funcs,&snes->max_funcs,PETSC_NULL);CHKERRQ(ierr); 490 ierr = PetscOptionsInt("-snes_max_fail","Maximum failures","SNESSetTolerances",snes->maxFailures,&snes->maxFailures,PETSC_NULL);CHKERRQ(ierr); 491 ierr = PetscOptionsInt("-snes_max_linear_solve_fail","Maximum failures in linear solves allowed","SNESSetMaxLinearSolveFailures",snes->maxLinearSolveFailures,&snes->maxLinearSolveFailures,PETSC_NULL);CHKERRQ(ierr); 492 ierr = PetscOptionsBool("-snes_error_if_not_converged","Generate error if solver does not converge","SNESSetErrorIfNotConverged",snes->errorifnotconverged,&snes->errorifnotconverged,PETSC_NULL);CHKERRQ(ierr); 493 494 ierr = PetscOptionsInt("-snes_lag_preconditioner","How often to rebuild preconditioner","SNESSetLagPreconditioner",snes->lagpreconditioner,&lag,&flg);CHKERRQ(ierr); 495 if (flg) { 496 ierr = SNESSetLagPreconditioner(snes,lag);CHKERRQ(ierr); 497 } 498 ierr = PetscOptionsInt("-snes_lag_jacobian","How often to rebuild Jacobian","SNESSetLagJacobian",snes->lagjacobian,&lag,&flg);CHKERRQ(ierr); 499 if (flg) { 500 ierr = SNESSetLagJacobian(snes,lag);CHKERRQ(ierr); 501 } 502 ierr = PetscOptionsInt("-snes_grid_sequence","Use grid sequencing to generate initial guess","SNESSetGridSequence",snes->gridsequence,&grids,&flg);CHKERRQ(ierr); 503 if (flg) { 504 ierr = SNESSetGridSequence(snes,grids);CHKERRQ(ierr); 505 } 506 507 ierr = PetscOptionsEList("-snes_convergence_test","Convergence test","SNESSetConvergenceTest",convtests,2,"default",&indx,&flg);CHKERRQ(ierr); 508 if (flg) { 509 switch (indx) { 510 case 0: ierr = SNESSetConvergenceTest(snes,SNESDefaultConverged,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); break; 511 case 1: ierr = SNESSetConvergenceTest(snes,SNESSkipConverged,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); break; 512 } 513 } 514 515 ierr = PetscOptionsBool("-snes_converged_reason","Print reason for converged or diverged","SNESSolve",snes->printreason,&snes->printreason,PETSC_NULL);CHKERRQ(ierr); 516 517 kctx = (SNESKSPEW *)snes->kspconvctx; 518 519 ierr = PetscOptionsBool("-snes_ksp_ew","Use Eisentat-Walker linear system convergence test","SNESKSPSetUseEW",snes->ksp_ewconv,&snes->ksp_ewconv,PETSC_NULL);CHKERRQ(ierr); 520 521 ierr = PetscOptionsInt("-snes_ksp_ew_version","Version 1, 2 or 3","SNESKSPSetParametersEW",kctx->version,&kctx->version,0);CHKERRQ(ierr); 522 ierr = PetscOptionsReal("-snes_ksp_ew_rtol0","0 <= rtol0 < 1","SNESKSPSetParametersEW",kctx->rtol_0,&kctx->rtol_0,0);CHKERRQ(ierr); 523 ierr = PetscOptionsReal("-snes_ksp_ew_rtolmax","0 <= rtolmax < 1","SNESKSPSetParametersEW",kctx->rtol_max,&kctx->rtol_max,0);CHKERRQ(ierr); 524 ierr = PetscOptionsReal("-snes_ksp_ew_gamma","0 <= gamma <= 1","SNESKSPSetParametersEW",kctx->gamma,&kctx->gamma,0);CHKERRQ(ierr); 525 ierr = PetscOptionsReal("-snes_ksp_ew_alpha","1 < alpha <= 2","SNESKSPSetParametersEW",kctx->alpha,&kctx->alpha,0);CHKERRQ(ierr); 526 ierr = PetscOptionsReal("-snes_ksp_ew_alpha2","alpha2","SNESKSPSetParametersEW",kctx->alpha2,&kctx->alpha2,0);CHKERRQ(ierr); 527 ierr = PetscOptionsReal("-snes_ksp_ew_threshold","0 < threshold < 1","SNESKSPSetParametersEW",kctx->threshold,&kctx->threshold,0);CHKERRQ(ierr); 528 529 flg = PETSC_FALSE; 530 ierr = PetscOptionsBool("-snes_monitor_cancel","Remove all monitors","SNESMonitorCancel",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 531 if (flg) {ierr = SNESMonitorCancel(snes);CHKERRQ(ierr);} 532 533 ierr = PetscOptionsString("-snes_monitor","Monitor norm of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 534 if (flg) { 535 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);CHKERRQ(ierr); 536 ierr = SNESMonitorSet(snes,SNESMonitorDefault,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 537 } 538 539 ierr = PetscOptionsString("-snes_monitor_range","Monitor range of elements of function","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 540 if (flg) { 541 ierr = SNESMonitorSet(snes,SNESMonitorRange,0,0);CHKERRQ(ierr); 542 } 543 544 ierr = PetscOptionsString("-snes_ratiomonitor","Monitor ratios of norms of function","SNESMonitorSetRatio","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 545 if (flg) { 546 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);CHKERRQ(ierr); 547 ierr = SNESMonitorSetRatio(snes,monviewer);CHKERRQ(ierr); 548 } 549 550 ierr = PetscOptionsString("-snes_monitor_short","Monitor norm of function (fewer digits)","SNESMonitorSet","stdout",monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 551 if (flg) { 552 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,monfilename,&monviewer);CHKERRQ(ierr); 553 ierr = SNESMonitorSet(snes,SNESMonitorDefaultShort,monviewer,(PetscErrorCode (*)(void**))PetscViewerDestroy);CHKERRQ(ierr); 554 } 555 556 ierr = PetscOptionsString("-snes_monitor_python","Use Python function","SNESMonitorSet",0,monfilename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 557 if (flg) {ierr = PetscPythonMonitorSet((PetscObject)snes,monfilename);CHKERRQ(ierr);} 558 559 flg = PETSC_FALSE; 560 ierr = PetscOptionsBool("-snes_monitor_solution","Plot solution at each iteration","SNESMonitorSolution",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 561 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorSolution,0,0);CHKERRQ(ierr);} 562 flg = PETSC_FALSE; 563 ierr = PetscOptionsBool("-snes_monitor_solution_update","Plot correction at each iteration","SNESMonitorSolutionUpdate",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 564 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorSolutionUpdate,0,0);CHKERRQ(ierr);} 565 flg = PETSC_FALSE; 566 ierr = PetscOptionsBool("-snes_monitor_residual","Plot residual at each iteration","SNESMonitorResidual",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 567 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorResidual,0,0);CHKERRQ(ierr);} 568 flg = PETSC_FALSE; 569 ierr = PetscOptionsBool("-snes_monitor_draw","Plot function norm at each iteration","SNESMonitorLG",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 570 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorLG,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);} 571 flg = PETSC_FALSE; 572 ierr = PetscOptionsBool("-snes_monitor_range_draw","Plot function range at each iteration","SNESMonitorLG",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 573 if (flg) {ierr = SNESMonitorSet(snes,SNESMonitorLGRange,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);} 574 575 flg = PETSC_FALSE; 576 ierr = PetscOptionsBool("-snes_fd","Use finite differences (slow) to compute Jacobian","SNESDefaultComputeJacobian",flg,&flg,PETSC_NULL);CHKERRQ(ierr); 577 if (flg) { 578 void *functx; 579 ierr = SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);CHKERRQ(ierr); 580 ierr = SNESSetJacobian(snes,snes->jacobian,snes->jacobian_pre,SNESDefaultComputeJacobian,functx);CHKERRQ(ierr); 581 ierr = PetscInfo(snes,"Setting default finite difference Jacobian matrix\n");CHKERRQ(ierr); 582 } 583 584 mf = mf_operator = PETSC_FALSE; 585 flg = PETSC_FALSE; 586 ierr = PetscOptionsBool("-snes_mf_operator","Use a Matrix-Free Jacobian with user-provided preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&mf_operator,&flg);CHKERRQ(ierr); 587 if (flg && mf_operator) { 588 snes->mf_operator = PETSC_TRUE; 589 mf = PETSC_TRUE; 590 } 591 flg = PETSC_FALSE; 592 ierr = PetscOptionsBool("-snes_mf","Use a Matrix-Free Jacobian with no preconditioner matrix","MatCreateSNESMF",PETSC_FALSE,&mf,&flg);CHKERRQ(ierr); 593 if (!flg && mf_operator) mf = PETSC_TRUE; 594 mf_version = 1; 595 ierr = PetscOptionsInt("-snes_mf_version","Matrix-Free routines version 1 or 2","None",mf_version,&mf_version,0);CHKERRQ(ierr); 596 597 598 /* GS Options */ 599 ierr = PetscOptionsInt("-snes_gs_sweeps","Number of sweeps of GS to apply","SNESComputeGS",snes->gssweeps,&snes->gssweeps,PETSC_NULL);CHKERRQ(ierr); 600 601 for(i = 0; i < numberofsetfromoptions; i++) { 602 ierr = (*othersetfromoptions[i])(snes);CHKERRQ(ierr); 603 } 604 605 if (snes->ops->setfromoptions) { 606 ierr = (*snes->ops->setfromoptions)(snes);CHKERRQ(ierr); 607 } 608 609 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 610 ierr = PetscObjectProcessOptionsHandlers((PetscObject)snes);CHKERRQ(ierr); 611 ierr = PetscOptionsEnd();CHKERRQ(ierr); 612 613 if (mf) { ierr = SNESSetUpMatrixFree_Private(snes, mf_operator, mf_version);CHKERRQ(ierr); } 614 615 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 616 ierr = KSPGetOperators(snes->ksp,PETSC_NULL,PETSC_NULL,&matflag); 617 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,matflag);CHKERRQ(ierr); 618 ierr = KSPSetFromOptions(snes->ksp);CHKERRQ(ierr); 619 620 if (!snes->linesearch) { 621 ierr = SNESGetPetscLineSearch(snes, &snes->linesearch);CHKERRQ(ierr); 622 } 623 ierr = PetscLineSearchSetFromOptions(snes->linesearch);CHKERRQ(ierr); 624 625 /* if someone has set the SNES PC type, create it. */ 626 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 627 ierr = PetscOptionsHasName(optionsprefix, "-npc_snes_type", &pcset);CHKERRQ(ierr); 628 if (pcset && (!snes->pc)) { 629 ierr = SNESGetPC(snes, &snes->pc);CHKERRQ(ierr); 630 } 631 if (snes->pc) { 632 ierr = SNESSetOptionsPrefix(snes->pc, optionsprefix);CHKERRQ(ierr); 633 ierr = SNESAppendOptionsPrefix(snes->pc, "npc_");CHKERRQ(ierr); 634 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 635 ierr = SNESSetFromOptions(snes->pc);CHKERRQ(ierr); 636 } 637 PetscFunctionReturn(0); 638 } 639 640 #undef __FUNCT__ 641 #define __FUNCT__ "SNESSetComputeApplicationContext" 642 /*@ 643 SNESSetComputeApplicationContext - Sets an optional function to compute a user-defined context for 644 the nonlinear solvers. 645 646 Logically Collective on SNES 647 648 Input Parameters: 649 + snes - the SNES context 650 . compute - function to compute the context 651 - destroy - function to destroy the context 652 653 Level: intermediate 654 655 .keywords: SNES, nonlinear, set, application, context 656 657 .seealso: SNESGetApplicationContext(), SNESSetComputeApplicationContext(), SNESGetApplicationContext() 658 @*/ 659 PetscErrorCode SNESSetComputeApplicationContext(SNES snes,PetscErrorCode (*compute)(SNES,void**),PetscErrorCode (*destroy)(void**)) 660 { 661 PetscFunctionBegin; 662 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 663 snes->ops->usercompute = compute; 664 snes->ops->userdestroy = destroy; 665 PetscFunctionReturn(0); 666 } 667 668 #undef __FUNCT__ 669 #define __FUNCT__ "SNESSetApplicationContext" 670 /*@ 671 SNESSetApplicationContext - Sets the optional user-defined context for 672 the nonlinear solvers. 673 674 Logically Collective on SNES 675 676 Input Parameters: 677 + snes - the SNES context 678 - usrP - optional user context 679 680 Level: intermediate 681 682 .keywords: SNES, nonlinear, set, application, context 683 684 .seealso: SNESGetApplicationContext(), SNESSetApplicationContext() 685 @*/ 686 PetscErrorCode SNESSetApplicationContext(SNES snes,void *usrP) 687 { 688 PetscErrorCode ierr; 689 KSP ksp; 690 691 PetscFunctionBegin; 692 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 693 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 694 ierr = KSPSetApplicationContext(ksp,usrP);CHKERRQ(ierr); 695 snes->user = usrP; 696 PetscFunctionReturn(0); 697 } 698 699 #undef __FUNCT__ 700 #define __FUNCT__ "SNESGetApplicationContext" 701 /*@ 702 SNESGetApplicationContext - Gets the user-defined context for the 703 nonlinear solvers. 704 705 Not Collective 706 707 Input Parameter: 708 . snes - SNES context 709 710 Output Parameter: 711 . usrP - user context 712 713 Level: intermediate 714 715 .keywords: SNES, nonlinear, get, application, context 716 717 .seealso: SNESSetApplicationContext() 718 @*/ 719 PetscErrorCode SNESGetApplicationContext(SNES snes,void *usrP) 720 { 721 PetscFunctionBegin; 722 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 723 *(void**)usrP = snes->user; 724 PetscFunctionReturn(0); 725 } 726 727 #undef __FUNCT__ 728 #define __FUNCT__ "SNESGetIterationNumber" 729 /*@ 730 SNESGetIterationNumber - Gets the number of nonlinear iterations completed 731 at this time. 732 733 Not Collective 734 735 Input Parameter: 736 . snes - SNES context 737 738 Output Parameter: 739 . iter - iteration number 740 741 Notes: 742 For example, during the computation of iteration 2 this would return 1. 743 744 This is useful for using lagged Jacobians (where one does not recompute the 745 Jacobian at each SNES iteration). For example, the code 746 .vb 747 ierr = SNESGetIterationNumber(snes,&it); 748 if (!(it % 2)) { 749 [compute Jacobian here] 750 } 751 .ve 752 can be used in your ComputeJacobian() function to cause the Jacobian to be 753 recomputed every second SNES iteration. 754 755 Level: intermediate 756 757 .keywords: SNES, nonlinear, get, iteration, number, 758 759 .seealso: SNESGetFunctionNorm(), SNESGetLinearSolveIterations() 760 @*/ 761 PetscErrorCode SNESGetIterationNumber(SNES snes,PetscInt* iter) 762 { 763 PetscFunctionBegin; 764 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 765 PetscValidIntPointer(iter,2); 766 *iter = snes->iter; 767 PetscFunctionReturn(0); 768 } 769 770 #undef __FUNCT__ 771 #define __FUNCT__ "SNESSetIterationNumber" 772 /*@ 773 SNESSetIterationNumber - Sets the current iteration number. 774 775 Not Collective 776 777 Input Parameter: 778 . snes - SNES context 779 . iter - iteration number 780 781 Level: developer 782 783 .keywords: SNES, nonlinear, set, iteration, number, 784 785 .seealso: SNESGetFunctionNorm(), SNESGetLinearSolveIterations() 786 @*/ 787 PetscErrorCode SNESSetIterationNumber(SNES snes,PetscInt iter) 788 { 789 PetscErrorCode ierr; 790 791 PetscFunctionBegin; 792 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 793 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 794 snes->iter = iter; 795 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 796 PetscFunctionReturn(0); 797 } 798 799 #undef __FUNCT__ 800 #define __FUNCT__ "SNESGetFunctionNorm" 801 /*@ 802 SNESGetFunctionNorm - Gets the norm of the current function that was set 803 with SNESSSetFunction(). 804 805 Collective on SNES 806 807 Input Parameter: 808 . snes - SNES context 809 810 Output Parameter: 811 . fnorm - 2-norm of function 812 813 Level: intermediate 814 815 .keywords: SNES, nonlinear, get, function, norm 816 817 .seealso: SNESGetFunction(), SNESGetIterationNumber(), SNESGetLinearSolveIterations() 818 @*/ 819 PetscErrorCode SNESGetFunctionNorm(SNES snes,PetscReal *fnorm) 820 { 821 PetscFunctionBegin; 822 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 823 PetscValidScalarPointer(fnorm,2); 824 *fnorm = snes->norm; 825 PetscFunctionReturn(0); 826 } 827 828 829 #undef __FUNCT__ 830 #define __FUNCT__ "SNESSetFunctionNorm" 831 /*@ 832 SNESSetFunctionNorm - Sets the 2-norm of the current function computed using VecNorm(). 833 834 Collective on SNES 835 836 Input Parameter: 837 . snes - SNES context 838 . fnorm - 2-norm of function 839 840 Level: developer 841 842 .keywords: SNES, nonlinear, set, function, norm 843 844 .seealso: SNESSetFunction(), SNESSetIterationNumber(), VecNorm(). 845 @*/ 846 PetscErrorCode SNESSetFunctionNorm(SNES snes,PetscReal fnorm) 847 { 848 849 PetscErrorCode ierr; 850 851 PetscFunctionBegin; 852 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 853 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 854 snes->norm = fnorm; 855 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 856 PetscFunctionReturn(0); 857 } 858 859 #undef __FUNCT__ 860 #define __FUNCT__ "SNESGetNonlinearStepFailures" 861 /*@ 862 SNESGetNonlinearStepFailures - Gets the number of unsuccessful steps 863 attempted by the nonlinear solver. 864 865 Not Collective 866 867 Input Parameter: 868 . snes - SNES context 869 870 Output Parameter: 871 . nfails - number of unsuccessful steps attempted 872 873 Notes: 874 This counter is reset to zero for each successive call to SNESSolve(). 875 876 Level: intermediate 877 878 .keywords: SNES, nonlinear, get, number, unsuccessful, steps 879 880 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(), 881 SNESSetMaxNonlinearStepFailures(), SNESGetMaxNonlinearStepFailures() 882 @*/ 883 PetscErrorCode SNESGetNonlinearStepFailures(SNES snes,PetscInt* nfails) 884 { 885 PetscFunctionBegin; 886 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 887 PetscValidIntPointer(nfails,2); 888 *nfails = snes->numFailures; 889 PetscFunctionReturn(0); 890 } 891 892 #undef __FUNCT__ 893 #define __FUNCT__ "SNESSetMaxNonlinearStepFailures" 894 /*@ 895 SNESSetMaxNonlinearStepFailures - Sets the maximum number of unsuccessful steps 896 attempted by the nonlinear solver before it gives up. 897 898 Not Collective 899 900 Input Parameters: 901 + snes - SNES context 902 - maxFails - maximum of unsuccessful steps 903 904 Level: intermediate 905 906 .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps 907 908 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(), 909 SNESGetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures() 910 @*/ 911 PetscErrorCode SNESSetMaxNonlinearStepFailures(SNES snes, PetscInt maxFails) 912 { 913 PetscFunctionBegin; 914 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 915 snes->maxFailures = maxFails; 916 PetscFunctionReturn(0); 917 } 918 919 #undef __FUNCT__ 920 #define __FUNCT__ "SNESGetMaxNonlinearStepFailures" 921 /*@ 922 SNESGetMaxNonlinearStepFailures - Gets the maximum number of unsuccessful steps 923 attempted by the nonlinear solver before it gives up. 924 925 Not Collective 926 927 Input Parameter: 928 . snes - SNES context 929 930 Output Parameter: 931 . maxFails - maximum of unsuccessful steps 932 933 Level: intermediate 934 935 .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps 936 937 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures(), 938 SNESSetMaxNonlinearStepFailures(), SNESGetNonlinearStepFailures() 939 940 @*/ 941 PetscErrorCode SNESGetMaxNonlinearStepFailures(SNES snes, PetscInt *maxFails) 942 { 943 PetscFunctionBegin; 944 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 945 PetscValidIntPointer(maxFails,2); 946 *maxFails = snes->maxFailures; 947 PetscFunctionReturn(0); 948 } 949 950 #undef __FUNCT__ 951 #define __FUNCT__ "SNESGetNumberFunctionEvals" 952 /*@ 953 SNESGetNumberFunctionEvals - Gets the number of user provided function evaluations 954 done by SNES. 955 956 Not Collective 957 958 Input Parameter: 959 . snes - SNES context 960 961 Output Parameter: 962 . nfuncs - number of evaluations 963 964 Level: intermediate 965 966 .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps 967 968 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), SNESGetLinearSolveFailures() 969 @*/ 970 PetscErrorCode SNESGetNumberFunctionEvals(SNES snes, PetscInt *nfuncs) 971 { 972 PetscFunctionBegin; 973 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 974 PetscValidIntPointer(nfuncs,2); 975 *nfuncs = snes->nfuncs; 976 PetscFunctionReturn(0); 977 } 978 979 #undef __FUNCT__ 980 #define __FUNCT__ "SNESGetLinearSolveFailures" 981 /*@ 982 SNESGetLinearSolveFailures - Gets the number of failed (non-converged) 983 linear solvers. 984 985 Not Collective 986 987 Input Parameter: 988 . snes - SNES context 989 990 Output Parameter: 991 . nfails - number of failed solves 992 993 Notes: 994 This counter is reset to zero for each successive call to SNESSolve(). 995 996 Level: intermediate 997 998 .keywords: SNES, nonlinear, get, number, unsuccessful, steps 999 1000 .seealso: SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures() 1001 @*/ 1002 PetscErrorCode SNESGetLinearSolveFailures(SNES snes,PetscInt* nfails) 1003 { 1004 PetscFunctionBegin; 1005 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1006 PetscValidIntPointer(nfails,2); 1007 *nfails = snes->numLinearSolveFailures; 1008 PetscFunctionReturn(0); 1009 } 1010 1011 #undef __FUNCT__ 1012 #define __FUNCT__ "SNESSetMaxLinearSolveFailures" 1013 /*@ 1014 SNESSetMaxLinearSolveFailures - the number of failed linear solve attempts 1015 allowed before SNES returns with a diverged reason of SNES_DIVERGED_LINEAR_SOLVE 1016 1017 Logically Collective on SNES 1018 1019 Input Parameters: 1020 + snes - SNES context 1021 - maxFails - maximum allowed linear solve failures 1022 1023 Level: intermediate 1024 1025 Notes: By default this is 0; that is SNES returns on the first failed linear solve 1026 1027 .keywords: SNES, nonlinear, set, maximum, unsuccessful, steps 1028 1029 .seealso: SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures(), SNESGetLinearSolveIterations() 1030 @*/ 1031 PetscErrorCode SNESSetMaxLinearSolveFailures(SNES snes, PetscInt maxFails) 1032 { 1033 PetscFunctionBegin; 1034 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1035 PetscValidLogicalCollectiveInt(snes,maxFails,2); 1036 snes->maxLinearSolveFailures = maxFails; 1037 PetscFunctionReturn(0); 1038 } 1039 1040 #undef __FUNCT__ 1041 #define __FUNCT__ "SNESGetMaxLinearSolveFailures" 1042 /*@ 1043 SNESGetMaxLinearSolveFailures - gets the maximum number of linear solve failures that 1044 are allowed before SNES terminates 1045 1046 Not Collective 1047 1048 Input Parameter: 1049 . snes - SNES context 1050 1051 Output Parameter: 1052 . maxFails - maximum of unsuccessful solves allowed 1053 1054 Level: intermediate 1055 1056 Notes: By default this is 1; that is SNES returns on the first failed linear solve 1057 1058 .keywords: SNES, nonlinear, get, maximum, unsuccessful, steps 1059 1060 .seealso: SNESGetLinearSolveFailures(), SNESGetLinearSolveIterations(), SNESSetMaxLinearSolveFailures(), 1061 @*/ 1062 PetscErrorCode SNESGetMaxLinearSolveFailures(SNES snes, PetscInt *maxFails) 1063 { 1064 PetscFunctionBegin; 1065 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1066 PetscValidIntPointer(maxFails,2); 1067 *maxFails = snes->maxLinearSolveFailures; 1068 PetscFunctionReturn(0); 1069 } 1070 1071 #undef __FUNCT__ 1072 #define __FUNCT__ "SNESGetLinearSolveIterations" 1073 /*@ 1074 SNESGetLinearSolveIterations - Gets the total number of linear iterations 1075 used by the nonlinear solver. 1076 1077 Not Collective 1078 1079 Input Parameter: 1080 . snes - SNES context 1081 1082 Output Parameter: 1083 . lits - number of linear iterations 1084 1085 Notes: 1086 This counter is reset to zero for each successive call to SNESSolve(). 1087 1088 Level: intermediate 1089 1090 .keywords: SNES, nonlinear, get, number, linear, iterations 1091 1092 .seealso: SNESGetIterationNumber(), SNESGetFunctionNorm(), SNESGetLinearSolveFailures(), SNESGetMaxLinearSolveFailures() 1093 @*/ 1094 PetscErrorCode SNESGetLinearSolveIterations(SNES snes,PetscInt* lits) 1095 { 1096 PetscFunctionBegin; 1097 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1098 PetscValidIntPointer(lits,2); 1099 *lits = snes->linear_its; 1100 PetscFunctionReturn(0); 1101 } 1102 1103 #undef __FUNCT__ 1104 #define __FUNCT__ "SNESGetKSP" 1105 /*@ 1106 SNESGetKSP - Returns the KSP context for a SNES solver. 1107 1108 Not Collective, but if SNES object is parallel, then KSP object is parallel 1109 1110 Input Parameter: 1111 . snes - the SNES context 1112 1113 Output Parameter: 1114 . ksp - the KSP context 1115 1116 Notes: 1117 The user can then directly manipulate the KSP context to set various 1118 options, etc. Likewise, the user can then extract and manipulate the 1119 PC contexts as well. 1120 1121 Level: beginner 1122 1123 .keywords: SNES, nonlinear, get, KSP, context 1124 1125 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 1126 @*/ 1127 PetscErrorCode SNESGetKSP(SNES snes,KSP *ksp) 1128 { 1129 PetscErrorCode ierr; 1130 1131 PetscFunctionBegin; 1132 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1133 PetscValidPointer(ksp,2); 1134 1135 if (!snes->ksp) { 1136 ierr = KSPCreate(((PetscObject)snes)->comm,&snes->ksp);CHKERRQ(ierr); 1137 ierr = PetscObjectIncrementTabLevel((PetscObject)snes->ksp,(PetscObject)snes,1);CHKERRQ(ierr); 1138 ierr = PetscLogObjectParent(snes,snes->ksp);CHKERRQ(ierr); 1139 } 1140 *ksp = snes->ksp; 1141 PetscFunctionReturn(0); 1142 } 1143 1144 #undef __FUNCT__ 1145 #define __FUNCT__ "SNESSetKSP" 1146 /*@ 1147 SNESSetKSP - Sets a KSP context for the SNES object to use 1148 1149 Not Collective, but the SNES and KSP objects must live on the same MPI_Comm 1150 1151 Input Parameters: 1152 + snes - the SNES context 1153 - ksp - the KSP context 1154 1155 Notes: 1156 The SNES object already has its KSP object, you can obtain with SNESGetKSP() 1157 so this routine is rarely needed. 1158 1159 The KSP object that is already in the SNES object has its reference count 1160 decreased by one. 1161 1162 Level: developer 1163 1164 .keywords: SNES, nonlinear, get, KSP, context 1165 1166 .seealso: KSPGetPC(), SNESCreate(), KSPCreate(), SNESSetKSP() 1167 @*/ 1168 PetscErrorCode SNESSetKSP(SNES snes,KSP ksp) 1169 { 1170 PetscErrorCode ierr; 1171 1172 PetscFunctionBegin; 1173 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1174 PetscValidHeaderSpecific(ksp,KSP_CLASSID,2); 1175 PetscCheckSameComm(snes,1,ksp,2); 1176 ierr = PetscObjectReference((PetscObject)ksp);CHKERRQ(ierr); 1177 if (snes->ksp) {ierr = PetscObjectDereference((PetscObject)snes->ksp);CHKERRQ(ierr);} 1178 snes->ksp = ksp; 1179 PetscFunctionReturn(0); 1180 } 1181 1182 #if 0 1183 #undef __FUNCT__ 1184 #define __FUNCT__ "SNESPublish_Petsc" 1185 static PetscErrorCode SNESPublish_Petsc(PetscObject obj) 1186 { 1187 PetscFunctionBegin; 1188 PetscFunctionReturn(0); 1189 } 1190 #endif 1191 1192 /* -----------------------------------------------------------*/ 1193 #undef __FUNCT__ 1194 #define __FUNCT__ "SNESCreate" 1195 /*@ 1196 SNESCreate - Creates a nonlinear solver context. 1197 1198 Collective on MPI_Comm 1199 1200 Input Parameters: 1201 . comm - MPI communicator 1202 1203 Output Parameter: 1204 . outsnes - the new SNES context 1205 1206 Options Database Keys: 1207 + -snes_mf - Activates default matrix-free Jacobian-vector products, 1208 and no preconditioning matrix 1209 . -snes_mf_operator - Activates default matrix-free Jacobian-vector 1210 products, and a user-provided preconditioning matrix 1211 as set by SNESSetJacobian() 1212 - -snes_fd - Uses (slow!) finite differences to compute Jacobian 1213 1214 Level: beginner 1215 1216 .keywords: SNES, nonlinear, create, context 1217 1218 .seealso: SNESSolve(), SNESDestroy(), SNES, SNESSetLagPreconditioner() 1219 1220 @*/ 1221 PetscErrorCode SNESCreate(MPI_Comm comm,SNES *outsnes) 1222 { 1223 PetscErrorCode ierr; 1224 SNES snes; 1225 SNESKSPEW *kctx; 1226 1227 PetscFunctionBegin; 1228 PetscValidPointer(outsnes,2); 1229 *outsnes = PETSC_NULL; 1230 #ifndef PETSC_USE_DYNAMIC_LIBRARIES 1231 ierr = SNESInitializePackage(PETSC_NULL);CHKERRQ(ierr); 1232 #endif 1233 1234 ierr = PetscHeaderCreate(snes,_p_SNES,struct _SNESOps,SNES_CLASSID,0,"SNES","Nonlinear solver","SNES",comm,SNESDestroy,SNESView);CHKERRQ(ierr); 1235 1236 snes->ops->converged = SNESDefaultConverged; 1237 snes->usesksp = PETSC_TRUE; 1238 snes->max_its = 50; 1239 snes->max_funcs = 10000; 1240 snes->norm = 0.0; 1241 snes->rtol = 1.e-8; 1242 snes->ttol = 0.0; 1243 snes->abstol = 1.e-50; 1244 snes->xtol = 1.e-8; 1245 snes->deltatol = 1.e-12; 1246 snes->nfuncs = 0; 1247 snes->numFailures = 0; 1248 snes->maxFailures = 1; 1249 snes->linear_its = 0; 1250 snes->lagjacobian = 1; 1251 snes->lagpreconditioner = 1; 1252 snes->numbermonitors = 0; 1253 snes->data = 0; 1254 snes->setupcalled = PETSC_FALSE; 1255 snes->ksp_ewconv = PETSC_FALSE; 1256 snes->nwork = 0; 1257 snes->work = 0; 1258 snes->nvwork = 0; 1259 snes->vwork = 0; 1260 snes->conv_hist_len = 0; 1261 snes->conv_hist_max = 0; 1262 snes->conv_hist = PETSC_NULL; 1263 snes->conv_hist_its = PETSC_NULL; 1264 snes->conv_hist_reset = PETSC_TRUE; 1265 snes->reason = SNES_CONVERGED_ITERATING; 1266 snes->gssweeps = 1; 1267 1268 snes->ops->linesearch = PETSC_NULL; 1269 snes->precheck = PETSC_NULL; 1270 snes->ops->precheckstep = PETSC_NULL; 1271 snes->postcheck = PETSC_NULL; 1272 snes->ops->postcheckstep= PETSC_NULL; 1273 1274 snes->numLinearSolveFailures = 0; 1275 snes->maxLinearSolveFailures = 1; 1276 1277 /* Create context to compute Eisenstat-Walker relative tolerance for KSP */ 1278 ierr = PetscNewLog(snes,SNESKSPEW,&kctx);CHKERRQ(ierr); 1279 snes->kspconvctx = (void*)kctx; 1280 kctx->version = 2; 1281 kctx->rtol_0 = .3; /* Eisenstat and Walker suggest rtol_0=.5, but 1282 this was too large for some test cases */ 1283 kctx->rtol_last = 0.0; 1284 kctx->rtol_max = .9; 1285 kctx->gamma = 1.0; 1286 kctx->alpha = .5*(1.0 + PetscSqrtReal(5.0)); 1287 kctx->alpha2 = kctx->alpha; 1288 kctx->threshold = .1; 1289 kctx->lresid_last = 0.0; 1290 kctx->norm_last = 0.0; 1291 1292 *outsnes = snes; 1293 PetscFunctionReturn(0); 1294 } 1295 1296 #undef __FUNCT__ 1297 #define __FUNCT__ "SNESSetFunction" 1298 /*@C 1299 SNESSetFunction - Sets the function evaluation routine and function 1300 vector for use by the SNES routines in solving systems of nonlinear 1301 equations. 1302 1303 Logically Collective on SNES 1304 1305 Input Parameters: 1306 + snes - the SNES context 1307 . r - vector to store function value 1308 . func - function evaluation routine 1309 - ctx - [optional] user-defined context for private data for the 1310 function evaluation routine (may be PETSC_NULL) 1311 1312 Calling sequence of func: 1313 $ func (SNES snes,Vec x,Vec f,void *ctx); 1314 1315 + snes - the SNES context 1316 . x - state at which to evaluate residual 1317 . f - vector to put residual 1318 - ctx - optional user-defined function context 1319 1320 Notes: 1321 The Newton-like methods typically solve linear systems of the form 1322 $ f'(x) x = -f(x), 1323 where f'(x) denotes the Jacobian matrix and f(x) is the function. 1324 1325 Level: beginner 1326 1327 .keywords: SNES, nonlinear, set, function 1328 1329 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetPicard() 1330 @*/ 1331 PetscErrorCode SNESSetFunction(SNES snes,Vec r,PetscErrorCode (*func)(SNES,Vec,Vec,void*),void *ctx) 1332 { 1333 PetscErrorCode ierr; 1334 DM dm; 1335 1336 PetscFunctionBegin; 1337 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1338 if (r) { 1339 PetscValidHeaderSpecific(r,VEC_CLASSID,2); 1340 PetscCheckSameComm(snes,1,r,2); 1341 ierr = PetscObjectReference((PetscObject)r);CHKERRQ(ierr); 1342 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 1343 snes->vec_func = r; 1344 } 1345 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1346 ierr = DMSNESSetFunction(dm,func,ctx);CHKERRQ(ierr); 1347 PetscFunctionReturn(0); 1348 } 1349 1350 1351 #undef __FUNCT__ 1352 #define __FUNCT__ "SNESSetGS" 1353 /*@C 1354 SNESSetGS - Sets the user nonlinear Gauss-Seidel routine for 1355 use with composed nonlinear solvers. 1356 1357 Input Parameters: 1358 + snes - the SNES context 1359 . gsfunc - function evaluation routine 1360 - ctx - [optional] user-defined context for private data for the 1361 smoother evaluation routine (may be PETSC_NULL) 1362 1363 Calling sequence of func: 1364 $ func (SNES snes,Vec x,Vec b,void *ctx); 1365 1366 + X - solution vector 1367 . B - RHS vector 1368 - ctx - optional user-defined Gauss-Seidel context 1369 1370 Notes: 1371 The GS routines are used by the composed nonlinear solver to generate 1372 a problem appropriate update to the solution, particularly FAS. 1373 1374 Level: intermediate 1375 1376 .keywords: SNES, nonlinear, set, Gauss-Seidel 1377 1378 .seealso: SNESGetFunction(), SNESComputeGS() 1379 @*/ 1380 PetscErrorCode SNESSetGS(SNES snes,PetscErrorCode (*gsfunc)(SNES,Vec,Vec,void*),void *ctx) 1381 { 1382 PetscErrorCode ierr; 1383 DM dm; 1384 1385 PetscFunctionBegin; 1386 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1387 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1388 ierr = DMSNESSetGS(dm,gsfunc,ctx);CHKERRQ(ierr); 1389 PetscFunctionReturn(0); 1390 } 1391 1392 #undef __FUNCT__ 1393 #define __FUNCT__ "SNESSetGSSweeps" 1394 /*@ 1395 SNESSetGSSweeps - Sets the number of sweeps of GS to use. 1396 1397 Input Parameters: 1398 + snes - the SNES context 1399 - sweeps - the number of sweeps of GS to perform. 1400 1401 Level: intermediate 1402 1403 .keywords: SNES, nonlinear, set, Gauss-Siedel 1404 1405 .seealso: SNESSetGS(), SNESGetGS(), SNESSetPC(), SNESGetGSSweeps() 1406 @*/ 1407 1408 PetscErrorCode SNESSetGSSweeps(SNES snes, PetscInt sweeps) { 1409 PetscFunctionBegin; 1410 snes->gssweeps = sweeps; 1411 PetscFunctionReturn(0); 1412 } 1413 1414 1415 #undef __FUNCT__ 1416 #define __FUNCT__ "SNESGetGSSweeps" 1417 /*@ 1418 SNESGetGSSweeps - Gets the number of sweeps GS will use. 1419 1420 Input Parameters: 1421 . snes - the SNES context 1422 1423 Output Parameters: 1424 . sweeps - the number of sweeps of GS to perform. 1425 1426 Level: intermediate 1427 1428 .keywords: SNES, nonlinear, set, Gauss-Siedel 1429 1430 .seealso: SNESSetGS(), SNESGetGS(), SNESSetPC(), SNESSetGSSweeps() 1431 @*/ 1432 PetscErrorCode SNESGetGSSweeps(SNES snes, PetscInt * sweeps) { 1433 PetscFunctionBegin; 1434 *sweeps = snes->gssweeps; 1435 PetscFunctionReturn(0); 1436 } 1437 1438 #undef __FUNCT__ 1439 #define __FUNCT__ "SNESPicardComputeFunction" 1440 PetscErrorCode SNESPicardComputeFunction(SNES snes,Vec x,Vec f,void *ctx) 1441 { 1442 PetscErrorCode ierr; 1443 void *functx,*jacctx; 1444 1445 PetscFunctionBegin; 1446 ierr = SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);CHKERRQ(ierr); 1447 ierr = SNESGetJacobian(snes,PETSC_NULL,PETSC_NULL,PETSC_NULL,&jacctx);CHKERRQ(ierr); 1448 /* A(x)*x - b(x) */ 1449 ierr = (*snes->ops->computepjacobian)(snes,x,&snes->jacobian,&snes->jacobian_pre,&snes->matstruct,jacctx);CHKERRQ(ierr); 1450 ierr = (*snes->ops->computepfunction)(snes,x,f,functx);CHKERRQ(ierr); 1451 ierr = VecView(x,PETSC_VIEWER_BINARY_WORLD);CHKERRQ(ierr); 1452 ierr = VecView(f,PETSC_VIEWER_BINARY_WORLD);CHKERRQ(ierr); 1453 ierr = VecScale(f,-1.0);CHKERRQ(ierr); 1454 ierr = MatMultAdd(snes->jacobian_pre,x,f,f);CHKERRQ(ierr); 1455 PetscFunctionReturn(0); 1456 } 1457 1458 #undef __FUNCT__ 1459 #define __FUNCT__ "SNESPicardComputeJacobian" 1460 PetscErrorCode SNESPicardComputeJacobian(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx) 1461 { 1462 PetscFunctionBegin; 1463 *flag = snes->matstruct; 1464 PetscFunctionReturn(0); 1465 } 1466 1467 #undef __FUNCT__ 1468 #define __FUNCT__ "SNESSetPicard" 1469 /*@C 1470 SNESSetPicard - Use SNES to solve the semilinear-system A(x) x = b(x) via a Picard type iteration (Picard linearization) 1471 1472 Logically Collective on SNES 1473 1474 Input Parameters: 1475 + snes - the SNES context 1476 . r - vector to store function value 1477 . func - function evaluation routine 1478 . jmat - normally the same as mat but you can pass another matrix for which you compute the Jacobian of A(x) x - b(x) (see jmat below) 1479 . mat - matrix to store A 1480 . mfunc - function to compute matrix value 1481 - ctx - [optional] user-defined context for private data for the 1482 function evaluation routine (may be PETSC_NULL) 1483 1484 Calling sequence of func: 1485 $ func (SNES snes,Vec x,Vec f,void *ctx); 1486 1487 + f - function vector 1488 - ctx - optional user-defined function context 1489 1490 Calling sequence of mfunc: 1491 $ mfunc (SNES snes,Vec x,Mat *jmat,Mat *mat,int *flag,void *ctx); 1492 1493 + x - input vector 1494 . jmat - Form Jacobian matrix of A(x) x - b(x) if available, not there is really no reason to use it in this way since then you can just use SNESSetJacobian(), 1495 normally just pass mat in this location 1496 . mat - form A(x) matrix 1497 . flag - flag indicating information about the preconditioner matrix 1498 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 1499 - ctx - [optional] user-defined Jacobian context 1500 1501 Notes: 1502 One can call SNESSetPicard() or SNESSetFunction() (and possibly SNESSetJacobian()) but cannot call both 1503 1504 $ Solves the equation A(x) x = b(x) via the defect correction algorithm A(x^{n}) (x^{n+1} - x^{n}) = b(x^{n}) - A(x^{n})x^{n} 1505 $ Note that when an exact solver is used this corresponds to the "classic" Picard A(x^{n}) x^{n+1} = b(x^{n}) iteration. 1506 1507 Run with -snes_mf_operator to solve the system with Newton's method using A(x^{n}) to construct the preconditioner. 1508 1509 We implement the defect correction form of the Picard iteration because it converges much more generally when inexact linear solvers are used then 1510 the direct Picard iteration A(x^n) x^{n+1} = b(x^n) 1511 1512 There is some controversity over the definition of a Picard iteration for nonlinear systems but almost everyone agrees that it involves a linear solve and some 1513 believe it is the iteration A(x^{n}) x^{n+1} = b(x^{n}) hence we use the name Picard. If anyone has an authoritative reference that defines the Picard iteration 1514 different please contact us at petsc-dev@mcs.anl.gov and we'll have an entirely new argument :-). 1515 1516 Level: beginner 1517 1518 .keywords: SNES, nonlinear, set, function 1519 1520 .seealso: SNESGetFunction(), SNESSetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESGetPicard(), SNESLineSearchPreCheckPicard() 1521 @*/ 1522 PetscErrorCode SNESSetPicard(SNES snes,Vec r,PetscErrorCode (*func)(SNES,Vec,Vec,void*),Mat jmat, Mat mat, PetscErrorCode (*mfunc)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 1523 { 1524 PetscErrorCode ierr; 1525 PetscFunctionBegin; 1526 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1527 snes->ops->computepfunction = func; 1528 snes->ops->computepjacobian = mfunc; 1529 ierr = SNESSetFunction(snes,r,SNESPicardComputeFunction,ctx);CHKERRQ(ierr); 1530 ierr = SNESSetJacobian(snes,jmat,mat,SNESPicardComputeJacobian,ctx);CHKERRQ(ierr); 1531 PetscFunctionReturn(0); 1532 } 1533 1534 #undef __FUNCT__ 1535 #define __FUNCT__ "SNESSetComputeInitialGuess" 1536 /*@C 1537 SNESSetComputeInitialGuess - Sets a routine used to compute an initial guess for the problem 1538 1539 Logically Collective on SNES 1540 1541 Input Parameters: 1542 + snes - the SNES context 1543 . func - function evaluation routine 1544 - ctx - [optional] user-defined context for private data for the 1545 function evaluation routine (may be PETSC_NULL) 1546 1547 Calling sequence of func: 1548 $ func (SNES snes,Vec x,void *ctx); 1549 1550 . f - function vector 1551 - ctx - optional user-defined function context 1552 1553 Level: intermediate 1554 1555 .keywords: SNES, nonlinear, set, function 1556 1557 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian() 1558 @*/ 1559 PetscErrorCode SNESSetComputeInitialGuess(SNES snes,PetscErrorCode (*func)(SNES,Vec,void*),void *ctx) 1560 { 1561 PetscFunctionBegin; 1562 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1563 if (func) snes->ops->computeinitialguess = func; 1564 if (ctx) snes->initialguessP = ctx; 1565 PetscFunctionReturn(0); 1566 } 1567 1568 /* --------------------------------------------------------------- */ 1569 #undef __FUNCT__ 1570 #define __FUNCT__ "SNESGetRhs" 1571 /*@C 1572 SNESGetRhs - Gets the vector for solving F(x) = rhs. If rhs is not set 1573 it assumes a zero right hand side. 1574 1575 Logically Collective on SNES 1576 1577 Input Parameter: 1578 . snes - the SNES context 1579 1580 Output Parameter: 1581 . rhs - the right hand side vector or PETSC_NULL if the right hand side vector is null 1582 1583 Level: intermediate 1584 1585 .keywords: SNES, nonlinear, get, function, right hand side 1586 1587 .seealso: SNESGetSolution(), SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 1588 @*/ 1589 PetscErrorCode SNESGetRhs(SNES snes,Vec *rhs) 1590 { 1591 PetscFunctionBegin; 1592 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1593 PetscValidPointer(rhs,2); 1594 *rhs = snes->vec_rhs; 1595 PetscFunctionReturn(0); 1596 } 1597 1598 #undef __FUNCT__ 1599 #define __FUNCT__ "SNESComputeFunction" 1600 /*@ 1601 SNESComputeFunction - Calls the function that has been set with 1602 SNESSetFunction(). 1603 1604 Collective on SNES 1605 1606 Input Parameters: 1607 + snes - the SNES context 1608 - x - input vector 1609 1610 Output Parameter: 1611 . y - function vector, as set by SNESSetFunction() 1612 1613 Notes: 1614 SNESComputeFunction() is typically used within nonlinear solvers 1615 implementations, so most users would not generally call this routine 1616 themselves. 1617 1618 Level: developer 1619 1620 .keywords: SNES, nonlinear, compute, function 1621 1622 .seealso: SNESSetFunction(), SNESGetFunction() 1623 @*/ 1624 PetscErrorCode SNESComputeFunction(SNES snes,Vec x,Vec y) 1625 { 1626 PetscErrorCode ierr; 1627 DM dm; 1628 SNESDM sdm; 1629 1630 PetscFunctionBegin; 1631 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1632 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 1633 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 1634 PetscCheckSameComm(snes,1,x,2); 1635 PetscCheckSameComm(snes,1,y,3); 1636 VecValidValues(x,2,PETSC_TRUE); 1637 1638 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1639 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 1640 ierr = PetscLogEventBegin(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 1641 if (sdm->computefunction) { 1642 PetscStackPush("SNES user function"); 1643 ierr = (*sdm->computefunction)(snes,x,y,sdm->functionctx);CHKERRQ(ierr); 1644 PetscStackPop; 1645 } else if (snes->dm) { 1646 ierr = DMComputeFunction(snes->dm,x,y);CHKERRQ(ierr); 1647 } else if (snes->vec_rhs) { 1648 ierr = MatMult(snes->jacobian, x, y);CHKERRQ(ierr); 1649 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetFunction() or SNESSetDM() before SNESComputeFunction(), likely called from SNESSolve()."); 1650 if (snes->vec_rhs) { 1651 ierr = VecAXPY(y,-1.0,snes->vec_rhs);CHKERRQ(ierr); 1652 } 1653 snes->nfuncs++; 1654 ierr = PetscLogEventEnd(SNES_FunctionEval,snes,x,y,0);CHKERRQ(ierr); 1655 VecValidValues(y,3,PETSC_FALSE); 1656 PetscFunctionReturn(0); 1657 } 1658 1659 #undef __FUNCT__ 1660 #define __FUNCT__ "SNESComputeGS" 1661 /*@ 1662 SNESComputeGS - Calls the Gauss-Seidel function that has been set with 1663 SNESSetGS(). 1664 1665 Collective on SNES 1666 1667 Input Parameters: 1668 + snes - the SNES context 1669 . x - input vector 1670 - b - rhs vector 1671 1672 Output Parameter: 1673 . x - new solution vector 1674 1675 Notes: 1676 SNESComputeGS() is typically used within composed nonlinear solver 1677 implementations, so most users would not generally call this routine 1678 themselves. 1679 1680 Level: developer 1681 1682 .keywords: SNES, nonlinear, compute, function 1683 1684 .seealso: SNESSetGS(), SNESComputeFunction() 1685 @*/ 1686 PetscErrorCode SNESComputeGS(SNES snes,Vec b,Vec x) 1687 { 1688 PetscErrorCode ierr; 1689 PetscInt i; 1690 DM dm; 1691 SNESDM sdm; 1692 1693 PetscFunctionBegin; 1694 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1695 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 1696 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,3); 1697 PetscCheckSameComm(snes,1,x,2); 1698 if(b) PetscCheckSameComm(snes,1,b,3); 1699 if (b) VecValidValues(b,2,PETSC_TRUE); 1700 ierr = PetscLogEventBegin(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 1701 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1702 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 1703 if (sdm->computegs) { 1704 for (i = 0; i < snes->gssweeps; i++) { 1705 PetscStackPush("SNES user GS"); 1706 ierr = (*sdm->computegs)(snes,x,b,sdm->gsctx);CHKERRQ(ierr); 1707 PetscStackPop; 1708 } 1709 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Must call SNESSetGS() before SNESComputeGS(), likely called from SNESSolve()."); 1710 ierr = PetscLogEventEnd(SNES_GSEval,snes,x,b,0);CHKERRQ(ierr); 1711 VecValidValues(x,3,PETSC_FALSE); 1712 PetscFunctionReturn(0); 1713 } 1714 1715 1716 #undef __FUNCT__ 1717 #define __FUNCT__ "SNESComputeJacobian" 1718 /*@ 1719 SNESComputeJacobian - Computes the Jacobian matrix that has been 1720 set with SNESSetJacobian(). 1721 1722 Collective on SNES and Mat 1723 1724 Input Parameters: 1725 + snes - the SNES context 1726 - x - input vector 1727 1728 Output Parameters: 1729 + A - Jacobian matrix 1730 . B - optional preconditioning matrix 1731 - flag - flag indicating matrix structure (one of, SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER) 1732 1733 Options Database Keys: 1734 + -snes_lag_preconditioner <lag> 1735 . -snes_lag_jacobian <lag> 1736 . -snes_compare_explicit - Compare the computed Jacobian to the finite difference Jacobian and output the differences 1737 . -snes_compare_explicit_draw - Compare the computed Jacobian to the finite difference Jacobian and draw the result 1738 . -snes_compare_explicit_contour - Compare the computed Jacobian to the finite difference Jacobian and draw a contour plot with the result 1739 . -snes_compare_operator - Make the comparison options above use the operator instead of the preconditioning matrix 1740 . -snes_compare_coloring - Compute the finite differece Jacobian using coloring and display norms of difference 1741 . -snes_compare_coloring_display - Compute the finite differece Jacobian using coloring and display verbose differences 1742 . -snes_compare_coloring_threshold - Display only those matrix entries that differ by more than a given threshold 1743 . -snes_compare_coloring_threshold_atol - Absolute tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 1744 . -snes_compare_coloring_threshold_rtol - Relative tolerance for difference in matrix entries to be displayed by -snes_compare_coloring_threshold 1745 . -snes_compare_coloring_draw - Compute the finite differece Jacobian using coloring and draw differences 1746 - -snes_compare_coloring_draw_contour - Compute the finite differece Jacobian using coloring and show contours of matrices and differences 1747 1748 1749 Notes: 1750 Most users should not need to explicitly call this routine, as it 1751 is used internally within the nonlinear solvers. 1752 1753 See KSPSetOperators() for important information about setting the 1754 flag parameter. 1755 1756 Level: developer 1757 1758 .keywords: SNES, compute, Jacobian, matrix 1759 1760 .seealso: SNESSetJacobian(), KSPSetOperators(), MatStructure, SNESSetLagPreconditioner(), SNESSetLagJacobian() 1761 @*/ 1762 PetscErrorCode SNESComputeJacobian(SNES snes,Vec X,Mat *A,Mat *B,MatStructure *flg) 1763 { 1764 PetscErrorCode ierr; 1765 PetscBool flag; 1766 DM dm; 1767 SNESDM sdm; 1768 1769 PetscFunctionBegin; 1770 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 1771 PetscValidHeaderSpecific(X,VEC_CLASSID,2); 1772 PetscValidPointer(flg,5); 1773 PetscCheckSameComm(snes,1,X,2); 1774 VecValidValues(X,2,PETSC_TRUE); 1775 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 1776 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 1777 if (!sdm->computejacobian) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_USER,"Must call SNESSetJacobian(), DMSNESSetJacobian(), DMDASNESSetJacobianLocal(), etc"); 1778 1779 /* make sure that MatAssemblyBegin/End() is called on A matrix if it is matrix free */ 1780 1781 if (snes->lagjacobian == -2) { 1782 snes->lagjacobian = -1; 1783 ierr = PetscInfo(snes,"Recomputing Jacobian/preconditioner because lag is -2 (means compute Jacobian, but then never again) \n");CHKERRQ(ierr); 1784 } else if (snes->lagjacobian == -1) { 1785 *flg = SAME_PRECONDITIONER; 1786 ierr = PetscInfo(snes,"Reusing Jacobian/preconditioner because lag is -1\n");CHKERRQ(ierr); 1787 ierr = PetscTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 1788 if (flag) { 1789 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1790 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1791 } 1792 PetscFunctionReturn(0); 1793 } else if (snes->lagjacobian > 1 && snes->iter % snes->lagjacobian) { 1794 *flg = SAME_PRECONDITIONER; 1795 ierr = PetscInfo2(snes,"Reusing Jacobian/preconditioner because lag is %D and SNES iteration is %D\n",snes->lagjacobian,snes->iter);CHKERRQ(ierr); 1796 ierr = PetscTypeCompare((PetscObject)*A,MATMFFD,&flag);CHKERRQ(ierr); 1797 if (flag) { 1798 ierr = MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1799 ierr = MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1800 } 1801 PetscFunctionReturn(0); 1802 } 1803 1804 *flg = DIFFERENT_NONZERO_PATTERN; 1805 ierr = PetscLogEventBegin(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 1806 PetscStackPush("SNES user Jacobian function"); 1807 ierr = (*sdm->computejacobian)(snes,X,A,B,flg,sdm->jacobianctx);CHKERRQ(ierr); 1808 PetscStackPop; 1809 ierr = PetscLogEventEnd(SNES_JacobianEval,snes,X,*A,*B);CHKERRQ(ierr); 1810 1811 if (snes->lagpreconditioner == -2) { 1812 ierr = PetscInfo(snes,"Rebuilding preconditioner exactly once since lag is -2\n");CHKERRQ(ierr); 1813 snes->lagpreconditioner = -1; 1814 } else if (snes->lagpreconditioner == -1) { 1815 *flg = SAME_PRECONDITIONER; 1816 ierr = PetscInfo(snes,"Reusing preconditioner because lag is -1\n");CHKERRQ(ierr); 1817 } else if (snes->lagpreconditioner > 1 && snes->iter % snes->lagpreconditioner) { 1818 *flg = SAME_PRECONDITIONER; 1819 ierr = PetscInfo2(snes,"Reusing preconditioner because lag is %D and SNES iteration is %D\n",snes->lagpreconditioner,snes->iter);CHKERRQ(ierr); 1820 } 1821 1822 /* make sure user returned a correct Jacobian and preconditioner */ 1823 /* PetscValidHeaderSpecific(*A,MAT_CLASSID,3); 1824 PetscValidHeaderSpecific(*B,MAT_CLASSID,4); */ 1825 { 1826 PetscBool flag = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_operator = PETSC_FALSE; 1827 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit",&flag,PETSC_NULL);CHKERRQ(ierr); 1828 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw",&flag_draw,PETSC_NULL);CHKERRQ(ierr); 1829 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_explicit_draw_contour",&flag_contour,PETSC_NULL);CHKERRQ(ierr); 1830 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_operator",&flag_operator,PETSC_NULL);CHKERRQ(ierr); 1831 if (flag || flag_draw || flag_contour) { 1832 Mat Bexp_mine = PETSC_NULL,Bexp,FDexp; 1833 MatStructure mstruct; 1834 PetscViewer vdraw,vstdout; 1835 PetscBool flg; 1836 if (flag_operator) { 1837 ierr = MatComputeExplicitOperator(*A,&Bexp_mine);CHKERRQ(ierr); 1838 Bexp = Bexp_mine; 1839 } else { 1840 /* See if the preconditioning matrix can be viewed and added directly */ 1841 ierr = PetscTypeCompareAny((PetscObject)*B,&flg,MATSEQAIJ,MATMPIAIJ,MATSEQDENSE,MATMPIDENSE,MATSEQBAIJ,MATMPIBAIJ,MATSEQSBAIJ,MATMPIBAIJ,"");CHKERRQ(ierr); 1842 if (flg) Bexp = *B; 1843 else { 1844 /* If the "preconditioning" matrix is itself MATSHELL or some other type without direct support */ 1845 ierr = MatComputeExplicitOperator(*B,&Bexp_mine);CHKERRQ(ierr); 1846 Bexp = Bexp_mine; 1847 } 1848 } 1849 ierr = MatConvert(Bexp,MATSAME,MAT_INITIAL_MATRIX,&FDexp);CHKERRQ(ierr); 1850 ierr = SNESDefaultComputeJacobian(snes,X,&FDexp,&FDexp,&mstruct,NULL);CHKERRQ(ierr); 1851 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);CHKERRQ(ierr); 1852 if (flag_draw || flag_contour) { 1853 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Explicit Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 1854 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 1855 } else vdraw = PETSC_NULL; 1856 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit %s\n",flag_operator?"Jacobian":"preconditioning Jacobian");CHKERRQ(ierr); 1857 if (flag) {ierr = MatView(Bexp,vstdout);CHKERRQ(ierr);} 1858 if (vdraw) {ierr = MatView(Bexp,vdraw);CHKERRQ(ierr);} 1859 ierr = PetscViewerASCIIPrintf(vstdout,"Finite difference Jacobian\n");CHKERRQ(ierr); 1860 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 1861 if (vdraw) {ierr = MatView(FDexp,vdraw);CHKERRQ(ierr);} 1862 ierr = MatAYPX(FDexp,-1.0,Bexp,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 1863 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian\n");CHKERRQ(ierr); 1864 if (flag) {ierr = MatView(FDexp,vstdout);CHKERRQ(ierr);} 1865 if (vdraw) { /* Always use contour for the difference */ 1866 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 1867 ierr = MatView(FDexp,vdraw);CHKERRQ(ierr); 1868 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 1869 } 1870 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 1871 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 1872 ierr = MatDestroy(&Bexp_mine);CHKERRQ(ierr); 1873 ierr = MatDestroy(&FDexp);CHKERRQ(ierr); 1874 } 1875 } 1876 { 1877 PetscBool flag = PETSC_FALSE,flag_display = PETSC_FALSE,flag_draw = PETSC_FALSE,flag_contour = PETSC_FALSE,flag_threshold = PETSC_FALSE; 1878 PetscReal threshold_atol = PETSC_SQRT_MACHINE_EPSILON,threshold_rtol = 10*PETSC_SQRT_MACHINE_EPSILON; 1879 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring",&flag,PETSC_NULL);CHKERRQ(ierr); 1880 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_display",&flag_display,PETSC_NULL);CHKERRQ(ierr); 1881 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw",&flag_draw,PETSC_NULL);CHKERRQ(ierr); 1882 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_draw_contour",&flag_contour,PETSC_NULL);CHKERRQ(ierr); 1883 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold",&flag_threshold,PETSC_NULL);CHKERRQ(ierr); 1884 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_rtol",&threshold_rtol,PETSC_NULL);CHKERRQ(ierr); 1885 ierr = PetscOptionsGetReal(((PetscObject)snes)->prefix,"-snes_compare_coloring_threshold_atol",&threshold_atol,PETSC_NULL);CHKERRQ(ierr); 1886 if (flag || flag_display || flag_draw || flag_contour || flag_threshold) { 1887 Mat Bfd; 1888 MatStructure mstruct; 1889 PetscViewer vdraw,vstdout; 1890 ISColoring iscoloring; 1891 MatFDColoring matfdcoloring; 1892 PetscErrorCode (*func)(SNES,Vec,Vec,void*); 1893 void *funcctx; 1894 PetscReal norm1,norm2,normmax; 1895 1896 ierr = MatDuplicate(*B,MAT_DO_NOT_COPY_VALUES,&Bfd);CHKERRQ(ierr); 1897 ierr = MatGetColoring(Bfd,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr); 1898 ierr = MatFDColoringCreate(Bfd,iscoloring,&matfdcoloring);CHKERRQ(ierr); 1899 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 1900 1901 /* This method of getting the function is currently unreliable since it doesn't work for DM local functions. */ 1902 ierr = SNESGetFunction(snes,PETSC_NULL,&func,&funcctx);CHKERRQ(ierr); 1903 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode(*)(void))func,funcctx);CHKERRQ(ierr); 1904 ierr = PetscObjectSetOptionsPrefix((PetscObject)matfdcoloring,((PetscObject)snes)->prefix);CHKERRQ(ierr); 1905 ierr = PetscObjectAppendOptionsPrefix((PetscObject)matfdcoloring,"coloring_");CHKERRQ(ierr); 1906 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 1907 ierr = MatFDColoringApply(Bfd,matfdcoloring,X,&mstruct,snes);CHKERRQ(ierr); 1908 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 1909 1910 ierr = PetscViewerASCIIGetStdout(((PetscObject)snes)->comm,&vstdout);CHKERRQ(ierr); 1911 if (flag_draw || flag_contour) { 1912 ierr = PetscViewerDrawOpen(((PetscObject)snes)->comm,0,"Colored Jacobians",PETSC_DECIDE,PETSC_DECIDE,300,300,&vdraw);CHKERRQ(ierr); 1913 if (flag_contour) {ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr);} 1914 } else vdraw = PETSC_NULL; 1915 ierr = PetscViewerASCIIPrintf(vstdout,"Explicit preconditioning Jacobian\n");CHKERRQ(ierr); 1916 if (flag_display) {ierr = MatView(*B,vstdout);CHKERRQ(ierr);} 1917 if (vdraw) {ierr = MatView(*B,vdraw);CHKERRQ(ierr);} 1918 ierr = PetscViewerASCIIPrintf(vstdout,"Colored Finite difference Jacobian\n");CHKERRQ(ierr); 1919 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 1920 if (vdraw) {ierr = MatView(Bfd,vdraw);CHKERRQ(ierr);} 1921 ierr = MatAYPX(Bfd,-1.0,*B,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 1922 ierr = MatNorm(Bfd,NORM_1,&norm1);CHKERRQ(ierr); 1923 ierr = MatNorm(Bfd,NORM_FROBENIUS,&norm2);CHKERRQ(ierr); 1924 ierr = MatNorm(Bfd,NORM_MAX,&normmax);CHKERRQ(ierr); 1925 ierr = PetscViewerASCIIPrintf(vstdout,"User-provided matrix minus finite difference Jacobian, norm1=%G normFrob=%G normmax=%G\n",norm1,norm2,normmax);CHKERRQ(ierr); 1926 if (flag_display) {ierr = MatView(Bfd,vstdout);CHKERRQ(ierr);} 1927 if (vdraw) { /* Always use contour for the difference */ 1928 ierr = PetscViewerPushFormat(vdraw,PETSC_VIEWER_DRAW_CONTOUR);CHKERRQ(ierr); 1929 ierr = MatView(Bfd,vdraw);CHKERRQ(ierr); 1930 ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr); 1931 } 1932 if (flag_contour) {ierr = PetscViewerPopFormat(vdraw);CHKERRQ(ierr);} 1933 1934 if (flag_threshold) { 1935 PetscInt bs,rstart,rend,i; 1936 ierr = MatGetBlockSize(*B,&bs);CHKERRQ(ierr); 1937 ierr = MatGetOwnershipRange(*B,&rstart,&rend);CHKERRQ(ierr); 1938 for (i=rstart; i<rend; i++) { 1939 const PetscScalar *ba,*ca; 1940 const PetscInt *bj,*cj; 1941 PetscInt bn,cn,j,maxentrycol = -1,maxdiffcol = -1,maxrdiffcol = -1; 1942 PetscReal maxentry = 0,maxdiff = 0,maxrdiff = 0; 1943 ierr = MatGetRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 1944 ierr = MatGetRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 1945 if (bn != cn) SETERRQ(((PetscObject)*A)->comm,PETSC_ERR_PLIB,"Unexpected different nonzero pattern in -snes_compare_coloring_threshold"); 1946 for (j=0; j<bn; j++) { 1947 PetscReal rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 1948 if (PetscAbsScalar(ba[j]) > PetscAbs(maxentry)) { 1949 maxentrycol = bj[j]; 1950 maxentry = PetscRealPart(ba[j]); 1951 } 1952 if (PetscAbsScalar(ca[j]) > PetscAbs(maxdiff)) { 1953 maxdiffcol = bj[j]; 1954 maxdiff = PetscRealPart(ca[j]); 1955 } 1956 if (rdiff > maxrdiff) { 1957 maxrdiffcol = bj[j]; 1958 maxrdiff = rdiff; 1959 } 1960 } 1961 if (maxrdiff > 1) { 1962 ierr = PetscViewerASCIIPrintf(vstdout,"row %D (maxentry=%G at %D, maxdiff=%G at %D, maxrdiff=%G at %D):",i,maxentry,maxentrycol,maxdiff,maxdiffcol,maxrdiff,maxrdiffcol);CHKERRQ(ierr); 1963 for (j=0; j<bn; j++) { 1964 PetscReal rdiff; 1965 rdiff = PetscAbsScalar(ca[j]) / (threshold_atol + threshold_rtol*PetscAbsScalar(ba[j])); 1966 if (rdiff > 1) { 1967 ierr = PetscViewerASCIIPrintf(vstdout," (%D,%G:%G)",bj[j],PetscRealPart(ba[j]),PetscRealPart(ca[j]));CHKERRQ(ierr); 1968 } 1969 } 1970 ierr = PetscViewerASCIIPrintf(vstdout,"\n",i,maxentry,maxdiff,maxrdiff);CHKERRQ(ierr); 1971 } 1972 ierr = MatRestoreRow(*B,i,&bn,&bj,&ba);CHKERRQ(ierr); 1973 ierr = MatRestoreRow(Bfd,i,&cn,&cj,&ca);CHKERRQ(ierr); 1974 } 1975 } 1976 ierr = PetscViewerDestroy(&vdraw);CHKERRQ(ierr); 1977 ierr = MatDestroy(&Bfd);CHKERRQ(ierr); 1978 } 1979 } 1980 PetscFunctionReturn(0); 1981 } 1982 1983 #undef __FUNCT__ 1984 #define __FUNCT__ "SNESSetJacobian" 1985 /*@C 1986 SNESSetJacobian - Sets the function to compute Jacobian as well as the 1987 location to store the matrix. 1988 1989 Logically Collective on SNES and Mat 1990 1991 Input Parameters: 1992 + snes - the SNES context 1993 . A - Jacobian matrix 1994 . B - preconditioner matrix (usually same as the Jacobian) 1995 . func - Jacobian evaluation routine (if PETSC_NULL then SNES retains any previously set value) 1996 - ctx - [optional] user-defined context for private data for the 1997 Jacobian evaluation routine (may be PETSC_NULL) (if PETSC_NULL then SNES retains any previously set value) 1998 1999 Calling sequence of func: 2000 $ func (SNES snes,Vec x,Mat *A,Mat *B,int *flag,void *ctx); 2001 2002 + x - input vector 2003 . A - Jacobian matrix 2004 . B - preconditioner matrix, usually the same as A 2005 . flag - flag indicating information about the preconditioner matrix 2006 structure (same as flag in KSPSetOperators()), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER 2007 - ctx - [optional] user-defined Jacobian context 2008 2009 Notes: 2010 See KSPSetOperators() for important information about setting the flag 2011 output parameter in the routine func(). Be sure to read this information! 2012 2013 The routine func() takes Mat * as the matrix arguments rather than Mat. 2014 This allows the Jacobian evaluation routine to replace A and/or B with a 2015 completely new new matrix structure (not just different matrix elements) 2016 when appropriate, for instance, if the nonzero structure is changing 2017 throughout the global iterations. 2018 2019 If the A matrix and B matrix are different you must call MatAssemblyBegin/End() on 2020 each matrix. 2021 2022 If using SNESDefaultComputeJacobianColor() to assemble a Jacobian, the ctx argument 2023 must be a MatFDColoring. 2024 2025 Other defect-correction schemes can be used by computing a different matrix in place of the Jacobian. One common 2026 example is to use the "Picard linearization" which only differentiates through the highest order parts of each term. 2027 2028 Level: beginner 2029 2030 .keywords: SNES, nonlinear, set, Jacobian, matrix 2031 2032 .seealso: KSPSetOperators(), SNESSetFunction(), MatMFFDComputeJacobian(), SNESDefaultComputeJacobianColor(), MatStructure 2033 @*/ 2034 PetscErrorCode SNESSetJacobian(SNES snes,Mat A,Mat B,PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void *ctx) 2035 { 2036 PetscErrorCode ierr; 2037 DM dm; 2038 2039 PetscFunctionBegin; 2040 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2041 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2); 2042 if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3); 2043 if (A) PetscCheckSameComm(snes,1,A,2); 2044 if (B) PetscCheckSameComm(snes,1,B,3); 2045 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2046 ierr = DMSNESSetJacobian(dm,func,ctx);CHKERRQ(ierr); 2047 if (A) { 2048 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2049 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2050 snes->jacobian = A; 2051 } 2052 if (B) { 2053 ierr = PetscObjectReference((PetscObject)B);CHKERRQ(ierr); 2054 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2055 snes->jacobian_pre = B; 2056 } 2057 PetscFunctionReturn(0); 2058 } 2059 2060 #undef __FUNCT__ 2061 #define __FUNCT__ "SNESGetJacobian" 2062 /*@C 2063 SNESGetJacobian - Returns the Jacobian matrix and optionally the user 2064 provided context for evaluating the Jacobian. 2065 2066 Not Collective, but Mat object will be parallel if SNES object is 2067 2068 Input Parameter: 2069 . snes - the nonlinear solver context 2070 2071 Output Parameters: 2072 + A - location to stash Jacobian matrix (or PETSC_NULL) 2073 . B - location to stash preconditioner matrix (or PETSC_NULL) 2074 . func - location to put Jacobian function (or PETSC_NULL) 2075 - ctx - location to stash Jacobian ctx (or PETSC_NULL) 2076 2077 Level: advanced 2078 2079 .seealso: SNESSetJacobian(), SNESComputeJacobian() 2080 @*/ 2081 PetscErrorCode SNESGetJacobian(SNES snes,Mat *A,Mat *B,PetscErrorCode (**func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*),void **ctx) 2082 { 2083 PetscErrorCode ierr; 2084 DM dm; 2085 SNESDM sdm; 2086 2087 PetscFunctionBegin; 2088 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2089 if (A) *A = snes->jacobian; 2090 if (B) *B = snes->jacobian_pre; 2091 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2092 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 2093 if (func) *func = sdm->computejacobian; 2094 if (ctx) *ctx = sdm->jacobianctx; 2095 PetscFunctionReturn(0); 2096 } 2097 2098 /* ----- Routines to initialize and destroy a nonlinear solver ---- */ 2099 2100 #undef __FUNCT__ 2101 #define __FUNCT__ "SNESSetUp" 2102 /*@ 2103 SNESSetUp - Sets up the internal data structures for the later use 2104 of a nonlinear solver. 2105 2106 Collective on SNES 2107 2108 Input Parameters: 2109 . snes - the SNES context 2110 2111 Notes: 2112 For basic use of the SNES solvers the user need not explicitly call 2113 SNESSetUp(), since these actions will automatically occur during 2114 the call to SNESSolve(). However, if one wishes to control this 2115 phase separately, SNESSetUp() should be called after SNESCreate() 2116 and optional routines of the form SNESSetXXX(), but before SNESSolve(). 2117 2118 Level: advanced 2119 2120 .keywords: SNES, nonlinear, setup 2121 2122 .seealso: SNESCreate(), SNESSolve(), SNESDestroy() 2123 @*/ 2124 PetscErrorCode SNESSetUp(SNES snes) 2125 { 2126 PetscErrorCode ierr; 2127 DM dm; 2128 SNESDM sdm; 2129 2130 PetscFunctionBegin; 2131 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2132 if (snes->setupcalled) PetscFunctionReturn(0); 2133 2134 if (!((PetscObject)snes)->type_name) { 2135 ierr = SNESSetType(snes,SNESLS);CHKERRQ(ierr); 2136 } 2137 2138 ierr = SNESGetFunction(snes,&snes->vec_func,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr); 2139 if (snes->vec_func == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be function vector"); 2140 if (snes->vec_rhs == snes->vec_sol) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_IDN,"Solution vector cannot be right hand side vector"); 2141 2142 if (!snes->vec_sol_update /* && snes->vec_sol */) { 2143 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 2144 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 2145 } 2146 2147 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 2148 ierr = DMSNESSetUpLegacy(dm);CHKERRQ(ierr); /* To be removed when function routines are taken out of the DM package */ 2149 ierr = DMSNESGetContext(dm,&sdm);CHKERRQ(ierr); 2150 if (!sdm->computefunction) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must provide a residual function with SNESSetFunction(), DMSNESSetFunction(), DMDASNESSetFunctionLocal(), etc"); 2151 if (!snes->vec_func) { 2152 ierr = DMCreateGlobalVector(dm,&snes->vec_func);CHKERRQ(ierr); 2153 } 2154 2155 if (!snes->ksp) {ierr = SNESGetKSP(snes, &snes->ksp);CHKERRQ(ierr);} 2156 2157 if (!snes->linesearch) {ierr = SNESGetPetscLineSearch(snes, &snes->linesearch);} 2158 2159 if (snes->ops->usercompute && !snes->user) { 2160 ierr = (*snes->ops->usercompute)(snes,(void**)&snes->user);CHKERRQ(ierr); 2161 } 2162 2163 if (snes->ops->setup) { 2164 ierr = (*snes->ops->setup)(snes);CHKERRQ(ierr); 2165 } 2166 2167 snes->setupcalled = PETSC_TRUE; 2168 PetscFunctionReturn(0); 2169 } 2170 2171 #undef __FUNCT__ 2172 #define __FUNCT__ "SNESReset" 2173 /*@ 2174 SNESReset - Resets a SNES context to the snessetupcalled = 0 state and removes any allocated Vecs and Mats 2175 2176 Collective on SNES 2177 2178 Input Parameter: 2179 . snes - iterative context obtained from SNESCreate() 2180 2181 Level: intermediate 2182 2183 Notes: Also calls the application context destroy routine set with SNESSetComputeApplicationContext() 2184 2185 .keywords: SNES, destroy 2186 2187 .seealso: SNESCreate(), SNESSetUp(), SNESSolve() 2188 @*/ 2189 PetscErrorCode SNESReset(SNES snes) 2190 { 2191 PetscErrorCode ierr; 2192 2193 PetscFunctionBegin; 2194 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2195 if (snes->ops->userdestroy && snes->user) { 2196 ierr = (*snes->ops->userdestroy)((void**)&snes->user);CHKERRQ(ierr); 2197 snes->user = PETSC_NULL; 2198 } 2199 if (snes->pc) { 2200 ierr = SNESReset(snes->pc);CHKERRQ(ierr); 2201 } 2202 2203 if (snes->ops->reset) { 2204 ierr = (*snes->ops->reset)(snes);CHKERRQ(ierr); 2205 } 2206 if (snes->ksp) { 2207 ierr = KSPReset(snes->ksp);CHKERRQ(ierr); 2208 } 2209 2210 if (snes->linesearch) { 2211 ierr = PetscLineSearchReset(snes->linesearch);CHKERRQ(ierr); 2212 } 2213 2214 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 2215 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 2216 ierr = VecDestroy(&snes->vec_sol_update);CHKERRQ(ierr); 2217 ierr = VecDestroy(&snes->vec_func);CHKERRQ(ierr); 2218 ierr = MatDestroy(&snes->jacobian);CHKERRQ(ierr); 2219 ierr = MatDestroy(&snes->jacobian_pre);CHKERRQ(ierr); 2220 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 2221 ierr = VecDestroyVecs(snes->nvwork,&snes->vwork);CHKERRQ(ierr); 2222 snes->nwork = snes->nvwork = 0; 2223 snes->setupcalled = PETSC_FALSE; 2224 PetscFunctionReturn(0); 2225 } 2226 2227 #undef __FUNCT__ 2228 #define __FUNCT__ "SNESDestroy" 2229 /*@ 2230 SNESDestroy - Destroys the nonlinear solver context that was created 2231 with SNESCreate(). 2232 2233 Collective on SNES 2234 2235 Input Parameter: 2236 . snes - the SNES context 2237 2238 Level: beginner 2239 2240 .keywords: SNES, nonlinear, destroy 2241 2242 .seealso: SNESCreate(), SNESSolve() 2243 @*/ 2244 PetscErrorCode SNESDestroy(SNES *snes) 2245 { 2246 PetscErrorCode ierr; 2247 2248 PetscFunctionBegin; 2249 if (!*snes) PetscFunctionReturn(0); 2250 PetscValidHeaderSpecific((*snes),SNES_CLASSID,1); 2251 if (--((PetscObject)(*snes))->refct > 0) {*snes = 0; PetscFunctionReturn(0);} 2252 2253 ierr = SNESReset((*snes));CHKERRQ(ierr); 2254 ierr = SNESDestroy(&(*snes)->pc);CHKERRQ(ierr); 2255 2256 /* if memory was published with AMS then destroy it */ 2257 ierr = PetscObjectDepublish((*snes));CHKERRQ(ierr); 2258 if ((*snes)->ops->destroy) {ierr = (*((*snes))->ops->destroy)((*snes));CHKERRQ(ierr);} 2259 2260 ierr = DMDestroy(&(*snes)->dm);CHKERRQ(ierr); 2261 ierr = KSPDestroy(&(*snes)->ksp);CHKERRQ(ierr); 2262 ierr = PetscLineSearchDestroy(&(*snes)->linesearch);CHKERRQ(ierr); 2263 2264 ierr = PetscFree((*snes)->kspconvctx);CHKERRQ(ierr); 2265 if ((*snes)->ops->convergeddestroy) { 2266 ierr = (*(*snes)->ops->convergeddestroy)((*snes)->cnvP);CHKERRQ(ierr); 2267 } 2268 if ((*snes)->conv_malloc) { 2269 ierr = PetscFree((*snes)->conv_hist);CHKERRQ(ierr); 2270 ierr = PetscFree((*snes)->conv_hist_its);CHKERRQ(ierr); 2271 } 2272 ierr = PetscViewerDestroy(&(*snes)->ls_monitor);CHKERRQ(ierr); 2273 ierr = SNESMonitorCancel((*snes));CHKERRQ(ierr); 2274 ierr = PetscHeaderDestroy(snes);CHKERRQ(ierr); 2275 PetscFunctionReturn(0); 2276 } 2277 2278 /* ----------- Routines to set solver parameters ---------- */ 2279 2280 #undef __FUNCT__ 2281 #define __FUNCT__ "SNESSetLagPreconditioner" 2282 /*@ 2283 SNESSetLagPreconditioner - Determines when the preconditioner is rebuilt in the nonlinear solve. 2284 2285 Logically Collective on SNES 2286 2287 Input Parameters: 2288 + snes - the SNES context 2289 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2290 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2291 2292 Options Database Keys: 2293 . -snes_lag_preconditioner <lag> 2294 2295 Notes: 2296 The default is 1 2297 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2298 If -1 is used before the very first nonlinear solve the preconditioner is still built because there is no previous preconditioner to use 2299 2300 Level: intermediate 2301 2302 .keywords: SNES, nonlinear, set, convergence, tolerances 2303 2304 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2305 2306 @*/ 2307 PetscErrorCode SNESSetLagPreconditioner(SNES snes,PetscInt lag) 2308 { 2309 PetscFunctionBegin; 2310 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2311 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2312 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2313 PetscValidLogicalCollectiveInt(snes,lag,2); 2314 snes->lagpreconditioner = lag; 2315 PetscFunctionReturn(0); 2316 } 2317 2318 #undef __FUNCT__ 2319 #define __FUNCT__ "SNESSetGridSequence" 2320 /*@ 2321 SNESSetGridSequence - sets the number of steps of grid sequencing that SNES does 2322 2323 Logically Collective on SNES 2324 2325 Input Parameters: 2326 + snes - the SNES context 2327 - steps - the number of refinements to do, defaults to 0 2328 2329 Options Database Keys: 2330 . -snes_grid_sequence <steps> 2331 2332 Level: intermediate 2333 2334 Notes: 2335 Use SNESGetSolution() to extract the fine grid solution after grid sequencing. 2336 2337 .keywords: SNES, nonlinear, set, convergence, tolerances 2338 2339 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagJacobian(), SNESGetLagJacobian() 2340 2341 @*/ 2342 PetscErrorCode SNESSetGridSequence(SNES snes,PetscInt steps) 2343 { 2344 PetscFunctionBegin; 2345 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2346 PetscValidLogicalCollectiveInt(snes,steps,2); 2347 snes->gridsequence = steps; 2348 PetscFunctionReturn(0); 2349 } 2350 2351 #undef __FUNCT__ 2352 #define __FUNCT__ "SNESGetLagPreconditioner" 2353 /*@ 2354 SNESGetLagPreconditioner - Indicates how often the preconditioner is rebuilt 2355 2356 Not Collective 2357 2358 Input Parameter: 2359 . snes - the SNES context 2360 2361 Output Parameter: 2362 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2363 the Jacobian is built etc. -2 indicates rebuild preconditioner at next chance but then never rebuild after that 2364 2365 Options Database Keys: 2366 . -snes_lag_preconditioner <lag> 2367 2368 Notes: 2369 The default is 1 2370 The preconditioner is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2371 2372 Level: intermediate 2373 2374 .keywords: SNES, nonlinear, set, convergence, tolerances 2375 2376 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagPreconditioner() 2377 2378 @*/ 2379 PetscErrorCode SNESGetLagPreconditioner(SNES snes,PetscInt *lag) 2380 { 2381 PetscFunctionBegin; 2382 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2383 *lag = snes->lagpreconditioner; 2384 PetscFunctionReturn(0); 2385 } 2386 2387 #undef __FUNCT__ 2388 #define __FUNCT__ "SNESSetLagJacobian" 2389 /*@ 2390 SNESSetLagJacobian - Determines when the Jacobian is rebuilt in the nonlinear solve. See SNESSetLagPreconditioner() for determining how 2391 often the preconditioner is rebuilt. 2392 2393 Logically Collective on SNES 2394 2395 Input Parameters: 2396 + snes - the SNES context 2397 - lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2398 the Jacobian is built etc. -2 means rebuild at next chance but then never again 2399 2400 Options Database Keys: 2401 . -snes_lag_jacobian <lag> 2402 2403 Notes: 2404 The default is 1 2405 The Jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2406 If -1 is used before the very first nonlinear solve the CODE WILL FAIL! because no Jacobian is used, use -2 to indicate you want it recomputed 2407 at the next Newton step but never again (unless it is reset to another value) 2408 2409 Level: intermediate 2410 2411 .keywords: SNES, nonlinear, set, convergence, tolerances 2412 2413 .seealso: SNESSetTrustRegionTolerance(), SNESGetLagPreconditioner(), SNESSetLagPreconditioner(), SNESGetLagJacobian() 2414 2415 @*/ 2416 PetscErrorCode SNESSetLagJacobian(SNES snes,PetscInt lag) 2417 { 2418 PetscFunctionBegin; 2419 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2420 if (lag < -2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be -2, -1, 1 or greater"); 2421 if (!lag) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag cannot be 0"); 2422 PetscValidLogicalCollectiveInt(snes,lag,2); 2423 snes->lagjacobian = lag; 2424 PetscFunctionReturn(0); 2425 } 2426 2427 #undef __FUNCT__ 2428 #define __FUNCT__ "SNESGetLagJacobian" 2429 /*@ 2430 SNESGetLagJacobian - Indicates how often the Jacobian is rebuilt. See SNESGetLagPreconditioner() to determine when the preconditioner is rebuilt 2431 2432 Not Collective 2433 2434 Input Parameter: 2435 . snes - the SNES context 2436 2437 Output Parameter: 2438 . lag - -1 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is computed within a single nonlinear solve, 2 means every second time 2439 the Jacobian is built etc. 2440 2441 Options Database Keys: 2442 . -snes_lag_jacobian <lag> 2443 2444 Notes: 2445 The default is 1 2446 The jacobian is ALWAYS built in the first iteration of a nonlinear solve unless lag is -1 2447 2448 Level: intermediate 2449 2450 .keywords: SNES, nonlinear, set, convergence, tolerances 2451 2452 .seealso: SNESSetTrustRegionTolerance(), SNESSetLagJacobian(), SNESSetLagPreconditioner(), SNESGetLagPreconditioner() 2453 2454 @*/ 2455 PetscErrorCode SNESGetLagJacobian(SNES snes,PetscInt *lag) 2456 { 2457 PetscFunctionBegin; 2458 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2459 *lag = snes->lagjacobian; 2460 PetscFunctionReturn(0); 2461 } 2462 2463 #undef __FUNCT__ 2464 #define __FUNCT__ "SNESSetTolerances" 2465 /*@ 2466 SNESSetTolerances - Sets various parameters used in convergence tests. 2467 2468 Logically Collective on SNES 2469 2470 Input Parameters: 2471 + snes - the SNES context 2472 . abstol - absolute convergence tolerance 2473 . rtol - relative convergence tolerance 2474 . stol - convergence tolerance in terms of the norm 2475 of the change in the solution between steps 2476 . maxit - maximum number of iterations 2477 - maxf - maximum number of function evaluations 2478 2479 Options Database Keys: 2480 + -snes_atol <abstol> - Sets abstol 2481 . -snes_rtol <rtol> - Sets rtol 2482 . -snes_stol <stol> - Sets stol 2483 . -snes_max_it <maxit> - Sets maxit 2484 - -snes_max_funcs <maxf> - Sets maxf 2485 2486 Notes: 2487 The default maximum number of iterations is 50. 2488 The default maximum number of function evaluations is 1000. 2489 2490 Level: intermediate 2491 2492 .keywords: SNES, nonlinear, set, convergence, tolerances 2493 2494 .seealso: SNESSetTrustRegionTolerance() 2495 @*/ 2496 PetscErrorCode SNESSetTolerances(SNES snes,PetscReal abstol,PetscReal rtol,PetscReal stol,PetscInt maxit,PetscInt maxf) 2497 { 2498 PetscFunctionBegin; 2499 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2500 PetscValidLogicalCollectiveReal(snes,abstol,2); 2501 PetscValidLogicalCollectiveReal(snes,rtol,3); 2502 PetscValidLogicalCollectiveReal(snes,stol,4); 2503 PetscValidLogicalCollectiveInt(snes,maxit,5); 2504 PetscValidLogicalCollectiveInt(snes,maxf,6); 2505 2506 if (abstol != PETSC_DEFAULT) { 2507 if (abstol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol); 2508 snes->abstol = abstol; 2509 } 2510 if (rtol != PETSC_DEFAULT) { 2511 if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %G must be non-negative and less than 1.0",rtol); 2512 snes->rtol = rtol; 2513 } 2514 if (stol != PETSC_DEFAULT) { 2515 if (stol < 0.0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Step tolerance %G must be non-negative",stol); 2516 snes->xtol = stol; 2517 } 2518 if (maxit != PETSC_DEFAULT) { 2519 if (maxit < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxit); 2520 snes->max_its = maxit; 2521 } 2522 if (maxf != PETSC_DEFAULT) { 2523 if (maxf < 0) SETERRQ1(((PetscObject)snes)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of function evaluations %D must be non-negative",maxf); 2524 snes->max_funcs = maxf; 2525 } 2526 PetscFunctionReturn(0); 2527 } 2528 2529 #undef __FUNCT__ 2530 #define __FUNCT__ "SNESGetTolerances" 2531 /*@ 2532 SNESGetTolerances - Gets various parameters used in convergence tests. 2533 2534 Not Collective 2535 2536 Input Parameters: 2537 + snes - the SNES context 2538 . atol - absolute convergence tolerance 2539 . rtol - relative convergence tolerance 2540 . stol - convergence tolerance in terms of the norm 2541 of the change in the solution between steps 2542 . maxit - maximum number of iterations 2543 - maxf - maximum number of function evaluations 2544 2545 Notes: 2546 The user can specify PETSC_NULL for any parameter that is not needed. 2547 2548 Level: intermediate 2549 2550 .keywords: SNES, nonlinear, get, convergence, tolerances 2551 2552 .seealso: SNESSetTolerances() 2553 @*/ 2554 PetscErrorCode SNESGetTolerances(SNES snes,PetscReal *atol,PetscReal *rtol,PetscReal *stol,PetscInt *maxit,PetscInt *maxf) 2555 { 2556 PetscFunctionBegin; 2557 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2558 if (atol) *atol = snes->abstol; 2559 if (rtol) *rtol = snes->rtol; 2560 if (stol) *stol = snes->xtol; 2561 if (maxit) *maxit = snes->max_its; 2562 if (maxf) *maxf = snes->max_funcs; 2563 PetscFunctionReturn(0); 2564 } 2565 2566 #undef __FUNCT__ 2567 #define __FUNCT__ "SNESSetTrustRegionTolerance" 2568 /*@ 2569 SNESSetTrustRegionTolerance - Sets the trust region parameter tolerance. 2570 2571 Logically Collective on SNES 2572 2573 Input Parameters: 2574 + snes - the SNES context 2575 - tol - tolerance 2576 2577 Options Database Key: 2578 . -snes_trtol <tol> - Sets tol 2579 2580 Level: intermediate 2581 2582 .keywords: SNES, nonlinear, set, trust region, tolerance 2583 2584 .seealso: SNESSetTolerances() 2585 @*/ 2586 PetscErrorCode SNESSetTrustRegionTolerance(SNES snes,PetscReal tol) 2587 { 2588 PetscFunctionBegin; 2589 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2590 PetscValidLogicalCollectiveReal(snes,tol,2); 2591 snes->deltatol = tol; 2592 PetscFunctionReturn(0); 2593 } 2594 2595 /* 2596 Duplicate the lg monitors for SNES from KSP; for some reason with 2597 dynamic libraries things don't work under Sun4 if we just use 2598 macros instead of functions 2599 */ 2600 #undef __FUNCT__ 2601 #define __FUNCT__ "SNESMonitorLG" 2602 PetscErrorCode SNESMonitorLG(SNES snes,PetscInt it,PetscReal norm,void *ctx) 2603 { 2604 PetscErrorCode ierr; 2605 2606 PetscFunctionBegin; 2607 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2608 ierr = KSPMonitorLG((KSP)snes,it,norm,ctx);CHKERRQ(ierr); 2609 PetscFunctionReturn(0); 2610 } 2611 2612 #undef __FUNCT__ 2613 #define __FUNCT__ "SNESMonitorLGCreate" 2614 PetscErrorCode SNESMonitorLGCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 2615 { 2616 PetscErrorCode ierr; 2617 2618 PetscFunctionBegin; 2619 ierr = KSPMonitorLGCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 2620 PetscFunctionReturn(0); 2621 } 2622 2623 #undef __FUNCT__ 2624 #define __FUNCT__ "SNESMonitorLGDestroy" 2625 PetscErrorCode SNESMonitorLGDestroy(PetscDrawLG *draw) 2626 { 2627 PetscErrorCode ierr; 2628 2629 PetscFunctionBegin; 2630 ierr = KSPMonitorLGDestroy(draw);CHKERRQ(ierr); 2631 PetscFunctionReturn(0); 2632 } 2633 2634 extern PetscErrorCode SNESMonitorRange_Private(SNES,PetscInt,PetscReal*); 2635 #undef __FUNCT__ 2636 #define __FUNCT__ "SNESMonitorLGRange" 2637 PetscErrorCode SNESMonitorLGRange(SNES snes,PetscInt n,PetscReal rnorm,void *monctx) 2638 { 2639 PetscDrawLG lg; 2640 PetscErrorCode ierr; 2641 PetscReal x,y,per; 2642 PetscViewer v = (PetscViewer)monctx; 2643 static PetscReal prev; /* should be in the context */ 2644 PetscDraw draw; 2645 PetscFunctionBegin; 2646 if (!monctx) { 2647 MPI_Comm comm; 2648 2649 ierr = PetscObjectGetComm((PetscObject)snes,&comm);CHKERRQ(ierr); 2650 v = PETSC_VIEWER_DRAW_(comm); 2651 } 2652 ierr = PetscViewerDrawGetDrawLG(v,0,&lg);CHKERRQ(ierr); 2653 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 2654 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 2655 ierr = PetscDrawSetTitle(draw,"Residual norm");CHKERRQ(ierr); 2656 x = (PetscReal) n; 2657 if (rnorm > 0.0) y = log10(rnorm); else y = -15.0; 2658 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 2659 if (n < 20 || !(n % 5)) { 2660 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 2661 } 2662 2663 ierr = PetscViewerDrawGetDrawLG(v,1,&lg);CHKERRQ(ierr); 2664 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 2665 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 2666 ierr = PetscDrawSetTitle(draw,"% elemts > .2*max elemt");CHKERRQ(ierr); 2667 ierr = SNESMonitorRange_Private(snes,n,&per);CHKERRQ(ierr); 2668 x = (PetscReal) n; 2669 y = 100.0*per; 2670 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 2671 if (n < 20 || !(n % 5)) { 2672 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 2673 } 2674 2675 ierr = PetscViewerDrawGetDrawLG(v,2,&lg);CHKERRQ(ierr); 2676 if (!n) {prev = rnorm;ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 2677 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 2678 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm");CHKERRQ(ierr); 2679 x = (PetscReal) n; 2680 y = (prev - rnorm)/prev; 2681 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 2682 if (n < 20 || !(n % 5)) { 2683 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 2684 } 2685 2686 ierr = PetscViewerDrawGetDrawLG(v,3,&lg);CHKERRQ(ierr); 2687 if (!n) {ierr = PetscDrawLGReset(lg);CHKERRQ(ierr);} 2688 ierr = PetscDrawLGGetDraw(lg,&draw);CHKERRQ(ierr); 2689 ierr = PetscDrawSetTitle(draw,"(norm -oldnorm)/oldnorm*(% > .2 max)");CHKERRQ(ierr); 2690 x = (PetscReal) n; 2691 y = (prev - rnorm)/(prev*per); 2692 if (n > 2) { /*skip initial crazy value */ 2693 ierr = PetscDrawLGAddPoint(lg,&x,&y);CHKERRQ(ierr); 2694 } 2695 if (n < 20 || !(n % 5)) { 2696 ierr = PetscDrawLGDraw(lg);CHKERRQ(ierr); 2697 } 2698 prev = rnorm; 2699 PetscFunctionReturn(0); 2700 } 2701 2702 #undef __FUNCT__ 2703 #define __FUNCT__ "SNESMonitorLGRangeCreate" 2704 PetscErrorCode SNESMonitorLGRangeCreate(const char host[],const char label[],int x,int y,int m,int n,PetscDrawLG *draw) 2705 { 2706 PetscErrorCode ierr; 2707 2708 PetscFunctionBegin; 2709 ierr = KSPMonitorLGCreate(host,label,x,y,m,n,draw);CHKERRQ(ierr); 2710 PetscFunctionReturn(0); 2711 } 2712 2713 #undef __FUNCT__ 2714 #define __FUNCT__ "SNESMonitorLGRangeDestroy" 2715 PetscErrorCode SNESMonitorLGRangeDestroy(PetscDrawLG *draw) 2716 { 2717 PetscErrorCode ierr; 2718 2719 PetscFunctionBegin; 2720 ierr = KSPMonitorLGDestroy(draw);CHKERRQ(ierr); 2721 PetscFunctionReturn(0); 2722 } 2723 2724 #undef __FUNCT__ 2725 #define __FUNCT__ "SNESMonitor" 2726 /*@ 2727 SNESMonitor - runs the user provided monitor routines, if they exist 2728 2729 Collective on SNES 2730 2731 Input Parameters: 2732 + snes - nonlinear solver context obtained from SNESCreate() 2733 . iter - iteration number 2734 - rnorm - relative norm of the residual 2735 2736 Notes: 2737 This routine is called by the SNES implementations. 2738 It does not typically need to be called by the user. 2739 2740 Level: developer 2741 2742 .seealso: SNESMonitorSet() 2743 @*/ 2744 PetscErrorCode SNESMonitor(SNES snes,PetscInt iter,PetscReal rnorm) 2745 { 2746 PetscErrorCode ierr; 2747 PetscInt i,n = snes->numbermonitors; 2748 2749 PetscFunctionBegin; 2750 for (i=0; i<n; i++) { 2751 ierr = (*snes->monitor[i])(snes,iter,rnorm,snes->monitorcontext[i]);CHKERRQ(ierr); 2752 } 2753 PetscFunctionReturn(0); 2754 } 2755 2756 /* ------------ Routines to set performance monitoring options ----------- */ 2757 2758 #undef __FUNCT__ 2759 #define __FUNCT__ "SNESMonitorSet" 2760 /*@C 2761 SNESMonitorSet - Sets an ADDITIONAL function that is to be used at every 2762 iteration of the nonlinear solver to display the iteration's 2763 progress. 2764 2765 Logically Collective on SNES 2766 2767 Input Parameters: 2768 + snes - the SNES context 2769 . func - monitoring routine 2770 . mctx - [optional] user-defined context for private data for the 2771 monitor routine (use PETSC_NULL if no context is desired) 2772 - monitordestroy - [optional] routine that frees monitor context 2773 (may be PETSC_NULL) 2774 2775 Calling sequence of func: 2776 $ int func(SNES snes,PetscInt its, PetscReal norm,void *mctx) 2777 2778 + snes - the SNES context 2779 . its - iteration number 2780 . norm - 2-norm function value (may be estimated) 2781 - mctx - [optional] monitoring context 2782 2783 Options Database Keys: 2784 + -snes_monitor - sets SNESMonitorDefault() 2785 . -snes_monitor_draw - sets line graph monitor, 2786 uses SNESMonitorLGCreate() 2787 - -snes_monitor_cancel - cancels all monitors that have 2788 been hardwired into a code by 2789 calls to SNESMonitorSet(), but 2790 does not cancel those set via 2791 the options database. 2792 2793 Notes: 2794 Several different monitoring routines may be set by calling 2795 SNESMonitorSet() multiple times; all will be called in the 2796 order in which they were set. 2797 2798 Fortran notes: Only a single monitor function can be set for each SNES object 2799 2800 Level: intermediate 2801 2802 .keywords: SNES, nonlinear, set, monitor 2803 2804 .seealso: SNESMonitorDefault(), SNESMonitorCancel() 2805 @*/ 2806 PetscErrorCode SNESMonitorSet(SNES snes,PetscErrorCode (*monitor)(SNES,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**)) 2807 { 2808 PetscInt i; 2809 PetscErrorCode ierr; 2810 2811 PetscFunctionBegin; 2812 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2813 if (snes->numbermonitors >= MAXSNESMONITORS) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many monitors set"); 2814 for (i=0; i<snes->numbermonitors;i++) { 2815 if (monitor == snes->monitor[i] && monitordestroy == snes->monitordestroy[i] && mctx == snes->monitorcontext[i]) { 2816 if (monitordestroy) { 2817 ierr = (*monitordestroy)(&mctx);CHKERRQ(ierr); 2818 } 2819 PetscFunctionReturn(0); 2820 } 2821 } 2822 snes->monitor[snes->numbermonitors] = monitor; 2823 snes->monitordestroy[snes->numbermonitors] = monitordestroy; 2824 snes->monitorcontext[snes->numbermonitors++] = (void*)mctx; 2825 PetscFunctionReturn(0); 2826 } 2827 2828 #undef __FUNCT__ 2829 #define __FUNCT__ "SNESMonitorCancel" 2830 /*@C 2831 SNESMonitorCancel - Clears all the monitor functions for a SNES object. 2832 2833 Logically Collective on SNES 2834 2835 Input Parameters: 2836 . snes - the SNES context 2837 2838 Options Database Key: 2839 . -snes_monitor_cancel - cancels all monitors that have been hardwired 2840 into a code by calls to SNESMonitorSet(), but does not cancel those 2841 set via the options database 2842 2843 Notes: 2844 There is no way to clear one specific monitor from a SNES object. 2845 2846 Level: intermediate 2847 2848 .keywords: SNES, nonlinear, set, monitor 2849 2850 .seealso: SNESMonitorDefault(), SNESMonitorSet() 2851 @*/ 2852 PetscErrorCode SNESMonitorCancel(SNES snes) 2853 { 2854 PetscErrorCode ierr; 2855 PetscInt i; 2856 2857 PetscFunctionBegin; 2858 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2859 for (i=0; i<snes->numbermonitors; i++) { 2860 if (snes->monitordestroy[i]) { 2861 ierr = (*snes->monitordestroy[i])(&snes->monitorcontext[i]);CHKERRQ(ierr); 2862 } 2863 } 2864 snes->numbermonitors = 0; 2865 PetscFunctionReturn(0); 2866 } 2867 2868 #undef __FUNCT__ 2869 #define __FUNCT__ "SNESSetConvergenceTest" 2870 /*@C 2871 SNESSetConvergenceTest - Sets the function that is to be used 2872 to test for convergence of the nonlinear iterative solution. 2873 2874 Logically Collective on SNES 2875 2876 Input Parameters: 2877 + snes - the SNES context 2878 . func - routine to test for convergence 2879 . cctx - [optional] context for private data for the convergence routine (may be PETSC_NULL) 2880 - destroy - [optional] destructor for the context (may be PETSC_NULL; PETSC_NULL_FUNCTION in Fortran) 2881 2882 Calling sequence of func: 2883 $ PetscErrorCode func (SNES snes,PetscInt it,PetscReal xnorm,PetscReal gnorm,PetscReal f,SNESConvergedReason *reason,void *cctx) 2884 2885 + snes - the SNES context 2886 . it - current iteration (0 is the first and is before any Newton step) 2887 . cctx - [optional] convergence context 2888 . reason - reason for convergence/divergence 2889 . xnorm - 2-norm of current iterate 2890 . gnorm - 2-norm of current step 2891 - f - 2-norm of function 2892 2893 Level: advanced 2894 2895 .keywords: SNES, nonlinear, set, convergence, test 2896 2897 .seealso: SNESDefaultConverged(), SNESSkipConverged() 2898 @*/ 2899 PetscErrorCode SNESSetConvergenceTest(SNES snes,PetscErrorCode (*func)(SNES,PetscInt,PetscReal,PetscReal,PetscReal,SNESConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*)) 2900 { 2901 PetscErrorCode ierr; 2902 2903 PetscFunctionBegin; 2904 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2905 if (!func) func = SNESSkipConverged; 2906 if (snes->ops->convergeddestroy) { 2907 ierr = (*snes->ops->convergeddestroy)(snes->cnvP);CHKERRQ(ierr); 2908 } 2909 snes->ops->converged = func; 2910 snes->ops->convergeddestroy = destroy; 2911 snes->cnvP = cctx; 2912 PetscFunctionReturn(0); 2913 } 2914 2915 #undef __FUNCT__ 2916 #define __FUNCT__ "SNESGetConvergedReason" 2917 /*@ 2918 SNESGetConvergedReason - Gets the reason the SNES iteration was stopped. 2919 2920 Not Collective 2921 2922 Input Parameter: 2923 . snes - the SNES context 2924 2925 Output Parameter: 2926 . reason - negative value indicates diverged, positive value converged, see SNESConvergedReason or the 2927 manual pages for the individual convergence tests for complete lists 2928 2929 Level: intermediate 2930 2931 Notes: Can only be called after the call the SNESSolve() is complete. 2932 2933 .keywords: SNES, nonlinear, set, convergence, test 2934 2935 .seealso: SNESSetConvergenceTest(), SNESConvergedReason 2936 @*/ 2937 PetscErrorCode SNESGetConvergedReason(SNES snes,SNESConvergedReason *reason) 2938 { 2939 PetscFunctionBegin; 2940 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2941 PetscValidPointer(reason,2); 2942 *reason = snes->reason; 2943 PetscFunctionReturn(0); 2944 } 2945 2946 #undef __FUNCT__ 2947 #define __FUNCT__ "SNESSetConvergenceHistory" 2948 /*@ 2949 SNESSetConvergenceHistory - Sets the array used to hold the convergence history. 2950 2951 Logically Collective on SNES 2952 2953 Input Parameters: 2954 + snes - iterative context obtained from SNESCreate() 2955 . a - array to hold history, this array will contain the function norms computed at each step 2956 . its - integer array holds the number of linear iterations for each solve. 2957 . na - size of a and its 2958 - reset - PETSC_TRUE indicates each new nonlinear solve resets the history counter to zero, 2959 else it continues storing new values for new nonlinear solves after the old ones 2960 2961 Notes: 2962 If 'a' and 'its' are PETSC_NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a 2963 default array of length 10000 is allocated. 2964 2965 This routine is useful, e.g., when running a code for purposes 2966 of accurate performance monitoring, when no I/O should be done 2967 during the section of code that is being timed. 2968 2969 Level: intermediate 2970 2971 .keywords: SNES, set, convergence, history 2972 2973 .seealso: SNESGetConvergenceHistory() 2974 2975 @*/ 2976 PetscErrorCode SNESSetConvergenceHistory(SNES snes,PetscReal a[],PetscInt its[],PetscInt na,PetscBool reset) 2977 { 2978 PetscErrorCode ierr; 2979 2980 PetscFunctionBegin; 2981 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 2982 if (na) PetscValidScalarPointer(a,2); 2983 if (its) PetscValidIntPointer(its,3); 2984 if (na == PETSC_DECIDE || na == PETSC_DEFAULT || !a) { 2985 if (na == PETSC_DECIDE || na == PETSC_DEFAULT) na = 1000; 2986 ierr = PetscMalloc(na*sizeof(PetscReal),&a);CHKERRQ(ierr); 2987 ierr = PetscMalloc(na*sizeof(PetscInt),&its);CHKERRQ(ierr); 2988 snes->conv_malloc = PETSC_TRUE; 2989 } 2990 snes->conv_hist = a; 2991 snes->conv_hist_its = its; 2992 snes->conv_hist_max = na; 2993 snes->conv_hist_len = 0; 2994 snes->conv_hist_reset = reset; 2995 PetscFunctionReturn(0); 2996 } 2997 2998 #if defined(PETSC_HAVE_MATLAB_ENGINE) 2999 #include <engine.h> /* MATLAB include file */ 3000 #include <mex.h> /* MATLAB include file */ 3001 EXTERN_C_BEGIN 3002 #undef __FUNCT__ 3003 #define __FUNCT__ "SNESGetConvergenceHistoryMatlab" 3004 mxArray *SNESGetConvergenceHistoryMatlab(SNES snes) 3005 { 3006 mxArray *mat; 3007 PetscInt i; 3008 PetscReal *ar; 3009 3010 PetscFunctionBegin; 3011 mat = mxCreateDoubleMatrix(snes->conv_hist_len,1,mxREAL); 3012 ar = (PetscReal*) mxGetData(mat); 3013 for (i=0; i<snes->conv_hist_len; i++) { 3014 ar[i] = snes->conv_hist[i]; 3015 } 3016 PetscFunctionReturn(mat); 3017 } 3018 EXTERN_C_END 3019 #endif 3020 3021 3022 #undef __FUNCT__ 3023 #define __FUNCT__ "SNESGetConvergenceHistory" 3024 /*@C 3025 SNESGetConvergenceHistory - Gets the array used to hold the convergence history. 3026 3027 Not Collective 3028 3029 Input Parameter: 3030 . snes - iterative context obtained from SNESCreate() 3031 3032 Output Parameters: 3033 . a - array to hold history 3034 . its - integer array holds the number of linear iterations (or 3035 negative if not converged) for each solve. 3036 - na - size of a and its 3037 3038 Notes: 3039 The calling sequence for this routine in Fortran is 3040 $ call SNESGetConvergenceHistory(SNES snes, integer na, integer ierr) 3041 3042 This routine is useful, e.g., when running a code for purposes 3043 of accurate performance monitoring, when no I/O should be done 3044 during the section of code that is being timed. 3045 3046 Level: intermediate 3047 3048 .keywords: SNES, get, convergence, history 3049 3050 .seealso: SNESSetConvergencHistory() 3051 3052 @*/ 3053 PetscErrorCode SNESGetConvergenceHistory(SNES snes,PetscReal *a[],PetscInt *its[],PetscInt *na) 3054 { 3055 PetscFunctionBegin; 3056 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3057 if (a) *a = snes->conv_hist; 3058 if (its) *its = snes->conv_hist_its; 3059 if (na) *na = snes->conv_hist_len; 3060 PetscFunctionReturn(0); 3061 } 3062 3063 #undef __FUNCT__ 3064 #define __FUNCT__ "SNESSetUpdate" 3065 /*@C 3066 SNESSetUpdate - Sets the general-purpose update function called 3067 at the beginning of every iteration of the nonlinear solve. Specifically 3068 it is called just before the Jacobian is "evaluated". 3069 3070 Logically Collective on SNES 3071 3072 Input Parameters: 3073 . snes - The nonlinear solver context 3074 . func - The function 3075 3076 Calling sequence of func: 3077 . func (SNES snes, PetscInt step); 3078 3079 . step - The current step of the iteration 3080 3081 Level: advanced 3082 3083 Note: This is NOT what one uses to update the ghost points before a function evaluation, that should be done at the beginning of your FormFunction() 3084 This is not used by most users. 3085 3086 .keywords: SNES, update 3087 3088 .seealso SNESDefaultUpdate(), SNESSetJacobian(), SNESSolve() 3089 @*/ 3090 PetscErrorCode SNESSetUpdate(SNES snes, PetscErrorCode (*func)(SNES, PetscInt)) 3091 { 3092 PetscFunctionBegin; 3093 PetscValidHeaderSpecific(snes, SNES_CLASSID,1); 3094 snes->ops->update = func; 3095 PetscFunctionReturn(0); 3096 } 3097 3098 #undef __FUNCT__ 3099 #define __FUNCT__ "SNESDefaultUpdate" 3100 /*@ 3101 SNESDefaultUpdate - The default update function which does nothing. 3102 3103 Not collective 3104 3105 Input Parameters: 3106 . snes - The nonlinear solver context 3107 . step - The current step of the iteration 3108 3109 Level: intermediate 3110 3111 .keywords: SNES, update 3112 .seealso SNESSetUpdate(), SNESDefaultRhsBC(), SNESDefaultShortolutionBC() 3113 @*/ 3114 PetscErrorCode SNESDefaultUpdate(SNES snes, PetscInt step) 3115 { 3116 PetscFunctionBegin; 3117 PetscFunctionReturn(0); 3118 } 3119 3120 #undef __FUNCT__ 3121 #define __FUNCT__ "SNESScaleStep_Private" 3122 /* 3123 SNESScaleStep_Private - Scales a step so that its length is less than the 3124 positive parameter delta. 3125 3126 Input Parameters: 3127 + snes - the SNES context 3128 . y - approximate solution of linear system 3129 . fnorm - 2-norm of current function 3130 - delta - trust region size 3131 3132 Output Parameters: 3133 + gpnorm - predicted function norm at the new point, assuming local 3134 linearization. The value is zero if the step lies within the trust 3135 region, and exceeds zero otherwise. 3136 - ynorm - 2-norm of the step 3137 3138 Note: 3139 For non-trust region methods such as SNESLS, the parameter delta 3140 is set to be the maximum allowable step size. 3141 3142 .keywords: SNES, nonlinear, scale, step 3143 */ 3144 PetscErrorCode SNESScaleStep_Private(SNES snes,Vec y,PetscReal *fnorm,PetscReal *delta,PetscReal *gpnorm,PetscReal *ynorm) 3145 { 3146 PetscReal nrm; 3147 PetscScalar cnorm; 3148 PetscErrorCode ierr; 3149 3150 PetscFunctionBegin; 3151 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3152 PetscValidHeaderSpecific(y,VEC_CLASSID,2); 3153 PetscCheckSameComm(snes,1,y,2); 3154 3155 ierr = VecNorm(y,NORM_2,&nrm);CHKERRQ(ierr); 3156 if (nrm > *delta) { 3157 nrm = *delta/nrm; 3158 *gpnorm = (1.0 - nrm)*(*fnorm); 3159 cnorm = nrm; 3160 ierr = VecScale(y,cnorm);CHKERRQ(ierr); 3161 *ynorm = *delta; 3162 } else { 3163 *gpnorm = 0.0; 3164 *ynorm = nrm; 3165 } 3166 PetscFunctionReturn(0); 3167 } 3168 3169 #undef __FUNCT__ 3170 #define __FUNCT__ "SNESSolve" 3171 /*@C 3172 SNESSolve - Solves a nonlinear system F(x) = b. 3173 Call SNESSolve() after calling SNESCreate() and optional routines of the form SNESSetXXX(). 3174 3175 Collective on SNES 3176 3177 Input Parameters: 3178 + snes - the SNES context 3179 . b - the constant part of the equation F(x) = b, or PETSC_NULL to use zero. 3180 - x - the solution vector. 3181 3182 Notes: 3183 The user should initialize the vector,x, with the initial guess 3184 for the nonlinear solve prior to calling SNESSolve. In particular, 3185 to employ an initial guess of zero, the user should explicitly set 3186 this vector to zero by calling VecSet(). 3187 3188 Level: beginner 3189 3190 .keywords: SNES, nonlinear, solve 3191 3192 .seealso: SNESCreate(), SNESDestroy(), SNESSetFunction(), SNESSetJacobian(), SNESSetGridSequence(), SNESGetSolution() 3193 @*/ 3194 PetscErrorCode SNESSolve(SNES snes,Vec b,Vec x) 3195 { 3196 PetscErrorCode ierr; 3197 PetscBool flg; 3198 char filename[PETSC_MAX_PATH_LEN]; 3199 PetscViewer viewer; 3200 PetscInt grid; 3201 Vec xcreated = PETSC_NULL; 3202 3203 PetscFunctionBegin; 3204 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3205 if (x) PetscValidHeaderSpecific(x,VEC_CLASSID,3); 3206 if (x) PetscCheckSameComm(snes,1,x,3); 3207 if (b) PetscValidHeaderSpecific(b,VEC_CLASSID,2); 3208 if (b) PetscCheckSameComm(snes,1,b,2); 3209 3210 if (!x && snes->dm) { 3211 ierr = DMCreateGlobalVector(snes->dm,&xcreated);CHKERRQ(ierr); 3212 x = xcreated; 3213 } 3214 3215 for (grid=0; grid<snes->gridsequence; grid++) {ierr = PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr);} 3216 for (grid=0; grid<snes->gridsequence+1; grid++) { 3217 3218 /* set solution vector */ 3219 if (!grid) {ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr);} 3220 ierr = VecDestroy(&snes->vec_sol);CHKERRQ(ierr); 3221 snes->vec_sol = x; 3222 /* set afine vector if provided */ 3223 if (b) { ierr = PetscObjectReference((PetscObject)b);CHKERRQ(ierr); } 3224 ierr = VecDestroy(&snes->vec_rhs);CHKERRQ(ierr); 3225 snes->vec_rhs = b; 3226 3227 ierr = SNESSetUp(snes);CHKERRQ(ierr); 3228 3229 if (!grid) { 3230 if (snes->ops->computeinitialguess) { 3231 ierr = (*snes->ops->computeinitialguess)(snes,snes->vec_sol,snes->initialguessP);CHKERRQ(ierr); 3232 } else if (snes->dm) { 3233 PetscBool ig; 3234 ierr = DMHasInitialGuess(snes->dm,&ig);CHKERRQ(ierr); 3235 if (ig) { 3236 ierr = DMComputeInitialGuess(snes->dm,snes->vec_sol);CHKERRQ(ierr); 3237 } 3238 } 3239 } 3240 3241 if (snes->conv_hist_reset) snes->conv_hist_len = 0; 3242 snes->nfuncs = 0; snes->linear_its = 0; snes->numFailures = 0; 3243 3244 ierr = PetscLogEventBegin(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3245 ierr = (*snes->ops->solve)(snes);CHKERRQ(ierr); 3246 ierr = PetscLogEventEnd(SNES_Solve,snes,0,0,0);CHKERRQ(ierr); 3247 if (snes->domainerror){ 3248 snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN; 3249 snes->domainerror = PETSC_FALSE; 3250 } 3251 if (!snes->reason) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason"); 3252 3253 ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 3254 if (flg && !PetscPreLoadingOn) { 3255 ierr = PetscViewerASCIIOpen(((PetscObject)snes)->comm,filename,&viewer);CHKERRQ(ierr); 3256 ierr = SNESView(snes,viewer);CHKERRQ(ierr); 3257 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3258 } 3259 3260 flg = PETSC_FALSE; 3261 ierr = PetscOptionsGetBool(((PetscObject)snes)->prefix,"-snes_test_local_min",&flg,PETSC_NULL);CHKERRQ(ierr); 3262 if (flg && !PetscPreLoadingOn) { ierr = SNESTestLocalMin(snes);CHKERRQ(ierr); } 3263 if (snes->printreason) { 3264 ierr = PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3265 if (snes->reason > 0) { 3266 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve converged due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3267 } else { 3268 ierr = PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),"Nonlinear solve did not converge due to %s iterations %D\n",SNESConvergedReasons[snes->reason],snes->iter);CHKERRQ(ierr); 3269 } 3270 ierr = PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm),((PetscObject)snes)->tablevel);CHKERRQ(ierr); 3271 } 3272 3273 flg = PETSC_FALSE; 3274 ierr = PetscOptionsGetString(((PetscObject)snes)->prefix,"-snes_view_solution_vtk",filename,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr); 3275 if (flg) { 3276 PetscViewer viewer; 3277 ierr = PetscViewerCreate(((PetscObject)snes)->comm,&viewer);CHKERRQ(ierr); 3278 ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr); 3279 ierr = PetscViewerFileSetName(viewer,filename);CHKERRQ(ierr); 3280 ierr = VecView(snes->vec_sol,viewer);CHKERRQ(ierr); 3281 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 3282 } 3283 3284 if (snes->errorifnotconverged && snes->reason < 0) SETERRQ(((PetscObject)snes)->comm,PETSC_ERR_NOT_CONVERGED,"SNESSolve has not converged"); 3285 if (grid < snes->gridsequence) { 3286 DM fine; 3287 Vec xnew; 3288 Mat interp; 3289 3290 ierr = DMRefine(snes->dm,((PetscObject)snes)->comm,&fine);CHKERRQ(ierr); 3291 ierr = DMCreateInterpolation(snes->dm,fine,&interp,PETSC_NULL);CHKERRQ(ierr); 3292 ierr = DMCreateGlobalVector(fine,&xnew);CHKERRQ(ierr); 3293 ierr = MatInterpolate(interp,x,xnew);CHKERRQ(ierr); 3294 ierr = MatDestroy(&interp);CHKERRQ(ierr); 3295 x = xnew; 3296 3297 ierr = SNESReset(snes);CHKERRQ(ierr); 3298 ierr = SNESSetDM(snes,fine);CHKERRQ(ierr); 3299 ierr = DMDestroy(&fine);CHKERRQ(ierr); 3300 ierr = PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm));CHKERRQ(ierr); 3301 } 3302 } 3303 ierr = VecDestroy(&xcreated);CHKERRQ(ierr); 3304 PetscFunctionReturn(0); 3305 } 3306 3307 /* --------- Internal routines for SNES Package --------- */ 3308 3309 #undef __FUNCT__ 3310 #define __FUNCT__ "SNESSetType" 3311 /*@C 3312 SNESSetType - Sets the method for the nonlinear solver. 3313 3314 Collective on SNES 3315 3316 Input Parameters: 3317 + snes - the SNES context 3318 - type - a known method 3319 3320 Options Database Key: 3321 . -snes_type <type> - Sets the method; use -help for a list 3322 of available methods (for instance, ls or tr) 3323 3324 Notes: 3325 See "petsc/include/petscsnes.h" for available methods (for instance) 3326 + SNESLS - Newton's method with line search 3327 (systems of nonlinear equations) 3328 . SNESTR - Newton's method with trust region 3329 (systems of nonlinear equations) 3330 3331 Normally, it is best to use the SNESSetFromOptions() command and then 3332 set the SNES solver type from the options database rather than by using 3333 this routine. Using the options database provides the user with 3334 maximum flexibility in evaluating the many nonlinear solvers. 3335 The SNESSetType() routine is provided for those situations where it 3336 is necessary to set the nonlinear solver independently of the command 3337 line or options database. This might be the case, for example, when 3338 the choice of solver changes during the execution of the program, 3339 and the user's application is taking responsibility for choosing the 3340 appropriate method. 3341 3342 Level: intermediate 3343 3344 .keywords: SNES, set, type 3345 3346 .seealso: SNESType, SNESCreate() 3347 3348 @*/ 3349 PetscErrorCode SNESSetType(SNES snes,const SNESType type) 3350 { 3351 PetscErrorCode ierr,(*r)(SNES); 3352 PetscBool match; 3353 3354 PetscFunctionBegin; 3355 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3356 PetscValidCharPointer(type,2); 3357 3358 ierr = PetscTypeCompare((PetscObject)snes,type,&match);CHKERRQ(ierr); 3359 if (match) PetscFunctionReturn(0); 3360 3361 ierr = PetscFListFind(SNESList,((PetscObject)snes)->comm,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 3362 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested SNES type %s",type); 3363 /* Destroy the previous private SNES context */ 3364 if (snes->ops->destroy) { 3365 ierr = (*(snes)->ops->destroy)(snes);CHKERRQ(ierr); 3366 snes->ops->destroy = PETSC_NULL; 3367 } 3368 /* Reinitialize function pointers in SNESOps structure */ 3369 snes->ops->setup = 0; 3370 snes->ops->solve = 0; 3371 snes->ops->view = 0; 3372 snes->ops->setfromoptions = 0; 3373 snes->ops->destroy = 0; 3374 /* Call the SNESCreate_XXX routine for this particular Nonlinear solver */ 3375 snes->setupcalled = PETSC_FALSE; 3376 ierr = PetscObjectChangeTypeName((PetscObject)snes,type);CHKERRQ(ierr); 3377 ierr = (*r)(snes);CHKERRQ(ierr); 3378 #if defined(PETSC_HAVE_AMS) 3379 if (PetscAMSPublishAll) { 3380 ierr = PetscObjectAMSPublish((PetscObject)snes);CHKERRQ(ierr); 3381 } 3382 #endif 3383 PetscFunctionReturn(0); 3384 } 3385 3386 3387 /* --------------------------------------------------------------------- */ 3388 #undef __FUNCT__ 3389 #define __FUNCT__ "SNESRegisterDestroy" 3390 /*@ 3391 SNESRegisterDestroy - Frees the list of nonlinear solvers that were 3392 registered by SNESRegisterDynamic(). 3393 3394 Not Collective 3395 3396 Level: advanced 3397 3398 .keywords: SNES, nonlinear, register, destroy 3399 3400 .seealso: SNESRegisterAll(), SNESRegisterAll() 3401 @*/ 3402 PetscErrorCode SNESRegisterDestroy(void) 3403 { 3404 PetscErrorCode ierr; 3405 3406 PetscFunctionBegin; 3407 ierr = PetscFListDestroy(&SNESList);CHKERRQ(ierr); 3408 SNESRegisterAllCalled = PETSC_FALSE; 3409 PetscFunctionReturn(0); 3410 } 3411 3412 #undef __FUNCT__ 3413 #define __FUNCT__ "SNESGetType" 3414 /*@C 3415 SNESGetType - Gets the SNES method type and name (as a string). 3416 3417 Not Collective 3418 3419 Input Parameter: 3420 . snes - nonlinear solver context 3421 3422 Output Parameter: 3423 . type - SNES method (a character string) 3424 3425 Level: intermediate 3426 3427 .keywords: SNES, nonlinear, get, type, name 3428 @*/ 3429 PetscErrorCode SNESGetType(SNES snes,const SNESType *type) 3430 { 3431 PetscFunctionBegin; 3432 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3433 PetscValidPointer(type,2); 3434 *type = ((PetscObject)snes)->type_name; 3435 PetscFunctionReturn(0); 3436 } 3437 3438 #undef __FUNCT__ 3439 #define __FUNCT__ "SNESGetSolution" 3440 /*@ 3441 SNESGetSolution - Returns the vector where the approximate solution is 3442 stored. This is the fine grid solution when using SNESSetGridSequence(). 3443 3444 Not Collective, but Vec is parallel if SNES is parallel 3445 3446 Input Parameter: 3447 . snes - the SNES context 3448 3449 Output Parameter: 3450 . x - the solution 3451 3452 Level: intermediate 3453 3454 .keywords: SNES, nonlinear, get, solution 3455 3456 .seealso: SNESGetSolutionUpdate(), SNESGetFunction() 3457 @*/ 3458 PetscErrorCode SNESGetSolution(SNES snes,Vec *x) 3459 { 3460 PetscFunctionBegin; 3461 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3462 PetscValidPointer(x,2); 3463 *x = snes->vec_sol; 3464 PetscFunctionReturn(0); 3465 } 3466 3467 #undef __FUNCT__ 3468 #define __FUNCT__ "SNESGetSolutionUpdate" 3469 /*@ 3470 SNESGetSolutionUpdate - Returns the vector where the solution update is 3471 stored. 3472 3473 Not Collective, but Vec is parallel if SNES is parallel 3474 3475 Input Parameter: 3476 . snes - the SNES context 3477 3478 Output Parameter: 3479 . x - the solution update 3480 3481 Level: advanced 3482 3483 .keywords: SNES, nonlinear, get, solution, update 3484 3485 .seealso: SNESGetSolution(), SNESGetFunction() 3486 @*/ 3487 PetscErrorCode SNESGetSolutionUpdate(SNES snes,Vec *x) 3488 { 3489 PetscFunctionBegin; 3490 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3491 PetscValidPointer(x,2); 3492 *x = snes->vec_sol_update; 3493 PetscFunctionReturn(0); 3494 } 3495 3496 #undef __FUNCT__ 3497 #define __FUNCT__ "SNESGetFunction" 3498 /*@C 3499 SNESGetFunction - Returns the vector where the function is stored. 3500 3501 Not Collective, but Vec is parallel if SNES is parallel. Collective if Vec is requested, but has not been created yet. 3502 3503 Input Parameter: 3504 . snes - the SNES context 3505 3506 Output Parameter: 3507 + r - the function (or PETSC_NULL) 3508 . func - the function (or PETSC_NULL) 3509 - ctx - the function context (or PETSC_NULL) 3510 3511 Level: advanced 3512 3513 .keywords: SNES, nonlinear, get, function 3514 3515 .seealso: SNESSetFunction(), SNESGetSolution() 3516 @*/ 3517 PetscErrorCode SNESGetFunction(SNES snes,Vec *r,PetscErrorCode (**func)(SNES,Vec,Vec,void*),void **ctx) 3518 { 3519 PetscErrorCode ierr; 3520 DM dm; 3521 3522 PetscFunctionBegin; 3523 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3524 if (r) { 3525 if (!snes->vec_func) { 3526 if (snes->vec_rhs) { 3527 ierr = VecDuplicate(snes->vec_rhs,&snes->vec_func);CHKERRQ(ierr); 3528 } else if (snes->vec_sol) { 3529 ierr = VecDuplicate(snes->vec_sol,&snes->vec_func);CHKERRQ(ierr); 3530 } else if (snes->dm) { 3531 ierr = DMCreateGlobalVector(snes->dm,&snes->vec_func);CHKERRQ(ierr); 3532 } 3533 } 3534 *r = snes->vec_func; 3535 } 3536 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3537 ierr = DMSNESGetFunction(dm,func,ctx);CHKERRQ(ierr); 3538 PetscFunctionReturn(0); 3539 } 3540 3541 /*@C 3542 SNESGetGS - Returns the GS function and context. 3543 3544 Input Parameter: 3545 . snes - the SNES context 3546 3547 Output Parameter: 3548 + gsfunc - the function (or PETSC_NULL) 3549 - ctx - the function context (or PETSC_NULL) 3550 3551 Level: advanced 3552 3553 .keywords: SNES, nonlinear, get, function 3554 3555 .seealso: SNESSetGS(), SNESGetFunction() 3556 @*/ 3557 3558 #undef __FUNCT__ 3559 #define __FUNCT__ "SNESGetGS" 3560 PetscErrorCode SNESGetGS (SNES snes, PetscErrorCode(**func)(SNES, Vec, Vec, void*), void ** ctx) 3561 { 3562 PetscErrorCode ierr; 3563 DM dm; 3564 3565 PetscFunctionBegin; 3566 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3567 ierr = SNESGetDM(snes,&dm);CHKERRQ(ierr); 3568 ierr = DMSNESGetGS(dm,func,ctx);CHKERRQ(ierr); 3569 PetscFunctionReturn(0); 3570 } 3571 3572 #undef __FUNCT__ 3573 #define __FUNCT__ "SNESSetOptionsPrefix" 3574 /*@C 3575 SNESSetOptionsPrefix - Sets the prefix used for searching for all 3576 SNES options in the database. 3577 3578 Logically Collective on SNES 3579 3580 Input Parameter: 3581 + snes - the SNES context 3582 - prefix - the prefix to prepend to all option names 3583 3584 Notes: 3585 A hyphen (-) must NOT be given at the beginning of the prefix name. 3586 The first character of all runtime options is AUTOMATICALLY the hyphen. 3587 3588 Level: advanced 3589 3590 .keywords: SNES, set, options, prefix, database 3591 3592 .seealso: SNESSetFromOptions() 3593 @*/ 3594 PetscErrorCode SNESSetOptionsPrefix(SNES snes,const char prefix[]) 3595 { 3596 PetscErrorCode ierr; 3597 3598 PetscFunctionBegin; 3599 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3600 ierr = PetscObjectSetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3601 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3602 ierr = KSPSetOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 3603 PetscFunctionReturn(0); 3604 } 3605 3606 #undef __FUNCT__ 3607 #define __FUNCT__ "SNESAppendOptionsPrefix" 3608 /*@C 3609 SNESAppendOptionsPrefix - Appends to the prefix used for searching for all 3610 SNES options in the database. 3611 3612 Logically Collective on SNES 3613 3614 Input Parameters: 3615 + snes - the SNES context 3616 - prefix - the prefix to prepend to all option names 3617 3618 Notes: 3619 A hyphen (-) must NOT be given at the beginning of the prefix name. 3620 The first character of all runtime options is AUTOMATICALLY the hyphen. 3621 3622 Level: advanced 3623 3624 .keywords: SNES, append, options, prefix, database 3625 3626 .seealso: SNESGetOptionsPrefix() 3627 @*/ 3628 PetscErrorCode SNESAppendOptionsPrefix(SNES snes,const char prefix[]) 3629 { 3630 PetscErrorCode ierr; 3631 3632 PetscFunctionBegin; 3633 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3634 ierr = PetscObjectAppendOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3635 if (!snes->ksp) {ierr = SNESGetKSP(snes,&snes->ksp);CHKERRQ(ierr);} 3636 ierr = KSPAppendOptionsPrefix(snes->ksp,prefix);CHKERRQ(ierr); 3637 PetscFunctionReturn(0); 3638 } 3639 3640 #undef __FUNCT__ 3641 #define __FUNCT__ "SNESGetOptionsPrefix" 3642 /*@C 3643 SNESGetOptionsPrefix - Sets the prefix used for searching for all 3644 SNES options in the database. 3645 3646 Not Collective 3647 3648 Input Parameter: 3649 . snes - the SNES context 3650 3651 Output Parameter: 3652 . prefix - pointer to the prefix string used 3653 3654 Notes: On the fortran side, the user should pass in a string 'prefix' of 3655 sufficient length to hold the prefix. 3656 3657 Level: advanced 3658 3659 .keywords: SNES, get, options, prefix, database 3660 3661 .seealso: SNESAppendOptionsPrefix() 3662 @*/ 3663 PetscErrorCode SNESGetOptionsPrefix(SNES snes,const char *prefix[]) 3664 { 3665 PetscErrorCode ierr; 3666 3667 PetscFunctionBegin; 3668 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3669 ierr = PetscObjectGetOptionsPrefix((PetscObject)snes,prefix);CHKERRQ(ierr); 3670 PetscFunctionReturn(0); 3671 } 3672 3673 3674 #undef __FUNCT__ 3675 #define __FUNCT__ "SNESRegister" 3676 /*@C 3677 SNESRegister - See SNESRegisterDynamic() 3678 3679 Level: advanced 3680 @*/ 3681 PetscErrorCode SNESRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(SNES)) 3682 { 3683 char fullname[PETSC_MAX_PATH_LEN]; 3684 PetscErrorCode ierr; 3685 3686 PetscFunctionBegin; 3687 ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); 3688 ierr = PetscFListAdd(&SNESList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 3689 PetscFunctionReturn(0); 3690 } 3691 3692 #undef __FUNCT__ 3693 #define __FUNCT__ "SNESTestLocalMin" 3694 PetscErrorCode SNESTestLocalMin(SNES snes) 3695 { 3696 PetscErrorCode ierr; 3697 PetscInt N,i,j; 3698 Vec u,uh,fh; 3699 PetscScalar value; 3700 PetscReal norm; 3701 3702 PetscFunctionBegin; 3703 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 3704 ierr = VecDuplicate(u,&uh);CHKERRQ(ierr); 3705 ierr = VecDuplicate(u,&fh);CHKERRQ(ierr); 3706 3707 /* currently only works for sequential */ 3708 ierr = PetscPrintf(PETSC_COMM_WORLD,"Testing FormFunction() for local min\n"); 3709 ierr = VecGetSize(u,&N);CHKERRQ(ierr); 3710 for (i=0; i<N; i++) { 3711 ierr = VecCopy(u,uh);CHKERRQ(ierr); 3712 ierr = PetscPrintf(PETSC_COMM_WORLD,"i = %D\n",i);CHKERRQ(ierr); 3713 for (j=-10; j<11; j++) { 3714 value = PetscSign(j)*exp(PetscAbs(j)-10.0); 3715 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 3716 ierr = SNESComputeFunction(snes,uh,fh);CHKERRQ(ierr); 3717 ierr = VecNorm(fh,NORM_2,&norm);CHKERRQ(ierr); 3718 ierr = PetscPrintf(PETSC_COMM_WORLD," j norm %D %18.16e\n",j,norm);CHKERRQ(ierr); 3719 value = -value; 3720 ierr = VecSetValue(uh,i,value,ADD_VALUES);CHKERRQ(ierr); 3721 } 3722 } 3723 ierr = VecDestroy(&uh);CHKERRQ(ierr); 3724 ierr = VecDestroy(&fh);CHKERRQ(ierr); 3725 PetscFunctionReturn(0); 3726 } 3727 3728 #undef __FUNCT__ 3729 #define __FUNCT__ "SNESKSPSetUseEW" 3730 /*@ 3731 SNESKSPSetUseEW - Sets SNES use Eisenstat-Walker method for 3732 computing relative tolerance for linear solvers within an inexact 3733 Newton method. 3734 3735 Logically Collective on SNES 3736 3737 Input Parameters: 3738 + snes - SNES context 3739 - flag - PETSC_TRUE or PETSC_FALSE 3740 3741 Options Database: 3742 + -snes_ksp_ew - use Eisenstat-Walker method for determining linear system convergence 3743 . -snes_ksp_ew_version ver - version of Eisenstat-Walker method 3744 . -snes_ksp_ew_rtol0 <rtol0> - Sets rtol0 3745 . -snes_ksp_ew_rtolmax <rtolmax> - Sets rtolmax 3746 . -snes_ksp_ew_gamma <gamma> - Sets gamma 3747 . -snes_ksp_ew_alpha <alpha> - Sets alpha 3748 . -snes_ksp_ew_alpha2 <alpha2> - Sets alpha2 3749 - -snes_ksp_ew_threshold <threshold> - Sets threshold 3750 3751 Notes: 3752 Currently, the default is to use a constant relative tolerance for 3753 the inner linear solvers. Alternatively, one can use the 3754 Eisenstat-Walker method, where the relative convergence tolerance 3755 is reset at each Newton iteration according progress of the nonlinear 3756 solver. 3757 3758 Level: advanced 3759 3760 Reference: 3761 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 3762 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 3763 3764 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 3765 3766 .seealso: SNESKSPGetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 3767 @*/ 3768 PetscErrorCode SNESKSPSetUseEW(SNES snes,PetscBool flag) 3769 { 3770 PetscFunctionBegin; 3771 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3772 PetscValidLogicalCollectiveBool(snes,flag,2); 3773 snes->ksp_ewconv = flag; 3774 PetscFunctionReturn(0); 3775 } 3776 3777 #undef __FUNCT__ 3778 #define __FUNCT__ "SNESKSPGetUseEW" 3779 /*@ 3780 SNESKSPGetUseEW - Gets if SNES is using Eisenstat-Walker method 3781 for computing relative tolerance for linear solvers within an 3782 inexact Newton method. 3783 3784 Not Collective 3785 3786 Input Parameter: 3787 . snes - SNES context 3788 3789 Output Parameter: 3790 . flag - PETSC_TRUE or PETSC_FALSE 3791 3792 Notes: 3793 Currently, the default is to use a constant relative tolerance for 3794 the inner linear solvers. Alternatively, one can use the 3795 Eisenstat-Walker method, where the relative convergence tolerance 3796 is reset at each Newton iteration according progress of the nonlinear 3797 solver. 3798 3799 Level: advanced 3800 3801 Reference: 3802 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 3803 inexact Newton method", SISC 17 (1), pp.16-32, 1996. 3804 3805 .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 3806 3807 .seealso: SNESKSPSetUseEW(), SNESKSPGetParametersEW(), SNESKSPSetParametersEW() 3808 @*/ 3809 PetscErrorCode SNESKSPGetUseEW(SNES snes, PetscBool *flag) 3810 { 3811 PetscFunctionBegin; 3812 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3813 PetscValidPointer(flag,2); 3814 *flag = snes->ksp_ewconv; 3815 PetscFunctionReturn(0); 3816 } 3817 3818 #undef __FUNCT__ 3819 #define __FUNCT__ "SNESKSPSetParametersEW" 3820 /*@ 3821 SNESKSPSetParametersEW - Sets parameters for Eisenstat-Walker 3822 convergence criteria for the linear solvers within an inexact 3823 Newton method. 3824 3825 Logically Collective on SNES 3826 3827 Input Parameters: 3828 + snes - SNES context 3829 . version - version 1, 2 (default is 2) or 3 3830 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 3831 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 3832 . gamma - multiplicative factor for version 2 rtol computation 3833 (0 <= gamma2 <= 1) 3834 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 3835 . alpha2 - power for safeguard 3836 - threshold - threshold for imposing safeguard (0 < threshold < 1) 3837 3838 Note: 3839 Version 3 was contributed by Luis Chacon, June 2006. 3840 3841 Use PETSC_DEFAULT to retain the default for any of the parameters. 3842 3843 Level: advanced 3844 3845 Reference: 3846 S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 3847 inexact Newton method", Utah State University Math. Stat. Dept. Res. 3848 Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 3849 3850 .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 3851 3852 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPGetParametersEW() 3853 @*/ 3854 PetscErrorCode SNESKSPSetParametersEW(SNES snes,PetscInt version,PetscReal rtol_0,PetscReal rtol_max, 3855 PetscReal gamma,PetscReal alpha,PetscReal alpha2,PetscReal threshold) 3856 { 3857 SNESKSPEW *kctx; 3858 PetscFunctionBegin; 3859 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3860 kctx = (SNESKSPEW*)snes->kspconvctx; 3861 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 3862 PetscValidLogicalCollectiveInt(snes,version,2); 3863 PetscValidLogicalCollectiveReal(snes,rtol_0,3); 3864 PetscValidLogicalCollectiveReal(snes,rtol_max,4); 3865 PetscValidLogicalCollectiveReal(snes,gamma,5); 3866 PetscValidLogicalCollectiveReal(snes,alpha,6); 3867 PetscValidLogicalCollectiveReal(snes,alpha2,7); 3868 PetscValidLogicalCollectiveReal(snes,threshold,8); 3869 3870 if (version != PETSC_DEFAULT) kctx->version = version; 3871 if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 3872 if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 3873 if (gamma != PETSC_DEFAULT) kctx->gamma = gamma; 3874 if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 3875 if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 3876 if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 3877 3878 if (kctx->version < 1 || kctx->version > 3) { 3879 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 and 3 are supported: %D",kctx->version); 3880 } 3881 if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) { 3882 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_0 < 1.0: %G",kctx->rtol_0); 3883 } 3884 if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) { 3885 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= rtol_max (%G) < 1.0\n",kctx->rtol_max); 3886 } 3887 if (kctx->gamma < 0.0 || kctx->gamma > 1.0) { 3888 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 <= gamma (%G) <= 1.0\n",kctx->gamma); 3889 } 3890 if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) { 3891 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"1.0 < alpha (%G) <= 2.0\n",kctx->alpha); 3892 } 3893 if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) { 3894 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"0.0 < threshold (%G) < 1.0\n",kctx->threshold); 3895 } 3896 PetscFunctionReturn(0); 3897 } 3898 3899 #undef __FUNCT__ 3900 #define __FUNCT__ "SNESKSPGetParametersEW" 3901 /*@ 3902 SNESKSPGetParametersEW - Gets parameters for Eisenstat-Walker 3903 convergence criteria for the linear solvers within an inexact 3904 Newton method. 3905 3906 Not Collective 3907 3908 Input Parameters: 3909 snes - SNES context 3910 3911 Output Parameters: 3912 + version - version 1, 2 (default is 2) or 3 3913 . rtol_0 - initial relative tolerance (0 <= rtol_0 < 1) 3914 . rtol_max - maximum relative tolerance (0 <= rtol_max < 1) 3915 . gamma - multiplicative factor for version 2 rtol computation 3916 (0 <= gamma2 <= 1) 3917 . alpha - power for version 2 rtol computation (1 < alpha <= 2) 3918 . alpha2 - power for safeguard 3919 - threshold - threshold for imposing safeguard (0 < threshold < 1) 3920 3921 Level: advanced 3922 3923 .keywords: SNES, KSP, Eisenstat, Walker, get, parameters 3924 3925 .seealso: SNESKSPSetUseEW(), SNESKSPGetUseEW(), SNESKSPSetParametersEW() 3926 @*/ 3927 PetscErrorCode SNESKSPGetParametersEW(SNES snes,PetscInt *version,PetscReal *rtol_0,PetscReal *rtol_max, 3928 PetscReal *gamma,PetscReal *alpha,PetscReal *alpha2,PetscReal *threshold) 3929 { 3930 SNESKSPEW *kctx; 3931 PetscFunctionBegin; 3932 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 3933 kctx = (SNESKSPEW*)snes->kspconvctx; 3934 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context existing"); 3935 if(version) *version = kctx->version; 3936 if(rtol_0) *rtol_0 = kctx->rtol_0; 3937 if(rtol_max) *rtol_max = kctx->rtol_max; 3938 if(gamma) *gamma = kctx->gamma; 3939 if(alpha) *alpha = kctx->alpha; 3940 if(alpha2) *alpha2 = kctx->alpha2; 3941 if(threshold) *threshold = kctx->threshold; 3942 PetscFunctionReturn(0); 3943 } 3944 3945 #undef __FUNCT__ 3946 #define __FUNCT__ "SNESKSPEW_PreSolve" 3947 static PetscErrorCode SNESKSPEW_PreSolve(SNES snes, KSP ksp, Vec b, Vec x) 3948 { 3949 PetscErrorCode ierr; 3950 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 3951 PetscReal rtol=PETSC_DEFAULT,stol; 3952 3953 PetscFunctionBegin; 3954 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 3955 if (!snes->iter) { /* first time in, so use the original user rtol */ 3956 rtol = kctx->rtol_0; 3957 } else { 3958 if (kctx->version == 1) { 3959 rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 3960 if (rtol < 0.0) rtol = -rtol; 3961 stol = pow(kctx->rtol_last,kctx->alpha2); 3962 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 3963 } else if (kctx->version == 2) { 3964 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 3965 stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha); 3966 if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 3967 } else if (kctx->version == 3) {/* contributed by Luis Chacon, June 2006. */ 3968 rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 3969 /* safeguard: avoid sharp decrease of rtol */ 3970 stol = kctx->gamma*pow(kctx->rtol_last,kctx->alpha); 3971 stol = PetscMax(rtol,stol); 3972 rtol = PetscMin(kctx->rtol_0,stol); 3973 /* safeguard: avoid oversolving */ 3974 stol = kctx->gamma*(snes->ttol)/snes->norm; 3975 stol = PetscMax(rtol,stol); 3976 rtol = PetscMin(kctx->rtol_0,stol); 3977 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only versions 1, 2 or 3 are supported: %D",kctx->version); 3978 } 3979 /* safeguard: avoid rtol greater than one */ 3980 rtol = PetscMin(rtol,kctx->rtol_max); 3981 ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); 3982 ierr = PetscInfo3(snes,"iter %D, Eisenstat-Walker (version %D) KSP rtol=%G\n",snes->iter,kctx->version,rtol);CHKERRQ(ierr); 3983 PetscFunctionReturn(0); 3984 } 3985 3986 #undef __FUNCT__ 3987 #define __FUNCT__ "SNESKSPEW_PostSolve" 3988 static PetscErrorCode SNESKSPEW_PostSolve(SNES snes, KSP ksp, Vec b, Vec x) 3989 { 3990 PetscErrorCode ierr; 3991 SNESKSPEW *kctx = (SNESKSPEW*)snes->kspconvctx; 3992 PCSide pcside; 3993 Vec lres; 3994 3995 PetscFunctionBegin; 3996 if (!kctx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"No Eisenstat-Walker context exists"); 3997 ierr = KSPGetTolerances(ksp,&kctx->rtol_last,0,0,0);CHKERRQ(ierr); 3998 ierr = SNESGetFunctionNorm(snes,&kctx->norm_last);CHKERRQ(ierr); 3999 if (kctx->version == 1) { 4000 ierr = KSPGetPCSide(ksp,&pcside);CHKERRQ(ierr); 4001 if (pcside == PC_RIGHT) { /* XXX Should we also test KSP_UNPRECONDITIONED_NORM ? */ 4002 /* KSP residual is true linear residual */ 4003 ierr = KSPGetResidualNorm(ksp,&kctx->lresid_last);CHKERRQ(ierr); 4004 } else { 4005 /* KSP residual is preconditioned residual */ 4006 /* compute true linear residual norm */ 4007 ierr = VecDuplicate(b,&lres);CHKERRQ(ierr); 4008 ierr = MatMult(snes->jacobian,x,lres);CHKERRQ(ierr); 4009 ierr = VecAYPX(lres,-1.0,b);CHKERRQ(ierr); 4010 ierr = VecNorm(lres,NORM_2,&kctx->lresid_last);CHKERRQ(ierr); 4011 ierr = VecDestroy(&lres);CHKERRQ(ierr); 4012 } 4013 } 4014 PetscFunctionReturn(0); 4015 } 4016 4017 #undef __FUNCT__ 4018 #define __FUNCT__ "SNES_KSPSolve" 4019 PetscErrorCode SNES_KSPSolve(SNES snes, KSP ksp, Vec b, Vec x) 4020 { 4021 PetscErrorCode ierr; 4022 4023 PetscFunctionBegin; 4024 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PreSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4025 ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); 4026 if (snes->ksp_ewconv) { ierr = SNESKSPEW_PostSolve(snes,ksp,b,x);CHKERRQ(ierr); } 4027 PetscFunctionReturn(0); 4028 } 4029 4030 #undef __FUNCT__ 4031 #define __FUNCT__ "SNESSetDM" 4032 /*@ 4033 SNESSetDM - Sets the DM that may be used by some preconditioners 4034 4035 Logically Collective on SNES 4036 4037 Input Parameters: 4038 + snes - the preconditioner context 4039 - dm - the dm 4040 4041 Level: intermediate 4042 4043 4044 .seealso: SNESGetDM(), KSPSetDM(), KSPGetDM() 4045 @*/ 4046 PetscErrorCode SNESSetDM(SNES snes,DM dm) 4047 { 4048 PetscErrorCode ierr; 4049 KSP ksp; 4050 SNESDM sdm; 4051 4052 PetscFunctionBegin; 4053 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4054 if (dm) {ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr);} 4055 if (snes->dm) { /* Move the SNESDM context over to the new DM unless the new DM already has one */ 4056 PetscContainer oldcontainer,container; 4057 ierr = PetscObjectQuery((PetscObject)snes->dm,"SNESDM",(PetscObject*)&oldcontainer);CHKERRQ(ierr); 4058 ierr = PetscObjectQuery((PetscObject)dm,"SNESDM",(PetscObject*)&container);CHKERRQ(ierr); 4059 if (oldcontainer && !container) { 4060 ierr = DMSNESCopyContext(snes->dm,dm);CHKERRQ(ierr); 4061 ierr = DMSNESGetContext(snes->dm,&sdm);CHKERRQ(ierr); 4062 if (sdm->originaldm == snes->dm) { /* Grant write privileges to the replacement DM */ 4063 sdm->originaldm = dm; 4064 } 4065 } 4066 ierr = DMDestroy(&snes->dm);CHKERRQ(ierr); 4067 } 4068 snes->dm = dm; 4069 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 4070 ierr = KSPSetDM(ksp,dm);CHKERRQ(ierr); 4071 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr); 4072 if (snes->pc) { 4073 ierr = SNESSetDM(snes->pc, snes->dm);CHKERRQ(ierr); 4074 } 4075 PetscFunctionReturn(0); 4076 } 4077 4078 #undef __FUNCT__ 4079 #define __FUNCT__ "SNESGetDM" 4080 /*@ 4081 SNESGetDM - Gets the DM that may be used by some preconditioners 4082 4083 Not Collective but DM obtained is parallel on SNES 4084 4085 Input Parameter: 4086 . snes - the preconditioner context 4087 4088 Output Parameter: 4089 . dm - the dm 4090 4091 Level: intermediate 4092 4093 4094 .seealso: SNESSetDM(), KSPSetDM(), KSPGetDM() 4095 @*/ 4096 PetscErrorCode SNESGetDM(SNES snes,DM *dm) 4097 { 4098 PetscErrorCode ierr; 4099 4100 PetscFunctionBegin; 4101 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4102 if (!snes->dm) { 4103 ierr = DMShellCreate(((PetscObject)snes)->comm,&snes->dm);CHKERRQ(ierr); 4104 } 4105 *dm = snes->dm; 4106 PetscFunctionReturn(0); 4107 } 4108 4109 #undef __FUNCT__ 4110 #define __FUNCT__ "SNESSetPC" 4111 /*@ 4112 SNESSetPC - Sets the nonlinear preconditioner to be used. 4113 4114 Collective on SNES 4115 4116 Input Parameters: 4117 + snes - iterative context obtained from SNESCreate() 4118 - pc - the preconditioner object 4119 4120 Notes: 4121 Use SNESGetPC() to retrieve the preconditioner context (for example, 4122 to configure it using the API). 4123 4124 Level: developer 4125 4126 .keywords: SNES, set, precondition 4127 .seealso: SNESGetPC() 4128 @*/ 4129 PetscErrorCode SNESSetPC(SNES snes, SNES pc) 4130 { 4131 PetscErrorCode ierr; 4132 4133 PetscFunctionBegin; 4134 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4135 PetscValidHeaderSpecific(pc, SNES_CLASSID, 2); 4136 PetscCheckSameComm(snes, 1, pc, 2); 4137 ierr = PetscObjectReference((PetscObject) pc);CHKERRQ(ierr); 4138 ierr = SNESDestroy(&snes->pc);CHKERRQ(ierr); 4139 snes->pc = pc; 4140 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4141 PetscFunctionReturn(0); 4142 } 4143 4144 #undef __FUNCT__ 4145 #define __FUNCT__ "SNESGetPC" 4146 /*@ 4147 SNESGetPC - Returns a pointer to the nonlinear preconditioning context set with SNESSetPC(). 4148 4149 Not Collective 4150 4151 Input Parameter: 4152 . snes - iterative context obtained from SNESCreate() 4153 4154 Output Parameter: 4155 . pc - preconditioner context 4156 4157 Level: developer 4158 4159 .keywords: SNES, get, preconditioner 4160 .seealso: SNESSetPC() 4161 @*/ 4162 PetscErrorCode SNESGetPC(SNES snes, SNES *pc) 4163 { 4164 PetscErrorCode ierr; 4165 4166 PetscFunctionBegin; 4167 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4168 PetscValidPointer(pc, 2); 4169 if (!snes->pc) { 4170 ierr = SNESCreate(((PetscObject) snes)->comm, &snes->pc);CHKERRQ(ierr); 4171 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->pc, (PetscObject) snes, 1);CHKERRQ(ierr); 4172 ierr = PetscLogObjectParent(snes, snes->pc);CHKERRQ(ierr); 4173 } 4174 *pc = snes->pc; 4175 PetscFunctionReturn(0); 4176 } 4177 4178 #undef __FUNCT__ 4179 #define __FUNCT__ "SNESSetPetscLineSearch" 4180 /*@ 4181 SNESSetPetscLineSearch - Sets the linesearch. 4182 4183 Collective on SNES 4184 4185 Input Parameters: 4186 + snes - iterative context obtained from SNESCreate() 4187 - linesearch - the linesearch object 4188 4189 Notes: 4190 Use SNESGetPetscLineSearch() to retrieve the preconditioner context (for example, 4191 to configure it using the API). 4192 4193 Level: developer 4194 4195 .keywords: SNES, set, linesearch 4196 .seealso: SNESGetPetscLineSearch() 4197 @*/ 4198 PetscErrorCode SNESSetPetscLineSearch(SNES snes, PetscLineSearch linesearch) 4199 { 4200 PetscErrorCode ierr; 4201 4202 PetscFunctionBegin; 4203 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4204 PetscValidHeaderSpecific(linesearch, PETSCLINESEARCH_CLASSID, 2); 4205 PetscCheckSameComm(snes, 1, linesearch, 2); 4206 ierr = PetscObjectReference((PetscObject) linesearch);CHKERRQ(ierr); 4207 ierr = PetscLineSearchDestroy(&snes->linesearch);CHKERRQ(ierr); 4208 snes->linesearch = linesearch; 4209 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4210 PetscFunctionReturn(0); 4211 } 4212 4213 #undef __FUNCT__ 4214 #define __FUNCT__ "SNESGetPetscLineSearch" 4215 /*@C 4216 SNESGetPetscLineSearch - Returns a pointer to the line search context set with SNESSetLineSearch(). 4217 4218 Not Collective 4219 4220 Input Parameter: 4221 . snes - iterative context obtained from SNESCreate() 4222 4223 Output Parameter: 4224 . linesearch - linesearch context 4225 4226 Level: developer 4227 4228 .keywords: SNES, get, linesearch 4229 .seealso: SNESSetPetscLineSearch() 4230 @*/ 4231 PetscErrorCode SNESGetPetscLineSearch(SNES snes, PetscLineSearch *linesearch) 4232 { 4233 PetscErrorCode ierr; 4234 const char *optionsprefix; 4235 4236 PetscFunctionBegin; 4237 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1); 4238 PetscValidPointer(linesearch, 2); 4239 if (!snes->linesearch) { 4240 ierr = SNESGetOptionsPrefix(snes, &optionsprefix);CHKERRQ(ierr); 4241 ierr = PetscLineSearchCreate(((PetscObject) snes)->comm, &snes->linesearch);CHKERRQ(ierr); 4242 ierr = PetscLineSearchSetSNES(snes->linesearch, snes);CHKERRQ(ierr); 4243 ierr = PetscObjectIncrementTabLevel((PetscObject) snes->linesearch, (PetscObject) snes, 1);CHKERRQ(ierr); 4244 ierr = PetscLogObjectParent(snes, snes->linesearch);CHKERRQ(ierr); 4245 } 4246 *linesearch = snes->linesearch; 4247 PetscFunctionReturn(0); 4248 } 4249 4250 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4251 #include <mex.h> 4252 4253 typedef struct {char *funcname; mxArray *ctx;} SNESMatlabContext; 4254 4255 #undef __FUNCT__ 4256 #define __FUNCT__ "SNESComputeFunction_Matlab" 4257 /* 4258 SNESComputeFunction_Matlab - Calls the function that has been set with 4259 SNESSetFunctionMatlab(). 4260 4261 Collective on SNES 4262 4263 Input Parameters: 4264 + snes - the SNES context 4265 - x - input vector 4266 4267 Output Parameter: 4268 . y - function vector, as set by SNESSetFunction() 4269 4270 Notes: 4271 SNESComputeFunction() is typically used within nonlinear solvers 4272 implementations, so most users would not generally call this routine 4273 themselves. 4274 4275 Level: developer 4276 4277 .keywords: SNES, nonlinear, compute, function 4278 4279 .seealso: SNESSetFunction(), SNESGetFunction() 4280 */ 4281 PetscErrorCode SNESComputeFunction_Matlab(SNES snes,Vec x,Vec y, void *ctx) 4282 { 4283 PetscErrorCode ierr; 4284 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4285 int nlhs = 1,nrhs = 5; 4286 mxArray *plhs[1],*prhs[5]; 4287 long long int lx = 0,ly = 0,ls = 0; 4288 4289 PetscFunctionBegin; 4290 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4291 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4292 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 4293 PetscCheckSameComm(snes,1,x,2); 4294 PetscCheckSameComm(snes,1,y,3); 4295 4296 /* call Matlab function in ctx with arguments x and y */ 4297 4298 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4299 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4300 ierr = PetscMemcpy(&ly,&y,sizeof(x));CHKERRQ(ierr); 4301 prhs[0] = mxCreateDoubleScalar((double)ls); 4302 prhs[1] = mxCreateDoubleScalar((double)lx); 4303 prhs[2] = mxCreateDoubleScalar((double)ly); 4304 prhs[3] = mxCreateString(sctx->funcname); 4305 prhs[4] = sctx->ctx; 4306 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeFunctionInternal");CHKERRQ(ierr); 4307 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4308 mxDestroyArray(prhs[0]); 4309 mxDestroyArray(prhs[1]); 4310 mxDestroyArray(prhs[2]); 4311 mxDestroyArray(prhs[3]); 4312 mxDestroyArray(plhs[0]); 4313 PetscFunctionReturn(0); 4314 } 4315 4316 4317 #undef __FUNCT__ 4318 #define __FUNCT__ "SNESSetFunctionMatlab" 4319 /* 4320 SNESSetFunctionMatlab - Sets the function evaluation routine and function 4321 vector for use by the SNES routines in solving systems of nonlinear 4322 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4323 4324 Logically Collective on SNES 4325 4326 Input Parameters: 4327 + snes - the SNES context 4328 . r - vector to store function value 4329 - func - function evaluation routine 4330 4331 Calling sequence of func: 4332 $ func (SNES snes,Vec x,Vec f,void *ctx); 4333 4334 4335 Notes: 4336 The Newton-like methods typically solve linear systems of the form 4337 $ f'(x) x = -f(x), 4338 where f'(x) denotes the Jacobian matrix and f(x) is the function. 4339 4340 Level: beginner 4341 4342 .keywords: SNES, nonlinear, set, function 4343 4344 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4345 */ 4346 PetscErrorCode SNESSetFunctionMatlab(SNES snes,Vec r,const char *func,mxArray *ctx) 4347 { 4348 PetscErrorCode ierr; 4349 SNESMatlabContext *sctx; 4350 4351 PetscFunctionBegin; 4352 /* currently sctx is memory bleed */ 4353 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4354 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4355 /* 4356 This should work, but it doesn't 4357 sctx->ctx = ctx; 4358 mexMakeArrayPersistent(sctx->ctx); 4359 */ 4360 sctx->ctx = mxDuplicateArray(ctx); 4361 ierr = SNESSetFunction(snes,r,SNESComputeFunction_Matlab,sctx);CHKERRQ(ierr); 4362 PetscFunctionReturn(0); 4363 } 4364 4365 #undef __FUNCT__ 4366 #define __FUNCT__ "SNESComputeJacobian_Matlab" 4367 /* 4368 SNESComputeJacobian_Matlab - Calls the function that has been set with 4369 SNESSetJacobianMatlab(). 4370 4371 Collective on SNES 4372 4373 Input Parameters: 4374 + snes - the SNES context 4375 . x - input vector 4376 . A, B - the matrices 4377 - ctx - user context 4378 4379 Output Parameter: 4380 . flag - structure of the matrix 4381 4382 Level: developer 4383 4384 .keywords: SNES, nonlinear, compute, function 4385 4386 .seealso: SNESSetFunction(), SNESGetFunction() 4387 @*/ 4388 PetscErrorCode SNESComputeJacobian_Matlab(SNES snes,Vec x,Mat *A,Mat *B,MatStructure *flag, void *ctx) 4389 { 4390 PetscErrorCode ierr; 4391 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4392 int nlhs = 2,nrhs = 6; 4393 mxArray *plhs[2],*prhs[6]; 4394 long long int lx = 0,lA = 0,ls = 0, lB = 0; 4395 4396 PetscFunctionBegin; 4397 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4398 PetscValidHeaderSpecific(x,VEC_CLASSID,2); 4399 4400 /* call Matlab function in ctx with arguments x and y */ 4401 4402 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4403 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4404 ierr = PetscMemcpy(&lA,A,sizeof(x));CHKERRQ(ierr); 4405 ierr = PetscMemcpy(&lB,B,sizeof(x));CHKERRQ(ierr); 4406 prhs[0] = mxCreateDoubleScalar((double)ls); 4407 prhs[1] = mxCreateDoubleScalar((double)lx); 4408 prhs[2] = mxCreateDoubleScalar((double)lA); 4409 prhs[3] = mxCreateDoubleScalar((double)lB); 4410 prhs[4] = mxCreateString(sctx->funcname); 4411 prhs[5] = sctx->ctx; 4412 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESComputeJacobianInternal");CHKERRQ(ierr); 4413 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4414 *flag = (MatStructure) mxGetScalar(plhs[1]);CHKERRQ(ierr); 4415 mxDestroyArray(prhs[0]); 4416 mxDestroyArray(prhs[1]); 4417 mxDestroyArray(prhs[2]); 4418 mxDestroyArray(prhs[3]); 4419 mxDestroyArray(prhs[4]); 4420 mxDestroyArray(plhs[0]); 4421 mxDestroyArray(plhs[1]); 4422 PetscFunctionReturn(0); 4423 } 4424 4425 4426 #undef __FUNCT__ 4427 #define __FUNCT__ "SNESSetJacobianMatlab" 4428 /* 4429 SNESSetJacobianMatlab - Sets the Jacobian function evaluation routine and two empty Jacobian matrices 4430 vector for use by the SNES routines in solving systems of nonlinear 4431 equations from MATLAB. Here the function is a string containing the name of a MATLAB function 4432 4433 Logically Collective on SNES 4434 4435 Input Parameters: 4436 + snes - the SNES context 4437 . A,B - Jacobian matrices 4438 . func - function evaluation routine 4439 - ctx - user context 4440 4441 Calling sequence of func: 4442 $ flag = func (SNES snes,Vec x,Mat A,Mat B,void *ctx); 4443 4444 4445 Level: developer 4446 4447 .keywords: SNES, nonlinear, set, function 4448 4449 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4450 */ 4451 PetscErrorCode SNESSetJacobianMatlab(SNES snes,Mat A,Mat B,const char *func,mxArray *ctx) 4452 { 4453 PetscErrorCode ierr; 4454 SNESMatlabContext *sctx; 4455 4456 PetscFunctionBegin; 4457 /* currently sctx is memory bleed */ 4458 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4459 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4460 /* 4461 This should work, but it doesn't 4462 sctx->ctx = ctx; 4463 mexMakeArrayPersistent(sctx->ctx); 4464 */ 4465 sctx->ctx = mxDuplicateArray(ctx); 4466 ierr = SNESSetJacobian(snes,A,B,SNESComputeJacobian_Matlab,sctx);CHKERRQ(ierr); 4467 PetscFunctionReturn(0); 4468 } 4469 4470 #undef __FUNCT__ 4471 #define __FUNCT__ "SNESMonitor_Matlab" 4472 /* 4473 SNESMonitor_Matlab - Calls the function that has been set with SNESMonitorSetMatlab(). 4474 4475 Collective on SNES 4476 4477 .seealso: SNESSetFunction(), SNESGetFunction() 4478 @*/ 4479 PetscErrorCode SNESMonitor_Matlab(SNES snes,PetscInt it, PetscReal fnorm, void *ctx) 4480 { 4481 PetscErrorCode ierr; 4482 SNESMatlabContext *sctx = (SNESMatlabContext *)ctx; 4483 int nlhs = 1,nrhs = 6; 4484 mxArray *plhs[1],*prhs[6]; 4485 long long int lx = 0,ls = 0; 4486 Vec x=snes->vec_sol; 4487 4488 PetscFunctionBegin; 4489 PetscValidHeaderSpecific(snes,SNES_CLASSID,1); 4490 4491 ierr = PetscMemcpy(&ls,&snes,sizeof(snes));CHKERRQ(ierr); 4492 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 4493 prhs[0] = mxCreateDoubleScalar((double)ls); 4494 prhs[1] = mxCreateDoubleScalar((double)it); 4495 prhs[2] = mxCreateDoubleScalar((double)fnorm); 4496 prhs[3] = mxCreateDoubleScalar((double)lx); 4497 prhs[4] = mxCreateString(sctx->funcname); 4498 prhs[5] = sctx->ctx; 4499 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscSNESMonitorInternal");CHKERRQ(ierr); 4500 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 4501 mxDestroyArray(prhs[0]); 4502 mxDestroyArray(prhs[1]); 4503 mxDestroyArray(prhs[2]); 4504 mxDestroyArray(prhs[3]); 4505 mxDestroyArray(prhs[4]); 4506 mxDestroyArray(plhs[0]); 4507 PetscFunctionReturn(0); 4508 } 4509 4510 4511 #undef __FUNCT__ 4512 #define __FUNCT__ "SNESMonitorSetMatlab" 4513 /* 4514 SNESMonitorSetMatlab - Sets the monitor function from MATLAB 4515 4516 Level: developer 4517 4518 .keywords: SNES, nonlinear, set, function 4519 4520 .seealso: SNESGetFunction(), SNESComputeFunction(), SNESSetJacobian(), SNESSetFunction() 4521 */ 4522 PetscErrorCode SNESMonitorSetMatlab(SNES snes,const char *func,mxArray *ctx) 4523 { 4524 PetscErrorCode ierr; 4525 SNESMatlabContext *sctx; 4526 4527 PetscFunctionBegin; 4528 /* currently sctx is memory bleed */ 4529 ierr = PetscMalloc(sizeof(SNESMatlabContext),&sctx);CHKERRQ(ierr); 4530 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 4531 /* 4532 This should work, but it doesn't 4533 sctx->ctx = ctx; 4534 mexMakeArrayPersistent(sctx->ctx); 4535 */ 4536 sctx->ctx = mxDuplicateArray(ctx); 4537 ierr = SNESMonitorSet(snes,SNESMonitor_Matlab,sctx,PETSC_NULL);CHKERRQ(ierr); 4538 PetscFunctionReturn(0); 4539 } 4540 4541 #endif 4542