Actual source code: itcreate.c
1: /*$Id: itcreate.c,v 1.206 2001/08/06 21:16:38 bsmith Exp $*/
2: /*
3: The basic KSP routines, Create, View etc. are here.
4: */
5: #include src/sles/ksp/kspimpl.h
6: #include petscsys.h
8: /* Logging support */
9: int KSP_COOKIE;
10: int KSP_GMRESOrthogonalization;
12: EXTERN int SLESInitializePackage(char *);
14: PetscTruth KSPRegisterAllCalled = PETSC_FALSE;
16: #undef __FUNCT__
18: /*@C
19: KSPView - Prints the KSP data structure.
21: Collective on KSP
23: Input Parameters:
24: + ksp - the Krylov space context
25: - viewer - visualization context
27: Note:
28: The available visualization contexts include
29: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
30: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
31: output where only the first processor opens
32: the file. All other processors send their
33: data to the first processor to print.
35: The user can open an alternative visualization context with
36: PetscViewerASCIIOpen() - output to a specified file.
38: Level: developer
40: .keywords: KSP, view
42: .seealso: PCView(), PetscViewerASCIIOpen()
43: @*/
44: int KSPView(KSP ksp,PetscViewer viewer)
45: {
46: char *type;
47: int ierr;
48: PetscTruth isascii;
52: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(ksp->comm);
56: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
57: if (isascii) {
58: KSPGetType(ksp,&type);
59: if (ksp->prefix) {
60: PetscViewerASCIIPrintf(viewer,"KSP Object:(%s)n",ksp->prefix);
61: } else {
62: PetscViewerASCIIPrintf(viewer,"KSP Object:n");
63: }
64: if (type) {
65: PetscViewerASCIIPrintf(viewer," type: %sn",type);
66: } else {
67: PetscViewerASCIIPrintf(viewer," type: not yet setn");
68: }
69: if (ksp->ops->view) {
70: PetscViewerASCIIPushTab(viewer);
71: (*ksp->ops->view)(ksp,viewer);
72: PetscViewerASCIIPopTab(viewer);
73: }
74: if (ksp->guess_zero) {PetscViewerASCIIPrintf(viewer," maximum iterations=%d, initial guess is zeron",ksp->max_it);}
75: else {PetscViewerASCIIPrintf(viewer," maximum iterations=%dn", ksp->max_it);}
76: if (ksp->guess_knoll) {PetscViewerASCIIPrintf(viewer," using preconditioner applied to right hand side for initial guessn");}
77: PetscViewerASCIIPrintf(viewer," tolerances: relative=%g, absolute=%g, divergence=%gn",ksp->rtol,ksp->atol,ksp->divtol);
78: if (ksp->pc_side == PC_RIGHT) {PetscViewerASCIIPrintf(viewer," right preconditioningn");}
79: else if (ksp->pc_side == PC_SYMMETRIC) {PetscViewerASCIIPrintf(viewer," symmetric preconditioningn");}
80: else {PetscViewerASCIIPrintf(viewer," left preconditioningn");}
81: } else {
82: if (ksp->ops->view) {
83: (*ksp->ops->view)(ksp,viewer);
84: }
85: }
86: return(0);
87: }
89: /*
90: Contains the list of registered KSP routines
91: */
92: PetscFList KSPList = 0;
94: #undef __FUNCT__
96: /*@C
97: KSPSetNormType - Sets the norm that is used for convergence testing.
99: Collective on KSP
101: Input Parameter:
102: + ksp - Krylov solver context
103: - normtype - one of
104: $ KSP_NO_NORM - skips computing the norm, this should only be used if you are using
105: $ the Krylov method as a smoother with a fixed small number of iterations.
106: $ You must also call KSPSetConvergenceTest(ksp,KSPSkipConverged,PETSC_NULL);
107: $ supported only by CG, Richardson, Bi-CG-stab, CR, and CGS methods.
108: $ KSP_PRECONDITIONED_NORM - the default for left preconditioned solves, uses the l2 norm
109: $ of the preconditioned residual
110: $ KSP_UNPRECONDITIONED_NORM - uses the l2 norm of the true b - Ax residual, supported only by
111: $ CG, CHEBYCHEV, and RICHARDSON
112: $ KSP_NATURAL_NORM - supported by cg, cr, and cgs
115: Options Database Key:
116: . -ksp_norm_type <none,preconditioned,unpreconditioned,natural>
118: Notes:
119: Currently only works with the CG, Richardson, Bi-CG-stab, CR, and CGS methods.
121: Level: advanced
123: .keywords: KSP, create, context, norms
125: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPSkipConverged()
126: @*/
127: int KSPSetNormType(KSP ksp,KSPNormType normtype)
128: {
132: ksp->normtype = normtype;
133: if (normtype == KSP_NO_NORM) {
134: PetscLogInfo(ksp,"KSPSetNormType:Warning seting KSPNormType to skip computing the normn
135: make sure you set the KSP convergence test to KSPSkipConvergencen");
136: }
137: return(0);
138: }
140: #undef __FUNCT__
142: static int KSPPublish_Petsc(PetscObject obj)
143: {
144: #if defined(PETSC_HAVE_AMS)
145: KSP v = (KSP) obj;
146: int ierr;
147: #endif
151: #if defined(PETSC_HAVE_AMS)
152: /* if it is already published then return */
153: if (v->amem >=0) return(0);
155: PetscObjectPublishBaseBegin(obj);
156: AMS_Memory_add_field((AMS_Memory)v->amem,"Iteration",&v->its,1,AMS_INT,AMS_READ,
157: AMS_COMMON,AMS_REDUCT_UNDEF);
158: AMS_Memory_add_field((AMS_Memory)v->amem,"Residual",&v->rnorm,1,AMS_DOUBLE,AMS_READ,
159: AMS_COMMON,AMS_REDUCT_UNDEF);
161: if (v->res_hist_max > 0) {
162: AMS_Memory_add_field((AMS_Memory)v->amem,"ResidualNormsCount",&v->res_hist_len,1,AMS_INT,AMS_READ,
163: AMS_COMMON,AMS_REDUCT_UNDEF);
164: AMS_Memory_add_field((AMS_Memory)v->amem,"ResidualNormsCountMax",&v->res_hist_max,1,AMS_INT,AMS_READ,
165: AMS_COMMON,AMS_REDUCT_UNDEF);
166: AMS_Memory_add_field((AMS_Memory)v->amem,"ResidualNorms",v->res_hist,v->res_hist_max,AMS_DOUBLE,AMS_READ,
167: AMS_COMMON,AMS_REDUCT_UNDEF);
168: }
170: PetscObjectPublishBaseEnd(obj);
171: #endif
173: return(0);
174: }
177: #undef __FUNCT__
179: /*@C
180: KSPCreate - Creates the default KSP context.
182: Collective on MPI_Comm
184: Input Parameter:
185: . comm - MPI communicator
187: Output Parameter:
188: . ksp - location to put the KSP context
190: Notes:
191: The default KSP type is GMRES with a restart of 30, using modified Gram-Schmidt
192: orthogonalization.
194: Level: developer
196: .keywords: KSP, create, context
198: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSP
199: @*/
200: int KSPCreate(MPI_Comm comm,KSP *inksp)
201: {
202: KSP ksp;
207: *inksp = 0;
208: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
209: SLESInitializePackage(PETSC_NULL);
210: #endif
212: PetscHeaderCreate(ksp,_p_KSP,struct _KSPOps,KSP_COOKIE,-1,"KSP",comm,KSPDestroy,KSPView);
213: PetscLogObjectCreate(ksp);
214: *inksp = ksp;
215: ksp->bops->publish = KSPPublish_Petsc;
217: ksp->type = -1;
218: ksp->max_it = 10000;
219: ksp->pc_side = PC_LEFT;
220: ksp->rtol = 1.e-5;
221: ksp->atol = 1.e-50;
222: ksp->divtol = 1.e4;
224: ksp->normtype = KSP_PRECONDITIONED_NORM;
225: ksp->rnorm = 0.0;
226: ksp->its = 0;
227: ksp->guess_zero = PETSC_TRUE;
228: ksp->calc_sings = PETSC_FALSE;
229: ksp->res_hist = PETSC_NULL;
230: ksp->res_hist_len = 0;
231: ksp->res_hist_max = 0;
232: ksp->res_hist_reset = PETSC_TRUE;
233: ksp->numbermonitors = 0;
234: ksp->converged = KSPDefaultConverged;
235: ksp->ops->buildsolution = KSPDefaultBuildSolution;
236: ksp->ops->buildresidual = KSPDefaultBuildResidual;
238: ksp->ops->setfromoptions = 0;
240: ksp->vec_sol = 0;
241: ksp->vec_rhs = 0;
242: ksp->B = 0;
244: ksp->ops->solve = 0;
245: ksp->ops->setup = 0;
246: ksp->ops->destroy = 0;
248: ksp->data = 0;
249: ksp->nwork = 0;
250: ksp->work = 0;
252: ksp->cnvP = 0;
254: ksp->reason = KSP_CONVERGED_ITERATING;
256: ksp->setupcalled = 0;
257: PetscPublishAll(ksp);
258: return(0);
259: }
260:
261: #undef __FUNCT__
263: /*@C
264: KSPSetType - Builds KSP for a particular solver.
266: Collective on KSP
268: Input Parameters:
269: + ksp - the Krylov space context
270: - type - a known method
272: Options Database Key:
273: . -ksp_type <method> - Sets the method; use -help for a list
274: of available methods (for instance, cg or gmres)
276: Notes:
277: See "petsc/include/petscksp.h" for available methods (for instance,
278: KSPCG or KSPGMRES).
280: Normally, it is best to use the SLESSetFromOptions() command and
281: then set the KSP type from the options database rather than by using
282: this routine. Using the options database provides the user with
283: maximum flexibility in evaluating the many different Krylov methods.
284: The KSPSetType() routine is provided for those situations where it
285: is necessary to set the iterative solver independently of the command
286: line or options database. This might be the case, for example, when
287: the choice of iterative solver changes during the execution of the
288: program, and the user's application is taking responsibility for
289: choosing the appropriate method. In other words, this routine is
290: not for beginners.
292: Level: intermediate
294: .keywords: KSP, set, method
296: .seealso: PCSetType(), KSPType
298: @*/
299: int KSPSetType(KSP ksp,KSPType type)
300: {
301: int ierr,(*r)(KSP);
302: PetscTruth match;
308: PetscTypeCompare((PetscObject)ksp,type,&match);
309: if (match) return(0);
311: if (ksp->data) {
312: /* destroy the old private KSP context */
313: (*ksp->ops->destroy)(ksp);
314: ksp->data = 0;
315: }
316: /* Get the function pointers for the iterative method requested */
317: if (!KSPRegisterAllCalled) {KSPRegisterAll(PETSC_NULL);}
319: PetscFListFind(ksp->comm,KSPList,type,(void (**)(void)) &r);
321: if (!r) SETERRQ1(1,"Unknown KSP type given: %s",type);
323: ksp->setupcalled = 0;
324: (*r)(ksp);
326: PetscObjectChangeTypeName((PetscObject)ksp,type);
327: return(0);
328: }
330: #undef __FUNCT__
332: /*@C
333: KSPRegisterDestroy - Frees the list of KSP methods that were
334: registered by KSPRegisterDynamic().
336: Not Collective
338: Level: advanced
340: .keywords: KSP, register, destroy
342: .seealso: KSPRegisterDynamic(), KSPRegisterAll()
343: @*/
344: int KSPRegisterDestroy(void)
345: {
349: if (KSPList) {
350: PetscFListDestroy(&KSPList);
351: KSPList = 0;
352: }
353: KSPRegisterAllCalled = PETSC_FALSE;
354: return(0);
355: }
357: #undef __FUNCT__
359: /*@C
360: KSPGetType - Gets the KSP type as a string from the KSP object.
362: Not Collective
364: Input Parameter:
365: . ksp - Krylov context
367: Output Parameter:
368: . name - name of KSP method
370: Level: intermediate
372: .keywords: KSP, get, method, name
374: .seealso: KSPSetType()
375: @*/
376: int KSPGetType(KSP ksp,KSPType *type)
377: {
380: *type = ksp->type_name;
381: return(0);
382: }
384: #undef __FUNCT__
386: /*@
387: KSPSetFromOptions - Sets KSP options from the options database.
388: This routine must be called before KSPSetUp() if the user is to be
389: allowed to set the Krylov type.
391: Collective on KSP
393: Input Parameters:
394: . ksp - the Krylov space context
396: Options Database Keys:
397: + -ksp_max_it - maximum number of linear iterations
398: . -ksp_rtol rtol - relative tolerance used in default determination of convergence, i.e.
399: if residual norm decreases by this factor than convergence is declared
400: . -ksp_atol atol - absolute tolerance used in default convergence test, i.e. if residual
401: norm is less than this then convergence is declared
402: . -ksp_divtol tol - if residual norm increases by this factor than divergence is declared
403: . -ksp_norm_type - none - skip norms used in convergence tests (useful only when not using
404: $ convergence test (say you always want to run with 5 iterations) to
405: $ save on communication overhead
406: $ preconditioned - default for left preconditioning
407: $ unpreconditioned - see KSPSetNormType()
408: $ natural - see KSPSetNormType()
409: . -ksp_knoll - compute initial guess by applying the preconditioner to the right hand side
410: . -ksp_cancelmonitors - cancel all previous convergene monitor routines set
411: . -ksp_monitor - print residual norm at each iteration
412: . -ksp_xmonitor - plot residual norm at each iteration
413: . -ksp_vecmonitor - plot solution at each iteration
414: - -ksp_singmonitor - monitor extremem singular values at each iteration
416: Notes:
417: To see all options, run your program with the -help option
418: or consult the users manual.
420: Level: developer
422: .keywords: KSP, set, from, options, database
424: .seealso:
425: @*/
426: int KSPSetFromOptions(KSP ksp)
427: {
428: int ierr;
429: char type[256],*stype[] = {"none","preconditioned","unpreconditioned","natural"};
430: PetscTruth flg;
434: if (!KSPRegisterAllCalled) {KSPRegisterAll(PETSC_NULL);}
435: PetscOptionsBegin(ksp->comm,ksp->prefix,"Krylov Method (KSP) Options","KSP");
436: PetscOptionsList("-ksp_type","Krylov method","KSPSetType",KSPList,(char*)(ksp->type_name?ksp->type_name:KSPGMRES),type,256,&flg);
437: if (flg) {
438: KSPSetType(ksp,type);
439: }
440: /*
441: Set the type if it was never set.
442: */
443: if (!ksp->type_name) {
444: KSPSetType(ksp,KSPGMRES);
445: }
447: PetscOptionsInt("-ksp_max_it","Maximum number of iterations","KSPSetTolerances",ksp->max_it,&ksp->max_it,PETSC_NULL);
448: PetscOptionsReal("-ksp_rtol","Relative decrease in residual norm","KSPSetTolerances",ksp->rtol,&ksp->rtol,PETSC_NULL);
449: PetscOptionsReal("-ksp_atol","Absolute value of residual norm","KSPSetTolerances",ksp->atol,&ksp->atol,PETSC_NULL);
450: PetscOptionsReal("-ksp_divtol","Residual norm increase cause divergence","KSPSetTolerances",ksp->divtol,&ksp->divtol,PETSC_NULL);
451: PetscOptionsLogical("-ksp_knoll","Use preconditioner applied to b for initial guess","KSPSetInitialGuessKnoll",ksp->guess_knoll,
452: &ksp->guess_knoll,PETSC_NULL);
454: PetscOptionsEList("-ksp_norm_type","KSP Norm type","KSPSetNormType",stype,4,"preconditioned",type,256,&flg);
455: if (flg) {
456: PetscTruth isnone,ispreconditioned,isunpreconditioned,isnatural;
458: PetscStrcmp(type,stype[0],&isnone);
459: PetscStrcmp(type,stype[1],&ispreconditioned);
460: PetscStrcmp(type,stype[2],&isunpreconditioned);
461: PetscStrcmp(type,stype[3],&isnatural);
463: if (isnone) {
464: KSPSetNormType(ksp,KSP_NO_NORM);
465: KSPSetConvergenceTest(ksp,KSPSkipConverged,0);
466: } else if (ispreconditioned) {
467: KSPSetNormType(ksp,KSP_PRECONDITIONED_NORM);
468: } else if (isunpreconditioned) {
469: KSPSetNormType(ksp,KSP_UNPRECONDITIONED_NORM);
470: } else if (isnatural) {
471: KSPSetNormType(ksp,KSP_NATURAL_NORM);
472: } else {
473: SETERRQ1(1,"Unknown KSP normtype %s",type);
474: }
475: }
477: PetscOptionsName("-ksp_cancelmonitors","Remove any hardwired monitor routines","KSPClearMonitor",&flg);
478: /* -----------------------------------------------------------------------*/
479: /*
480: Cancels all monitors hardwired into code before call to KSPSetFromOptions()
481: */
482: if (flg) {
483: KSPClearMonitor(ksp);
484: }
485: /*
486: Prints preconditioned residual norm at each iteration
487: */
488: PetscOptionsName("-ksp_monitor","Monitor preconditioned residual norm","KSPSetMonitor",&flg);
489: if (flg) {
490: KSPSetMonitor(ksp,KSPDefaultMonitor,PETSC_NULL,PETSC_NULL);
491: }
492: /*
493: Plots the vector solution
494: */
495: PetscOptionsName("-ksp_vecmonitor","Monitor solution graphically","KSPSetMonitor",&flg);
496: if (flg) {
497: KSPSetMonitor(ksp,KSPVecViewMonitor,PETSC_NULL,PETSC_NULL);
498: }
499: /*
500: Prints preconditioned and true residual norm at each iteration
501: */
502: PetscOptionsName("-ksp_truemonitor","Monitor true (unpreconditioned) residual norm","KSPSetMonitor",&flg);
503: if (flg) {
504: KSPSetMonitor(ksp,KSPTrueMonitor,PETSC_NULL,PETSC_NULL);
505: }
506: /*
507: Prints extreme eigenvalue estimates at each iteration
508: */
509: PetscOptionsName("-ksp_singmonitor","Monitor singular values","KSPSetMonitor",&flg);
510: if (flg) {
511: KSPSetComputeSingularValues(ksp,PETSC_TRUE);
512: KSPSetMonitor(ksp,KSPSingularValueMonitor,PETSC_NULL,PETSC_NULL);
513: }
514: /*
515: Prints preconditioned residual norm with fewer digits
516: */
517: PetscOptionsName("-ksp_smonitor","Monitor preconditioned residual norm with fewer digitis","KSPSetMonitor",&flg);
518: if (flg) {
519: KSPSetMonitor(ksp,KSPDefaultSMonitor,PETSC_NULL,PETSC_NULL);
520: }
521: /*
522: Graphically plots preconditioned residual norm
523: */
524: PetscOptionsName("-ksp_xmonitor","Monitor graphically preconditioned residual norm","KSPSetMonitor",&flg);
525: if (flg) {
526: KSPSetMonitor(ksp,KSPLGMonitor,PETSC_NULL,PETSC_NULL);
527: }
528: /*
529: Graphically plots preconditioned and true residual norm
530: */
531: PetscOptionsName("-ksp_xtruemonitor","Monitor graphically true residual norm","KSPSetMonitor",&flg);
532: if (flg){
533: KSPSetMonitor(ksp,KSPLGTrueMonitor,PETSC_NULL,PETSC_NULL);
534: }
536: /* -----------------------------------------------------------------------*/
538: PetscOptionsLogicalGroupBegin("-ksp_left_pc","Use left preconditioning","KSPSetPreconditionerSide",&flg);
539: if (flg) { KSPSetPreconditionerSide(ksp,PC_LEFT); }
540: PetscOptionsLogicalGroup("-ksp_right_pc","Use right preconditioning","KSPSetPreconditionerSide",&flg);
541: if (flg) { KSPSetPreconditionerSide(ksp,PC_RIGHT);}
542: PetscOptionsLogicalGroupEnd("-ksp_symmetric_pc","Use symmetric (factorized) preconditioning","KSPSetPreconditionerSide",&flg);
543: if (flg) { KSPSetPreconditionerSide(ksp,PC_SYMMETRIC);}
545: PetscOptionsName("-ksp_compute_singularvalues","Compute singular values of preconditioned operator","KSPSetComputeSingularValues",&flg);
546: if (flg) { KSPSetComputeSingularValues(ksp,PETSC_TRUE); }
547: PetscOptionsName("-ksp_compute_eigenvalues","Compute eigenvalues of preconditioned operator","KSPSetComputeSingularValues",&flg);
548: if (flg) { KSPSetComputeSingularValues(ksp,PETSC_TRUE); }
549: PetscOptionsName("-ksp_plot_eigenvalues","Scatter plot extreme eigenvalues","KSPSetComputeSingularValues",&flg);
550: if (flg) { KSPSetComputeSingularValues(ksp,PETSC_TRUE); }
552: if (ksp->ops->setfromoptions) {
553: (*ksp->ops->setfromoptions)(ksp);
554: }
555: PetscOptionsEnd();
558: return(0);
559: }
561: /*MC
562: KSPRegisterDynamic - Adds a method to the Krylov subspace solver package.
564: Synopsis:
565: int KSPRegisterDynamic(char *name_solver,char *path,char *name_create,int (*routine_create)(KSP))
567: Not Collective
569: Input Parameters:
570: + name_solver - name of a new user-defined solver
571: . path - path (either absolute or relative) the library containing this solver
572: . name_create - name of routine to create method context
573: - routine_create - routine to create method context
575: Notes:
576: KSPRegisterDynamic() may be called multiple times to add several user-defined solvers.
578: If dynamic libraries are used, then the fourth input argument (routine_create)
579: is ignored.
581: Sample usage:
582: .vb
583: KSPRegisterDynamic("my_solver",/home/username/my_lib/lib/libO/solaris/mylib.a,
584: "MySolverCreate",MySolverCreate);
585: .ve
587: Then, your solver can be chosen with the procedural interface via
588: $ KSPSetType(ksp,"my_solver")
589: or at runtime via the option
590: $ -ksp_type my_solver
592: Level: advanced
594: Notes: Environmental variables such as ${PETSC_ARCH}, ${PETSC_DIR}, ${PETSC_LIB_DIR}, ${BOPT},
595: and others of the form ${any_environmental_variable} occuring in pathname will be
596: replaced with appropriate values.
597: If your function is not being put into a shared library then use KSPRegister() instead
599: .keywords: KSP, register
601: .seealso: KSPRegisterAll(), KSPRegisterDestroy()
603: M*/
605: #undef __FUNCT__
607: /*@C
608: KSPRegister - See KSPRegisterDynamic()
610: Level: advanced
611: @*/
612: int KSPRegister(char *sname,char *path,char *name,int (*function)(KSP))
613: {
614: int ierr;
615: char fullname[256];
618: PetscFListConcat(path,name,fullname);
619: PetscFListAdd(&KSPList,sname,fullname,(void (*)(void))function);
620: return(0);
621: }