Actual source code: specest.c

petsc-3.5.0 2014-06-30
Report Typos and Errors
  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,&real,its,&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",(double)spec->min,(double)spec->max,(double)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",(double)spec->min,(double)spec->max,(double)spec->radius);
 99:     PetscViewerASCIIPrintf(viewer,"  Using scaling factors min=%g max=%g rich=%g\n",(double)spec->minfactor,(double)spec->maxfactor,(double)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,&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,KSPConvergedSkip,0,0);
231:   KSPSetNormType(spec->kspcheap,KSP_NORM_NONE);
232:   KSPSetTolerances(spec->kspcheap,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,5);
233:   return(0);
234: }