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