1
2 #include <../src/ksp/ksp/impls/cheby/chebyshevimpl.h> /*I "petscksp.h" I*/
3
KSPReset_Chebyshev(KSP ksp)4 static PetscErrorCode KSPReset_Chebyshev(KSP ksp)
5 {
6 KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
7 PetscErrorCode ierr;
8
9 PetscFunctionBegin;
10 if (cheb->kspest) {
11 ierr = KSPReset(cheb->kspest);CHKERRQ(ierr);
12 }
13 PetscFunctionReturn(0);
14 }
15
chebyhash(PetscInt xx)16 PETSC_STATIC_INLINE PetscScalar chebyhash(PetscInt xx)
17 {
18 unsigned int x = xx;
19 x = ((x >> 16) ^ x) * 0x45d9f3b;
20 x = ((x >> 16) ^ x) * 0x45d9f3b;
21 x = ((x >> 16) ^ x);
22 return (PetscScalar)((PetscInt64)x-2147483648)*5.e-10; /* center around zero, scaled about -1. to 1.*/
23 }
24
25 /*
26 * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
27 */
KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal * emin,PetscReal * emax)28 static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
29 {
30 PetscErrorCode ierr;
31 PetscInt n,neig;
32 PetscReal *re,*im,min,max;
33
34 PetscFunctionBegin;
35 ierr = KSPGetIterationNumber(kspest,&n);CHKERRQ(ierr);
36 ierr = PetscMalloc2(n,&re,n,&im);CHKERRQ(ierr);
37 ierr = KSPComputeEigenvalues(kspest,n,re,im,&neig);CHKERRQ(ierr);
38 min = PETSC_MAX_REAL;
39 max = PETSC_MIN_REAL;
40 for (n=0; n<neig; n++) {
41 min = PetscMin(min,re[n]);
42 max = PetscMax(max,re[n]);
43 }
44 ierr = PetscFree2(re,im);CHKERRQ(ierr);
45 *emax = max;
46 *emin = min;
47 PetscFunctionReturn(0);
48 }
49
KSPSetUp_Chebyshev(KSP ksp)50 static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
51 {
52 KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
53 PetscErrorCode ierr;
54 PetscBool flg;
55 Mat Pmat,Amat;
56 PetscObjectId amatid, pmatid;
57 PetscObjectState amatstate, pmatstate;
58 PetscFunctionBegin;
59 ierr = KSPSetWorkVecs(ksp,3);CHKERRQ(ierr);
60 if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) { /* We need to estimate eigenvalues */
61 ierr = KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
62 }
63 if (cheb->kspest) {
64 ierr = KSPGetOperators(ksp,&Amat,&Pmat);CHKERRQ(ierr);
65 ierr = MatGetOption(Pmat, MAT_SPD, &flg);CHKERRQ(ierr);
66 if (flg) {
67 const char *prefix;
68 ierr = KSPGetOptionsPrefix(cheb->kspest,&prefix);CHKERRQ(ierr);
69 ierr = PetscOptionsHasName(NULL,prefix,"-ksp_type",&flg);CHKERRQ(ierr);
70 if (!flg) {
71 ierr = KSPSetType(cheb->kspest, KSPCG);CHKERRQ(ierr);
72 }
73 }
74 ierr = PetscObjectGetId((PetscObject)Amat,&amatid);CHKERRQ(ierr);
75 ierr = PetscObjectGetId((PetscObject)Pmat,&pmatid);CHKERRQ(ierr);
76 ierr = PetscObjectStateGet((PetscObject)Amat,&amatstate);CHKERRQ(ierr);
77 ierr = PetscObjectStateGet((PetscObject)Pmat,&pmatstate);CHKERRQ(ierr);
78 if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
79 PetscReal max=0.0,min=0.0;
80 Vec B;
81 KSPConvergedReason reason;
82 ierr = KSPSetPC(cheb->kspest,ksp->pc);CHKERRQ(ierr);
83 if (cheb->usenoisy) {
84 PetscInt n,i,istart;
85 PetscScalar *xx;
86
87 B = ksp->work[1];
88 ierr = VecGetOwnershipRange(B,&istart,NULL);CHKERRQ(ierr);
89 ierr = VecGetLocalSize(B,&n);CHKERRQ(ierr);
90 ierr = VecGetArrayWrite(B,&xx);CHKERRQ(ierr);
91 for (i=0; i<n; i++) xx[i] = chebyhash(i+istart);
92 ierr = VecRestoreArrayWrite(B,&xx);CHKERRQ(ierr);
93 } else {
94 PetscBool change;
95
96 ierr = PCPreSolveChangeRHS(ksp->pc,&change);CHKERRQ(ierr);
97 if (change) {
98 B = ksp->work[1];
99 ierr = VecCopy(ksp->vec_rhs,B);CHKERRQ(ierr);
100 } else B = ksp->vec_rhs;
101 }
102 ierr = KSPSolve(cheb->kspest,B,ksp->work[0]);CHKERRQ(ierr);
103 ierr = KSPGetConvergedReason(cheb->kspest,&reason);CHKERRQ(ierr);
104 if (reason == KSP_DIVERGED_ITS) {
105 ierr = PetscInfo(ksp,"Eigen estimator ran for prescribed number of iterations\n");CHKERRQ(ierr);
106 } else if (reason == KSP_DIVERGED_PC_FAILED) {
107 PetscInt its;
108 PCFailedReason pcreason;
109
110 ierr = KSPGetIterationNumber(cheb->kspest,&its);CHKERRQ(ierr);
111 if (ksp->normtype == KSP_NORM_NONE) {
112 PetscInt sendbuf,recvbuf;
113 ierr = PCGetFailedReasonRank(ksp->pc,&pcreason);CHKERRQ(ierr);
114 sendbuf = (PetscInt)pcreason;
115 ierr = MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr);
116 ierr = PCSetFailedReason(ksp->pc,(PCFailedReason) recvbuf);CHKERRQ(ierr);
117 }
118 ierr = PCGetFailedReason(ksp->pc,&pcreason);CHKERRQ(ierr);
119 ksp->reason = KSP_DIVERGED_PC_FAILED;
120 ierr = VecSetInf(ksp->vec_sol);CHKERRQ(ierr);
121 ierr = PetscInfo3(ksp,"Eigen estimator failed: %s %s at iteration %D",KSPConvergedReasons[reason],PCFailedReasons[pcreason],its);CHKERRQ(ierr);
122 PetscFunctionReturn(0);
123 } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
124 ierr = PetscInfo(ksp,"Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n");CHKERRQ(ierr);
125 } else if (reason < 0) {
126 ierr = PetscInfo1(ksp,"Eigen estimator failed %s, using estimates anyway\n",KSPConvergedReasons[reason]);CHKERRQ(ierr);
127 }
128
129 ierr = KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);CHKERRQ(ierr);
130 ierr = KSPSetPC(cheb->kspest,NULL);CHKERRQ(ierr);
131
132 cheb->emin_computed = min;
133 cheb->emax_computed = max;
134 cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
135 cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;
136
137 cheb->amatid = amatid;
138 cheb->pmatid = pmatid;
139 cheb->amatstate = amatstate;
140 cheb->pmatstate = pmatstate;
141 }
142 }
143 PetscFunctionReturn(0);
144 }
145
KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp,PetscReal emax,PetscReal emin)146 static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp,PetscReal emax,PetscReal emin)
147 {
148 KSP_Chebyshev *chebyshevP = (KSP_Chebyshev*)ksp->data;
149 PetscErrorCode ierr;
150
151 PetscFunctionBegin;
152 if (emax <= emin) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Maximum eigenvalue must be larger than minimum: max %g min %g",(double)emax,(double)emin);
153 if (emax*emin <= 0.0) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Both eigenvalues must be of the same sign: max %g min %g",(double)emax,(double)emin);
154 chebyshevP->emax = emax;
155 chebyshevP->emin = emin;
156
157 ierr = KSPChebyshevEstEigSet(ksp,0.,0.,0.,0.);CHKERRQ(ierr); /* Destroy any estimation setup */
158 PetscFunctionReturn(0);
159 }
160
KSPChebyshevEstEigSet_Chebyshev(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)161 static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
162 {
163 KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
164 PetscErrorCode ierr;
165
166 PetscFunctionBegin;
167 if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
168 if (!cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
169 ierr = KSPCreate(PetscObjectComm((PetscObject)ksp),&cheb->kspest);CHKERRQ(ierr);
170 ierr = KSPSetErrorIfNotConverged(cheb->kspest,ksp->errorifnotconverged);CHKERRQ(ierr);
171 ierr = PetscObjectIncrementTabLevel((PetscObject)cheb->kspest,(PetscObject)ksp,1);CHKERRQ(ierr);
172 /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
173 ierr = PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest,((PetscObject)ksp)->prefix);CHKERRQ(ierr);
174 ierr = PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest,"esteig_");CHKERRQ(ierr);
175 ierr = KSPSetSkipPCSetFromOptions(cheb->kspest,PETSC_TRUE);CHKERRQ(ierr);
176
177 ierr = KSPSetComputeEigenvalues(cheb->kspest,PETSC_TRUE);CHKERRQ(ierr);
178
179 /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
180 ierr = KSPSetTolerances(cheb->kspest,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,cheb->eststeps);CHKERRQ(ierr);
181 }
182 if (a >= 0) cheb->tform[0] = a;
183 if (b >= 0) cheb->tform[1] = b;
184 if (c >= 0) cheb->tform[2] = c;
185 if (d >= 0) cheb->tform[3] = d;
186 cheb->amatid = 0;
187 cheb->pmatid = 0;
188 cheb->amatstate = -1;
189 cheb->pmatstate = -1;
190 } else {
191 ierr = KSPDestroy(&cheb->kspest);CHKERRQ(ierr);
192 }
193 PetscFunctionReturn(0);
194 }
195
KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp,PetscBool use)196 static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp,PetscBool use)
197 {
198 KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
199
200 PetscFunctionBegin;
201 cheb->usenoisy = use;
202 PetscFunctionReturn(0);
203 }
204
205 /*@
206 KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
207 of the preconditioned problem.
208
209 Logically Collective on ksp
210
211 Input Parameters:
212 + ksp - the Krylov space context
213 - emax, emin - the eigenvalue estimates
214
215 Options Database:
216 . -ksp_chebyshev_eigenvalues emin,emax
217
218 Note: Call KSPChebyshevEstEigSet() or use the option -ksp_chebyshev_esteig a,b,c,d to have the KSP
219 estimate the eigenvalues and use these estimated values automatically
220
221 Level: intermediate
222
223 @*/
KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)224 PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
225 {
226 PetscErrorCode ierr;
227
228 PetscFunctionBegin;
229 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
230 PetscValidLogicalCollectiveReal(ksp,emax,2);
231 PetscValidLogicalCollectiveReal(ksp,emin,3);
232 ierr = PetscTryMethod(ksp,"KSPChebyshevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));CHKERRQ(ierr);
233 PetscFunctionReturn(0);
234 }
235
236 /*@
237 KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev
238
239 Logically Collective on ksp
240
241 Input Parameters:
242 + ksp - the Krylov space context
243 . a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
244 . b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
245 . c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
246 - d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
247
248 Options Database:
249 . -ksp_chebyshev_esteig a,b,c,d
250
251 Notes:
252 The Chebyshev bounds are set using
253 .vb
254 minbound = a*minest + b*maxest
255 maxbound = c*minest + d*maxest
256 .ve
257 The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
258 The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.
259
260 If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.
261
262 The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.
263
264 The eigenvalues are estimated using the Lanczo (KSPCG) or Arnoldi (KSPGMRES) process using a noisy right hand side vector.
265
266 Level: intermediate
267
268 @*/
KSPChebyshevEstEigSet(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)269 PetscErrorCode KSPChebyshevEstEigSet(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
270 {
271 PetscErrorCode ierr;
272
273 PetscFunctionBegin;
274 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
275 PetscValidLogicalCollectiveReal(ksp,a,2);
276 PetscValidLogicalCollectiveReal(ksp,b,3);
277 PetscValidLogicalCollectiveReal(ksp,c,4);
278 PetscValidLogicalCollectiveReal(ksp,d,5);
279 ierr = PetscTryMethod(ksp,"KSPChebyshevEstEigSet_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));CHKERRQ(ierr);
280 PetscFunctionReturn(0);
281 }
282
283 /*@
284 KSPChebyshevEstEigSetUseNoisy - use a noisy right hand side in order to do the estimate instead of the given right hand side
285
286 Logically Collective
287
288 Input Arguments:
289 + ksp - linear solver context
290 - use - PETSC_TRUE to use noisy
291
292 Options Database:
293 . -ksp_chebyshev_esteig_noisy <true,false>
294
295 Notes:
296 This alledgely works better for multigrid smoothers
297
298 Level: intermediate
299
300 .seealso: KSPChebyshevEstEigSet()
301 @*/
KSPChebyshevEstEigSetUseNoisy(KSP ksp,PetscBool use)302 PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp,PetscBool use)
303 {
304 PetscErrorCode ierr;
305
306 PetscFunctionBegin;
307 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
308 ierr = PetscTryMethod(ksp,"KSPChebyshevEstEigSetUseNoisy_C",(KSP,PetscBool),(ksp,use));CHKERRQ(ierr);
309 PetscFunctionReturn(0);
310 }
311
312 /*@
313 KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate eigenvalues for the Chebyshev method. If
314 a Krylov method is not being used for this purpose, NULL is returned. The reference count of the returned KSP is
315 not incremented: it should not be destroyed by the user.
316
317 Input Parameters:
318 . ksp - the Krylov space context
319
320 Output Parameters:
321 . kspest the eigenvalue estimation Krylov space context
322
323 Level: intermediate
324
325 .seealso: KSPChebyshevEstEigSet()
326 @*/
KSPChebyshevEstEigGetKSP(KSP ksp,KSP * kspest)327 PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
328 {
329 PetscErrorCode ierr;
330
331 PetscFunctionBegin;
332 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1);
333 PetscValidPointer(kspest,2);
334 *kspest = NULL;
335 ierr = PetscTryMethod(ksp,"KSPChebyshevEstEigGetKSP_C",(KSP,KSP*),(ksp,kspest));CHKERRQ(ierr);
336 PetscFunctionReturn(0);
337 }
338
KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp,KSP * kspest)339 static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
340 {
341 KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
342
343 PetscFunctionBegin;
344 *kspest = cheb->kspest;
345 PetscFunctionReturn(0);
346 }
347
KSPSetFromOptions_Chebyshev(PetscOptionItems * PetscOptionsObject,KSP ksp)348 static PetscErrorCode KSPSetFromOptions_Chebyshev(PetscOptionItems *PetscOptionsObject,KSP ksp)
349 {
350 KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
351 PetscErrorCode ierr;
352 PetscInt neigarg = 2, nestarg = 4;
353 PetscReal eminmax[2] = {0., 0.};
354 PetscReal tform[4] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
355 PetscBool flgeig, flgest;
356
357 PetscFunctionBegin;
358 ierr = PetscOptionsHead(PetscOptionsObject,"KSP Chebyshev Options");CHKERRQ(ierr);
359 ierr = PetscOptionsInt("-ksp_chebyshev_esteig_steps","Number of est steps in Chebyshev","",cheb->eststeps,&cheb->eststeps,NULL);CHKERRQ(ierr);
360 ierr = PetscOptionsRealArray("-ksp_chebyshev_eigenvalues","extreme eigenvalues","KSPChebyshevSetEigenvalues",eminmax,&neigarg,&flgeig);CHKERRQ(ierr);
361 if (flgeig) {
362 if (neigarg != 2) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
363 ierr = KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]);CHKERRQ(ierr);
364 }
365 ierr = PetscOptionsRealArray("-ksp_chebyshev_esteig","estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds","KSPChebyshevEstEigSet",tform,&nestarg,&flgest);CHKERRQ(ierr);
366 if (flgest) {
367 switch (nestarg) {
368 case 0:
369 ierr = KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
370 break;
371 case 2: /* Base everything on the max eigenvalues */
372 ierr = KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,tform[0],PETSC_DECIDE,tform[1]);CHKERRQ(ierr);
373 break;
374 case 4: /* Use the full 2x2 linear transformation */
375 ierr = KSPChebyshevEstEigSet(ksp,tform[0],tform[1],tform[2],tform[3]);CHKERRQ(ierr);
376 break;
377 default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
378 }
379 }
380
381 /* We need to estimate eigenvalues; need to set this here so that KSPSetFromOptions() is called on the estimator */
382 if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) {
383 ierr = KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
384 }
385
386 if (cheb->kspest) {
387 ierr = PetscOptionsBool("-ksp_chebyshev_esteig_noisy","Use noisy right hand side for estimate","KSPChebyshevEstEigSetUseNoisy",cheb->usenoisy,&cheb->usenoisy,NULL);CHKERRQ(ierr);
388 ierr = KSPSetFromOptions(cheb->kspest);CHKERRQ(ierr);
389 }
390 ierr = PetscOptionsTail();CHKERRQ(ierr);
391 PetscFunctionReturn(0);
392 }
393
KSPSolve_Chebyshev(KSP ksp)394 static PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
395 {
396 KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
397 PetscErrorCode ierr;
398 PetscInt k,kp1,km1,ktmp,i;
399 PetscScalar alpha,omegaprod,mu,omega,Gamma,c[3],scale;
400 PetscReal rnorm = 0.0;
401 Vec sol_orig,b,p[3],r;
402 Mat Amat,Pmat;
403 PetscBool diagonalscale;
404
405 PetscFunctionBegin;
406 ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
407 if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
408
409 ierr = PCGetOperators(ksp->pc,&Amat,&Pmat);CHKERRQ(ierr);
410 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
411 ksp->its = 0;
412 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
413 /* These three point to the three active solutions, we
414 rotate these three at each solution update */
415 km1 = 0; k = 1; kp1 = 2;
416 sol_orig = ksp->vec_sol; /* ksp->vec_sol will be asigned to rotating vector p[k], thus save its address */
417 b = ksp->vec_rhs;
418 p[km1] = sol_orig;
419 p[k] = ksp->work[0];
420 p[kp1] = ksp->work[1];
421 r = ksp->work[2];
422
423 /* use scale*B as our preconditioner */
424 scale = 2.0/(cheb->emax + cheb->emin);
425
426 /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
427 alpha = 1.0 - scale*(cheb->emin);
428 Gamma = 1.0;
429 mu = 1.0/alpha;
430 omegaprod = 2.0/alpha;
431
432 c[km1] = 1.0;
433 c[k] = mu;
434
435 if (!ksp->guess_zero) {
436 ierr = KSP_MatMult(ksp,Amat,sol_orig,r);CHKERRQ(ierr); /* r = b - A*p[km1] */
437 ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr);
438 } else {
439 ierr = VecCopy(b,r);CHKERRQ(ierr);
440 }
441
442 /* calculate residual norm if requested, we have done one iteration */
443 if (ksp->normtype) {
444 switch (ksp->normtype) {
445 case KSP_NORM_PRECONDITIONED:
446 ierr = KSP_PCApply(ksp,r,p[k]);CHKERRQ(ierr); /* p[k] = B^{-1}r */
447 ierr = VecNorm(p[k],NORM_2,&rnorm);CHKERRQ(ierr);
448 break;
449 case KSP_NORM_UNPRECONDITIONED:
450 case KSP_NORM_NATURAL:
451 ierr = VecNorm(r,NORM_2,&rnorm);CHKERRQ(ierr);
452 break;
453 default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
454 }
455 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
456 ksp->rnorm = rnorm;
457 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
458 ierr = KSPLogResidualHistory(ksp,rnorm);CHKERRQ(ierr);
459 ierr = KSPMonitor(ksp,0,rnorm);CHKERRQ(ierr);
460 ierr = (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
461 } else ksp->reason = KSP_CONVERGED_ITERATING;
462 if (ksp->reason || ksp->max_it==0) {
463 if (ksp->max_it==0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
464 PetscFunctionReturn(0);
465 }
466 if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
467 ierr = KSP_PCApply(ksp,r,p[k]);CHKERRQ(ierr); /* p[k] = B^{-1}r */
468 }
469 ierr = VecAYPX(p[k],scale,p[km1]);CHKERRQ(ierr); /* p[k] = scale B^{-1}r + p[km1] */
470 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
471 ksp->its = 1;
472 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
473
474 for (i=1; i<ksp->max_it; i++) {
475 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
476 ksp->its++;
477 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
478
479 ierr = KSP_MatMult(ksp,Amat,p[k],r);CHKERRQ(ierr); /* r = b - Ap[k] */
480 ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr);
481 /* calculate residual norm if requested */
482 if (ksp->normtype) {
483 switch (ksp->normtype) {
484 case KSP_NORM_PRECONDITIONED:
485 ierr = KSP_PCApply(ksp,r,p[kp1]);CHKERRQ(ierr); /* p[kp1] = B^{-1}r */
486 ierr = VecNorm(p[kp1],NORM_2,&rnorm);CHKERRQ(ierr);
487 break;
488 case KSP_NORM_UNPRECONDITIONED:
489 case KSP_NORM_NATURAL:
490 ierr = VecNorm(r,NORM_2,&rnorm);CHKERRQ(ierr);
491 break;
492 default:
493 rnorm = 0.0;
494 break;
495 }
496 KSPCheckNorm(ksp,rnorm);
497 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
498 ksp->rnorm = rnorm;
499 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
500 ierr = KSPLogResidualHistory(ksp,rnorm);CHKERRQ(ierr);
501 ierr = KSPMonitor(ksp,i,rnorm);CHKERRQ(ierr);
502 ierr = (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
503 if (ksp->reason) break;
504 if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
505 ierr = KSP_PCApply(ksp,r,p[kp1]);CHKERRQ(ierr); /* p[kp1] = B^{-1}r */
506 }
507 } else {
508 ierr = KSP_PCApply(ksp,r,p[kp1]);CHKERRQ(ierr); /* p[kp1] = B^{-1}r */
509 }
510 ksp->vec_sol = p[k];
511
512 c[kp1] = 2.0*mu*c[k] - c[km1];
513 omega = omegaprod*c[k]/c[kp1];
514
515 /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
516 ierr = VecAXPBYPCZ(p[kp1],1.0-omega,omega,omega*Gamma*scale,p[km1],p[k]);CHKERRQ(ierr);
517
518 ktmp = km1;
519 km1 = k;
520 k = kp1;
521 kp1 = ktmp;
522 }
523 if (!ksp->reason) {
524 if (ksp->normtype) {
525 ierr = KSP_MatMult(ksp,Amat,p[k],r);CHKERRQ(ierr); /* r = b - Ap[k] */
526 ierr = VecAYPX(r,-1.0,b);CHKERRQ(ierr);
527 switch (ksp->normtype) {
528 case KSP_NORM_PRECONDITIONED:
529 ierr = KSP_PCApply(ksp,r,p[kp1]);CHKERRQ(ierr); /* p[kp1] = B^{-1}r */
530 ierr = VecNorm(p[kp1],NORM_2,&rnorm);CHKERRQ(ierr);
531 break;
532 case KSP_NORM_UNPRECONDITIONED:
533 case KSP_NORM_NATURAL:
534 ierr = VecNorm(r,NORM_2,&rnorm);CHKERRQ(ierr);
535 break;
536 default:
537 rnorm = 0.0;
538 break;
539 }
540 KSPCheckNorm(ksp,rnorm);
541 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr);
542 ksp->rnorm = rnorm;
543 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr);
544 ierr = KSPLogResidualHistory(ksp,rnorm);CHKERRQ(ierr);
545 ierr = KSPMonitor(ksp,i,rnorm);CHKERRQ(ierr);
546 }
547 if (ksp->its >= ksp->max_it) {
548 if (ksp->normtype != KSP_NORM_NONE) {
549 ierr = (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
550 if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
551 } else ksp->reason = KSP_CONVERGED_ITS;
552 }
553 }
554
555 /* make sure solution is in vector x */
556 ksp->vec_sol = sol_orig;
557 if (k) {
558 ierr = VecCopy(p[k],sol_orig);CHKERRQ(ierr);
559 }
560 PetscFunctionReturn(0);
561 }
562
KSPView_Chebyshev(KSP ksp,PetscViewer viewer)563 static PetscErrorCode KSPView_Chebyshev(KSP ksp,PetscViewer viewer)
564 {
565 KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
566 PetscErrorCode ierr;
567 PetscBool iascii;
568
569 PetscFunctionBegin;
570 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
571 if (iascii) {
572 ierr = PetscViewerASCIIPrintf(viewer," eigenvalue estimates used: min = %g, max = %g\n",(double)cheb->emin,(double)cheb->emax);CHKERRQ(ierr);
573 if (cheb->kspest) {
574 ierr = PetscViewerASCIIPrintf(viewer," eigenvalues estimate via %s min %g, max %g\n",((PetscObject)(cheb->kspest))->type_name,(double)cheb->emin_computed,(double)cheb->emax_computed);CHKERRQ(ierr);
575 ierr = PetscViewerASCIIPrintf(viewer," eigenvalues estimated using %s with translations [%g %g; %g %g]\n",((PetscObject) cheb->kspest)->type_name,(double)cheb->tform[0],(double)cheb->tform[1],(double)cheb->tform[2],(double)cheb->tform[3]);CHKERRQ(ierr);
576 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
577 ierr = KSPView(cheb->kspest,viewer);CHKERRQ(ierr);
578 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
579 if (cheb->usenoisy) {
580 ierr = PetscViewerASCIIPrintf(viewer," estimating eigenvalues using noisy right hand side\n");CHKERRQ(ierr);
581 }
582 }
583 }
584 PetscFunctionReturn(0);
585 }
586
KSPDestroy_Chebyshev(KSP ksp)587 static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
588 {
589 KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
590 PetscErrorCode ierr;
591
592 PetscFunctionBegin;
593 ierr = KSPDestroy(&cheb->kspest);CHKERRQ(ierr);
594 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",NULL);CHKERRQ(ierr);
595 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",NULL);CHKERRQ(ierr);
596 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",NULL);CHKERRQ(ierr);
597 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",NULL);CHKERRQ(ierr);
598 ierr = KSPDestroyDefault(ksp);CHKERRQ(ierr);
599 PetscFunctionReturn(0);
600 }
601
602 /*MC
603 KSPCHEBYSHEV - The preconditioned Chebyshev iterative method
604
605 Options Database Keys:
606 + -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
607 of the preconditioned operator. If these are accurate you will get much faster convergence.
608 . -ksp_chebyshev_esteig <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
609 transform for Chebyshev eigenvalue bounds (KSPChebyshevEstEigSet())
610 . -ksp_chebyshev_esteig_steps - number of estimation steps
611 - -ksp_chebyshev_esteig_noisy - use noisy number generator to create right hand side for eigenvalue estimator
612
613 Level: beginner
614
615 Notes:
616 The Chebyshev method requires both the matrix and preconditioner to
617 be symmetric positive (semi) definite.
618 Only support for left preconditioning.
619
620 Chebyshev is configured as a smoother by default, targetting the "upper" part of the spectrum.
621 The user should call KSPChebyshevSetEigenvalues() if they have eigenvalue estimates.
622
623 .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
624 KSPChebyshevSetEigenvalues(), KSPChebyshevEstEigSet(), KSPChebyshevEstEigSetUseNoisy()
625 KSPRICHARDSON, KSPCG, PCMG
626
627 M*/
628
KSPCreate_Chebyshev(KSP ksp)629 PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
630 {
631 PetscErrorCode ierr;
632 KSP_Chebyshev *chebyshevP;
633
634 PetscFunctionBegin;
635 ierr = PetscNewLog(ksp,&chebyshevP);CHKERRQ(ierr);
636
637 ksp->data = (void*)chebyshevP;
638 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);CHKERRQ(ierr);
639 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr);
640 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);CHKERRQ(ierr);
641 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);CHKERRQ(ierr);
642
643 chebyshevP->emin = 0.;
644 chebyshevP->emax = 0.;
645
646 chebyshevP->tform[0] = 0.0;
647 chebyshevP->tform[1] = 0.1;
648 chebyshevP->tform[2] = 0;
649 chebyshevP->tform[3] = 1.1;
650 chebyshevP->eststeps = 10;
651 chebyshevP->usenoisy = PETSC_TRUE;
652 ksp->setupnewmatrix = PETSC_TRUE;
653
654 ksp->ops->setup = KSPSetUp_Chebyshev;
655 ksp->ops->solve = KSPSolve_Chebyshev;
656 ksp->ops->destroy = KSPDestroy_Chebyshev;
657 ksp->ops->buildsolution = KSPBuildSolutionDefault;
658 ksp->ops->buildresidual = KSPBuildResidualDefault;
659 ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
660 ksp->ops->view = KSPView_Chebyshev;
661 ksp->ops->reset = KSPReset_Chebyshev;
662
663 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",KSPChebyshevSetEigenvalues_Chebyshev);CHKERRQ(ierr);
664 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",KSPChebyshevEstEigSet_Chebyshev);CHKERRQ(ierr);
665 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",KSPChebyshevEstEigSetUseNoisy_Chebyshev);CHKERRQ(ierr);
666 ierr = PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",KSPChebyshevEstEigGetKSP_Chebyshev);CHKERRQ(ierr);
667 PetscFunctionReturn(0);
668 }
669