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