Actual source code: ex15.c

petsc-master 2017-06-27
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*/

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

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

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

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

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

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

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

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

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

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

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

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

126:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127:                 Create the linear solver and set various options
128:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

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

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

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

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

177:   } else {
178:     PCSetType(pc,PCJACOBI);
179:   }

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

190:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
191:                       Solve the linear system
192:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

194:   KSPSolve(ksp,b,x);

196:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197:                       Check solution and clean up
198:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

216:   PetscFinalize();
217:   return ierr;

219: }

221: /***********************************************************************/
222: /*          Routines for a user-defined shell preconditioner           */
223: /***********************************************************************/

225: /*
226:    SampleShellPCCreate - This routine creates a user-defined
227:    preconditioner context.

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

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

247:    Input Parameters:
248: .  pc    - preconditioner object
249: .  pmat  - preconditioner matrix
250: .  x     - vector

252:    Output Parameter:
253: .  shell - fully set up user-defined preconditioner context

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

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

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

280:    Input Parameters:
281: +  pc - preconditioner object
282: -  x - input vector

284:    Output Parameter:
285: .  y - preconditioned vector

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

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

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

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

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

319:   return 0;
320: }