Actual source code: specest.c
petsc-3.4.5 2014-06-29
2: #include <petsc-private/kspimpl.h>
4: typedef struct {
5: KSP kspest; /* KSP capable of estimating eigenvalues */
6: KSP kspcheap; /* Cheap smoother (should have few dot products) */
7: PC pcnone; /* Dummy PC to drop in so PCSetFromOptions doesn't get called extra times */
8: PetscReal min,max; /* Singular value estimates */
9: PetscReal radius; /* Spectral radius of 1-B where B is the preconditioned operator */
10: PetscBool current; /* Eigenvalue estimates are current */
11: PetscReal minfactor,maxfactor;
12: PetscReal richfactor;
13: } KSP_SpecEst;
17: static PetscErrorCode KSPSetUp_SpecEst(KSP ksp)
18: {
19: KSP_SpecEst *spec = (KSP_SpecEst*)ksp->data;
21: PetscBool nonzero;
24: KSPSetPC(spec->kspest,ksp->pc);
25: KSPSetPC(spec->kspcheap,ksp->pc);
26: KSPGetInitialGuessNonzero(ksp,&nonzero);
27: KSPSetInitialGuessNonzero(spec->kspest,nonzero);
28: KSPSetInitialGuessNonzero(spec->kspcheap,nonzero);
29: KSPSetComputeSingularValues(spec->kspest,PETSC_TRUE);
30: KSPSetUp(spec->kspest);
32: spec->current = PETSC_FALSE;
33: return(0);
34: }
38: static PetscErrorCode KSPSpecEstPropagateUp(KSP ksp,KSP subksp)
39: {
43: KSPGetConvergedReason(subksp,&ksp->reason);
44: KSPGetIterationNumber(subksp,&ksp->its);
45: return(0);
46: }
50: static PetscErrorCode KSPSolve_SpecEst(KSP ksp)
51: {
53: KSP_SpecEst *spec = (KSP_SpecEst*)ksp->data;
56: if (spec->current) {
57: KSPSolve(spec->kspcheap,ksp->vec_rhs,ksp->vec_sol);
58: KSPSpecEstPropagateUp(ksp,spec->kspcheap);
59: } else {
60: PetscInt i,its,neig;
61: PetscReal *real,*imag,rad = 0;
62: KSPSolve(spec->kspest,ksp->vec_rhs,ksp->vec_sol);
63: KSPSpecEstPropagateUp(ksp,spec->kspest);
64: KSPComputeExtremeSingularValues(spec->kspest,&spec->max,&spec->min);
66: KSPGetIterationNumber(spec->kspest,&its);
67: PetscMalloc2(its,PetscReal,&real,its,PetscReal,&imag);
68: KSPComputeEigenvalues(spec->kspest,its,real,imag,&neig);
69: for (i=0; i<neig; i++) {
70: /* We would really like to compute w (nominally 1/radius) to minimize |1-wB|. Empirically it
71: is better to compute rad = |1-B| than rad = |B|. There must be a cheap way to do better. */
72: rad = PetscMax(rad,PetscRealPart(PetscSqrtScalar((PetscScalar)(PetscSqr(real[i]-1.) + PetscSqr(imag[i])))));
73: }
74: PetscFree2(real,imag);
76: spec->radius = rad;
78: KSPChebyshevSetEigenvalues(spec->kspcheap,spec->max*spec->maxfactor,spec->min*spec->minfactor);
79: KSPRichardsonSetScale(spec->kspcheap,spec->richfactor/spec->radius);
80: PetscInfo3(ksp,"Estimated singular value min=%G max=%G, spectral radius=%G",spec->min,spec->max,spec->radius);
82: spec->current = PETSC_TRUE;
83: }
84: return(0);
85: }
89: static PetscErrorCode KSPView_SpecEst(KSP ksp,PetscViewer viewer)
90: {
91: KSP_SpecEst *spec = (KSP_SpecEst*)ksp->data;
93: PetscBool iascii;
96: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
97: if (iascii) {
98: PetscViewerASCIIPrintf(viewer," SpecEst: last singular value estimate min=%G max=%G rad=%G\n",spec->min,spec->max,spec->radius);
99: PetscViewerASCIIPrintf(viewer," Using scaling factors min=%G max=%G rich=%G\n",spec->minfactor,spec->maxfactor,spec->richfactor);
100: PetscViewerASCIIPrintf(viewer," Sub KSP used for estimating spectrum:\n");
101: PetscViewerASCIIPushTab(viewer);
102: KSPView(spec->kspest,viewer);
103: PetscViewerASCIIPopTab(viewer);
104: PetscViewerASCIIPrintf(viewer," Sub KSP used for subsequent smoothing steps:\n");
105: PetscViewerASCIIPushTab(viewer);
106: KSPView(spec->kspcheap,viewer);
107: PetscViewerASCIIPopTab(viewer);
108: }
109: return(0);
110: }
114: static PetscErrorCode KSPSetFromOptions_SpecEst(KSP ksp)
115: {
117: KSP_SpecEst *spec = (KSP_SpecEst*)ksp->data;
118: char prefix[256];
121: PetscOptionsHead("KSP SpecEst Options");
122: PetscOptionsReal("-ksp_specest_minfactor","Multiplier on the minimum eigen/singular value","None",spec->minfactor,&spec->minfactor,NULL);
123: PetscOptionsReal("-ksp_specest_maxfactor","Multiplier on the maximum eigen/singular value","None",spec->maxfactor,&spec->maxfactor,NULL);
124: PetscOptionsReal("-ksp_specest_richfactor","Multiplier on the richimum eigen/singular value","None",spec->richfactor,&spec->richfactor,NULL);
125: PetscOptionsTail();
127: /* Mask the PC so that PCSetFromOptions does not do anything */
128: KSPSetPC(spec->kspest,spec->pcnone);
129: KSPSetPC(spec->kspcheap,spec->pcnone);
131: PetscSNPrintf(prefix,sizeof(prefix),"%sspecest_",((PetscObject)ksp)->prefix ? ((PetscObject)ksp)->prefix : "");
132: KSPSetOptionsPrefix(spec->kspest,prefix);
133: PetscSNPrintf(prefix,sizeof(prefix),"%sspeccheap_",((PetscObject)ksp)->prefix ? ((PetscObject)ksp)->prefix : "");
134: KSPSetOptionsPrefix(spec->kspcheap,prefix);
136: if (!((PetscObject)spec->kspest)->type_name) {
137: KSPSetType(spec->kspest,KSPGMRES);
138: }
139: if (!((PetscObject)spec->kspcheap)->type_name) {
140: KSPSetType(spec->kspcheap,KSPCHEBYSHEV);
141: }
142: KSPSetFromOptions(spec->kspest);
143: KSPSetFromOptions(spec->kspcheap);
145: /* Unmask the PC */
146: KSPSetPC(spec->kspest,ksp->pc);
147: KSPSetPC(spec->kspcheap,ksp->pc);
148: return(0);
149: }
154: static PetscErrorCode KSPDestroy_SpecEst(KSP ksp)
155: {
157: KSP_SpecEst *spec = (KSP_SpecEst*)ksp->data;
160: KSPDestroy(&spec->kspest);
161: KSPDestroy(&spec->kspcheap);
162: PCDestroy(&spec->pcnone);
163: PetscFree(ksp->data);
164: return(0);
165: }
169: /*MC
170: KSPSPECEST - Estimate the spectrum on the first KSPSolve, then use cheaper smoother for subsequent solves.
172: Options Database Keys:
173: + -ksp_specest_minfactor <0.9> - Multiplier on the minimum eigen/singular value
174: . -ksp_specest_maxfactor <1.1> - Multiplier on the maximum eigen/singular value
175: . -ksp_specest_richfactor <1> - Multiplier on the richimum eigen/singular value
176: . -specest_ksp_type <type> - KSP used to estimate the spectrum (usually CG or GMRES)
177: . -speccheap_ksp_type <type> - KSP used as a cheap smoother once the spectrum has been estimated (usually Chebyshev or Richardson)
178: - see KSPSolve() for more
180: Notes:
181: This KSP estimates the extremal singular values on the first pass, then uses them to configure a smoother that
182: uses fewer dot products. It is intended for use on the levels of multigrid, especially at high process counts,
183: where dot products are very expensive.
185: The same PC is used for both the estimator and the cheap smoother, it is only set up once. There are no options
186: keys for -specest_pc_ or speccheap_pc_ since it is the same object as -pc_.
188: Level: intermediate
190: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGMRES, KSPCG, KSPCHEBYSHEV, KSPRICHARDSON
191: M*/
192: PETSC_EXTERN PetscErrorCode KSPCreate_SpecEst(KSP ksp)
193: {
194: KSP_SpecEst *spec;
198: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
199: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_RIGHT,1);
200: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);
201: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,1);
203: PetscNewLog(ksp,KSP_SpecEst,&spec);
205: ksp->data = (void*)spec;
206: ksp->ops->setup = KSPSetUp_SpecEst;
207: ksp->ops->solve = KSPSolve_SpecEst;
208: ksp->ops->destroy = KSPDestroy_SpecEst;
209: ksp->ops->buildsolution = KSPBuildSolutionDefault;
210: ksp->ops->buildresidual = KSPBuildResidualDefault;
211: ksp->ops->setfromoptions = KSPSetFromOptions_SpecEst;
212: ksp->ops->view = KSPView_SpecEst;
214: spec->minfactor = 0.9;
215: spec->maxfactor = 1.1;
216: spec->richfactor = 1.0;
218: KSPCreate(PetscObjectComm((PetscObject)ksp),&spec->kspest);
219: KSPCreate(PetscObjectComm((PetscObject)ksp),&spec->kspcheap);
221: /* Hold an empty PC */
222: KSPGetPC(spec->kspest,&spec->pcnone);
223: PetscObjectReference((PetscObject)spec->pcnone);
224: PCSetType(spec->pcnone,PCNONE);
225: KSPSetPC(spec->kspcheap,spec->pcnone);
227: KSPSetTolerances(spec->kspest,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,5);
229: /* Make the "cheap" preconditioner cheap by default */
230: KSPSetConvergenceTest(spec->kspcheap,KSPSkipConverged,0,0);
231: KSPSetNormType(spec->kspcheap,KSP_NORM_NONE);
232: KSPSetTolerances(spec->kspcheap,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,5);
233: return(0);
234: }