1
2 /*
3 The PC (preconditioner) interface routines, callable by users.
4 */
5 #include <petsc/private/pcimpl.h> /*I "petscksp.h" I*/
6 #include <petscdm.h>
7
8 /* Logging support */
9 PetscClassId PC_CLASSID;
10 PetscLogEvent PC_SetUp, PC_SetUpOnBlocks, PC_Apply, PC_MatApply, PC_ApplyCoarse, PC_ApplyMultiple, PC_ApplySymmetricLeft;
11 PetscLogEvent PC_ApplySymmetricRight, PC_ModifySubMatrices, PC_ApplyOnBlocks, PC_ApplyTransposeOnBlocks;
12 PetscInt PetscMGLevelId;
13
PCGetDefaultType_Private(PC pc,const char * type[])14 PetscErrorCode PCGetDefaultType_Private(PC pc,const char *type[])
15 {
16 PetscErrorCode ierr;
17 PetscMPIInt size;
18 PetscBool hasop,flg1,flg2,set,flg3;
19
20 PetscFunctionBegin;
21 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr);
22 if (pc->pmat) {
23 ierr = MatHasOperation(pc->pmat,MATOP_GET_DIAGONAL_BLOCK,&hasop);CHKERRQ(ierr);
24 if (size == 1) {
25 ierr = MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ICC,&flg1);CHKERRQ(ierr);
26 ierr = MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ILU,&flg2);CHKERRQ(ierr);
27 ierr = MatIsSymmetricKnown(pc->pmat,&set,&flg3);CHKERRQ(ierr);
28 if (flg1 && (!flg2 || (set && flg3))) {
29 *type = PCICC;
30 } else if (flg2) {
31 *type = PCILU;
32 } else if (hasop) { /* likely is a parallel matrix run on one processor */
33 *type = PCBJACOBI;
34 } else {
35 *type = PCNONE;
36 }
37 } else {
38 if (hasop) {
39 *type = PCBJACOBI;
40 } else {
41 *type = PCNONE;
42 }
43 }
44 } else {
45 if (size == 1) {
46 *type = PCILU;
47 } else {
48 *type = PCBJACOBI;
49 }
50 }
51 PetscFunctionReturn(0);
52 }
53
54 /*@
55 PCReset - Resets a PC context to the pcsetupcalled = 0 state and removes any allocated Vecs and Mats
56
57 Collective on PC
58
59 Input Parameter:
60 . pc - the preconditioner context
61
62 Level: developer
63
64 Notes:
65 This allows a PC to be reused for a different sized linear system but using the same options that have been previously set in the PC
66
67 .seealso: PCCreate(), PCSetUp()
68 @*/
PCReset(PC pc)69 PetscErrorCode PCReset(PC pc)
70 {
71 PetscErrorCode ierr;
72
73 PetscFunctionBegin;
74 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
75 if (pc->ops->reset) {
76 ierr = (*pc->ops->reset)(pc);CHKERRQ(ierr);
77 }
78 ierr = VecDestroy(&pc->diagonalscaleright);CHKERRQ(ierr);
79 ierr = VecDestroy(&pc->diagonalscaleleft);CHKERRQ(ierr);
80 ierr = MatDestroy(&pc->pmat);CHKERRQ(ierr);
81 ierr = MatDestroy(&pc->mat);CHKERRQ(ierr);
82
83 pc->setupcalled = 0;
84 PetscFunctionReturn(0);
85 }
86
87 /*@
88 PCDestroy - Destroys PC context that was created with PCCreate().
89
90 Collective on PC
91
92 Input Parameter:
93 . pc - the preconditioner context
94
95 Level: developer
96
97 .seealso: PCCreate(), PCSetUp()
98 @*/
PCDestroy(PC * pc)99 PetscErrorCode PCDestroy(PC *pc)
100 {
101 PetscErrorCode ierr;
102
103 PetscFunctionBegin;
104 if (!*pc) PetscFunctionReturn(0);
105 PetscValidHeaderSpecific((*pc),PC_CLASSID,1);
106 if (--((PetscObject)(*pc))->refct > 0) {*pc = NULL; PetscFunctionReturn(0);}
107
108 ierr = PCReset(*pc);CHKERRQ(ierr);
109
110 /* if memory was published with SAWs then destroy it */
111 ierr = PetscObjectSAWsViewOff((PetscObject)*pc);CHKERRQ(ierr);
112 if ((*pc)->ops->destroy) {ierr = (*(*pc)->ops->destroy)((*pc));CHKERRQ(ierr);}
113 ierr = DMDestroy(&(*pc)->dm);CHKERRQ(ierr);
114 ierr = PetscHeaderDestroy(pc);CHKERRQ(ierr);
115 PetscFunctionReturn(0);
116 }
117
118 /*@C
119 PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
120 scaling as needed by certain time-stepping codes.
121
122 Logically Collective on PC
123
124 Input Parameter:
125 . pc - the preconditioner context
126
127 Output Parameter:
128 . flag - PETSC_TRUE if it applies the scaling
129
130 Level: developer
131
132 Notes:
133 If this returns PETSC_TRUE then the system solved via the Krylov method is
134 $ D M A D^{-1} y = D M b for left preconditioning or
135 $ D A M D^{-1} z = D b for right preconditioning
136
137 .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCSetDiagonalScale()
138 @*/
PCGetDiagonalScale(PC pc,PetscBool * flag)139 PetscErrorCode PCGetDiagonalScale(PC pc,PetscBool *flag)
140 {
141 PetscFunctionBegin;
142 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
143 PetscValidBoolPointer(flag,2);
144 *flag = pc->diagonalscale;
145 PetscFunctionReturn(0);
146 }
147
148 /*@
149 PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
150 scaling as needed by certain time-stepping codes.
151
152 Logically Collective on PC
153
154 Input Parameters:
155 + pc - the preconditioner context
156 - s - scaling vector
157
158 Level: intermediate
159
160 Notes:
161 The system solved via the Krylov method is
162 $ D M A D^{-1} y = D M b for left preconditioning or
163 $ D A M D^{-1} z = D b for right preconditioning
164
165 PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
166
167 .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCGetDiagonalScale()
168 @*/
PCSetDiagonalScale(PC pc,Vec s)169 PetscErrorCode PCSetDiagonalScale(PC pc,Vec s)
170 {
171 PetscErrorCode ierr;
172
173 PetscFunctionBegin;
174 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
175 PetscValidHeaderSpecific(s,VEC_CLASSID,2);
176 pc->diagonalscale = PETSC_TRUE;
177
178 ierr = PetscObjectReference((PetscObject)s);CHKERRQ(ierr);
179 ierr = VecDestroy(&pc->diagonalscaleleft);CHKERRQ(ierr);
180
181 pc->diagonalscaleleft = s;
182
183 ierr = VecDuplicate(s,&pc->diagonalscaleright);CHKERRQ(ierr);
184 ierr = VecCopy(s,pc->diagonalscaleright);CHKERRQ(ierr);
185 ierr = VecReciprocal(pc->diagonalscaleright);CHKERRQ(ierr);
186 PetscFunctionReturn(0);
187 }
188
189 /*@
190 PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
191
192 Logically Collective on PC
193
194 Input Parameters:
195 + pc - the preconditioner context
196 . in - input vector
197 - out - scaled vector (maybe the same as in)
198
199 Level: intermediate
200
201 Notes:
202 The system solved via the Krylov method is
203 $ D M A D^{-1} y = D M b for left preconditioning or
204 $ D A M D^{-1} z = D b for right preconditioning
205
206 PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
207
208 If diagonal scaling is turned off and in is not out then in is copied to out
209
210 .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
211 @*/
PCDiagonalScaleLeft(PC pc,Vec in,Vec out)212 PetscErrorCode PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
213 {
214 PetscErrorCode ierr;
215
216 PetscFunctionBegin;
217 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
218 PetscValidHeaderSpecific(in,VEC_CLASSID,2);
219 PetscValidHeaderSpecific(out,VEC_CLASSID,3);
220 if (pc->diagonalscale) {
221 ierr = VecPointwiseMult(out,pc->diagonalscaleleft,in);CHKERRQ(ierr);
222 } else if (in != out) {
223 ierr = VecCopy(in,out);CHKERRQ(ierr);
224 }
225 PetscFunctionReturn(0);
226 }
227
228 /*@
229 PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
230
231 Logically Collective on PC
232
233 Input Parameters:
234 + pc - the preconditioner context
235 . in - input vector
236 - out - scaled vector (maybe the same as in)
237
238 Level: intermediate
239
240 Notes:
241 The system solved via the Krylov method is
242 $ D M A D^{-1} y = D M b for left preconditioning or
243 $ D A M D^{-1} z = D b for right preconditioning
244
245 PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
246
247 If diagonal scaling is turned off and in is not out then in is copied to out
248
249 .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
250 @*/
PCDiagonalScaleRight(PC pc,Vec in,Vec out)251 PetscErrorCode PCDiagonalScaleRight(PC pc,Vec in,Vec out)
252 {
253 PetscErrorCode ierr;
254
255 PetscFunctionBegin;
256 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
257 PetscValidHeaderSpecific(in,VEC_CLASSID,2);
258 PetscValidHeaderSpecific(out,VEC_CLASSID,3);
259 if (pc->diagonalscale) {
260 ierr = VecPointwiseMult(out,pc->diagonalscaleright,in);CHKERRQ(ierr);
261 } else if (in != out) {
262 ierr = VecCopy(in,out);CHKERRQ(ierr);
263 }
264 PetscFunctionReturn(0);
265 }
266
267 /*@
268 PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
269 operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
270 TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
271
272 Logically Collective on PC
273
274 Input Parameters:
275 + pc - the preconditioner context
276 - flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
277
278 Options Database Key:
279 . -pc_use_amat <true,false>
280
281 Notes:
282 For the common case in which the linear system matrix and the matrix used to construct the
283 preconditioner are identical, this routine is does nothing.
284
285 Level: intermediate
286
287 .seealso: PCGetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
288 @*/
PCSetUseAmat(PC pc,PetscBool flg)289 PetscErrorCode PCSetUseAmat(PC pc,PetscBool flg)
290 {
291 PetscFunctionBegin;
292 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
293 pc->useAmat = flg;
294 PetscFunctionReturn(0);
295 }
296
297 /*@
298 PCSetErrorIfFailure - Causes PC to generate an error if a FPE, for example a zero pivot, is detected.
299
300 Logically Collective on PC
301
302 Input Parameters:
303 + pc - iterative context obtained from PCCreate()
304 - flg - PETSC_TRUE indicates you want the error generated
305
306 Level: advanced
307
308 Notes:
309 Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call KSPGetConvergedReason() after a KSPSolve()
310 to determine if it has converged or failed. Or use -ksp_error_if_not_converged to cause the program to terminate as soon as lack of convergence is
311 detected.
312
313 This is propagated into KSPs used by this PC, which then propagate it into PCs used by those KSPs
314
315 .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
316 @*/
PCSetErrorIfFailure(PC pc,PetscBool flg)317 PetscErrorCode PCSetErrorIfFailure(PC pc,PetscBool flg)
318 {
319 PetscFunctionBegin;
320 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
321 PetscValidLogicalCollectiveBool(pc,flg,2);
322 pc->erroriffailure = flg;
323 PetscFunctionReturn(0);
324 }
325
326 /*@
327 PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
328 operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
329 TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
330
331 Logically Collective on PC
332
333 Input Parameter:
334 . pc - the preconditioner context
335
336 Output Parameter:
337 . flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
338
339 Notes:
340 For the common case in which the linear system matrix and the matrix used to construct the
341 preconditioner are identical, this routine is does nothing.
342
343 Level: intermediate
344
345 .seealso: PCSetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
346 @*/
PCGetUseAmat(PC pc,PetscBool * flg)347 PetscErrorCode PCGetUseAmat(PC pc,PetscBool *flg)
348 {
349 PetscFunctionBegin;
350 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
351 *flg = pc->useAmat;
352 PetscFunctionReturn(0);
353 }
354
355 /*@
356 PCCreate - Creates a preconditioner context.
357
358 Collective
359
360 Input Parameter:
361 . comm - MPI communicator
362
363 Output Parameter:
364 . pc - location to put the preconditioner context
365
366 Notes:
367 The default preconditioner for sparse matrices is PCILU or PCICC with 0 fill on one process and block Jacobi with PCILU or PCICC
368 in parallel. For dense matrices it is always PCNONE.
369
370 Level: developer
371
372 .seealso: PCSetUp(), PCApply(), PCDestroy()
373 @*/
PCCreate(MPI_Comm comm,PC * newpc)374 PetscErrorCode PCCreate(MPI_Comm comm,PC *newpc)
375 {
376 PC pc;
377 PetscErrorCode ierr;
378
379 PetscFunctionBegin;
380 PetscValidPointer(newpc,1);
381 *newpc = NULL;
382 ierr = PCInitializePackage();CHKERRQ(ierr);
383
384 ierr = PetscHeaderCreate(pc,PC_CLASSID,"PC","Preconditioner","PC",comm,PCDestroy,PCView);CHKERRQ(ierr);
385
386 pc->mat = NULL;
387 pc->pmat = NULL;
388 pc->setupcalled = 0;
389 pc->setfromoptionscalled = 0;
390 pc->data = NULL;
391 pc->diagonalscale = PETSC_FALSE;
392 pc->diagonalscaleleft = NULL;
393 pc->diagonalscaleright = NULL;
394
395 pc->modifysubmatrices = NULL;
396 pc->modifysubmatricesP = NULL;
397
398 *newpc = pc;
399 PetscFunctionReturn(0);
400
401 }
402
403 /* -------------------------------------------------------------------------------*/
404
405 /*@
406 PCApply - Applies the preconditioner to a vector.
407
408 Collective on PC
409
410 Input Parameters:
411 + pc - the preconditioner context
412 - x - input vector
413
414 Output Parameter:
415 . y - output vector
416
417 Level: developer
418
419 .seealso: PCApplyTranspose(), PCApplyBAorAB()
420 @*/
PCApply(PC pc,Vec x,Vec y)421 PetscErrorCode PCApply(PC pc,Vec x,Vec y)
422 {
423 PetscErrorCode ierr;
424 PetscInt m,n,mv,nv;
425
426 PetscFunctionBegin;
427 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
428 PetscValidHeaderSpecific(x,VEC_CLASSID,2);
429 PetscValidHeaderSpecific(y,VEC_CLASSID,3);
430 if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
431 if (pc->erroriffailure) {ierr = VecValidValues(x,2,PETSC_TRUE);CHKERRQ(ierr);}
432 /* use pmat to check vector sizes since for KSPLQR the pmat may be of a different size than mat */
433 ierr = MatGetLocalSize(pc->pmat,&m,&n);CHKERRQ(ierr);
434 ierr = VecGetLocalSize(x,&nv);CHKERRQ(ierr);
435 ierr = VecGetLocalSize(y,&mv);CHKERRQ(ierr);
436 if (mv != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Preconditioner number of local rows %D does not equal resulting vector number of rows %D",m,mv);
437 if (nv != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Preconditioner number of local columns %D does not equal resulting vector number of rows %D",n,nv);
438 ierr = VecSetErrorIfLocked(y,3);CHKERRQ(ierr);
439
440 ierr = PCSetUp(pc);CHKERRQ(ierr);
441 if (!pc->ops->apply) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply");
442 ierr = VecLockReadPush(x);CHKERRQ(ierr);
443 ierr = PetscLogEventBegin(PC_Apply,pc,x,y,0);CHKERRQ(ierr);
444 ierr = (*pc->ops->apply)(pc,x,y);CHKERRQ(ierr);
445 ierr = PetscLogEventEnd(PC_Apply,pc,x,y,0);CHKERRQ(ierr);
446 if (pc->erroriffailure) {ierr = VecValidValues(y,3,PETSC_FALSE);CHKERRQ(ierr);}
447 ierr = VecLockReadPop(x);CHKERRQ(ierr);
448 PetscFunctionReturn(0);
449 }
450
451 /*@
452 PCMatApply - Applies the preconditioner to multiple vectors stored as a MATDENSE. Like PCApply(), Y and X must be different matrices.
453
454 Collective on PC
455
456 Input Parameters:
457 + pc - the preconditioner context
458 - X - block of input vectors
459
460 Output Parameter:
461 . Y - block of output vectors
462
463 Level: developer
464
465 .seealso: PCApply(), KSPMatSolve()
466 @*/
PCMatApply(PC pc,Mat X,Mat Y)467 PetscErrorCode PCMatApply(PC pc,Mat X,Mat Y)
468 {
469 Mat A;
470 Vec cy, cx;
471 PetscInt m1, M1, m2, M2, n1, N1, n2, N2;
472 PetscBool match;
473 PetscErrorCode ierr;
474
475 PetscFunctionBegin;
476 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
477 PetscValidHeaderSpecific(Y, MAT_CLASSID, 2);
478 PetscValidHeaderSpecific(X, MAT_CLASSID, 3);
479 PetscCheckSameComm(pc, 1, Y, 2);
480 PetscCheckSameComm(pc, 1, X, 3);
481 if (Y == X) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_IDN, "Y and X must be different matrices");
482 ierr = PCGetOperators(pc, NULL, &A);CHKERRQ(ierr);
483 ierr = MatGetLocalSize(A, &m1, NULL);CHKERRQ(ierr);
484 ierr = MatGetLocalSize(Y, &m2, &n2);CHKERRQ(ierr);
485 ierr = MatGetSize(A, &M1, NULL);CHKERRQ(ierr);
486 ierr = MatGetSize(X, &M2, &N2);CHKERRQ(ierr);
487 if (m1 != m2 || M1 != M2) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot use a block of input vectors with (m2,M2) = (%D,%D) for a preconditioner with (m1,M1) = (%D,%D)", m2, M2, m1, M1);
488 ierr = MatGetLocalSize(Y, &m1, &n1);CHKERRQ(ierr);
489 ierr = MatGetSize(Y, &M1, &N1);CHKERRQ(ierr);
490 if (m1 != m2 || M1 != M2 || n1 != n2 || N1 != N2) SETERRQ8(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible block of input vectors (m2,M2)x(n2,N2) = (%D,%D)x(%D,%D) and output vectors (m1,M1)x(n1,N1) = (%D,%D)x(%D,%D)", m2, M2, n2, N2, m1, M1, n1, N1);
491 ierr = PetscObjectBaseTypeCompareAny((PetscObject)Y, &match, MATSEQDENSE, MATMPIDENSE, "");CHKERRQ(ierr);
492 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of output vectors not stored in a dense Mat");
493 ierr = PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, "");CHKERRQ(ierr);
494 if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of input vectors not stored in a dense Mat");
495 ierr = PCSetUp(pc);CHKERRQ(ierr);
496 if (pc->ops->matapply) {
497 ierr = PetscLogEventBegin(PC_MatApply, pc, X, Y, 0);CHKERRQ(ierr);
498 ierr = (*pc->ops->matapply)(pc, X, Y);CHKERRQ(ierr);
499 ierr = PetscLogEventEnd(PC_MatApply, pc, X, Y, 0);CHKERRQ(ierr);
500 } else {
501 ierr = PetscInfo1(pc, "PC type %s applying column by column\n", ((PetscObject)pc)->type_name);CHKERRQ(ierr);
502 for (n2 = 0; n2 < N2; ++n2) {
503 ierr = MatDenseGetColumnVecRead(X, n2, &cx);CHKERRQ(ierr);
504 ierr = MatDenseGetColumnVecWrite(Y, n2, &cy);CHKERRQ(ierr);
505 ierr = PCApply(pc, cx, cy);CHKERRQ(ierr);
506 ierr = MatDenseRestoreColumnVecWrite(Y, n2, &cy);CHKERRQ(ierr);
507 ierr = MatDenseRestoreColumnVecRead(X, n2, &cx);CHKERRQ(ierr);
508 }
509 }
510 PetscFunctionReturn(0);
511 }
512
513 /*@
514 PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
515
516 Collective on PC
517
518 Input Parameters:
519 + pc - the preconditioner context
520 - x - input vector
521
522 Output Parameter:
523 . y - output vector
524
525 Notes:
526 Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
527
528 Level: developer
529
530 .seealso: PCApply(), PCApplySymmetricRight()
531 @*/
PCApplySymmetricLeft(PC pc,Vec x,Vec y)532 PetscErrorCode PCApplySymmetricLeft(PC pc,Vec x,Vec y)
533 {
534 PetscErrorCode ierr;
535
536 PetscFunctionBegin;
537 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
538 PetscValidHeaderSpecific(x,VEC_CLASSID,2);
539 PetscValidHeaderSpecific(y,VEC_CLASSID,3);
540 if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
541 if (pc->erroriffailure) {ierr = VecValidValues(x,2,PETSC_TRUE);CHKERRQ(ierr);}
542 ierr = PCSetUp(pc);CHKERRQ(ierr);
543 if (!pc->ops->applysymmetricleft) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
544 ierr = VecLockReadPush(x);CHKERRQ(ierr);
545 ierr = PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);CHKERRQ(ierr);
546 ierr = (*pc->ops->applysymmetricleft)(pc,x,y);CHKERRQ(ierr);
547 ierr = PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);CHKERRQ(ierr);
548 ierr = VecLockReadPop(x);CHKERRQ(ierr);
549 if (pc->erroriffailure) {ierr = VecValidValues(y,3,PETSC_FALSE);CHKERRQ(ierr);}
550 PetscFunctionReturn(0);
551 }
552
553 /*@
554 PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
555
556 Collective on PC
557
558 Input Parameters:
559 + pc - the preconditioner context
560 - x - input vector
561
562 Output Parameter:
563 . y - output vector
564
565 Level: developer
566
567 Notes:
568 Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
569
570 .seealso: PCApply(), PCApplySymmetricLeft()
571 @*/
PCApplySymmetricRight(PC pc,Vec x,Vec y)572 PetscErrorCode PCApplySymmetricRight(PC pc,Vec x,Vec y)
573 {
574 PetscErrorCode ierr;
575
576 PetscFunctionBegin;
577 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
578 PetscValidHeaderSpecific(x,VEC_CLASSID,2);
579 PetscValidHeaderSpecific(y,VEC_CLASSID,3);
580 if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
581 if (pc->erroriffailure) {ierr = VecValidValues(x,2,PETSC_TRUE);CHKERRQ(ierr);}
582 ierr = PCSetUp(pc);CHKERRQ(ierr);
583 if (!pc->ops->applysymmetricright) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
584 ierr = VecLockReadPush(x);CHKERRQ(ierr);
585 ierr = PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);CHKERRQ(ierr);
586 ierr = (*pc->ops->applysymmetricright)(pc,x,y);CHKERRQ(ierr);
587 ierr = PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);CHKERRQ(ierr);
588 ierr = VecLockReadPop(x);CHKERRQ(ierr);
589 if (pc->erroriffailure) {ierr = VecValidValues(y,3,PETSC_FALSE);CHKERRQ(ierr);}
590 PetscFunctionReturn(0);
591 }
592
593 /*@
594 PCApplyTranspose - Applies the transpose of preconditioner to a vector.
595
596 Collective on PC
597
598 Input Parameters:
599 + pc - the preconditioner context
600 - x - input vector
601
602 Output Parameter:
603 . y - output vector
604
605 Notes:
606 For complex numbers this applies the non-Hermitian transpose.
607
608 Developer Notes:
609 We need to implement a PCApplyHermitianTranspose()
610
611 Level: developer
612
613 .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
614 @*/
PCApplyTranspose(PC pc,Vec x,Vec y)615 PetscErrorCode PCApplyTranspose(PC pc,Vec x,Vec y)
616 {
617 PetscErrorCode ierr;
618
619 PetscFunctionBegin;
620 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
621 PetscValidHeaderSpecific(x,VEC_CLASSID,2);
622 PetscValidHeaderSpecific(y,VEC_CLASSID,3);
623 if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
624 if (pc->erroriffailure) {ierr = VecValidValues(x,2,PETSC_TRUE);CHKERRQ(ierr);}
625 ierr = PCSetUp(pc);CHKERRQ(ierr);
626 if (!pc->ops->applytranspose) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply transpose");
627 ierr = VecLockReadPush(x);CHKERRQ(ierr);
628 ierr = PetscLogEventBegin(PC_Apply,pc,x,y,0);CHKERRQ(ierr);
629 ierr = (*pc->ops->applytranspose)(pc,x,y);CHKERRQ(ierr);
630 ierr = PetscLogEventEnd(PC_Apply,pc,x,y,0);CHKERRQ(ierr);
631 ierr = VecLockReadPop(x);CHKERRQ(ierr);
632 if (pc->erroriffailure) {ierr = VecValidValues(y,3,PETSC_FALSE);CHKERRQ(ierr);}
633 PetscFunctionReturn(0);
634 }
635
636 /*@
637 PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
638
639 Collective on PC
640
641 Input Parameters:
642 . pc - the preconditioner context
643
644 Output Parameter:
645 . flg - PETSC_TRUE if a transpose operation is defined
646
647 Level: developer
648
649 .seealso: PCApplyTranspose()
650 @*/
PCApplyTransposeExists(PC pc,PetscBool * flg)651 PetscErrorCode PCApplyTransposeExists(PC pc,PetscBool *flg)
652 {
653 PetscFunctionBegin;
654 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
655 PetscValidBoolPointer(flg,2);
656 if (pc->ops->applytranspose) *flg = PETSC_TRUE;
657 else *flg = PETSC_FALSE;
658 PetscFunctionReturn(0);
659 }
660
661 /*@
662 PCApplyBAorAB - Applies the preconditioner and operator to a vector. y = B*A*x or y = A*B*x.
663
664 Collective on PC
665
666 Input Parameters:
667 + pc - the preconditioner context
668 . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
669 . x - input vector
670 - work - work vector
671
672 Output Parameter:
673 . y - output vector
674
675 Level: developer
676
677 Notes:
678 If the PC has had PCSetDiagonalScale() set then D M A D^{-1} for left preconditioning or D A M D^{-1} is actually applied. Note that the
679 specific KSPSolve() method must also be written to handle the post-solve "correction" for the diagonal scaling.
680
681 .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
682 @*/
PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)683 PetscErrorCode PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
684 {
685 PetscErrorCode ierr;
686
687 PetscFunctionBegin;
688 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
689 PetscValidLogicalCollectiveEnum(pc,side,2);
690 PetscValidHeaderSpecific(x,VEC_CLASSID,3);
691 PetscValidHeaderSpecific(y,VEC_CLASSID,4);
692 PetscValidHeaderSpecific(work,VEC_CLASSID,5);
693 PetscCheckSameComm(pc,1,x,3);
694 PetscCheckSameComm(pc,1,y,4);
695 PetscCheckSameComm(pc,1,work,5);
696 if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
697 if (side != PC_LEFT && side != PC_SYMMETRIC && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right, left, or symmetric");
698 if (pc->diagonalscale && side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot include diagonal scaling with symmetric preconditioner application");
699 if (pc->erroriffailure) {ierr = VecValidValues(x,3,PETSC_TRUE);CHKERRQ(ierr);}
700
701 ierr = PCSetUp(pc);CHKERRQ(ierr);
702 if (pc->diagonalscale) {
703 if (pc->ops->applyBA) {
704 Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
705 ierr = VecDuplicate(x,&work2);CHKERRQ(ierr);
706 ierr = PCDiagonalScaleRight(pc,x,work2);CHKERRQ(ierr);
707 ierr = (*pc->ops->applyBA)(pc,side,work2,y,work);CHKERRQ(ierr);
708 ierr = PCDiagonalScaleLeft(pc,y,y);CHKERRQ(ierr);
709 ierr = VecDestroy(&work2);CHKERRQ(ierr);
710 } else if (side == PC_RIGHT) {
711 ierr = PCDiagonalScaleRight(pc,x,y);CHKERRQ(ierr);
712 ierr = PCApply(pc,y,work);CHKERRQ(ierr);
713 ierr = MatMult(pc->mat,work,y);CHKERRQ(ierr);
714 ierr = PCDiagonalScaleLeft(pc,y,y);CHKERRQ(ierr);
715 } else if (side == PC_LEFT) {
716 ierr = PCDiagonalScaleRight(pc,x,y);CHKERRQ(ierr);
717 ierr = MatMult(pc->mat,y,work);CHKERRQ(ierr);
718 ierr = PCApply(pc,work,y);CHKERRQ(ierr);
719 ierr = PCDiagonalScaleLeft(pc,y,y);CHKERRQ(ierr);
720 } else if (side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric application of preconditioner");
721 } else {
722 if (pc->ops->applyBA) {
723 ierr = (*pc->ops->applyBA)(pc,side,x,y,work);CHKERRQ(ierr);
724 } else if (side == PC_RIGHT) {
725 ierr = PCApply(pc,x,work);CHKERRQ(ierr);
726 ierr = MatMult(pc->mat,work,y);CHKERRQ(ierr);
727 } else if (side == PC_LEFT) {
728 ierr = MatMult(pc->mat,x,work);CHKERRQ(ierr);
729 ierr = PCApply(pc,work,y);CHKERRQ(ierr);
730 } else if (side == PC_SYMMETRIC) {
731 /* There's an extra copy here; maybe should provide 2 work vectors instead? */
732 ierr = PCApplySymmetricRight(pc,x,work);CHKERRQ(ierr);
733 ierr = MatMult(pc->mat,work,y);CHKERRQ(ierr);
734 ierr = VecCopy(y,work);CHKERRQ(ierr);
735 ierr = PCApplySymmetricLeft(pc,work,y);CHKERRQ(ierr);
736 }
737 }
738 if (pc->erroriffailure) {ierr = VecValidValues(y,4,PETSC_FALSE);CHKERRQ(ierr);}
739 PetscFunctionReturn(0);
740 }
741
742 /*@
743 PCApplyBAorABTranspose - Applies the transpose of the preconditioner
744 and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
745 NOT tr(B*A) = tr(A)*tr(B).
746
747 Collective on PC
748
749 Input Parameters:
750 + pc - the preconditioner context
751 . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
752 . x - input vector
753 - work - work vector
754
755 Output Parameter:
756 . y - output vector
757
758
759 Notes:
760 this routine is used internally so that the same Krylov code can be used to solve A x = b and A' x = b, with a preconditioner
761 defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)
762
763 Level: developer
764
765 .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
766 @*/
PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)767 PetscErrorCode PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
768 {
769 PetscErrorCode ierr;
770
771 PetscFunctionBegin;
772 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
773 PetscValidHeaderSpecific(x,VEC_CLASSID,3);
774 PetscValidHeaderSpecific(y,VEC_CLASSID,4);
775 PetscValidHeaderSpecific(work,VEC_CLASSID,5);
776 if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
777 if (pc->erroriffailure) {ierr = VecValidValues(x,3,PETSC_TRUE);CHKERRQ(ierr);}
778 if (pc->ops->applyBAtranspose) {
779 ierr = (*pc->ops->applyBAtranspose)(pc,side,x,y,work);CHKERRQ(ierr);
780 if (pc->erroriffailure) {ierr = VecValidValues(y,4,PETSC_FALSE);CHKERRQ(ierr);}
781 PetscFunctionReturn(0);
782 }
783 if (side != PC_LEFT && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");
784
785 ierr = PCSetUp(pc);CHKERRQ(ierr);
786 if (side == PC_RIGHT) {
787 ierr = PCApplyTranspose(pc,x,work);CHKERRQ(ierr);
788 ierr = MatMultTranspose(pc->mat,work,y);CHKERRQ(ierr);
789 } else if (side == PC_LEFT) {
790 ierr = MatMultTranspose(pc->mat,x,work);CHKERRQ(ierr);
791 ierr = PCApplyTranspose(pc,work,y);CHKERRQ(ierr);
792 }
793 /* add support for PC_SYMMETRIC */
794 if (pc->erroriffailure) {ierr = VecValidValues(y,4,PETSC_FALSE);CHKERRQ(ierr);}
795 PetscFunctionReturn(0);
796 }
797
798 /* -------------------------------------------------------------------------------*/
799
800 /*@
801 PCApplyRichardsonExists - Determines whether a particular preconditioner has a
802 built-in fast application of Richardson's method.
803
804 Not Collective
805
806 Input Parameter:
807 . pc - the preconditioner
808
809 Output Parameter:
810 . exists - PETSC_TRUE or PETSC_FALSE
811
812 Level: developer
813
814 .seealso: PCApplyRichardson()
815 @*/
PCApplyRichardsonExists(PC pc,PetscBool * exists)816 PetscErrorCode PCApplyRichardsonExists(PC pc,PetscBool *exists)
817 {
818 PetscFunctionBegin;
819 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
820 PetscValidPointer(exists,2);
821 if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
822 else *exists = PETSC_FALSE;
823 PetscFunctionReturn(0);
824 }
825
826 /*@
827 PCApplyRichardson - Applies several steps of Richardson iteration with
828 the particular preconditioner. This routine is usually used by the
829 Krylov solvers and not the application code directly.
830
831 Collective on PC
832
833 Input Parameters:
834 + pc - the preconditioner context
835 . b - the right hand side
836 . w - one work vector
837 . rtol - relative decrease in residual norm convergence criteria
838 . abstol - absolute residual norm convergence criteria
839 . dtol - divergence residual norm increase criteria
840 . its - the number of iterations to apply.
841 - guesszero - if the input x contains nonzero initial guess
842
843 Output Parameter:
844 + outits - number of iterations actually used (for SOR this always equals its)
845 . reason - the reason the apply terminated
846 - y - the solution (also contains initial guess if guesszero is PETSC_FALSE
847
848 Notes:
849 Most preconditioners do not support this function. Use the command
850 PCApplyRichardsonExists() to determine if one does.
851
852 Except for the multigrid PC this routine ignores the convergence tolerances
853 and always runs for the number of iterations
854
855 Level: developer
856
857 .seealso: PCApplyRichardsonExists()
858 @*/
PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt * outits,PCRichardsonConvergedReason * reason)859 PetscErrorCode PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
860 {
861 PetscErrorCode ierr;
862
863 PetscFunctionBegin;
864 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
865 PetscValidHeaderSpecific(b,VEC_CLASSID,2);
866 PetscValidHeaderSpecific(y,VEC_CLASSID,3);
867 PetscValidHeaderSpecific(w,VEC_CLASSID,4);
868 if (b == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"b and y must be different vectors");
869 ierr = PCSetUp(pc);CHKERRQ(ierr);
870 if (!pc->ops->applyrichardson) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply richardson");
871 ierr = (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);CHKERRQ(ierr);
872 PetscFunctionReturn(0);
873 }
874
875 /*@
876 PCSetFailedReason - Sets the reason a PCSetUp() failed or PC_NOERROR if it did not fail
877
878 Logically Collective on PC
879
880 Input Parameter:
881 + pc - the preconditioner context
882 - reason - the reason it failedx
883
884 Level: advanced
885
886 .seealso: PCCreate(), PCApply(), PCDestroy(), PCFailedReason
887 @*/
PCSetFailedReason(PC pc,PCFailedReason reason)888 PetscErrorCode PCSetFailedReason(PC pc,PCFailedReason reason)
889 {
890 PetscFunctionBegin;
891 pc->failedreason = reason;
892 PetscFunctionReturn(0);
893 }
894
895 /*@
896 PCGetFailedReason - Gets the reason a PCSetUp() failed or PC_NOERROR if it did not fail
897
898 Logically Collective on PC
899
900 Input Parameter:
901 . pc - the preconditioner context
902
903 Output Parameter:
904 . reason - the reason it failed
905
906 Level: advanced
907
908 Notes: This is the maximum over reason over all ranks in the PC communicator. It is only valid after
909 a call KSPCheckDot() or KSPCheckNorm() inside a KSPSolve(). It is not valid immediately after a PCSetUp()
910 or PCApply(), then use PCGetFailedReasonRank()
911
912 .seealso: PCCreate(), PCApply(), PCDestroy(), PCGetFailedReasonRank(), PCSetFailedReason()
913 @*/
PCGetFailedReason(PC pc,PCFailedReason * reason)914 PetscErrorCode PCGetFailedReason(PC pc,PCFailedReason *reason)
915 {
916 PetscFunctionBegin;
917 if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
918 else *reason = pc->failedreason;
919 PetscFunctionReturn(0);
920 }
921
922 /*@
923 PCGetFailedReasonRank - Gets the reason a PCSetUp() failed or PC_NOERROR if it did not fail on this MPI rank
924
925 Not Collective on PC
926
927 Input Parameter:
928 . pc - the preconditioner context
929
930 Output Parameter:
931 . reason - the reason it failed
932
933 Notes:
934 Different ranks may have different reasons or no reason, see PCGetFailedReason()
935
936 Level: advanced
937
938 .seealso: PCCreate(), PCApply(), PCDestroy(), PCGetFailedReason(), PCSetFailedReason()
939 @*/
PCGetFailedReasonRank(PC pc,PCFailedReason * reason)940 PetscErrorCode PCGetFailedReasonRank(PC pc,PCFailedReason *reason)
941 {
942 PetscFunctionBegin;
943 if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
944 else *reason = pc->failedreason;
945 PetscFunctionReturn(0);
946 }
947
948 /*
949 a setupcall of 0 indicates never setup,
950 1 indicates has been previously setup
951 -1 indicates a PCSetUp() was attempted and failed
952 */
953 /*@
954 PCSetUp - Prepares for the use of a preconditioner.
955
956 Collective on PC
957
958 Input Parameter:
959 . pc - the preconditioner context
960
961 Level: developer
962
963 .seealso: PCCreate(), PCApply(), PCDestroy()
964 @*/
PCSetUp(PC pc)965 PetscErrorCode PCSetUp(PC pc)
966 {
967 PetscErrorCode ierr;
968 const char *def;
969 PetscObjectState matstate, matnonzerostate;
970
971 PetscFunctionBegin;
972 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
973 if (!pc->mat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");
974
975 if (pc->setupcalled && pc->reusepreconditioner) {
976 ierr = PetscInfo(pc,"Leaving PC with identical preconditioner since reuse preconditioner is set\n");CHKERRQ(ierr);
977 PetscFunctionReturn(0);
978 }
979
980 ierr = PetscObjectStateGet((PetscObject)pc->pmat,&matstate);CHKERRQ(ierr);
981 ierr = MatGetNonzeroState(pc->pmat,&matnonzerostate);CHKERRQ(ierr);
982 if (!pc->setupcalled) {
983 ierr = PetscInfo(pc,"Setting up PC for first time\n");CHKERRQ(ierr);
984 pc->flag = DIFFERENT_NONZERO_PATTERN;
985 } else if (matstate == pc->matstate) {
986 ierr = PetscInfo(pc,"Leaving PC with identical preconditioner since operator is unchanged\n");CHKERRQ(ierr);
987 PetscFunctionReturn(0);
988 } else {
989 if (matnonzerostate > pc->matnonzerostate) {
990 ierr = PetscInfo(pc,"Setting up PC with different nonzero pattern\n");CHKERRQ(ierr);
991 pc->flag = DIFFERENT_NONZERO_PATTERN;
992 } else {
993 ierr = PetscInfo(pc,"Setting up PC with same nonzero pattern\n");CHKERRQ(ierr);
994 pc->flag = SAME_NONZERO_PATTERN;
995 }
996 }
997 pc->matstate = matstate;
998 pc->matnonzerostate = matnonzerostate;
999
1000 if (!((PetscObject)pc)->type_name) {
1001 ierr = PCGetDefaultType_Private(pc,&def);CHKERRQ(ierr);
1002 ierr = PCSetType(pc,def);CHKERRQ(ierr);
1003 }
1004
1005 ierr = MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);CHKERRQ(ierr);
1006 ierr = MatSetErrorIfFailure(pc->mat,pc->erroriffailure);CHKERRQ(ierr);
1007 ierr = PetscLogEventBegin(PC_SetUp,pc,0,0,0);CHKERRQ(ierr);
1008 if (pc->ops->setup) {
1009 ierr = (*pc->ops->setup)(pc);CHKERRQ(ierr);
1010 }
1011 ierr = PetscLogEventEnd(PC_SetUp,pc,0,0,0);CHKERRQ(ierr);
1012 if (!pc->setupcalled) pc->setupcalled = 1;
1013 PetscFunctionReturn(0);
1014 }
1015
1016 /*@
1017 PCSetUpOnBlocks - Sets up the preconditioner for each block in
1018 the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
1019 methods.
1020
1021 Collective on PC
1022
1023 Input Parameters:
1024 . pc - the preconditioner context
1025
1026 Level: developer
1027
1028 .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
1029 @*/
PCSetUpOnBlocks(PC pc)1030 PetscErrorCode PCSetUpOnBlocks(PC pc)
1031 {
1032 PetscErrorCode ierr;
1033
1034 PetscFunctionBegin;
1035 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1036 if (!pc->ops->setuponblocks) PetscFunctionReturn(0);
1037 ierr = PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);CHKERRQ(ierr);
1038 ierr = (*pc->ops->setuponblocks)(pc);CHKERRQ(ierr);
1039 ierr = PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);CHKERRQ(ierr);
1040 PetscFunctionReturn(0);
1041 }
1042
1043 /*@C
1044 PCSetModifySubMatrices - Sets a user-defined routine for modifying the
1045 submatrices that arise within certain subdomain-based preconditioners.
1046 The basic submatrices are extracted from the preconditioner matrix as
1047 usual; the user can then alter these (for example, to set different boundary
1048 conditions for each submatrix) before they are used for the local solves.
1049
1050 Logically Collective on PC
1051
1052 Input Parameters:
1053 + pc - the preconditioner context
1054 . func - routine for modifying the submatrices
1055 - ctx - optional user-defined context (may be null)
1056
1057 Calling sequence of func:
1058 $ func (PC pc,PetscInt nsub,IS *row,IS *col,Mat *submat,void *ctx);
1059
1060 + row - an array of index sets that contain the global row numbers
1061 that comprise each local submatrix
1062 . col - an array of index sets that contain the global column numbers
1063 that comprise each local submatrix
1064 . submat - array of local submatrices
1065 - ctx - optional user-defined context for private data for the
1066 user-defined func routine (may be null)
1067
1068 Notes:
1069 PCSetModifySubMatrices() MUST be called before KSPSetUp() and
1070 KSPSolve().
1071
1072 A routine set by PCSetModifySubMatrices() is currently called within
1073 the block Jacobi (PCBJACOBI) and additive Schwarz (PCASM)
1074 preconditioners. All other preconditioners ignore this routine.
1075
1076 Level: advanced
1077
1078 .seealso: PCModifySubMatrices()
1079 @*/
PCSetModifySubMatrices(PC pc,PetscErrorCode (* func)(PC,PetscInt,const IS[],const IS[],Mat[],void *),void * ctx)1080 PetscErrorCode PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
1081 {
1082 PetscFunctionBegin;
1083 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1084 pc->modifysubmatrices = func;
1085 pc->modifysubmatricesP = ctx;
1086 PetscFunctionReturn(0);
1087 }
1088
1089 /*@C
1090 PCModifySubMatrices - Calls an optional user-defined routine within
1091 certain preconditioners if one has been set with PCSetModifySubMatrices().
1092
1093 Collective on PC
1094
1095 Input Parameters:
1096 + pc - the preconditioner context
1097 . nsub - the number of local submatrices
1098 . row - an array of index sets that contain the global row numbers
1099 that comprise each local submatrix
1100 . col - an array of index sets that contain the global column numbers
1101 that comprise each local submatrix
1102 . submat - array of local submatrices
1103 - ctx - optional user-defined context for private data for the
1104 user-defined routine (may be null)
1105
1106 Output Parameter:
1107 . submat - array of local submatrices (the entries of which may
1108 have been modified)
1109
1110 Notes:
1111 The user should NOT generally call this routine, as it will
1112 automatically be called within certain preconditioners (currently
1113 block Jacobi, additive Schwarz) if set.
1114
1115 The basic submatrices are extracted from the preconditioner matrix
1116 as usual; the user can then alter these (for example, to set different
1117 boundary conditions for each submatrix) before they are used for the
1118 local solves.
1119
1120 Level: developer
1121
1122 .seealso: PCSetModifySubMatrices()
1123 @*/
PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void * ctx)1124 PetscErrorCode PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1125 {
1126 PetscErrorCode ierr;
1127
1128 PetscFunctionBegin;
1129 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1130 if (!pc->modifysubmatrices) PetscFunctionReturn(0);
1131 ierr = PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);CHKERRQ(ierr);
1132 ierr = (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);CHKERRQ(ierr);
1133 ierr = PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);CHKERRQ(ierr);
1134 PetscFunctionReturn(0);
1135 }
1136
1137 /*@
1138 PCSetOperators - Sets the matrix associated with the linear system and
1139 a (possibly) different one associated with the preconditioner.
1140
1141 Logically Collective on PC
1142
1143 Input Parameters:
1144 + pc - the preconditioner context
1145 . Amat - the matrix that defines the linear system
1146 - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1147
1148 Notes:
1149 Passing a NULL for Amat or Pmat removes the matrix that is currently used.
1150
1151 If you wish to replace either Amat or Pmat but leave the other one untouched then
1152 first call KSPGetOperators() to get the one you wish to keep, call PetscObjectReference()
1153 on it and then pass it back in in your call to KSPSetOperators().
1154
1155 More Notes about Repeated Solution of Linear Systems:
1156 PETSc does NOT reset the matrix entries of either Amat or Pmat
1157 to zero after a linear solve; the user is completely responsible for
1158 matrix assembly. See the routine MatZeroEntries() if desiring to
1159 zero all elements of a matrix.
1160
1161 Level: intermediate
1162
1163 .seealso: PCGetOperators(), MatZeroEntries()
1164 @*/
PCSetOperators(PC pc,Mat Amat,Mat Pmat)1165 PetscErrorCode PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1166 {
1167 PetscErrorCode ierr;
1168 PetscInt m1,n1,m2,n2;
1169
1170 PetscFunctionBegin;
1171 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1172 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2);
1173 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3);
1174 if (Amat) PetscCheckSameComm(pc,1,Amat,2);
1175 if (Pmat) PetscCheckSameComm(pc,1,Pmat,3);
1176 if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1177 ierr = MatGetLocalSize(Amat,&m1,&n1);CHKERRQ(ierr);
1178 ierr = MatGetLocalSize(pc->mat,&m2,&n2);CHKERRQ(ierr);
1179 if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot change local size of Amat after use old sizes %D %D new sizes %D %D",m2,n2,m1,n1);
1180 ierr = MatGetLocalSize(Pmat,&m1,&n1);CHKERRQ(ierr);
1181 ierr = MatGetLocalSize(pc->pmat,&m2,&n2);CHKERRQ(ierr);
1182 if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot change local size of Pmat after use old sizes %D %D new sizes %D %D",m2,n2,m1,n1);
1183 }
1184
1185 if (Pmat != pc->pmat) {
1186 /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1187 pc->matnonzerostate = -1;
1188 pc->matstate = -1;
1189 }
1190
1191 /* reference first in case the matrices are the same */
1192 if (Amat) {ierr = PetscObjectReference((PetscObject)Amat);CHKERRQ(ierr);}
1193 ierr = MatDestroy(&pc->mat);CHKERRQ(ierr);
1194 if (Pmat) {ierr = PetscObjectReference((PetscObject)Pmat);CHKERRQ(ierr);}
1195 ierr = MatDestroy(&pc->pmat);CHKERRQ(ierr);
1196 pc->mat = Amat;
1197 pc->pmat = Pmat;
1198 PetscFunctionReturn(0);
1199 }
1200
1201 /*@
1202 PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.
1203
1204 Logically Collective on PC
1205
1206 Input Parameters:
1207 + pc - the preconditioner context
1208 - flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1209
1210 Level: intermediate
1211
1212 .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1213 @*/
PCSetReusePreconditioner(PC pc,PetscBool flag)1214 PetscErrorCode PCSetReusePreconditioner(PC pc,PetscBool flag)
1215 {
1216 PetscFunctionBegin;
1217 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1218 PetscValidLogicalCollectiveBool(pc,flag,2);
1219 pc->reusepreconditioner = flag;
1220 PetscFunctionReturn(0);
1221 }
1222
1223 /*@
1224 PCGetReusePreconditioner - Determines if the PC reuses the current preconditioner even if the operator in the preconditioner has changed.
1225
1226 Not Collective
1227
1228 Input Parameter:
1229 . pc - the preconditioner context
1230
1231 Output Parameter:
1232 . flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1233
1234 Level: intermediate
1235
1236 .seealso: PCGetOperators(), MatZeroEntries(), PCSetReusePreconditioner()
1237 @*/
PCGetReusePreconditioner(PC pc,PetscBool * flag)1238 PetscErrorCode PCGetReusePreconditioner(PC pc,PetscBool *flag)
1239 {
1240 PetscFunctionBegin;
1241 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1242 PetscValidPointer(flag,2);
1243 *flag = pc->reusepreconditioner;
1244 PetscFunctionReturn(0);
1245 }
1246
1247 /*@
1248 PCGetOperators - Gets the matrix associated with the linear system and
1249 possibly a different one associated with the preconditioner.
1250
1251 Not collective, though parallel Mats are returned if the PC is parallel
1252
1253 Input Parameter:
1254 . pc - the preconditioner context
1255
1256 Output Parameters:
1257 + Amat - the matrix defining the linear system
1258 - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1259
1260 Level: intermediate
1261
1262 Notes:
1263 Does not increase the reference count of the matrices, so you should not destroy them
1264
1265 Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
1266 are created in PC and returned to the user. In this case, if both operators
1267 mat and pmat are requested, two DIFFERENT operators will be returned. If
1268 only one is requested both operators in the PC will be the same (i.e. as
1269 if one had called KSP/PCSetOperators() with the same argument for both Mats).
1270 The user must set the sizes of the returned matrices and their type etc just
1271 as if the user created them with MatCreate(). For example,
1272
1273 $ KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1274 $ set size, type, etc of Amat
1275
1276 $ MatCreate(comm,&mat);
1277 $ KSP/PCSetOperators(ksp/pc,Amat,Amat);
1278 $ PetscObjectDereference((PetscObject)mat);
1279 $ set size, type, etc of Amat
1280
1281 and
1282
1283 $ KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1284 $ set size, type, etc of Amat and Pmat
1285
1286 $ MatCreate(comm,&Amat);
1287 $ MatCreate(comm,&Pmat);
1288 $ KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1289 $ PetscObjectDereference((PetscObject)Amat);
1290 $ PetscObjectDereference((PetscObject)Pmat);
1291 $ set size, type, etc of Amat and Pmat
1292
1293 The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
1294 of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
1295 managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
1296 at this is when you create a SNES you do not NEED to create a KSP and attach it to
1297 the SNES object (the SNES object manages it for you). Similarly when you create a KSP
1298 you do not need to attach a PC to it (the KSP object manages the PC object for you).
1299 Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
1300 it can be created for you?
1301
1302
1303 .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1304 @*/
PCGetOperators(PC pc,Mat * Amat,Mat * Pmat)1305 PetscErrorCode PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1306 {
1307 PetscErrorCode ierr;
1308
1309 PetscFunctionBegin;
1310 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1311 if (Amat) {
1312 if (!pc->mat) {
1313 if (pc->pmat && !Pmat) { /* Apmat has been set, but user did not request it, so use for Amat */
1314 pc->mat = pc->pmat;
1315 ierr = PetscObjectReference((PetscObject)pc->mat);CHKERRQ(ierr);
1316 } else { /* both Amat and Pmat are empty */
1317 ierr = MatCreate(PetscObjectComm((PetscObject)pc),&pc->mat);CHKERRQ(ierr);
1318 if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1319 pc->pmat = pc->mat;
1320 ierr = PetscObjectReference((PetscObject)pc->pmat);CHKERRQ(ierr);
1321 }
1322 }
1323 }
1324 *Amat = pc->mat;
1325 }
1326 if (Pmat) {
1327 if (!pc->pmat) {
1328 if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1329 pc->pmat = pc->mat;
1330 ierr = PetscObjectReference((PetscObject)pc->pmat);CHKERRQ(ierr);
1331 } else {
1332 ierr = MatCreate(PetscObjectComm((PetscObject)pc),&pc->pmat);CHKERRQ(ierr);
1333 if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1334 pc->mat = pc->pmat;
1335 ierr = PetscObjectReference((PetscObject)pc->mat);CHKERRQ(ierr);
1336 }
1337 }
1338 }
1339 *Pmat = pc->pmat;
1340 }
1341 PetscFunctionReturn(0);
1342 }
1343
1344 /*@C
1345 PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1346 possibly a different one associated with the preconditioner have been set in the PC.
1347
1348 Not collective, though the results on all processes should be the same
1349
1350 Input Parameter:
1351 . pc - the preconditioner context
1352
1353 Output Parameters:
1354 + mat - the matrix associated with the linear system was set
1355 - pmat - matrix associated with the preconditioner was set, usually the same
1356
1357 Level: intermediate
1358
1359 .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1360 @*/
PCGetOperatorsSet(PC pc,PetscBool * mat,PetscBool * pmat)1361 PetscErrorCode PCGetOperatorsSet(PC pc,PetscBool *mat,PetscBool *pmat)
1362 {
1363 PetscFunctionBegin;
1364 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1365 if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1366 if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1367 PetscFunctionReturn(0);
1368 }
1369
1370 /*@
1371 PCFactorGetMatrix - Gets the factored matrix from the
1372 preconditioner context. This routine is valid only for the LU,
1373 incomplete LU, Cholesky, and incomplete Cholesky methods.
1374
1375 Not Collective on PC though Mat is parallel if PC is parallel
1376
1377 Input Parameters:
1378 . pc - the preconditioner context
1379
1380 Output parameters:
1381 . mat - the factored matrix
1382
1383 Level: advanced
1384
1385 Notes:
1386 Does not increase the reference count for the matrix so DO NOT destroy it
1387
1388 @*/
PCFactorGetMatrix(PC pc,Mat * mat)1389 PetscErrorCode PCFactorGetMatrix(PC pc,Mat *mat)
1390 {
1391 PetscErrorCode ierr;
1392
1393 PetscFunctionBegin;
1394 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1395 PetscValidPointer(mat,2);
1396 if (pc->ops->getfactoredmatrix) {
1397 ierr = (*pc->ops->getfactoredmatrix)(pc,mat);CHKERRQ(ierr);
1398 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1399 PetscFunctionReturn(0);
1400 }
1401
1402 /*@C
1403 PCSetOptionsPrefix - Sets the prefix used for searching for all
1404 PC options in the database.
1405
1406 Logically Collective on PC
1407
1408 Input Parameters:
1409 + pc - the preconditioner context
1410 - prefix - the prefix string to prepend to all PC option requests
1411
1412 Notes:
1413 A hyphen (-) must NOT be given at the beginning of the prefix name.
1414 The first character of all runtime options is AUTOMATICALLY the
1415 hyphen.
1416
1417 Level: advanced
1418
1419 .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1420 @*/
PCSetOptionsPrefix(PC pc,const char prefix[])1421 PetscErrorCode PCSetOptionsPrefix(PC pc,const char prefix[])
1422 {
1423 PetscErrorCode ierr;
1424
1425 PetscFunctionBegin;
1426 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1427 ierr = PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);CHKERRQ(ierr);
1428 PetscFunctionReturn(0);
1429 }
1430
1431 /*@C
1432 PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1433 PC options in the database.
1434
1435 Logically Collective on PC
1436
1437 Input Parameters:
1438 + pc - the preconditioner context
1439 - prefix - the prefix string to prepend to all PC option requests
1440
1441 Notes:
1442 A hyphen (-) must NOT be given at the beginning of the prefix name.
1443 The first character of all runtime options is AUTOMATICALLY the
1444 hyphen.
1445
1446 Level: advanced
1447
1448 .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1449 @*/
PCAppendOptionsPrefix(PC pc,const char prefix[])1450 PetscErrorCode PCAppendOptionsPrefix(PC pc,const char prefix[])
1451 {
1452 PetscErrorCode ierr;
1453
1454 PetscFunctionBegin;
1455 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1456 ierr = PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);CHKERRQ(ierr);
1457 PetscFunctionReturn(0);
1458 }
1459
1460 /*@C
1461 PCGetOptionsPrefix - Gets the prefix used for searching for all
1462 PC options in the database.
1463
1464 Not Collective
1465
1466 Input Parameters:
1467 . pc - the preconditioner context
1468
1469 Output Parameters:
1470 . prefix - pointer to the prefix string used, is returned
1471
1472 Notes:
1473 On the fortran side, the user should pass in a string 'prifix' of
1474 sufficient length to hold the prefix.
1475
1476 Level: advanced
1477
1478 .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1479 @*/
PCGetOptionsPrefix(PC pc,const char * prefix[])1480 PetscErrorCode PCGetOptionsPrefix(PC pc,const char *prefix[])
1481 {
1482 PetscErrorCode ierr;
1483
1484 PetscFunctionBegin;
1485 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1486 PetscValidPointer(prefix,2);
1487 ierr = PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);CHKERRQ(ierr);
1488 PetscFunctionReturn(0);
1489 }
1490
1491 /*
1492 Indicates the right hand side will be changed by KSPSolve(), this occurs for a few
1493 preconditioners including BDDC and Eisentat that transform the equations before applying
1494 the Krylov methods
1495 */
PCPreSolveChangeRHS(PC pc,PetscBool * change)1496 PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc,PetscBool *change)
1497 {
1498 PetscErrorCode ierr;
1499
1500 PetscFunctionBegin;
1501 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1502 PetscValidPointer(change,2);
1503 *change = PETSC_FALSE;
1504 ierr = PetscTryMethod(pc,"PCPreSolveChangeRHS_C",(PC,PetscBool*),(pc,change));CHKERRQ(ierr);
1505 PetscFunctionReturn(0);
1506 }
1507
1508 /*@
1509 PCPreSolve - Optional pre-solve phase, intended for any
1510 preconditioner-specific actions that must be performed before
1511 the iterative solve itself.
1512
1513 Collective on PC
1514
1515 Input Parameters:
1516 + pc - the preconditioner context
1517 - ksp - the Krylov subspace context
1518
1519 Level: developer
1520
1521 Sample of Usage:
1522 .vb
1523 PCPreSolve(pc,ksp);
1524 KSPSolve(ksp,b,x);
1525 PCPostSolve(pc,ksp);
1526 .ve
1527
1528 Notes:
1529 The pre-solve phase is distinct from the PCSetUp() phase.
1530
1531 KSPSolve() calls this directly, so is rarely called by the user.
1532
1533 .seealso: PCPostSolve()
1534 @*/
PCPreSolve(PC pc,KSP ksp)1535 PetscErrorCode PCPreSolve(PC pc,KSP ksp)
1536 {
1537 PetscErrorCode ierr;
1538 Vec x,rhs;
1539
1540 PetscFunctionBegin;
1541 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1542 PetscValidHeaderSpecific(ksp,KSP_CLASSID,2);
1543 pc->presolvedone++;
1544 if (pc->presolvedone > 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot embed PCPreSolve() more than twice");
1545 ierr = KSPGetSolution(ksp,&x);CHKERRQ(ierr);
1546 ierr = KSPGetRhs(ksp,&rhs);CHKERRQ(ierr);
1547
1548 if (pc->ops->presolve) {
1549 ierr = (*pc->ops->presolve)(pc,ksp,rhs,x);CHKERRQ(ierr);
1550 }
1551 PetscFunctionReturn(0);
1552 }
1553
1554 /*@
1555 PCPostSolve - Optional post-solve phase, intended for any
1556 preconditioner-specific actions that must be performed after
1557 the iterative solve itself.
1558
1559 Collective on PC
1560
1561 Input Parameters:
1562 + pc - the preconditioner context
1563 - ksp - the Krylov subspace context
1564
1565 Sample of Usage:
1566 .vb
1567 PCPreSolve(pc,ksp);
1568 KSPSolve(ksp,b,x);
1569 PCPostSolve(pc,ksp);
1570 .ve
1571
1572 Note:
1573 KSPSolve() calls this routine directly, so it is rarely called by the user.
1574
1575 Level: developer
1576
1577 .seealso: PCPreSolve(), KSPSolve()
1578 @*/
PCPostSolve(PC pc,KSP ksp)1579 PetscErrorCode PCPostSolve(PC pc,KSP ksp)
1580 {
1581 PetscErrorCode ierr;
1582 Vec x,rhs;
1583
1584 PetscFunctionBegin;
1585 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1586 PetscValidHeaderSpecific(ksp,KSP_CLASSID,2);
1587 pc->presolvedone--;
1588 ierr = KSPGetSolution(ksp,&x);CHKERRQ(ierr);
1589 ierr = KSPGetRhs(ksp,&rhs);CHKERRQ(ierr);
1590 if (pc->ops->postsolve) {
1591 ierr = (*pc->ops->postsolve)(pc,ksp,rhs,x);CHKERRQ(ierr);
1592 }
1593 PetscFunctionReturn(0);
1594 }
1595
1596 /*@C
1597 PCLoad - Loads a PC that has been stored in binary with PCView().
1598
1599 Collective on PetscViewer
1600
1601 Input Parameters:
1602 + newdm - the newly loaded PC, this needs to have been created with PCCreate() or
1603 some related function before a call to PCLoad().
1604 - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
1605
1606 Level: intermediate
1607
1608 Notes:
1609 The type is determined by the data in the file, any type set into the PC before this call is ignored.
1610
1611 Notes for advanced users:
1612 Most users should not need to know the details of the binary storage
1613 format, since PCLoad() and PCView() completely hide these details.
1614 But for anyone who's interested, the standard binary matrix storage
1615 format is
1616 .vb
1617 has not yet been determined
1618 .ve
1619
1620 .seealso: PetscViewerBinaryOpen(), PCView(), MatLoad(), VecLoad()
1621 @*/
PCLoad(PC newdm,PetscViewer viewer)1622 PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1623 {
1624 PetscErrorCode ierr;
1625 PetscBool isbinary;
1626 PetscInt classid;
1627 char type[256];
1628
1629 PetscFunctionBegin;
1630 PetscValidHeaderSpecific(newdm,PC_CLASSID,1);
1631 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
1632 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
1633 if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1634
1635 ierr = PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);CHKERRQ(ierr);
1636 if (classid != PC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not PC next in file");
1637 ierr = PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);CHKERRQ(ierr);
1638 ierr = PCSetType(newdm, type);CHKERRQ(ierr);
1639 if (newdm->ops->load) {
1640 ierr = (*newdm->ops->load)(newdm,viewer);CHKERRQ(ierr);
1641 }
1642 PetscFunctionReturn(0);
1643 }
1644
1645 #include <petscdraw.h>
1646 #if defined(PETSC_HAVE_SAWS)
1647 #include <petscviewersaws.h>
1648 #endif
1649
1650 /*@C
1651 PCViewFromOptions - View from Options
1652
1653 Collective on PC
1654
1655 Input Parameters:
1656 + A - the PC context
1657 . obj - Optional object
1658 - name - command line option
1659
1660 Level: intermediate
1661 .seealso: PC, PCView, PetscObjectViewFromOptions(), PCCreate()
1662 @*/
PCViewFromOptions(PC A,PetscObject obj,const char name[])1663 PetscErrorCode PCViewFromOptions(PC A,PetscObject obj,const char name[])
1664 {
1665 PetscErrorCode ierr;
1666
1667 PetscFunctionBegin;
1668 PetscValidHeaderSpecific(A,PC_CLASSID,1);
1669 ierr = PetscObjectViewFromOptions((PetscObject)A,obj,name);CHKERRQ(ierr);
1670 PetscFunctionReturn(0);
1671 }
1672
1673 /*@C
1674 PCView - Prints the PC data structure.
1675
1676 Collective on PC
1677
1678 Input Parameters:
1679 + PC - the PC context
1680 - viewer - optional visualization context
1681
1682 Note:
1683 The available visualization contexts include
1684 + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1685 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1686 output where only the first processor opens
1687 the file. All other processors send their
1688 data to the first processor to print.
1689
1690 The user can open an alternative visualization contexts with
1691 PetscViewerASCIIOpen() (output to a specified file).
1692
1693 Level: developer
1694
1695 .seealso: KSPView(), PetscViewerASCIIOpen()
1696 @*/
PCView(PC pc,PetscViewer viewer)1697 PetscErrorCode PCView(PC pc,PetscViewer viewer)
1698 {
1699 PCType cstr;
1700 PetscErrorCode ierr;
1701 PetscBool iascii,isstring,isbinary,isdraw;
1702 #if defined(PETSC_HAVE_SAWS)
1703 PetscBool issaws;
1704 #endif
1705
1706 PetscFunctionBegin;
1707 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1708 if (!viewer) {
1709 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);CHKERRQ(ierr);
1710 }
1711 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
1712 PetscCheckSameComm(pc,1,viewer,2);
1713
1714 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
1715 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr);
1716 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
1717 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
1718 #if defined(PETSC_HAVE_SAWS)
1719 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);CHKERRQ(ierr);
1720 #endif
1721
1722 if (iascii) {
1723 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)pc,viewer);CHKERRQ(ierr);
1724 if (!pc->setupcalled) {
1725 ierr = PetscViewerASCIIPrintf(viewer," PC has not been set up so information may be incomplete\n");CHKERRQ(ierr);
1726 }
1727 if (pc->ops->view) {
1728 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
1729 ierr = (*pc->ops->view)(pc,viewer);CHKERRQ(ierr);
1730 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
1731 }
1732 if (pc->mat) {
1733 ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);CHKERRQ(ierr);
1734 if (pc->pmat == pc->mat) {
1735 ierr = PetscViewerASCIIPrintf(viewer," linear system matrix = precond matrix:\n");CHKERRQ(ierr);
1736 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
1737 ierr = MatView(pc->mat,viewer);CHKERRQ(ierr);
1738 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
1739 } else {
1740 if (pc->pmat) {
1741 ierr = PetscViewerASCIIPrintf(viewer," linear system matrix followed by preconditioner matrix:\n");CHKERRQ(ierr);
1742 } else {
1743 ierr = PetscViewerASCIIPrintf(viewer," linear system matrix:\n");CHKERRQ(ierr);
1744 }
1745 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
1746 ierr = MatView(pc->mat,viewer);CHKERRQ(ierr);
1747 if (pc->pmat) {ierr = MatView(pc->pmat,viewer);CHKERRQ(ierr);}
1748 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
1749 }
1750 ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
1751 }
1752 } else if (isstring) {
1753 ierr = PCGetType(pc,&cstr);CHKERRQ(ierr);
1754 ierr = PetscViewerStringSPrintf(viewer," PCType: %-7.7s",cstr);CHKERRQ(ierr);
1755 if (pc->ops->view) {ierr = (*pc->ops->view)(pc,viewer);CHKERRQ(ierr);}
1756 if (pc->mat) {ierr = MatView(pc->mat,viewer);CHKERRQ(ierr);}
1757 if (pc->pmat && pc->pmat != pc->mat) {ierr = MatView(pc->pmat,viewer);CHKERRQ(ierr);}
1758 } else if (isbinary) {
1759 PetscInt classid = PC_FILE_CLASSID;
1760 MPI_Comm comm;
1761 PetscMPIInt rank;
1762 char type[256];
1763
1764 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
1765 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
1766 if (!rank) {
1767 ierr = PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT);CHKERRQ(ierr);
1768 ierr = PetscStrncpy(type,((PetscObject)pc)->type_name,256);CHKERRQ(ierr);
1769 ierr = PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR);CHKERRQ(ierr);
1770 }
1771 if (pc->ops->view) {
1772 ierr = (*pc->ops->view)(pc,viewer);CHKERRQ(ierr);
1773 }
1774 } else if (isdraw) {
1775 PetscDraw draw;
1776 char str[25];
1777 PetscReal x,y,bottom,h;
1778 PetscInt n;
1779
1780 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
1781 ierr = PetscDrawGetCurrentPoint(draw,&x,&y);CHKERRQ(ierr);
1782 if (pc->mat) {
1783 ierr = MatGetSize(pc->mat,&n,NULL);CHKERRQ(ierr);
1784 ierr = PetscSNPrintf(str,25,"PC: %s (%D)",((PetscObject)pc)->type_name,n);CHKERRQ(ierr);
1785 } else {
1786 ierr = PetscSNPrintf(str,25,"PC: %s",((PetscObject)pc)->type_name);CHKERRQ(ierr);
1787 }
1788 ierr = PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);CHKERRQ(ierr);
1789 bottom = y - h;
1790 ierr = PetscDrawPushCurrentPoint(draw,x,bottom);CHKERRQ(ierr);
1791 if (pc->ops->view) {
1792 ierr = (*pc->ops->view)(pc,viewer);CHKERRQ(ierr);
1793 }
1794 ierr = PetscDrawPopCurrentPoint(draw);CHKERRQ(ierr);
1795 #if defined(PETSC_HAVE_SAWS)
1796 } else if (issaws) {
1797 PetscMPIInt rank;
1798
1799 ierr = PetscObjectName((PetscObject)pc);CHKERRQ(ierr);
1800 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
1801 if (!((PetscObject)pc)->amsmem && !rank) {
1802 ierr = PetscObjectViewSAWs((PetscObject)pc,viewer);CHKERRQ(ierr);
1803 }
1804 if (pc->mat) {ierr = MatView(pc->mat,viewer);CHKERRQ(ierr);}
1805 if (pc->pmat && pc->pmat != pc->mat) {ierr = MatView(pc->pmat,viewer);CHKERRQ(ierr);}
1806 #endif
1807 }
1808 PetscFunctionReturn(0);
1809 }
1810
1811 /*@C
1812 PCRegister - Adds a method to the preconditioner package.
1813
1814 Not collective
1815
1816 Input Parameters:
1817 + name_solver - name of a new user-defined solver
1818 - routine_create - routine to create method context
1819
1820 Notes:
1821 PCRegister() may be called multiple times to add several user-defined preconditioners.
1822
1823 Sample usage:
1824 .vb
1825 PCRegister("my_solver", MySolverCreate);
1826 .ve
1827
1828 Then, your solver can be chosen with the procedural interface via
1829 $ PCSetType(pc,"my_solver")
1830 or at runtime via the option
1831 $ -pc_type my_solver
1832
1833 Level: advanced
1834
1835 .seealso: PCRegisterAll()
1836 @*/
PCRegister(const char sname[],PetscErrorCode (* function)(PC))1837 PetscErrorCode PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1838 {
1839 PetscErrorCode ierr;
1840
1841 PetscFunctionBegin;
1842 ierr = PCInitializePackage();CHKERRQ(ierr);
1843 ierr = PetscFunctionListAdd(&PCList,sname,function);CHKERRQ(ierr);
1844 PetscFunctionReturn(0);
1845 }
1846
MatMult_PC(Mat A,Vec X,Vec Y)1847 static PetscErrorCode MatMult_PC(Mat A,Vec X,Vec Y)
1848 {
1849 PC pc;
1850 PetscErrorCode ierr;
1851
1852 PetscFunctionBegin;
1853 ierr = MatShellGetContext(A,&pc);CHKERRQ(ierr);
1854 ierr = PCApply(pc,X,Y);CHKERRQ(ierr);
1855 PetscFunctionReturn(0);
1856 }
1857
1858 /*@
1859 PCComputeOperator - Computes the explicit preconditioned operator.
1860
1861 Collective on PC
1862
1863 Input Parameter:
1864 + pc - the preconditioner object
1865 - mattype - the matrix type to be used for the operator
1866
1867 Output Parameter:
1868 . mat - the explict preconditioned operator
1869
1870 Notes:
1871 This computation is done by applying the operators to columns of the identity matrix.
1872 This routine is costly in general, and is recommended for use only with relatively small systems.
1873 Currently, this routine uses a dense matrix format when mattype == NULL
1874
1875 Level: advanced
1876
1877 .seealso: KSPComputeOperator(), MatType
1878
1879 @*/
PCComputeOperator(PC pc,MatType mattype,Mat * mat)1880 PetscErrorCode PCComputeOperator(PC pc,MatType mattype,Mat *mat)
1881 {
1882 PetscErrorCode ierr;
1883 PetscInt N,M,m,n;
1884 Mat A,Apc;
1885
1886 PetscFunctionBegin;
1887 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1888 PetscValidPointer(mat,3);
1889 ierr = PCGetOperators(pc,&A,NULL);CHKERRQ(ierr);
1890 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
1891 ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr);
1892 ierr = MatCreateShell(PetscObjectComm((PetscObject)pc),m,n,M,N,pc,&Apc);CHKERRQ(ierr);
1893 ierr = MatShellSetOperation(Apc,MATOP_MULT,(void (*)(void))MatMult_PC);CHKERRQ(ierr);
1894 ierr = MatComputeOperator(Apc,mattype,mat);CHKERRQ(ierr);
1895 ierr = MatDestroy(&Apc);CHKERRQ(ierr);
1896 PetscFunctionReturn(0);
1897 }
1898
1899 /*@
1900 PCSetCoordinates - sets the coordinates of all the nodes on the local process
1901
1902 Collective on PC
1903
1904 Input Parameters:
1905 + pc - the solver context
1906 . dim - the dimension of the coordinates 1, 2, or 3
1907 . nloc - the blocked size of the coordinates array
1908 - coords - the coordinates array
1909
1910 Level: intermediate
1911
1912 Notes:
1913 coords is an array of the dim coordinates for the nodes on
1914 the local processor, of size dim*nloc.
1915 If there are 108 equation on a processor
1916 for a displacement finite element discretization of elasticity (so
1917 that there are nloc = 36 = 108/3 nodes) then the array must have 108
1918 double precision values (ie, 3 * 36). These x y z coordinates
1919 should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1920 ... , N-1.z ].
1921
1922 .seealso: MatSetNearNullSpace()
1923 @*/
PCSetCoordinates(PC pc,PetscInt dim,PetscInt nloc,PetscReal coords[])1924 PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
1925 {
1926 PetscErrorCode ierr;
1927
1928 PetscFunctionBegin;
1929 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1930 PetscValidLogicalCollectiveInt(pc,dim,2);
1931 ierr = PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));CHKERRQ(ierr);
1932 PetscFunctionReturn(0);
1933 }
1934
1935 /*@
1936 PCGetInterpolations - Gets interpolation matrices for all levels (except level 0)
1937
1938 Logically Collective on PC
1939
1940 Input Parameters:
1941 + pc - the precondition context
1942
1943 Output Parameter:
1944 - num_levels - the number of levels
1945 . interpolations - the interpolation matrices (size of num_levels-1)
1946
1947 Level: advanced
1948
1949 .keywords: MG, GAMG, BoomerAMG, multigrid, interpolation, level
1950
1951 .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetInterpolation(), PCGetCoarseOperators()
1952 @*/
PCGetInterpolations(PC pc,PetscInt * num_levels,Mat * interpolations[])1953 PetscErrorCode PCGetInterpolations(PC pc,PetscInt *num_levels,Mat *interpolations[])
1954 {
1955 PetscErrorCode ierr;
1956
1957 PetscFunctionBegin;
1958 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1959 PetscValidPointer(num_levels,2);
1960 PetscValidPointer(interpolations,3);
1961 ierr = PetscUseMethod(pc,"PCGetInterpolations_C",(PC,PetscInt*,Mat*[]),(pc,num_levels,interpolations));CHKERRQ(ierr);
1962 PetscFunctionReturn(0);
1963 }
1964
1965 /*@
1966 PCGetCoarseOperators - Gets coarse operator matrices for all levels (except the finest level)
1967
1968 Logically Collective on PC
1969
1970 Input Parameters:
1971 + pc - the precondition context
1972
1973 Output Parameter:
1974 - num_levels - the number of levels
1975 . coarseOperators - the coarse operator matrices (size of num_levels-1)
1976
1977 Level: advanced
1978
1979 .keywords: MG, GAMG, BoomerAMG, get, multigrid, interpolation, level
1980
1981 .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetRScale(), PCMGGetInterpolation(), PCGetInterpolations()
1982 @*/
PCGetCoarseOperators(PC pc,PetscInt * num_levels,Mat * coarseOperators[])1983 PetscErrorCode PCGetCoarseOperators(PC pc,PetscInt *num_levels,Mat *coarseOperators[])
1984 {
1985 PetscErrorCode ierr;
1986
1987 PetscFunctionBegin;
1988 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1989 PetscValidPointer(num_levels,2);
1990 PetscValidPointer(coarseOperators,3);
1991 ierr = PetscUseMethod(pc,"PCGetCoarseOperators_C",(PC,PetscInt*,Mat*[]),(pc,num_levels,coarseOperators));CHKERRQ(ierr);
1992 PetscFunctionReturn(0);
1993 }
1994