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