Actual source code: ex15.c

petsc-3.14.6 2021-03-30
Report Typos and Errors

  2: static char help[] = "Solves a linear system in parallel with KSP.  Also\n\
  3: illustrates setting a user-defined shell preconditioner and using the\n\
  4: Input parameters include:\n\
  5:   -user_defined_pc : Activate a user-defined preconditioner\n\n";

  7: /*T
  8:    Concepts: KSP^basic parallel example
  9:    Concepts: PC^setting a user-defined shell preconditioner
 10:    Processors: n
 11: T*/



 15: /*
 16:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
 17:   automatically includes:
 18:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 19:      petscmat.h - matrices
 20:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 21:      petscviewer.h - viewers               petscpc.h  - preconditioners
 22: */
 23: #include <petscksp.h>

 25: /* Define context for user-provided preconditioner */
 26: typedef struct {
 27:   Vec diag;
 28: } SampleShellPC;

 30: /* Declare routines for user-provided preconditioner */
 31: extern PetscErrorCode SampleShellPCCreate(SampleShellPC**);
 32: extern PetscErrorCode SampleShellPCSetUp(PC,Mat,Vec);
 33: extern PetscErrorCode SampleShellPCApply(PC,Vec x,Vec y);
 34: extern PetscErrorCode SampleShellPCDestroy(PC);

 36: /*
 37:    User-defined routines.  Note that immediately before each routine below,
 38:    If defined, this macro is used in the PETSc error handlers to provide a
 39:    complete traceback of routine names.  All PETSc library routines use this
 40:    macro, and users can optionally employ it as well in their application
 41:    codes.  Note that users can get a traceback of PETSc errors regardless of
 42:    provides the added traceback detail of the application routine names.
 43: */

 45: int main(int argc,char **args)
 46: {
 47:   Vec            x,b,u;   /* approx solution, RHS, exact solution */
 48:   Mat            A;         /* linear system matrix */
 49:   KSP            ksp;      /* linear solver context */
 50:   PC             pc;        /* preconditioner context */
 51:   PetscReal      norm;      /* norm of solution error */
 52:   SampleShellPC  *shell;    /* user-defined preconditioner context */
 53:   PetscScalar    v,one = 1.0,none = -1.0;
 54:   PetscInt       i,j,Ii,J,Istart,Iend,m = 8,n = 7,its;
 56:   PetscBool      user_defined_pc = PETSC_FALSE;

 58:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 59:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 60:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);

 62:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 63:          Compute the matrix and right-hand-side vector that define
 64:          the linear system, Ax = b.
 65:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 66:   /*
 67:      Create parallel matrix, specifying only its global dimensions.
 68:      When using MatCreate(), the matrix format can be specified at
 69:      runtime. Also, the parallel partioning of the matrix is
 70:      determined by PETSc at runtime.
 71:   */
 72:   MatCreate(PETSC_COMM_WORLD,&A);
 73:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
 74:   MatSetFromOptions(A);
 75:   MatSetUp(A);

 77:   /*
 78:      Currently, all PETSc parallel matrix formats are partitioned by
 79:      contiguous chunks of rows across the processors.  Determine which
 80:      rows of the matrix are locally owned.
 81:   */
 82:   MatGetOwnershipRange(A,&Istart,&Iend);

 84:   /*
 85:      Set matrix elements for the 2-D, five-point stencil in parallel.
 86:       - Each processor needs to insert only elements that it owns
 87:         locally (but any non-local elements will be sent to the
 88:         appropriate processor during matrix assembly).
 89:       - Always specify global rows and columns of matrix entries.
 90:    */
 91:   for (Ii=Istart; Ii<Iend; Ii++) {
 92:     v = -1.0; i = Ii/n; j = Ii - i*n;
 93:     if (i>0)   {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 94:     if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 95:     if (j>0)   {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 96:     if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
 97:     v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
 98:   }

100:   /*
101:      Assemble matrix, using the 2-step process:
102:        MatAssemblyBegin(), MatAssemblyEnd()
103:      Computations can be done while messages are in transition
104:      by placing code between these two statements.
105:   */
106:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
107:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

109:   /*
110:      Create parallel vectors.
111:       - When using VecCreate() VecSetSizes() and VecSetFromOptions(),
112:         we specify only the vector's global
113:         dimension; the parallel partitioning is determined at runtime.
114:       - Note: We form 1 vector from scratch and then duplicate as needed.
115:   */
116:   VecCreate(PETSC_COMM_WORLD,&u);
117:   VecSetSizes(u,PETSC_DECIDE,m*n);
118:   VecSetFromOptions(u);
119:   VecDuplicate(u,&b);
120:   VecDuplicate(b,&x);

122:   /*
123:      Set exact solution; then compute right-hand-side vector.
124:   */
125:   VecSet(u,one);
126:   MatMult(A,u,b);

128:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129:                 Create the linear solver and set various options
130:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

132:   /*
133:      Create linear solver context
134:   */
135:   KSPCreate(PETSC_COMM_WORLD,&ksp);

137:   /*
138:      Set operators. Here the matrix that defines the linear system
139:      also serves as the preconditioning matrix.
140:   */
141:   KSPSetOperators(ksp,A,A);

143:   /*
144:      Set linear solver defaults for this problem (optional).
145:      - By extracting the KSP and PC contexts from the KSP context,
146:        we can then directly call any KSP and PC routines
147:        to set various options.
148:   */
149:   KSPGetPC(ksp,&pc);
150:   KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,
151:                           PETSC_DEFAULT);

153:   /*
154:      Set a user-defined "shell" preconditioner if desired
155:   */
156:   PetscOptionsGetBool(NULL,NULL,"-user_defined_pc",&user_defined_pc,NULL);
157:   if (user_defined_pc) {
158:     /* (Required) Indicate to PETSc that we're using a "shell" preconditioner */
159:     PCSetType(pc,PCSHELL);

161:     /* (Optional) Create a context for the user-defined preconditioner; this
162:        context can be used to contain any application-specific data. */
163:     SampleShellPCCreate(&shell);

165:     /* (Required) Set the user-defined routine for applying the preconditioner */
166:     PCShellSetApply(pc,SampleShellPCApply);
167:     PCShellSetContext(pc,shell);

169:     /* (Optional) Set user-defined function to free objects used by custom preconditioner */
170:     PCShellSetDestroy(pc,SampleShellPCDestroy);

172:     /* (Optional) Set a name for the preconditioner, used for PCView() */
173:     PCShellSetName(pc,"MyPreconditioner");

175:     /* (Optional) Do any setup required for the preconditioner */
176:     /* Note: This function could be set with PCShellSetSetUp and it would be called when necessary */
177:     SampleShellPCSetUp(pc,A,x);

179:   } else {
180:     PCSetType(pc,PCJACOBI);
181:   }

183:   /*
184:     Set runtime options, e.g.,
185:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
186:     These options will override those specified above as long as
187:     KSPSetFromOptions() is called _after_ any other customization
188:     routines.
189:   */
190:   KSPSetFromOptions(ksp);

192:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
193:                       Solve the linear system
194:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

196:   KSPSolve(ksp,b,x);

198:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199:                       Check solution and clean up
200:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

202:   /*
203:      Check the error
204:   */
205:   VecAXPY(x,none,u);
206:   VecNorm(x,NORM_2,&norm);
207:   KSPGetIterationNumber(ksp,&its);
208:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);

