Actual source code: ex15.c


  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: /*
  8:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
  9:   automatically includes:
 10:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 11:      petscmat.h - matrices
 12:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 13:      petscviewer.h - viewers               petscpc.h  - preconditioners
 14: */
 15: #include <petscksp.h>

 17: /* Define context for user-provided preconditioner */
 18: typedef struct {
 19:   Vec diag;
 20: } SampleShellPC;

 22: /* Declare routines for user-provided preconditioner */
 23: extern PetscErrorCode SampleShellPCCreate(SampleShellPC**);
 24: extern PetscErrorCode SampleShellPCSetUp(PC,Mat,Vec);
 25: extern PetscErrorCode SampleShellPCApply(PC,Vec x,Vec y);
 26: extern PetscErrorCode SampleShellPCDestroy(PC);

 28: /*
 29:    User-defined routines.  Note that immediately before each routine below,
 30:    If defined, this macro is used in the PETSc error handlers to provide a
 31:    complete traceback of routine names.  All PETSc library routines use this
 32:    macro, and users can optionally employ it as well in their application
 33:    codes.  Note that users can get a traceback of PETSc errors regardless of
 34:    provides the added traceback detail of the application routine names.
 35: */

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

 50:   PetscInitialize(&argc,&args,(char*)0,help);
 51:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 52:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);

 54:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 55:          Compute the matrix and right-hand-side vector that define
 56:          the linear system, Ax = b.
 57:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 58:   /*
 59:      Create parallel matrix, specifying only its global dimensions.
 60:      When using MatCreate(), the matrix format can be specified at
 61:      runtime. Also, the parallel partioning of the matrix is
 62:      determined by PETSc at runtime.
 63:   */
 64:   MatCreate(PETSC_COMM_WORLD,&A);
 65:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
 66:   MatSetFromOptions(A);
 67:   MatSetUp(A);

 69:   /*
 70:      Currently, all PETSc parallel matrix formats are partitioned by
 71:      contiguous chunks of rows across the processors.  Determine which
 72:      rows of the matrix are locally owned.
 73:   */
 74:   MatGetOwnershipRange(A,&Istart,&Iend);

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

 92:   /*
 93:      Assemble matrix, using the 2-step process:
 94:        MatAssemblyBegin(), MatAssemblyEnd()
 95:      Computations can be done while messages are in transition
 96:      by placing code between these two statements.
 97:   */
 98:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 99:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

101:   /*
102:      Create parallel vectors.
103:       - When using VecCreate() VecSetSizes() and VecSetFromOptions(),
104:         we specify only the vector's global
105:         dimension; the parallel partitioning is determined at runtime.
106:       - Note: We form 1 vector from scratch and then duplicate as needed.
107:   */
108:   VecCreate(PETSC_COMM_WORLD,&u);
109:   VecSetSizes(u,PETSC_DECIDE,m*n);
110:   VecSetFromOptions(u);
111:   VecDuplicate(u,&b);
112:   VecDuplicate(b,&x);

114:   /*
115:      Set exact solution; then compute right-hand-side vector.
116:   */
117:   VecSet(u,one);
118:   MatMult(A,u,b);

120:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121:                 Create the linear solver and set various options
122:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

124:   /*
125:      Create linear solver context
126:   */
127:   KSPCreate(PETSC_COMM_WORLD,&ksp);

129:   /*
130:      Set operators. Here the matrix that defines the linear system
131:      also serves as the preconditioning matrix.
132:   */
133:   KSPSetOperators(ksp,A,A);

135:   /*
136:      Set linear solver defaults for this problem (optional).
137:      - By extracting the KSP and PC contexts from the KSP context,
138:        we can then directly call any KSP and PC routines
139:        to set various options.
140:   */
141:   KSPGetPC(ksp,&pc);
142:   KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,
143:                           PETSC_DEFAULT);

145:   /*
146:      Set a user-defined "shell" preconditioner if desired
147:   */
148:   PetscOptionsGetBool(NULL,NULL,"-user_defined_pc",&user_defined_pc,NULL);
149:   if (user_defined_pc) {
150:     /* (Required) Indicate to PETSc that we're using a "shell" preconditioner */
151:     PCSetType(pc,PCSHELL);

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

157:     /* (Required) Set the user-defined routine for applying the preconditioner */
158:     PCShellSetApply(pc,SampleShellPCApply);
159:     PCShellSetContext(pc,shell);

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

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

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

171:   } else {
172:     PCSetType(pc,PCJACOBI);
173:   }

