xref: /petsc/src/ksp/ksp/interface/itcreate.c (revision d78dad28c753525b71e92287f5cd5b2dad6bd135)
1 
2 /*
3      The basic KSP routines, Create, View etc. are here.
4 */
5 #include <petsc-private/kspimpl.h>      /*I "petscksp.h" I*/
6 
7 /* Logging support */
8 PetscClassId  KSP_CLASSID;
9 PetscClassId  DMKSP_CLASSID;
10 PetscLogEvent  KSP_GMRESOrthogonalization, KSP_SetUp, KSP_Solve;
11 
12 /*
13    Contains the list of registered KSP routines
14 */
15 PetscFList KSPList = 0;
16 PetscBool  KSPRegisterAllCalled = PETSC_FALSE;
17 
18 #undef __FUNCT__
19 #define __FUNCT__ "KSPLoad"
20 /*@C
21   KSPLoad - Loads a KSP that has been stored in binary  with KSPView().
22 
23   Collective on PetscViewer
24 
25   Input Parameters:
26 + newdm - the newly loaded KSP, this needs to have been created with KSPCreate() or
27            some related function before a call to KSPLoad().
28 - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
29 
30    Level: intermediate
31 
32   Notes:
33    The type is determined by the data in the file, any type set into the KSP before this call is ignored.
34 
35   Notes for advanced users:
36   Most users should not need to know the details of the binary storage
37   format, since KSPLoad() and KSPView() completely hide these details.
38   But for anyone who's interested, the standard binary matrix storage
39   format is
40 .vb
41      has not yet been determined
42 .ve
43 
44 .seealso: PetscViewerBinaryOpen(), KSPView(), MatLoad(), VecLoad()
45 @*/
46 PetscErrorCode  KSPLoad(KSP newdm, PetscViewer viewer)
47 {
48   PetscErrorCode ierr;
49   PetscBool      isbinary;
50   PetscInt       classid;
51   char           type[256];
52   PC             pc;
53 
54   PetscFunctionBegin;
55   PetscValidHeaderSpecific(newdm,KSP_CLASSID,1);
56   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
57   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
58   if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
59 
60   ierr = PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);CHKERRQ(ierr);
61   if (classid != KSP_FILE_CLASSID) SETERRQ(((PetscObject)newdm)->comm,PETSC_ERR_ARG_WRONG,"Not KSP next in file");
62   ierr = PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);CHKERRQ(ierr);
63   ierr = KSPSetType(newdm, type);CHKERRQ(ierr);
64   if (newdm->ops->load) {
65     ierr = (*newdm->ops->load)(newdm,viewer);CHKERRQ(ierr);
66   }
67   ierr = KSPGetPC(newdm,&pc);CHKERRQ(ierr);
68   ierr = PCLoad(pc,viewer);CHKERRQ(ierr);
69   PetscFunctionReturn(0);
70 }
71 
72 #undef __FUNCT__
73 #define __FUNCT__ "KSPView"
74 /*@C
75    KSPView - Prints the KSP data structure.
76 
77    Collective on KSP
78 
79    Input Parameters:
80 +  ksp - the Krylov space context
81 -  viewer - visualization context
82 
83    Options Database Keys:
84 .  -ksp_view - print the ksp data structure at the end of a KSPSolve call
85 
86    Note:
87    The available visualization contexts include
88 +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
89 -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
90          output where only the first processor opens
91          the file.  All other processors send their
92          data to the first processor to print.
93 
94    The user can open an alternative visualization context with
95    PetscViewerASCIIOpen() - output to a specified file.
96 
97    Level: beginner
98 
99 .keywords: KSP, view
100 
101 .seealso: PCView(), PetscViewerASCIIOpen()
102 @*/
103 PetscErrorCode  KSPView(KSP ksp,PetscViewer viewer)
104 {
105   PetscErrorCode ierr;
106   PetscBool      iascii,isbinary;
107 
108   PetscFunctionBegin;
109   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
110   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm);
111   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
112   PetscCheckSameComm(ksp,1,viewer,2);
113 
114   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
115   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
116   if (iascii) {
117     ierr = PetscObjectPrintClassNamePrefixType((PetscObject)ksp,viewer,"KSP Object");CHKERRQ(ierr);
118     if (ksp->ops->view) {
119       ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
120       ierr = (*ksp->ops->view)(ksp,viewer);CHKERRQ(ierr);
121       ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
122     }
123     if (ksp->guess_zero) {ierr = PetscViewerASCIIPrintf(viewer,"  maximum iterations=%D, initial guess is zero\n",ksp->max_it);CHKERRQ(ierr);}
124     else                 {ierr = PetscViewerASCIIPrintf(viewer,"  maximum iterations=%D\n", ksp->max_it);CHKERRQ(ierr);}
125     if (ksp->guess_knoll) {ierr = PetscViewerASCIIPrintf(viewer,"  using preconditioner applied to right hand side for initial guess\n");CHKERRQ(ierr);}
126     ierr = PetscViewerASCIIPrintf(viewer,"  tolerances:  relative=%G, absolute=%G, divergence=%G\n",ksp->rtol,ksp->abstol,ksp->divtol);CHKERRQ(ierr);
127     if (ksp->pc_side == PC_RIGHT)          {ierr = PetscViewerASCIIPrintf(viewer,"  right preconditioning\n");CHKERRQ(ierr);}
128     else if (ksp->pc_side == PC_SYMMETRIC) {ierr = PetscViewerASCIIPrintf(viewer,"  symmetric preconditioning\n");CHKERRQ(ierr);}
129     else                                   {ierr = PetscViewerASCIIPrintf(viewer,"  left preconditioning\n");CHKERRQ(ierr);}
130     if (ksp->guess) {ierr = PetscViewerASCIIPrintf(viewer,"  using Fischers initial guess method %D with size %D\n",ksp->guess->method,ksp->guess->maxl);CHKERRQ(ierr);}
131     if (ksp->dscale) {ierr = PetscViewerASCIIPrintf(viewer,"  diagonally scaled system\n");CHKERRQ(ierr);}
132     if (ksp->nullsp) {ierr = PetscViewerASCIIPrintf(viewer,"  has attached null space\n");CHKERRQ(ierr);}
133     if (!ksp->guess_zero) {ierr = PetscViewerASCIIPrintf(viewer,"  using nonzero initial guess\n");CHKERRQ(ierr);}
134     ierr = PetscViewerASCIIPrintf(viewer,"  using %s norm type for convergence test\n",KSPNormTypes[ksp->normtype]);CHKERRQ(ierr);
135   } else if (isbinary) {
136     PetscInt         classid = KSP_FILE_CLASSID;
137     MPI_Comm         comm;
138     PetscMPIInt      rank;
139     char             type[256];
140 
141     ierr = PetscObjectGetComm((PetscObject)ksp,&comm);CHKERRQ(ierr);
142     ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
143     if (!rank) {
144       ierr = PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr);
145       ierr = PetscStrncpy(type,((PetscObject)ksp)->type_name,256);CHKERRQ(ierr);
146       ierr = PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);CHKERRQ(ierr);
147     }
148     if (ksp->ops->view) {
149       ierr = (*ksp->ops->view)(ksp,viewer);CHKERRQ(ierr);
150     }
151   } else {
152     if (ksp->ops->view) {
153       ierr = (*ksp->ops->view)(ksp,viewer);CHKERRQ(ierr);
154     }
155   }
156   if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
157   ierr = PCView(ksp->pc,viewer);CHKERRQ(ierr);
158   PetscFunctionReturn(0);
159 }
160 
161 
162 #undef __FUNCT__
163 #define __FUNCT__ "KSPSetNormType"
164 /*@
165    KSPSetNormType - Sets the norm that is used for convergence testing.
166 
167    Logically Collective on KSP
168 
169    Input Parameter:
170 +  ksp - Krylov solver context
171 -  normtype - one of
172 $   KSP_NORM_NONE - skips computing the norm, this should only be used if you are using
173 $                 the Krylov method as a smoother with a fixed small number of iterations.
174 $                 Implicitly sets KSPSkipConverged as KSP convergence test.
175 $                 Supported only by CG, Richardson, Bi-CG-stab, CR, and CGS methods.
176 $   KSP_NORM_PRECONDITIONED - the default for left preconditioned solves, uses the l2 norm
177 $                 of the preconditioned residual
178 $   KSP_NORM_UNPRECONDITIONED - uses the l2 norm of the true b - Ax residual, supported only by
179 $                 CG, CHEBYSHEV, and RICHARDSON, automatically true for right (see KSPSetPCSide())
180 $                 preconditioning..
181 $   KSP_NORM_NATURAL - supported  by KSPCG, KSPCR, KSPCGNE, KSPCGS
182 
183 
184    Options Database Key:
185 .   -ksp_norm_type <none,preconditioned,unpreconditioned,natural>
186 
187    Notes:
188    Currently only works with the CG, Richardson, Bi-CG-stab, CR, and CGS methods.
189 
190    Level: advanced
191 
192 .keywords: KSP, create, context, norms
193 
194 .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPSkipConverged(), KSPSetCheckNormIteration()
195 @*/
196 PetscErrorCode  KSPSetNormType(KSP ksp,KSPNormType normtype)
197 {
198   PetscErrorCode ierr;
199 
200   PetscFunctionBegin;
201   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
202   PetscValidLogicalCollectiveEnum(ksp,normtype,2);
203   ksp->normtype = normtype;
204   if (normtype == KSP_NORM_NONE) {
205     ierr = KSPSetConvergenceTest(ksp,KSPSkipConverged,0,0);CHKERRQ(ierr);
206     ierr = PetscInfo(ksp,"Warning: setting KSPNormType to skip computing the norm\n\
207  KSP convergence test is implicitly set to KSPSkipConverged\n");CHKERRQ(ierr);
208   }
209   PetscFunctionReturn(0);
210 }
211 
212 #undef __FUNCT__
213 #define __FUNCT__ "KSPSetCheckNormIteration"
214 /*@
215    KSPSetCheckNormIteration - Sets the first iteration at which the norm of the residual will be
216      computed and used in the convergence test.
217 
218    Logically Collective on KSP
219 
220    Input Parameter:
221 +  ksp - Krylov solver context
222 -  it  - use -1 to check at all iterations
223 
224    Notes:
225    Currently only works with KSPCG, KSPBCGS and KSPIBCGS
226 
227    Use KSPSetNormType(ksp,KSP_NORM_NONE) to never check the norm
228 
229    On steps where the norm is not computed, the previous norm is still in the variable, so if you run with, for example,
230     -ksp_monitor the residual norm will appear to be unchanged for several iterations (though it is not really unchanged).
231    Level: advanced
232 
233 .keywords: KSP, create, context, norms
234 
235 .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPSkipConverged(), KSPSetNormType()
236 @*/
237 PetscErrorCode  KSPSetCheckNormIteration(KSP ksp,PetscInt it)
238 {
239   PetscFunctionBegin;
240   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
241   PetscValidLogicalCollectiveInt(ksp,it,2);
242   ksp->chknorm = it;
243   PetscFunctionReturn(0);
244 }
245 
246 #undef __FUNCT__
247 #define __FUNCT__ "KSPSetLagNorm"
248 /*@
249    KSPSetLagNorm - Lags the residual norm calculation so that it is computed as part of the MPI_Allreduce() for
250    computing the inner products for the next iteration.  This can reduce communication costs at the expense of doing
251    one additional iteration.
252 
253 
254    Logically Collective on KSP
255 
256    Input Parameter:
257 +  ksp - Krylov solver context
258 -  flg - PETSC_TRUE or PETSC_FALSE
259 
260    Options Database Keys:
261 .  -ksp_lag_norm - lag the calculated residual norm
262 
263    Notes:
264    Currently only works with KSPIBCGS.
265 
266    Use KSPSetNormType(ksp,KSP_NORM_NONE) to never check the norm
267 
268    If you lag the norm and run with, for example, -ksp_monitor, the residual norm reported will be the lagged one.
269    Level: advanced
270 
271 .keywords: KSP, create, context, norms
272 
273 .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPSkipConverged(), KSPSetNormType(), KSPSetCheckNormIteration()
274 @*/
275 PetscErrorCode  KSPSetLagNorm(KSP ksp,PetscBool  flg)
276 {
277   PetscFunctionBegin;
278   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
279   PetscValidLogicalCollectiveBool(ksp,flg,2);
280   ksp->lagnorm = flg;
281   PetscFunctionReturn(0);
282 }
283 
284 #undef __FUNCT__
285 #define __FUNCT__ "KSPSetSupportedNorm"
286 /*@
287    KSPSetSupportedNorm - Sets a norm and preconditioner side supported by a KSP
288 
289    Logically Collective
290 
291    Input Arguments:
292 +  ksp - Krylov method
293 .  normtype - supported norm type
294 .  pcside - preconditioner side that can be used with this norm
295 -  preference - integer preference for this combination, larger values have higher priority
296 
297    Level: developer
298 
299    Notes:
300    This function should be called from the implementation files KSPCreate_XXX() to declare
301    which norms and preconditioner sides are supported. Users should not need to call this
302    function.
303 
304    KSP_NORM_NONE is supported by default with all KSP methods and any PC side. If a KSP explicitly does not support
305    KSP_NORM_NONE, it should set this by setting priority=0.
306 
307 .seealso: KSPSetNormType(), KSPSetPCSide()
308 @*/
309 PetscErrorCode KSPSetSupportedNorm(KSP ksp,KSPNormType normtype,PCSide pcside,PetscInt priority)
310 {
311 
312   PetscFunctionBegin;
313   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
314   ksp->normsupporttable[normtype][pcside] = priority;
315   PetscFunctionReturn(0);
316 }
317 
318 #undef __FUNCT__
319 #define __FUNCT__ "KSPNormSupportTableReset_Private"
320 PetscErrorCode KSPNormSupportTableReset_Private(KSP ksp)
321 {
322   PetscErrorCode ierr;
323 
324   PetscFunctionBegin;
325   ierr = PetscMemzero(ksp->normsupporttable,sizeof(ksp->normsupporttable));CHKERRQ(ierr);
326   ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);CHKERRQ(ierr);
327   ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);CHKERRQ(ierr);
328   PetscFunctionReturn(0);
329 }
330 
331 #undef __FUNCT__
332 #define __FUNCT__ "KSPSetUpNorms_Private"
333 PetscErrorCode KSPSetUpNorms_Private(KSP ksp,KSPNormType *normtype,PCSide *pcside)
334 {
335   PetscInt i,j,best,ibest = 0,jbest = 0;
336 
337   PetscFunctionBegin;
338   best = 0;
339   for (i=0; i<KSP_NORM_MAX; i++) {
340     for (j=0; j<PC_SIDE_MAX; j++) {
341       if ((ksp->normtype == KSP_NORM_DEFAULT || ksp->normtype == i)
342           && (ksp->pc_side == PC_SIDE_DEFAULT || ksp->pc_side == j)
343           && (ksp->normsupporttable[i][j] > best)) {
344         if (ksp->normtype == KSP_NORM_DEFAULT && i == KSP_NORM_NONE && ksp->normsupporttable[i][j] <= 1)
345           continue; /* Skip because we don't want to default to no norms unless set by the KSP (preonly). */
346         best = ksp->normsupporttable[i][j];
347         ibest = i;
348         jbest = j;
349       }
350     }
351   }
352   if (best < 1) {
353     if (ksp->normtype == KSP_NORM_DEFAULT && ksp->pc_side == PC_SIDE_DEFAULT) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_PLIB,"The %s KSP implementation did not call KSPSetSupportedNorm()",((PetscObject)ksp)->type_name);
354     if (ksp->normtype == KSP_NORM_DEFAULT) SETERRQ2(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"KSP %s does not support %s",((PetscObject)ksp)->type_name,PCSides[ksp->pc_side]);
355     if (ksp->pc_side == PC_SIDE_DEFAULT) SETERRQ2(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"KSP %s does not support %s",((PetscObject)ksp)->type_name,KSPNormTypes[ksp->normtype]);
356     SETERRQ3(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"KSP %s does not support %s with %s",((PetscObject)ksp)->type_name,KSPNormTypes[ksp->normtype],PCSides[ksp->pc_side]);
357   }
358   *normtype = (KSPNormType)ibest;
359   *pcside = (PCSide)jbest;
360   PetscFunctionReturn(0);
361 }
362 
363 #undef __FUNCT__
364 #define __FUNCT__ "KSPGetNormType"
365 /*@
366    KSPGetNormType - Gets the norm that is used for convergence testing.
367 
368    Not Collective
369 
370    Input Parameter:
371 .  ksp - Krylov solver context
372 
373    Output Parameter:
374 .  normtype - norm that is used for convergence testing
375 
376    Level: advanced
377 
378 .keywords: KSP, create, context, norms
379 
380 .seealso: KSPNormType, KSPSetNormType(), KSPSkipConverged()
381 @*/
382 PetscErrorCode  KSPGetNormType(KSP ksp, KSPNormType *normtype)
383 {
384   PetscErrorCode ierr;
385 
386   PetscFunctionBegin;
387   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
388   PetscValidPointer(normtype,2);
389   ierr = KSPSetUpNorms_Private(ksp,&ksp->normtype,&ksp->pc_side);CHKERRQ(ierr);
390   *normtype = ksp->normtype;
391   PetscFunctionReturn(0);
392 }
393 
394 #if defined(PETSC_HAVE_AMS)
395 #undef __FUNCT__
396 #define __FUNCT__ "KSPPublish_Petsc"
397 static PetscErrorCode KSPPublish_Petsc(PetscObject obj)
398 {
399   KSP            ksp = (KSP) obj;
400   PetscErrorCode ierr;
401 
402   PetscFunctionBegin;
403   ierr = AMS_Memory_add_field(obj->amem,"its",&ksp->its,1,AMS_INT,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF);CHKERRQ(ierr);
404   PetscFunctionReturn(0);
405 }
406 #endif
407 
408 #undef __FUNCT__
409 #define __FUNCT__ "KSPSetOperators"
410 /*@
411    KSPSetOperators - Sets the matrix associated with the linear system
412    and a (possibly) different one associated with the preconditioner.
413 
414    Collective on KSP and Mat
415 
416    Input Parameters:
417 +  ksp - the KSP context
418 .  Amat - the matrix associated with the linear system
419 .  Pmat - the matrix to be used in constructing the preconditioner, usually the
420           same as Amat.
421 -  flag - flag indicating information about the preconditioner matrix structure
422    during successive linear solves.  This flag is ignored the first time a
423    linear system is solved, and thus is irrelevant when solving just one linear
424    system.
425 
426    Notes:
427    The flag can be used to eliminate unnecessary work in the preconditioner
428    during the repeated solution of linear systems of the same size.  The
429    available options are
430 $    SAME_PRECONDITIONER -
431 $      Pmat is identical during successive linear solves.
432 $      This option is intended for folks who are using
433 $      different Amat and Pmat matrices and want to reuse the
434 $      same preconditioner matrix.  For example, this option
435 $      saves work by not recomputing incomplete factorization
436 $      for ILU/ICC preconditioners.
437 $    SAME_NONZERO_PATTERN -
438 $      Pmat has the same nonzero structure during
439 $      successive linear solves.
440 $    DIFFERENT_NONZERO_PATTERN -
441 $      Pmat does not have the same nonzero structure.
442 
443     All future calls to KSPSetOperators() must use the same size matrices!
444 
445     Passing a PETSC_NULL for Amat or Pmat removes the matrix that is currently used.
446 
447     If you wish to replace either Amat or Pmat but leave the other one untouched then
448     first call KSPGetOperators() to get the one you wish to keep, call PetscObjectReference()
449     on it and then pass it back in in your call to KSPSetOperators().
450 
451     Caution:
452     If you specify SAME_NONZERO_PATTERN, PETSc believes your assertion
453     and does not check the structure of the matrix.  If you erroneously
454     claim that the structure is the same when it actually is not, the new
455     preconditioner will not function correctly.  Thus, use this optimization
456     feature carefully!
457 
458     If in doubt about whether your preconditioner matrix has changed
459     structure or not, use the flag DIFFERENT_NONZERO_PATTERN.
460 
461     Level: beginner
462 
463    Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
464       are created in PC and returned to the user. In this case, if both operators
465       mat and pmat are requested, two DIFFERENT operators will be returned. If
466       only one is requested both operators in the PC will be the same (i.e. as
467       if one had called KSP/PCSetOperators() with the same argument for both Mats).
468       The user must set the sizes of the returned matrices and their type etc just
469       as if the user created them with MatCreate(). For example,
470 
471 $         KSP/PCGetOperators(ksp/pc,&mat,PETSC_NULL,PETSC_NULL); is equivalent to
472 $           set size, type, etc of mat
473 
474 $         MatCreate(comm,&mat);
475 $         KSP/PCSetOperators(ksp/pc,mat,mat,SAME_NONZERO_PATTERN);
476 $         PetscObjectDereference((PetscObject)mat);
477 $           set size, type, etc of mat
478 
479      and
480 
481 $         KSP/PCGetOperators(ksp/pc,&mat,&pmat,PETSC_NULL); is equivalent to
482 $           set size, type, etc of mat and pmat
483 
484 $         MatCreate(comm,&mat);
485 $         MatCreate(comm,&pmat);
486 $         KSP/PCSetOperators(ksp/pc,mat,pmat,SAME_NONZERO_PATTERN);
487 $         PetscObjectDereference((PetscObject)mat);
488 $         PetscObjectDereference((PetscObject)pmat);
489 $           set size, type, etc of mat and pmat
490 
491     The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
492     of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
493     managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
494     at this is when you create a SNES you do not NEED to create a KSP and attach it to
495     the SNES object (the SNES object manages it for you). Similarly when you create a KSP
496     you do not need to attach a PC to it (the KSP object manages the PC object for you).
497     Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
498     it can be created for you?
499 
500 .keywords: KSP, set, operators, matrix, preconditioner, linear system
501 
502 .seealso: KSPSolve(), KSPGetPC(), PCGetOperators(), PCSetOperators(), KSPGetOperators()
503 @*/
504 PetscErrorCode  KSPSetOperators(KSP ksp,Mat Amat,Mat Pmat,MatStructure flag)
505 {
506   MatNullSpace   nullsp;
507   PetscErrorCode ierr;
508 
509   PetscFunctionBegin;
510   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
511   if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2);
512   if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3);
513   if (Amat) PetscCheckSameComm(ksp,1,Amat,2);
514   if (Pmat) PetscCheckSameComm(ksp,1,Pmat,3);
515   if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
516   ierr = PCSetOperators(ksp->pc,Amat,Pmat,flag);CHKERRQ(ierr);
517   if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = KSP_SETUP_NEWMATRIX;  /* so that next solve call will call PCSetUp() on new matrix */
518   if (ksp->guess) {
519     ierr = KSPFischerGuessReset(ksp->guess);CHKERRQ(ierr);
520   }
521   if (Pmat) {
522     ierr = MatGetNullSpace(Pmat, &nullsp);CHKERRQ(ierr);
523     if (nullsp) {
524       ierr = KSPSetNullSpace(ksp, nullsp);CHKERRQ(ierr);
525     }
526   }
527   PetscFunctionReturn(0);
528 }
529 
530 #undef __FUNCT__
531 #define __FUNCT__ "KSPGetOperators"
532 /*@
533    KSPGetOperators - Gets the matrix associated with the linear system
534    and a (possibly) different one associated with the preconditioner.
535 
536    Collective on KSP and Mat
537 
538    Input Parameter:
539 .  ksp - the KSP context
540 
541    Output Parameters:
542 +  Amat - the matrix associated with the linear system
543 .  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
544 -  flag - flag indicating information about the preconditioner matrix structure
545    during successive linear solves.  This flag is ignored the first time a
546    linear system is solved, and thus is irrelevant when solving just one linear
547    system.
548 
549     Level: intermediate
550 
551    Notes: DOES NOT increase the reference counts of the matrix, so you should NOT destroy them.
552 
553 .keywords: KSP, set, get, operators, matrix, preconditioner, linear system
554 
555 .seealso: KSPSolve(), KSPGetPC(), PCGetOperators(), PCSetOperators(), KSPSetOperators(), KSPGetOperatorsSet()
556 @*/
557 PetscErrorCode  KSPGetOperators(KSP ksp,Mat *Amat,Mat *Pmat,MatStructure *flag)
558 {
559   PetscErrorCode ierr;
560 
561   PetscFunctionBegin;
562   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
563   if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
564   ierr = PCGetOperators(ksp->pc,Amat,Pmat,flag);CHKERRQ(ierr);
565   PetscFunctionReturn(0);
566 }
567 
568 #undef __FUNCT__
569 #define __FUNCT__ "KSPGetOperatorsSet"
570 /*@C
571    KSPGetOperatorsSet - Determines if the matrix associated with the linear system and
572    possibly a different one associated with the preconditioner have been set in the KSP.
573 
574    Not collective, though the results on all processes should be the same
575 
576    Input Parameter:
577 .  pc - the KSP context
578 
579    Output Parameters:
580 +  mat - the matrix associated with the linear system was set
581 -  pmat - matrix associated with the preconditioner was set, usually the same
582 
583    Level: intermediate
584 
585 .keywords: KSP, get, operators, matrix, linear system
586 
587 .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators(), PCGetOperatorsSet()
588 @*/
589 PetscErrorCode  KSPGetOperatorsSet(KSP ksp,PetscBool  *mat,PetscBool  *pmat)
590 {
591   PetscErrorCode ierr;
592 
593   PetscFunctionBegin;
594   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
595   if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);}
596   ierr = PCGetOperatorsSet(ksp->pc,mat,pmat);CHKERRQ(ierr);
597   PetscFunctionReturn(0);
598 }
599 
600 #undef __FUNCT__
601 #define __FUNCT__ "KSPCreate"
602 /*@
603    KSPCreate - Creates the default KSP context.
604 
605    Collective on MPI_Comm
606 
607    Input Parameter:
608 .  comm - MPI communicator
609 
610    Output Parameter:
611 .  ksp - location to put the KSP context
612 
613    Notes:
614    The default KSP type is GMRES with a restart of 30, using modified Gram-Schmidt
615    orthogonalization.
616 
617    Level: beginner
618 
619 .keywords: KSP, create, context
620 
621 .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSP
622 @*/
623 PetscErrorCode  KSPCreate(MPI_Comm comm,KSP *inksp)
624 {
625   KSP            ksp;
626   PetscErrorCode ierr;
627   void           *ctx;
628 
629   PetscFunctionBegin;
630   PetscValidPointer(inksp,2);
631   *inksp = 0;
632 #ifndef PETSC_USE_DYNAMIC_LIBRARIES
633   ierr = KSPInitializePackage(PETSC_NULL);CHKERRQ(ierr);
634 #endif
635 
636   ierr = PetscHeaderCreate(ksp,_p_KSP,struct _KSPOps,KSP_CLASSID,-1,"KSP","Krylov Method","KSP",comm,KSPDestroy,KSPView);CHKERRQ(ierr);
637 
638   ksp->max_it        = 10000;
639   ksp->pc_side       = PC_SIDE_DEFAULT;
640   ksp->rtol          = 1.e-5;
641   ksp->abstol        = 1.e-50;
642   ksp->divtol        = 1.e4;
643 
644   ksp->chknorm             = -1;
645   ksp->normtype            = KSP_NORM_DEFAULT;
646   ksp->rnorm               = 0.0;
647   ksp->its                 = 0;
648   ksp->guess_zero          = PETSC_TRUE;
649   ksp->calc_sings          = PETSC_FALSE;
650   ksp->res_hist            = PETSC_NULL;
651   ksp->res_hist_alloc      = PETSC_NULL;
652   ksp->res_hist_len        = 0;
653   ksp->res_hist_max        = 0;
654   ksp->res_hist_reset      = PETSC_TRUE;
655   ksp->numbermonitors      = 0;
656 
657   ierr = KSPDefaultConvergedCreate(&ctx);CHKERRQ(ierr);
658   ierr = KSPSetConvergenceTest(ksp,KSPDefaultConverged,ctx,KSPDefaultConvergedDestroy);CHKERRQ(ierr);
659   ksp->ops->buildsolution  = KSPDefaultBuildSolution;
660   ksp->ops->buildresidual  = KSPDefaultBuildResidual;
661 #if defined(PETSC_HAVE_AMS)
662   ((PetscObject)ksp)->bops->publish       = KSPPublish_Petsc;
663 #endif
664 
665   ksp->vec_sol         = 0;
666   ksp->vec_rhs         = 0;
667   ksp->pc              = 0;
668   ksp->data            = 0;
669   ksp->nwork           = 0;
670   ksp->work            = 0;
671   ksp->reason          = KSP_CONVERGED_ITERATING;
672   ksp->setupstage      = KSP_SETUP_NEW;
673 
674   ierr = KSPNormSupportTableReset_Private(ksp);CHKERRQ(ierr);
675 
676   *inksp = ksp;
677   PetscFunctionReturn(0);
678 }
679 
680 #undef __FUNCT__
681 #define __FUNCT__ "KSPSetType"
682 /*@C
683    KSPSetType - Builds KSP for a particular solver.
684 
685    Logically Collective on KSP
686 
687    Input Parameters:
688 +  ksp      - the Krylov space context
689 -  type - a known method
690 
691    Options Database Key:
692 .  -ksp_type  <method> - Sets the method; use -help for a list
693     of available methods (for instance, cg or gmres)
694 
695    Notes:
696    See "petsc/include/petscksp.h" for available methods (for instance,
697    KSPCG or KSPGMRES).
698 
699   Normally, it is best to use the KSPSetFromOptions() command and
700   then set the KSP type from the options database rather than by using
701   this routine.  Using the options database provides the user with
702   maximum flexibility in evaluating the many different Krylov methods.
703   The KSPSetType() routine is provided for those situations where it
704   is necessary to set the iterative solver independently of the command
705   line or options database.  This might be the case, for example, when
706   the choice of iterative solver changes during the execution of the
707   program, and the user's application is taking responsibility for
708   choosing the appropriate method.  In other words, this routine is
709   not for beginners.
710 
711   Level: intermediate
712 
713 .keywords: KSP, set, method
714 
715 .seealso: PCSetType(), KSPType
716 
717 @*/
718 PetscErrorCode  KSPSetType(KSP ksp, KSPType type)
719 {
720   PetscErrorCode ierr,(*r)(KSP);
721   PetscBool      match;
722 
723   PetscFunctionBegin;
724   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
725   PetscValidCharPointer(type,2);
726 
727   ierr = PetscObjectTypeCompare((PetscObject)ksp,type,&match);CHKERRQ(ierr);
728   if (match) PetscFunctionReturn(0);
729 
730   ierr =  PetscFListFind(KSPList,((PetscObject)ksp)->comm,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr);
731   if (!r) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested KSP type %s",type);
732   /* Destroy the previous private KSP context */
733   if (ksp->ops->destroy) {
734     ierr = (*ksp->ops->destroy)(ksp);CHKERRQ(ierr);
735     ksp->ops->destroy = PETSC_NULL;
736   }
737   /* Reinitialize function pointers in KSPOps structure */
738   ierr = PetscMemzero(ksp->ops,sizeof(struct _KSPOps));CHKERRQ(ierr);
739   ksp->ops->buildsolution = KSPDefaultBuildSolution;
740   ksp->ops->buildresidual = KSPDefaultBuildResidual;
741   ierr = KSPNormSupportTableReset_Private(ksp);CHKERRQ(ierr);
742   /* Call the KSPCreate_XXX routine for this particular Krylov solver */
743   ksp->setupstage = KSP_SETUP_NEW;
744   ierr = PetscObjectChangeTypeName((PetscObject)ksp,type);CHKERRQ(ierr);
745   ierr = (*r)(ksp);CHKERRQ(ierr);
746 #if defined(PETSC_HAVE_AMS)
747   if (PetscAMSPublishAll) {
748     ierr = PetscObjectAMSPublish((PetscObject)ksp);CHKERRQ(ierr);
749   }
750 #endif
751   PetscFunctionReturn(0);
752 }
753 
754 #undef __FUNCT__
755 #define __FUNCT__ "KSPRegisterDestroy"
756 /*@
757    KSPRegisterDestroy - Frees the list of KSP methods that were
758    registered by KSPRegisterDynamic().
759 
760    Not Collective
761 
762    Level: advanced
763 
764 .keywords: KSP, register, destroy
765 
766 .seealso: KSPRegisterDynamic(), KSPRegisterAll()
767 @*/
768 PetscErrorCode  KSPRegisterDestroy(void)
769 {
770   PetscErrorCode ierr;
771 
772   PetscFunctionBegin;
773   ierr = PetscFListDestroy(&KSPList);CHKERRQ(ierr);
774   KSPRegisterAllCalled = PETSC_FALSE;
775   PetscFunctionReturn(0);
776 }
777 
778 #undef __FUNCT__
779 #define __FUNCT__ "KSPGetType"
780 /*@C
781    KSPGetType - Gets the KSP type as a string from the KSP object.
782 
783    Not Collective
784 
785    Input Parameter:
786 .  ksp - Krylov context
787 
788    Output Parameter:
789 .  name - name of KSP method
790 
791    Level: intermediate
792 
793 .keywords: KSP, get, method, name
794 
795 .seealso: KSPSetType()
796 @*/
797 PetscErrorCode  KSPGetType(KSP ksp,KSPType *type)
798 {
799   PetscFunctionBegin;
800   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
801   PetscValidPointer(type,2);
802   *type = ((PetscObject)ksp)->type_name;
803   PetscFunctionReturn(0);
804 }
805 
806 #undef __FUNCT__
807 #define __FUNCT__ "KSPRegister"
808 /*@C
809   KSPRegister - See KSPRegisterDynamic()
810 
811   Level: advanced
812 @*/
813 PetscErrorCode  KSPRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(KSP))
814 {
815   PetscErrorCode ierr;
816   char           fullname[PETSC_MAX_PATH_LEN];
817 
818   PetscFunctionBegin;
819   ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr);
820   ierr = PetscFListAdd(&KSPList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr);
821   PetscFunctionReturn(0);
822 }
823 
824 #undef __FUNCT__
825 #define __FUNCT__ "KSPSetNullSpace"
826 /*@
827   KSPSetNullSpace - Sets the null space of the operator
828 
829   Logically Collective on KSP
830 
831   Input Parameters:
832 +  ksp - the Krylov space object
833 -  nullsp - the null space of the operator
834 
835   Level: advanced
836 
837 .seealso: KSPSetOperators(), MatNullSpaceCreate(), KSPGetNullSpace(), MatSetNullSpace()
838 @*/
839 PetscErrorCode  KSPSetNullSpace(KSP ksp,MatNullSpace nullsp)
840 {
841   PetscErrorCode ierr;
842 
843   PetscFunctionBegin;
844   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
845   PetscValidHeaderSpecific(nullsp,MAT_NULLSPACE_CLASSID,2);
846   ierr = PetscObjectReference((PetscObject)nullsp);CHKERRQ(ierr);
847   if (ksp->nullsp) { ierr = MatNullSpaceDestroy(&ksp->nullsp);CHKERRQ(ierr); }
848   ksp->nullsp = nullsp;
849   PetscFunctionReturn(0);
850 }
851 
852 #undef __FUNCT__
853 #define __FUNCT__ "KSPGetNullSpace"
854 /*@
855   KSPGetNullSpace - Gets the null space of the operator
856 
857   Not Collective
858 
859   Input Parameters:
860 +  ksp - the Krylov space object
861 -  nullsp - the null space of the operator
862 
863   Level: advanced
864 
865 .seealso: KSPSetOperators(), MatNullSpaceCreate(), KSPSetNullSpace()
866 @*/
867 PetscErrorCode  KSPGetNullSpace(KSP ksp,MatNullSpace *nullsp)
868 {
869   PetscFunctionBegin;
870   PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
871   PetscValidPointer(nullsp,2);
872   *nullsp = ksp->nullsp;
873   PetscFunctionReturn(0);
874 }
875 
876