210:   /*
211:      Free work space.  All PETSc objects should be destroyed when they
212:      are no longer needed.
213:   */
214:   KSPDestroy(&ksp);
215:   VecDestroy(&u);  VecDestroy(&x);
216:   VecDestroy(&b);  MatDestroy(&A);

218:   PetscFinalize();
219:   return ierr;

221: }

223: /***********************************************************************/
224: /*          Routines for a user-defined shell preconditioner           */
225: /***********************************************************************/

227: /*
228:    SampleShellPCCreate - This routine creates a user-defined
229:    preconditioner context.

231:    Output Parameter:
232: .  shell - user-defined preconditioner context
233: */
234: PetscErrorCode SampleShellPCCreate(SampleShellPC **shell)
235: {
236:   SampleShellPC  *newctx;

239:   PetscNew(&newctx);
240:   newctx->diag = 0;
241:   *shell       = newctx;
242:   return 0;
243: }
244: /* ------------------------------------------------------------------- */
245: /*
246:    SampleShellPCSetUp - This routine sets up a user-defined
247:    preconditioner context.

249:    Input Parameters:
250: .  pc    - preconditioner object
251: .  pmat  - preconditioner matrix
252: .  x     - vector

254:    Output Parameter:
255: .  shell - fully set up user-defined preconditioner context

257:    Notes:
258:    In this example, we define the shell preconditioner to be Jacobi's
259:    method.  Thus, here we create a work vector for storing the reciprocal
260:    of the diagonal of the preconditioner matrix; this vector is then
261:    used within the routine SampleShellPCApply().
262: */
263: PetscErrorCode SampleShellPCSetUp(PC pc,Mat pmat,Vec x)
264: {
265:   SampleShellPC  *shell;
266:   Vec            diag;

269:   PCShellGetContext(pc,(void**)&shell);
270:   VecDuplicate(x,&diag);
271:   MatGetDiagonal(pmat,diag);
272:   VecReciprocal(diag);

274:   shell->diag = diag;
275:   return 0;
276: }
277: /* ------------------------------------------------------------------- */
278: /*
279:    SampleShellPCApply - This routine demonstrates the use of a
280:    user-provided preconditioner.

282:    Input Parameters:
283: +  pc - preconditioner object
284: -  x - input vector

286:    Output Parameter:
287: .  y - preconditioned vector

289:    Notes:
290:    This code implements the Jacobi preconditioner, merely as an
291:    example of working with a PCSHELL.  Note that the Jacobi method
292:    is already provided within PETSc.
293: */
294: PetscErrorCode SampleShellPCApply(PC pc,Vec x,Vec y)
295: {
296:   SampleShellPC  *shell;

299:   PCShellGetContext(pc,(void**)&shell);
300:   VecPointwiseMult(y,x,shell->diag);

302:   return 0;
303: }
304: /* ------------------------------------------------------------------- */
305: /*
306:    SampleShellPCDestroy - This routine destroys a user-defined
307:    preconditioner context.

309:    Input Parameter:
310: .  shell - user-defined preconditioner context
311: */
312: PetscErrorCode SampleShellPCDestroy(PC pc)
313: {
314:   SampleShellPC  *shell;

317:   PCShellGetContext(pc,(void**)&shell);
318:   VecDestroy(&shell->diag);
319:   PetscFree(shell);

321:   return 0;
322: }


325: /*TEST

327:    build:
328:       requires: !complex !single

330:    test:
331:       nsize: 2
332:       args: -ksp_view -user_defined_pc -ksp_gmres_cgs_refinement_type refine_always

334:    test:
335:       suffix: tsirm
336:       args: -m 60 -n 60 -ksp_type tsirm -pc_type ksp -ksp_monitor_short -ksp_ksp_type fgmres -ksp_ksp_rtol 1e-10 -ksp_pc_type mg -ksp_ksp_max_it 30
337:       timeoutfactor: 4

339: TEST*/