175:   /*
176:     Set runtime options, e.g.,
177:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
178:     These options will override those specified above as long as
179:     KSPSetFromOptions() is called _after_ any other customization
180:     routines.
181:   */
182:   KSPSetFromOptions(ksp);

184:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185:                       Solve the linear system
186:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

188:   KSPSolve(ksp,b,x);

190:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
191:                       Check solution and clean up
192:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

194:   /*
195:      Check the error
196:   */
197:   VecAXPY(x,none,u);
198:   VecNorm(x,NORM_2,&norm);
199:   KSPGetIterationNumber(ksp,&its);
200:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);

202:   /*
203:      Free work space.  All PETSc objects should be destroyed when they
204:      are no longer needed.
205:   */
206:   KSPDestroy(&ksp);
207:   VecDestroy(&u));  PetscCall(VecDestroy(&x);
208:   VecDestroy(&b));  PetscCall(MatDestroy(&A);

210:   PetscFinalize();
211:   return 0;

213: }

215: /***********************************************************************/
216: /*          Routines for a user-defined shell preconditioner           */
217: /***********************************************************************/

219: /*
220:    SampleShellPCCreate - This routine creates a user-defined
221:    preconditioner context.

223:    Output Parameter:
224: .  shell - user-defined preconditioner context
225: */
226: PetscErrorCode SampleShellPCCreate(SampleShellPC **shell)
227: {
228:   SampleShellPC  *newctx;

230:   PetscNew(&newctx);
231:   newctx->diag = 0;
232:   *shell       = newctx;
233:   return 0;
234: }
235: /* ------------------------------------------------------------------- */
236: /*
237:    SampleShellPCSetUp - This routine sets up a user-defined
238:    preconditioner context.

240:    Input Parameters:
241: .  pc    - preconditioner object
242: .  pmat  - preconditioner matrix
243: .  x     - vector

245:    Output Parameter:
246: .  shell - fully set up user-defined preconditioner context

248:    Notes:
249:    In this example, we define the shell preconditioner to be Jacobi's
250:    method.  Thus, here we create a work vector for storing the reciprocal
251:    of the diagonal of the preconditioner matrix; this vector is then
252:    used within the routine SampleShellPCApply().
253: */
254: PetscErrorCode SampleShellPCSetUp(PC pc,Mat pmat,Vec x)
255: {
256:   SampleShellPC  *shell;
257:   Vec            diag;

259:   PCShellGetContext(pc,&shell);
260:   VecDuplicate(x,&diag);
261:   MatGetDiagonal(pmat,diag);
262:   VecReciprocal(diag);

264:   shell->diag = diag;
265:   return 0;
266: }
267: /* ------------------------------------------------------------------- */
268: /*
269:    SampleShellPCApply - This routine demonstrates the use of a
270:    user-provided preconditioner.

272:    Input Parameters:
273: +  pc - preconditioner object
274: -  x - input vector

276:    Output Parameter:
277: .  y - preconditioned vector

279:    Notes:
280:    This code implements the Jacobi preconditioner, merely as an
281:    example of working with a PCSHELL.  Note that the Jacobi method
282:    is already provided within PETSc.
283: */
284: PetscErrorCode SampleShellPCApply(PC pc,Vec x,Vec y)
285: {
286:   SampleShellPC  *shell;

288:   PCShellGetContext(pc,&shell);
289:   VecPointwiseMult(y,x,shell->diag);

291:   return 0;
292: }
293: /* ------------------------------------------------------------------- */
294: /*
295:    SampleShellPCDestroy - This routine destroys a user-defined
296:    preconditioner context.

298:    Input Parameter:
299: .  shell - user-defined preconditioner context
300: */
301: PetscErrorCode SampleShellPCDestroy(PC pc)
302: {
303:   SampleShellPC  *shell;

305:   PCShellGetContext(pc,&shell);
306:   VecDestroy(&shell->diag);
307:   PetscFree(shell);

309:   return 0;
310: }

312: /*TEST

314:    build:
315:       requires: !complex !single

317:    test:
318:       nsize: 2
319:       args: -ksp_view -user_defined_pc -ksp_gmres_cgs_refinement_type refine_always

321:    test:
322:       suffix: tsirm
323:       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
324:       timeoutfactor: 4

326: TEST*/