Actual source code: ex15.c

petsc-3.5.2 2014-09-08
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: macro __FUNCT__ to define routine names for use in error handling.\n\
  5: Input parameters include:\n\
  6:   -user_defined_pc : Activate a user-defined preconditioner\n\n";

  8: /*T
  9:    Concepts: KSP^basic parallel example
 10:    Concepts: PC^setting a user-defined shell preconditioner
 11:    Concepts: error handling^Using the macro __FUNCT__ to define routine names;
 12:    Processors: n
 13: 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:    we define the macro __FUNCT__ to be a string containing the routine name.
 39:    If defined, this macro is used in the PETSc error handlers to provide a
 40:    complete traceback of routine names.  All PETSc library routines use this
 41:    macro, and users can optionally employ it as well in their application
 42:    codes.  Note that users can get a traceback of PETSc errors regardless of
 43:    whether they define __FUNCT__ in application codes; this macro merely
 44:    provides the added traceback detail of the application routine names.
 45: */

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

 62:   PetscInitialize(&argc,&args,(char*)0,help);
 63:   PetscOptionsGetInt(NULL,"-m",&m,NULL);
 64:   PetscOptionsGetInt(NULL,"-n",&n,NULL);

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

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

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

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

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

126:   /*
127:      Set exact solution; then compute right-hand-side vector.
128:   */
129:   VecSet(u,one);
130:   MatMult(A,u,b);

132:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133:                 Create the linear solver and set various options
134:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

136:   /*
137:      Create linear solver context
138:   */
139:   KSPCreate(PETSC_COMM_WORLD,&ksp);

141:   /*
142:      Set operators. Here the matrix that defines the linear system
143:      also serves as the preconditioning matrix.
144:   */
145:   KSPSetOperators(ksp,A,A);

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

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

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

169:     /* (Required) Set the user-defined routine for applying the preconditioner */
170:     PCShellSetApply(pc,SampleShellPCApply);
171:     PCShellSetContext(pc,shell);

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

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

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

183:   } else {
184:     PCSetType(pc,PCJACOBI);
185:   }

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

196:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197:                       Solve the linear system
198:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

200:   KSPSolve(ksp,b,x);

202:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203:                       Check solution and clean up
204:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

222:   PetscFinalize();
223:   return 0;

225: }

227: /***********************************************************************/
228: /*          Routines for a user-defined shell preconditioner           */
229: /***********************************************************************/

233: /*
234:    SampleShellPCCreate - This routine creates a user-defined
235:    preconditioner context.

237:    Output Parameter:
238: .  shell - user-defined preconditioner context
239: */
240: PetscErrorCode SampleShellPCCreate(SampleShellPC **shell)
241: {
242:   SampleShellPC  *newctx;

245:   PetscNew(&newctx);
246:   newctx->diag = 0;
247:   *shell       = newctx;
248:   return 0;
249: }
250: /* ------------------------------------------------------------------- */
253: /*
254:    SampleShellPCSetUp - This routine sets up a user-defined
255:    preconditioner context.

257:    Input Parameters:
258: .  pc    - preconditioner object
259: .  pmat  - preconditioner matrix
260: .  x     - vector

262:    Output Parameter:
263: .  shell - fully set up user-defined preconditioner context

265:    Notes:
266:    In this example, we define the shell preconditioner to be Jacobi's
267:    method.  Thus, here we create a work vector for storing the reciprocal
268:    of the diagonal of the preconditioner matrix; this vector is then
269:    used within the routine SampleShellPCApply().
270: */
271: PetscErrorCode SampleShellPCSetUp(PC pc,Mat pmat,Vec x)
272: {
273:   SampleShellPC  *shell;
274:   Vec            diag;

277:   PCShellGetContext(pc,(void**)&shell);
278:   VecDuplicate(x,&diag);
279:   MatGetDiagonal(pmat,diag);
280:   VecReciprocal(diag);

282:   shell->diag = diag;
283:   return 0;
284: }
285: /* ------------------------------------------------------------------- */
288: /*
289:    SampleShellPCApply - This routine demonstrates the use of a
290:    user-provided preconditioner.

292:    Input Parameters:
293: +  pc - preconditioner object
294: -  x - input vector

296:    Output Parameter:
297: .  y - preconditioned vector

299:    Notes:
300:    This code implements the Jacobi preconditioner, merely as an
301:    example of working with a PCSHELL.  Note that the Jacobi method
302:    is already provided within PETSc.
303: */
304: PetscErrorCode SampleShellPCApply(PC pc,Vec x,Vec y)
305: {
306:   SampleShellPC  *shell;

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

312:   return 0;
313: }
314: /* ------------------------------------------------------------------- */
317: /*
318:    SampleShellPCDestroy - This routine destroys a user-defined
319:    preconditioner context.

321:    Input Parameter:
322: .  shell - user-defined preconditioner context
323: */
324: PetscErrorCode SampleShellPCDestroy(PC pc)
325: {
326:   SampleShellPC  *shell;

329:   PCShellGetContext(pc,(void**)&shell);
330:   VecDestroy(&shell->diag);
331:   PetscFree(shell);

333:   return 0;
334: }