1 static char help[] = "Newton methods to solve u'' + u^{2} = f in parallel.\n\
2 This example employs a user-defined monitoring routine and optionally a user-defined\n\
3 routine to check candidate iterates produced by line search routines.\n\
4 The command line options include:\n\
5 -pre_check_iterates : activate checking of iterates\n\
6 -post_check_iterates : activate checking of iterates\n\
7 -check_tol <tol>: set tolerance for iterate checking\n\
8 -user_precond : activate a (trivial) user-defined preconditioner\n\n";
9
10 /*
11 Include "petscdm.h" so that we can use data management objects (DMs)
12 Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
13 Include "petscsnes.h" so that we can use SNES solvers. Note that this
14 file automatically includes:
15 petscsys.h - base PETSc routines
16 petscvec.h - vectors
17 petscmat.h - matrices
18 petscis.h - index sets
19 petscksp.h - Krylov subspace methods
20 petscviewer.h - viewers
21 petscpc.h - preconditioners
22 petscksp.h - linear solvers
23 */
24
25 #include <petscdm.h>
26 #include <petscdmda.h>
27 #include <petscsnes.h>
28
29 /*
30 User-defined routines.
31 */
32 PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
33 PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
34 PetscErrorCode FormInitialGuess(Vec);
35 PetscErrorCode Monitor(SNES, PetscInt, PetscReal, void *);
36 PetscErrorCode PreCheck(SNESLineSearch, Vec, Vec, PetscBool *, void *);
37 PetscErrorCode PostCheck(SNESLineSearch, Vec, Vec, Vec, PetscBool *, PetscBool *, void *);
38 PetscErrorCode PostSetSubKSP(SNESLineSearch, Vec, Vec, Vec, PetscBool *, PetscBool *, void *);
39 PetscErrorCode MatrixFreePreconditioner(PC, Vec, Vec);
40
41 /*
42 User-defined application context
43 */
44 typedef struct {
45 DM da; /* distributed array */
46 Vec F; /* right-hand side of PDE */
47 PetscMPIInt rank; /* rank of processor */
48 PetscMPIInt size; /* size of communicator */
49 PetscReal h; /* mesh spacing */
50 PetscBool sjerr; /* if or not to test jacobian domain error */
51 } ApplicationCtx;
52
53 /*
54 User-defined context for monitoring
55 */
56 typedef struct {
57 PetscViewer viewer;
58 } MonitorCtx;
59
60 /*
61 User-defined context for checking candidate iterates that are
62 determined by line search methods
63 */
64 typedef struct {
65 Vec last_step; /* previous iterate */
66 PetscReal tolerance; /* tolerance for changes between successive iterates */
67 ApplicationCtx *user;
68 } StepCheckCtx;
69
70 typedef struct {
71 PetscInt its0; /* num of previous outer KSP iterations */
72 } SetSubKSPCtx;
73
main(int argc,char ** argv)74 int main(int argc, char **argv)
75 {
76 SNES snes; /* SNES context */
77 SNESLineSearch linesearch; /* SNESLineSearch context */
78 Mat J; /* Jacobian matrix */
79 ApplicationCtx ctx; /* user-defined context */
80 Vec x, r, U, F; /* vectors */
81 MonitorCtx monP; /* monitoring context */
82 StepCheckCtx checkP; /* step-checking context */
83 SetSubKSPCtx checkP1;
84 PetscBool pre_check, post_check, post_setsubksp; /* flag indicating whether we're checking candidate iterates */
85 PetscScalar xp, *FF, *UU, none = -1.0;
86 PetscInt its, N = 5, i, maxit, maxf, xs, xm;
87 PetscReal abstol, rtol, stol, norm;
88 PetscBool flg, viewinitial = PETSC_FALSE;
89
90 PetscFunctionBeginUser;
91 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
92 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank));
93 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size));
94 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &N, NULL));
95 ctx.h = 1.0 / (N - 1);
96 ctx.sjerr = PETSC_FALSE;
97 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_jacobian_domain_error", &ctx.sjerr, NULL));
98 PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_initial", &viewinitial, NULL));
99
100 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101 Create nonlinear solver context
102 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103
104 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
105
106 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107 Create vector data structures; set function evaluation routine
108 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109
110 /*
111 Create distributed array (DMDA) to manage parallel grid and vectors
112 */
113 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, N, 1, 1, NULL, &ctx.da));
114 PetscCall(DMSetFromOptions(ctx.da));
115 PetscCall(DMSetUp(ctx.da));
116
117 /*
118 Extract global and local vectors from DMDA; then duplicate for remaining
119 vectors that are the same types
120 */
121 PetscCall(DMCreateGlobalVector(ctx.da, &x));
122 PetscCall(VecDuplicate(x, &r));
123 PetscCall(VecDuplicate(x, &F));
124 ctx.F = F;
125 PetscCall(VecDuplicate(x, &U));
126
127 /*
128 Set function evaluation routine and vector. Whenever the nonlinear
129 solver needs to compute the nonlinear function, it will call this
130 routine.
131 - Note that the final routine argument is the user-defined
132 context that provides application-specific data for the
133 function evaluation routine.
134 */
135 PetscCall(SNESSetFunction(snes, r, FormFunction, &ctx));
136
137 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138 Create matrix data structure; set Jacobian evaluation routine
139 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140
141 PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
142 PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, N, N));
143 PetscCall(MatSetFromOptions(J));
144 PetscCall(MatSeqAIJSetPreallocation(J, 3, NULL));
145 PetscCall(MatMPIAIJSetPreallocation(J, 3, NULL, 3, NULL));
146
147 /*
148 Set Jacobian matrix data structure and default Jacobian evaluation
149 routine. Whenever the nonlinear solver needs to compute the
150 Jacobian matrix, it will call this routine.
151 - Note that the final routine argument is the user-defined
152 context that provides application-specific data for the
153 Jacobian evaluation routine.
154 */
155 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, &ctx));
156
157 /*
158 Optionally allow user-provided preconditioner
159 */
160 PetscCall(PetscOptionsHasName(NULL, NULL, "-user_precond", &flg));
161 if (flg) {
162 KSP ksp;
163 PC pc;
164 PetscCall(SNESGetKSP(snes, &ksp));
165 PetscCall(KSPGetPC(ksp, &pc));
166 PetscCall(PCSetType(pc, PCSHELL));
167 PetscCall(PCShellSetApply(pc, MatrixFreePreconditioner));
168 }
169
170 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171 Customize nonlinear solver; set runtime options
172 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
173
174 /*
175 Set an optional user-defined monitoring routine
176 */
177 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, 0, 0, 0, 400, 400, &monP.viewer));
178 PetscCall(SNESMonitorSet(snes, Monitor, &monP, 0));
179
180 /*
181 Set names for some vectors to facilitate monitoring (optional)
182 */
183 PetscCall(PetscObjectSetName((PetscObject)x, "Approximate Solution"));
184 PetscCall(PetscObjectSetName((PetscObject)U, "Exact Solution"));
185
186 /*
187 Set SNES/KSP/KSP/PC runtime options, e.g.,
188 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
189 */
190 PetscCall(SNESSetFromOptions(snes));
191
192 /*
193 Set an optional user-defined routine to check the validity of candidate
194 iterates that are determined by line search methods
195 */
196 PetscCall(SNESGetLineSearch(snes, &linesearch));
197 PetscCall(PetscOptionsHasName(NULL, NULL, "-post_check_iterates", &post_check));
198
199 if (post_check) {
200 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Activating post step checking routine\n"));
201 PetscCall(SNESLineSearchSetPostCheck(linesearch, PostCheck, &checkP));
202 PetscCall(VecDuplicate(x, &checkP.last_step));
203
204 checkP.tolerance = 1.0;
205 checkP.user = &ctx;
206
207 PetscCall(PetscOptionsGetReal(NULL, NULL, "-check_tol", &checkP.tolerance, NULL));
208 }
209
210 PetscCall(PetscOptionsHasName(NULL, NULL, "-post_setsubksp", &post_setsubksp));
211 if (post_setsubksp) {
212 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Activating post setsubksp\n"));
213 PetscCall(SNESLineSearchSetPostCheck(linesearch, PostSetSubKSP, &checkP1));
214 }
215
216 PetscCall(PetscOptionsHasName(NULL, NULL, "-pre_check_iterates", &pre_check));
217 if (pre_check) {
218 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Activating pre step checking routine\n"));
219 PetscCall(SNESLineSearchSetPreCheck(linesearch, PreCheck, &checkP));
220 }
221
222 /*
223 Print parameters used for convergence testing (optional) ... just
224 to demonstrate this routine; this information is also printed with
225 the option -snes_view
226 */
227 PetscCall(SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf));
228 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "atol=%g, rtol=%g, stol=%g, maxit=%" PetscInt_FMT ", maxf=%" PetscInt_FMT "\n", (double)abstol, (double)rtol, (double)stol, maxit, maxf));
229
230 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231 Initialize application:
232 Store right-hand side of PDE and exact solution
233 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
234
235 /*
236 Get local grid boundaries (for 1-dimensional DMDA):
237 xs, xm - starting grid index, width of local grid (no ghost points)
238 */
239 PetscCall(DMDAGetCorners(ctx.da, &xs, NULL, NULL, &xm, NULL, NULL));
240
241 /*
242 Get pointers to vector data
243 */
244 PetscCall(DMDAVecGetArray(ctx.da, F, &FF));
245 PetscCall(DMDAVecGetArray(ctx.da, U, &UU));
246
247 /*
248 Compute local vector entries
249 */
250 xp = ctx.h * xs;
251 for (i = xs; i < xs + xm; i++) {
252 FF[i] = 6.0 * xp + PetscPowScalar(xp + 1.e-12, 6.0); /* +1.e-12 is to prevent 0^6 */
253 UU[i] = xp * xp * xp;
254 xp += ctx.h;
255 }
256
257 /*
258 Restore vectors
259 */
260 PetscCall(DMDAVecRestoreArray(ctx.da, F, &FF));
261 PetscCall(DMDAVecRestoreArray(ctx.da, U, &UU));
262 if (viewinitial) {
263 PetscCall(VecView(U, 0));
264 PetscCall(VecView(F, 0));
265 }
266
267 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
268 Evaluate initial guess; then solve nonlinear system
269 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
270
271 /*
272 Note: The user should initialize the vector, x, with the initial guess
273 for the nonlinear solver prior to calling SNESSolve(). In particular,
274 to employ an initial guess of zero, the user should explicitly set
275 this vector to zero by calling VecSet().
276 */
277 PetscCall(FormInitialGuess(x));
278 PetscCall(SNESSolve(snes, NULL, x));
279 PetscCall(SNESGetIterationNumber(snes, &its));
280 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
281
282 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
283 Check solution and clean up
284 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
285
286 /*
287 Check the error
288 */
289 PetscCall(VecAXPY(x, none, U));
290 PetscCall(VecNorm(x, NORM_2, &norm));
291 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Iterations %" PetscInt_FMT "\n", (double)norm, its));
292 if (ctx.sjerr) {
293 SNESType snestype;
294 PetscCall(SNESGetType(snes, &snestype));
295 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SNES Type %s\n", snestype));
296 }
297
298 /*
299 Free work space. All PETSc objects should be destroyed when they
300 are no longer needed.
301 */
302 PetscCall(PetscViewerDestroy(&monP.viewer));
303 if (post_check) PetscCall(VecDestroy(&checkP.last_step));
304 PetscCall(VecDestroy(&x));
305 PetscCall(VecDestroy(&r));
306 PetscCall(VecDestroy(&U));
307 PetscCall(VecDestroy(&F));
308 PetscCall(MatDestroy(&J));
309 PetscCall(SNESDestroy(&snes));
310 PetscCall(DMDestroy(&ctx.da));
311 PetscCall(PetscFinalize());
312 return 0;
313 }
314
315 /* ------------------------------------------------------------------- */
316 /*
317 FormInitialGuess - Computes initial guess.
318
319 Input/Output Parameter:
320 . x - the solution vector
321 */
FormInitialGuess(Vec x)322 PetscErrorCode FormInitialGuess(Vec x)
323 {
324 PetscScalar pfive = .50;
325
326 PetscFunctionBeginUser;
327 PetscCall(VecSet(x, pfive));
328 PetscFunctionReturn(PETSC_SUCCESS);
329 }
330
331 /* ------------------------------------------------------------------- */
332 /*
333 FormFunction - Evaluates nonlinear function, F(x).
334
335 Input Parameters:
336 . snes - the SNES context
337 . x - input vector
338 . ctx - optional user-defined context, as set by SNESSetFunction()
339
340 Output Parameter:
341 . f - function vector
342
343 Note:
344 The user-defined context can contain any application-specific
345 data needed for the function evaluation.
346 */
FormFunction(SNES snes,Vec x,Vec f,PetscCtx ctx)347 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, PetscCtx ctx)
348 {
349 ApplicationCtx *user = (ApplicationCtx *)ctx;
350 DM da = user->da;
351 PetscScalar *ff, d;
352 const PetscScalar *xx, *FF;
353 PetscInt i, M, xs, xm;
354 Vec xlocal;
355
356 PetscFunctionBeginUser;
357 PetscCall(DMGetLocalVector(da, &xlocal));
358 /*
359 Scatter ghost points to local vector, using the 2-step process
360 DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
361 By placing code between these two statements, computations can
362 be done while messages are in transition.
363 */
364 PetscCall(DMGlobalToLocalBegin(da, x, INSERT_VALUES, xlocal));
365 PetscCall(DMGlobalToLocalEnd(da, x, INSERT_VALUES, xlocal));
366
367 /*
368 Get pointers to vector data.
369 - The vector xlocal includes ghost point; the vectors x and f do
370 NOT include ghost points.
371 - Using DMDAVecGetArray() allows accessing the values using global ordering
372 */
373 PetscCall(DMDAVecGetArrayRead(da, xlocal, (void *)&xx));
374 PetscCall(DMDAVecGetArray(da, f, &ff));
375 PetscCall(DMDAVecGetArrayRead(da, user->F, (void *)&FF));
376
377 /*
378 Get local grid boundaries (for 1-dimensional DMDA):
379 xs, xm - starting grid index, width of local grid (no ghost points)
380 */
381 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
382 PetscCall(DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
383
384 /*
385 Set function values for boundary points; define local interior grid point range:
386 xsi - starting interior grid index
387 xei - ending interior grid index
388 */
389 if (xs == 0) { /* left boundary */
390 ff[0] = xx[0];
391 xs++;
392 xm--;
393 }
394 if (xs + xm == M) { /* right boundary */
395 ff[xs + xm - 1] = xx[xs + xm - 1] - 1.0;
396 xm--;
397 }
398
399 /*
400 Compute function over locally owned part of the grid (interior points only)
401 */
402 d = 1.0 / (user->h * user->h);
403 for (i = xs; i < xs + xm; i++) ff[i] = d * (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) + xx[i] * xx[i] - FF[i];
404
405 /*
406 Restore vectors
407 */
408 PetscCall(DMDAVecRestoreArrayRead(da, xlocal, (void *)&xx));
409 PetscCall(DMDAVecRestoreArray(da, f, &ff));
410 PetscCall(DMDAVecRestoreArrayRead(da, user->F, (void *)&FF));
411 PetscCall(DMRestoreLocalVector(da, &xlocal));
412 PetscFunctionReturn(PETSC_SUCCESS);
413 }
414
415 /* ------------------------------------------------------------------- */
416 /*
417 FormJacobian - Evaluates Jacobian matrix.
418
419 Input Parameters:
420 . snes - the SNES context
421 . x - input vector
422 . dummy - optional user-defined context (not used here)
423
424 Output Parameters:
425 . jac - Jacobian matrix
426 . B - optionally different matrix used to construct the preconditioner
427
428 */
FormJacobian(SNES snes,Vec x,Mat jac,Mat B,PetscCtx ctx)429 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, PetscCtx ctx)
430 {
431 ApplicationCtx *user = (ApplicationCtx *)ctx;
432 PetscScalar *xx, d, A[3];
433 PetscInt i, j[3], M, xs, xm;
434 DM da = user->da;
435
436 PetscFunctionBeginUser;
437 if (user->sjerr) {
438 PetscCall(SNESSetJacobianDomainError(snes));
439 PetscFunctionReturn(PETSC_SUCCESS);
440 }
441 /*
442 Get pointer to vector data
443 */
444 PetscCall(DMDAVecGetArrayRead(da, x, &xx));
445 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
446
447 /*
448 Get range of locally owned matrix
449 */
450 PetscCall(DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
451
452 /*
453 Determine starting and ending local indices for interior grid points.
454 Set Jacobian entries for boundary points.
455 */
456
457 if (xs == 0) { /* left boundary */
458 i = 0;
459 A[0] = 1.0;
460
461 PetscCall(MatSetValues(jac, 1, &i, 1, &i, A, INSERT_VALUES));
462 xs++;
463 xm--;
464 }
465 if (xs + xm == M) { /* right boundary */
466 i = M - 1;
467 A[0] = 1.0;
468 PetscCall(MatSetValues(jac, 1, &i, 1, &i, A, INSERT_VALUES));
469 xm--;
470 }
471
472 /*
473 Interior grid points
474 - Note that in this case we set all elements for a particular
475 row at once.
476 */
477 d = 1.0 / (user->h * user->h);
478 for (i = xs; i < xs + xm; i++) {
479 j[0] = i - 1;
480 j[1] = i;
481 j[2] = i + 1;
482 A[0] = A[2] = d;
483 A[1] = -2.0 * d + 2.0 * xx[i];
484 PetscCall(MatSetValues(jac, 1, &i, 3, j, A, INSERT_VALUES));
485 }
486
487 /*
488 Assemble matrix, using the 2-step process:
489 MatAssemblyBegin(), MatAssemblyEnd().
490 By placing code between these two statements, computations can be
491 done while messages are in transition.
492
493 Also, restore vector.
494 */
495
496 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
497 PetscCall(DMDAVecRestoreArrayRead(da, x, &xx));
498 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
499 PetscFunctionReturn(PETSC_SUCCESS);
500 }
501
502 /* ------------------------------------------------------------------- */
503 /*
504 Monitor - Optional user-defined monitoring routine that views the
505 current iterate with an x-window plot. Set by SNESMonitorSet().
506
507 Input Parameters:
508 snes - the SNES context
509 its - iteration number
510 norm - 2-norm function value (may be estimated)
511 ctx - optional user-defined context for private data for the
512 monitor routine, as set by SNESMonitorSet()
513
514 Note:
515 See the manpage for PetscViewerDrawOpen() for useful runtime options,
516 such as -nox to deactivate all x-window output.
517 */
Monitor(SNES snes,PetscInt its,PetscReal fnorm,PetscCtx ctx)518 PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, PetscCtx ctx)
519 {
520 MonitorCtx *monP = (MonitorCtx *)ctx;
521 Vec x;
522
523 PetscFunctionBeginUser;
524 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "iter = %" PetscInt_FMT ",SNES Function norm %g\n", its, (double)fnorm));
525 PetscCall(SNESGetSolution(snes, &x));
526 PetscCall(VecView(x, monP->viewer));
527 PetscFunctionReturn(PETSC_SUCCESS);
528 }
529
530 /* ------------------------------------------------------------------- */
531 /*
532 PreCheck - Optional user-defined routine that checks the validity of
533 candidate steps of a line search method. Set by SNESLineSearchSetPreCheck().
534
535 Input Parameters:
536 snes - the SNES context
537 xcurrent - current solution
538 y - search direction and length
539
540 Output Parameters:
541 y - proposed step (search direction and length) (possibly changed)
542 changed_y - tells if the step has changed or not
543 */
PreCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y,PetscBool * changed_y,PetscCtx ctx)544 PetscErrorCode PreCheck(SNESLineSearch linesearch, Vec xcurrent, Vec y, PetscBool *changed_y, PetscCtx ctx)
545 {
546 PetscFunctionBeginUser;
547 *changed_y = PETSC_FALSE;
548 PetscFunctionReturn(PETSC_SUCCESS);
549 }
550
551 /* ------------------------------------------------------------------- */
552 /*
553 PostCheck - Optional user-defined routine that checks the validity of
554 candidate steps of a line search method. Set by SNESLineSearchSetPostCheck().
555
556 Input Parameters:
557 snes - the SNES context
558 ctx - optional user-defined context for private data for the
559 monitor routine, as set by SNESLineSearchSetPostCheck()
560 xcurrent - current solution
561 y - search direction and length
562 x - the new candidate iterate
563
564 Output Parameters:
565 y - proposed step (search direction and length) (possibly changed)
566 x - current iterate (possibly modified)
567
568 */
PostCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool * changed_y,PetscBool * changed_x,PetscCtx ctx)569 PetscErrorCode PostCheck(SNESLineSearch linesearch, Vec xcurrent, Vec y, Vec x, PetscBool *changed_y, PetscBool *changed_x, PetscCtx ctx)
570 {
571 PetscInt i, iter, xs, xm;
572 StepCheckCtx *check;
573 ApplicationCtx *user;
574 PetscScalar *xa, *xa_last, tmp;
575 PetscReal rdiff;
576 DM da;
577 SNES snes;
578
579 PetscFunctionBeginUser;
580 *changed_x = PETSC_FALSE;
581 *changed_y = PETSC_FALSE;
582
583 PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
584 check = (StepCheckCtx *)ctx;
585 user = check->user;
586 PetscCall(SNESGetIterationNumber(snes, &iter));
587
588 /* iteration 1 indicates we are working on the second iteration */
589 if (iter > 0) {
590 da = user->da;
591 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Checking candidate step at iteration %" PetscInt_FMT " with tolerance %g\n", iter, (double)check->tolerance));
592
593 /* Access local array data */
594 PetscCall(DMDAVecGetArray(da, check->last_step, &xa_last));
595 PetscCall(DMDAVecGetArray(da, x, &xa));
596 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
597
598 /*
599 If we fail the user-defined check for validity of the candidate iterate,
600 then modify the iterate as we like. (Note that the particular modification
601 below is intended simply to demonstrate how to manipulate this data, not
602 as a meaningful or appropriate choice.)
603 */
604 for (i = xs; i < xs + xm; i++) {
605 if (!PetscAbsScalar(xa[i])) rdiff = 2 * check->tolerance;
606 else rdiff = PetscAbsScalar((xa[i] - xa_last[i]) / xa[i]);
607 if (rdiff > check->tolerance) {
608 tmp = xa[i];
609 xa[i] = .5 * (xa[i] + xa_last[i]);
610 *changed_x = PETSC_TRUE;
611 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Altering entry %" PetscInt_FMT ": x=%g, x_last=%g, diff=%g, x_new=%g\n", i, (double)PetscAbsScalar(tmp), (double)PetscAbsScalar(xa_last[i]), (double)rdiff, (double)PetscAbsScalar(xa[i])));
612 }
613 }
614 PetscCall(DMDAVecRestoreArray(da, check->last_step, &xa_last));
615 PetscCall(DMDAVecRestoreArray(da, x, &xa));
616 }
617 PetscCall(VecCopy(x, check->last_step));
618 PetscFunctionReturn(PETSC_SUCCESS);
619 }
620
621 /* ------------------------------------------------------------------- */
622 /*
623 PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
624 e.g,
625 mpiexec -n 8 ./ex3 -nox -n 10000 -ksp_type fgmres -pc_type bjacobi -pc_bjacobi_blocks 4 -sub_ksp_type gmres -sub_ksp_max_it 3 -post_setsubksp -sub_ksp_rtol 1.e-16
626 Set by SNESLineSearchSetPostCheck().
627
628 Input Parameters:
629 linesearch - the LineSearch context
630 xcurrent - current solution
631 y - search direction and length
632 x - the new candidate iterate
633
634 Output Parameters:
635 y - proposed step (search direction and length) (possibly changed)
636 x - current iterate (possibly modified)
637
638 */
PostSetSubKSP(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool * changed_y,PetscBool * changed_x,PetscCtx ctx)639 PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch, Vec xcurrent, Vec y, Vec x, PetscBool *changed_y, PetscBool *changed_x, PetscCtx ctx)
640 {
641 SetSubKSPCtx *check;
642 PetscInt iter, its, sub_its, maxit;
643 KSP ksp, sub_ksp, *sub_ksps;
644 PC pc;
645 PetscReal ksp_ratio;
646 SNES snes;
647
648 PetscFunctionBeginUser;
649 PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
650 check = (SetSubKSPCtx *)ctx;
651 PetscCall(SNESGetIterationNumber(snes, &iter));
652 PetscCall(SNESGetKSP(snes, &ksp));
653 PetscCall(KSPGetPC(ksp, &pc));
654 PetscCall(PCBJacobiGetSubKSP(pc, NULL, NULL, &sub_ksps));
655 sub_ksp = sub_ksps[0];
656 PetscCall(KSPGetIterationNumber(ksp, &its)); /* outer KSP iteration number */
657 PetscCall(KSPGetIterationNumber(sub_ksp, &sub_its)); /* inner KSP iteration number */
658
659 if (iter) {
660 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " ...PostCheck snes iteration %" PetscInt_FMT ", ksp_it %" PetscInt_FMT " %" PetscInt_FMT ", subksp_it %" PetscInt_FMT "\n", iter, check->its0, its, sub_its));
661 ksp_ratio = (PetscReal)its / check->its0;
662 maxit = (PetscInt)(ksp_ratio * sub_its + 0.5);
663 if (maxit < 2) maxit = 2;
664 PetscCall(KSPSetTolerances(sub_ksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, maxit));
665 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " ...ksp_ratio %g, new maxit %" PetscInt_FMT "\n\n", (double)ksp_ratio, maxit));
666 }
667 check->its0 = its; /* save current outer KSP iteration number */
668 PetscFunctionReturn(PETSC_SUCCESS);
669 }
670
671 /* ------------------------------------------------------------------- */
672 /*
673 MatrixFreePreconditioner - This routine demonstrates the use of a
674 user-provided preconditioner. This code implements just the null
675 preconditioner, which of course is not recommended for general use.
676
677 Input Parameters:
678 + pc - preconditioner
679 - x - input vector
680
681 Output Parameter:
682 . y - preconditioned vector
683 */
MatrixFreePreconditioner(PC pc,Vec x,Vec y)684 PetscErrorCode MatrixFreePreconditioner(PC pc, Vec x, Vec y)
685 {
686 PetscFunctionBeginUser;
687 PetscCall(VecCopy(x, y));
688 PetscFunctionReturn(PETSC_SUCCESS);
689 }
690
691 /*TEST
692
693 test:
694 args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
695
696 test:
697 suffix: 2
698 nsize: 3
699 args: -nox -pc_type asm -mat_type mpiaij -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
700
701 test:
702 suffix: 3
703 nsize: 2
704 args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
705
706 test:
707 suffix: 4
708 args: -nox -pre_check_iterates -post_check_iterates
709
710 test:
711 suffix: 5
712 requires: double !complex !single
713 nsize: 2
714 args: -nox -snes_test_jacobian -snes_test_jacobian_view
715
716 test:
717 suffix: 6
718 requires: double !complex !single
719 nsize: 4
720 args: -test_jacobian_domain_error -snes_converged_reason -snes_check_jacobian_domain_error 1
721
722 test:
723 suffix: 7
724 requires: double !complex !single
725 nsize: 4
726 args: -test_jacobian_domain_error -snes_converged_reason -snes_type newtontr -snes_check_jacobian_domain_error 1
727
728 test:
729 suffix: 8
730 requires: double !complex !single
731 nsize: 4
732 args: -test_jacobian_domain_error -snes_converged_reason -snes_type vinewtonrsls -snes_check_jacobian_domain_error 1
733
734 test:
735 suffix: 9
736 requires: double !complex !single
737 nsize: 4
738 args: -test_jacobian_domain_error -snes_converged_reason -snes_type vinewtonssls -snes_check_jacobian_domain_error 1
739
740 test:
741 suffix: 10
742 requires: double !complex !single
743 nsize: 4
744 args: -test_jacobian_domain_error -snes_converged_reason -snes_type qn -snes_qn_scale_type jacobian -snes_check_jacobian_domain_error 1
745
746 test:
747 suffix: 11
748 requires: double !complex !single
749 nsize: 4
750 args: -test_jacobian_domain_error -snes_converged_reason -snes_type ms -snes_ms_type m62 -snes_ms_damping 0.9 -snes_check_jacobian_domain_error 1
751
752 test:
753 suffix: 12
754 args: -view_initial
755 filter: grep -v "type:"
756
757 test:
758 suffix: 13
759 requires: double !complex !single
760 nsize: 4
761 args: -test_jacobian_domain_error -snes_converged_reason -snes_type newtontrdc -snes_check_jacobian_domain_error 1
762
763 TEST*/
764