Actual source code: rich.c
1: /*
2: This implements Richardson Iteration.
3: */
4: #include <../src/ksp/ksp/impls/rich/richardsonimpl.h>
6: static PetscErrorCode KSPSetUp_Richardson(KSP ksp)
7: {
8: KSP_Richardson *richardsonP = (KSP_Richardson *)ksp->data;
10: PetscFunctionBegin;
11: if (richardsonP->selfscale) PetscCall(KSPSetWorkVecs(ksp, 4));
12: else PetscCall(KSPSetWorkVecs(ksp, 2));
13: PetscFunctionReturn(PETSC_SUCCESS);
14: }
16: static PetscErrorCode KSPSolve_Richardson(KSP ksp)
17: {
18: PetscReal rnorm = 0.0, abr;
19: PetscScalar scale, rdot;
20: Vec x, b, r, z, w = NULL, y = NULL;
21: PetscInt i, maxit;
22: Mat Amat, Pmat;
23: KSP_Richardson *richardsonP = (KSP_Richardson *)ksp->data;
24: PetscBool exists;
25: MatNullSpace nullsp;
27: PetscFunctionBegin;
28: ksp->its = 0;
29: PetscCheck(ksp->nwork == (richardsonP->selfscale ? 4 : 2), PetscObjectComm((PetscObject)ksp), PETSC_ERR_COR, "Unexpected number of work vectors %" PetscInt_FMT " != %d", ksp->nwork, richardsonP->selfscale ? 4 : 2);
31: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
32: x = ksp->vec_sol;
33: b = ksp->vec_rhs;
34: r = ksp->work[0];
35: z = ksp->work[1];
36: if (richardsonP->selfscale) {
37: w = ksp->work[2];
38: y = ksp->work[3];
39: }
40: maxit = ksp->max_it;
42: /* if user has provided fast Richardson code use that */
43: PetscCall(PCApplyRichardsonExists(ksp->pc, &exists));
44: PetscCall(MatGetNullSpace(Pmat, &nullsp));
45: if (exists && maxit > 0 && richardsonP->scale == 1.0 && (ksp->converged == KSPConvergedDefault || ksp->converged == KSPConvergedSkip) && !ksp->numbermonitors && !ksp->transpose_solve && !nullsp) {
46: PCRichardsonConvergedReason reason;
47: PetscCall(PCApplyRichardson(ksp->pc, b, x, r, ksp->rtol, ksp->abstol, ksp->divtol, maxit, ksp->guess_zero, &ksp->its, &reason));
48: ksp->reason = (KSPConvergedReason)reason;
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: if (!ksp->guess_zero) { /* r <- b - A x */
53: PetscCall(KSP_MatMult(ksp, Amat, x, r));
54: PetscCall(VecAYPX(r, -1.0, b));
55: } else {
56: PetscCall(VecCopy(b, r));
57: }
59: ksp->its = 0;
60: if (richardsonP->selfscale) {
61: PetscCall(KSP_PCApply(ksp, r, z)); /* z <- B r */
62: for (i = 0; i < maxit; i++) {
63: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
64: PetscCall(VecNorm(r, NORM_2, &rnorm)); /* rnorm <- r'*r */
65: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
66: PetscCall(VecNorm(z, NORM_2, &rnorm)); /* rnorm <- z'*z */
67: } else rnorm = 0.0;
69: KSPCheckNorm(ksp, rnorm);
70: ksp->rnorm = rnorm;
71: PetscCall(KSPMonitor(ksp, i, rnorm));
72: PetscCall(KSPLogResidualHistory(ksp, rnorm));
73: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
74: if (ksp->reason) break;
75: PetscCall(KSP_PCApplyBAorAB(ksp, z, y, w)); /* y = BAz = BABr */
76: PetscCall(VecDotNorm2(z, y, &rdot, &abr)); /* rdot = (Br)^T(BABR); abr = (BABr)^T (BABr) */
77: scale = rdot / abr;
78: PetscCall(PetscInfo(ksp, "Self-scale factor %g\n", (double)PetscRealPart(scale)));
79: PetscCall(VecAXPY(x, scale, z)); /* x <- x + scale z */
80: PetscCall(VecAXPY(r, -scale, w)); /* r <- r - scale*Az */
81: PetscCall(VecAXPY(z, -scale, y)); /* z <- z - scale*y */
82: ksp->its++;
83: }
84: } else {
85: for (i = 0; i < maxit; i++) {
86: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
87: PetscCall(VecNorm(r, NORM_2, &rnorm)); /* rnorm <- r'*r */
88: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
89: PetscCall(KSP_PCApply(ksp, r, z)); /* z <- B r */
90: PetscCall(VecNorm(z, NORM_2, &rnorm)); /* rnorm <- z'*z */
91: } else rnorm = 0.0;
92: ksp->rnorm = rnorm;
93: PetscCall(KSPMonitor(ksp, i, rnorm));
94: PetscCall(KSPLogResidualHistory(ksp, rnorm));
95: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
96: if (ksp->reason) break;
97: if (ksp->normtype != KSP_NORM_PRECONDITIONED) PetscCall(KSP_PCApply(ksp, r, z)); /* z <- B r */
99: PetscCall(VecAXPY(x, richardsonP->scale, z)); /* x <- x + scale z */
100: ksp->its++;
102: if (i + 1 < maxit || ksp->normtype != KSP_NORM_NONE) {
103: PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r <- b - Ax */
104: PetscCall(VecAYPX(r, -1.0, b));
105: }
106: }
107: }
108: if (!ksp->reason) {
109: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
110: PetscCall(VecNorm(r, NORM_2, &rnorm)); /* rnorm <- r'*r */
111: } else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
112: PetscCall(KSP_PCApply(ksp, r, z)); /* z <- B r */
113: PetscCall(VecNorm(z, NORM_2, &rnorm)); /* rnorm <- z'*z */
114: } else rnorm = 0.0;
116: KSPCheckNorm(ksp, rnorm);
117: ksp->rnorm = rnorm;
118: PetscCall(KSPLogResidualHistory(ksp, rnorm));
119: PetscCall(KSPMonitor(ksp, i, rnorm));
120: if (ksp->its >= ksp->max_it) {
121: if (ksp->normtype != KSP_NORM_NONE) {
122: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
123: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
124: } else {
125: ksp->reason = KSP_CONVERGED_ITS;
126: }
127: }
128: }
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: /*
133: Replacement for KSPCheckNorm() in KSPMatSolve_Richardson(), where ksp->vec_sol is not available since the solution is a block of vectors
134: */
135: static PetscErrorCode KSPMatSolveCheckNorm_Richardson(KSP ksp, PetscReal rnorm, Mat X)
136: {
137: PCFailedReason pcreason;
139: PetscFunctionBegin;
140: if (!PetscIsInfOrNanReal(rnorm)) PetscFunctionReturn(PETSC_SUCCESS);
141: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPMatSolve%s() has not converged due to infinity or NaN norm", ksp->transpose.solve_requested ? "Transpose" : "");
142: PetscCall(PCReduceFailedReason(ksp->pc));
143: PetscCall(PCGetFailedReason(ksp->pc, &pcreason));
144: /* as with VecFlag() in KSPCheckNorm(), the state of the block of solutions increases whether or not it is flagged, so that an outer solver detects the failure */
145: PetscCall(MatFlag(X, pcreason));
146: ksp->reason = pcreason ? KSP_DIVERGED_PC_FAILED : KSP_DIVERGED_NANORINF;
147: ksp->rnorm = rnorm;
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: /*
152: Block analog of KSP_PCApply(), the null space of Amat is removed from each column of the block of preconditioned vectors as in KSP_RemoveNullSpaceMat()
153: */
154: static PetscErrorCode KSPMatSolvePCMatApply_Richardson(KSP ksp, Mat R, Mat Z)
155: {
156: PetscFunctionBegin;
157: PetscCall(KSP_PCMatApply(ksp, R, Z));
158: PetscCall(KSP_RemoveNullSpaceMat(ksp, Z));
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: /*
163: Releases the work blocks of vectors cached by the block iteration of KSPMatSolve_Richardson(), together with the product cached in R, but not the work vector
164: of its PCApplyRichardson() fallback, which the two paths do not share
165: */
166: static PetscErrorCode KSPMatSolveResetBlocks_Richardson(KSP ksp)
167: {
168: KSP_Richardson *richardsonP = (KSP_Richardson *)ksp->data;
170: PetscFunctionBegin;
171: PetscCall(MatDestroy(&richardsonP->R));
172: PetscCall(MatDestroy(&richardsonP->Z));
173: richardsonP->transpose = PETSC_FALSE;
174: richardsonP->id = 0;
175: richardsonP->state = 0;
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: static PetscErrorCode KSPReset_Richardson(KSP ksp)
180: {
181: KSP_Richardson *richardsonP = (KSP_Richardson *)ksp->data;
183: PetscFunctionBegin;
184: PetscCall(KSPMatSolveResetBlocks_Richardson(ksp));
185: PetscCall(VecDestroy(&richardsonP->w));
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: static PetscErrorCode KSPMatSolve_Richardson(KSP ksp, Mat B, Mat X)
190: {
191: PetscReal rnorm = 0.0;
192: Mat Amat, Pmat, R, Z;
193: Vec cb, cx;
194: PetscInt i, maxit, m, mn, n, nn, N, NN;
195: KSP_Richardson *richardsonP = (KSP_Richardson *)ksp->data;
196: PetscBool exists, matexists, match = PETSC_FALSE, reuse = PETSC_FALSE;
197: PetscObjectId id;
198: PetscObjectState state;
199: MatNullSpace nullsp;
201: PetscFunctionBegin;
202: ksp->its = 0;
203: ksp->reason = KSP_CONVERGED_ITERATING;
204: maxit = ksp->max_it;
205: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
207: /* if user has provided fast Richardson code use that, with the same conditions as in KSPSolve_Richardson() except for the monitors, which are not called during KSPMatSolve() */
208: PetscCall(PCApplyRichardsonExists(ksp->pc, &exists));
209: PetscCall(PCMatApplyRichardsonExists(ksp->pc, &matexists));
210: PetscCall(MatGetNullSpace(Pmat, &nullsp));
211: if ((exists || matexists) && maxit > 0 && richardsonP->scale == 1.0 && (ksp->converged == KSPConvergedDefault || ksp->converged == KSPConvergedSkip) && !ksp->transpose_solve && !nullsp) {
212: PCRichardsonConvergedReason reason;
214: /* the block iteration is not run in this call, so the work blocks of vectors it may have cached in an earlier one are released */
215: PetscCall(KSPMatSolveResetBlocks_Richardson(ksp));
216: if (matexists) {
217: PetscCall(PetscInfo(ksp, "Using PCMatApplyRichardson() on each batch of right-hand sides, by default the whole block\n"));
218: PetscCall(PCMatApplyRichardson(ksp->pc, B, X, NULL, ksp->rtol, ksp->abstol, ksp->divtol, maxit, ksp->guess_zero, &ksp->its, &reason));
219: ksp->reason = (KSPConvergedReason)reason;
220: } else {
221: /* the block iteration below does not compute the same thing as PCApplyRichardson(), so the right-hand sides are solved one at a time to match KSPSolve() */
222: PetscCall(PetscInfo(ksp, "Using PCApplyRichardson() on one right-hand side at a time\n"));
223: PetscCall(MatGetSize(B, NULL, &N));
224: PetscCall(MatGetLocalSize(B, &m, NULL));
225: /* ksp->work may have the size of an earlier operator since setup is not run again, so a work vector is cached here instead, it is rebuilt only when the
226: layout or the type of the columns of the block of right-hand sides changes, which a column of that block reports directly */
227: if (N > 0) {
228: PetscCall(MatDenseGetColumnVecRead(B, 0, &cb));
229: if (richardsonP->w) { /* a process without a cached work vector contributes PETSC_FALSE to the reduction below */
230: PetscCall(VecGetLocalSize(richardsonP->w, &n));
231: if (n == m) PetscCall(PetscObjectTypeCompare((PetscObject)richardsonP->w, ((PetscObject)cb)->type_name, &match));
232: }
233: /* the local size above may match on some processes only, while VecDestroy() and VecDuplicate() below are collective, so the verdict must be the same on all of them */
234: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &match, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)ksp)));
235: if (!match) PetscCall(VecDestroy(&richardsonP->w));
236: if (!richardsonP->w) PetscCall(VecDuplicate(cb, &richardsonP->w));
237: PetscCall(MatDenseRestoreColumnVecRead(B, 0, &cb));
238: }
239: /* as with the column-by-column fallback of KSPMatSolve_Private(), ksp->reason and ksp->its reflect the last right-hand side */
240: for (i = 0; i < N; i++) {
241: PetscCall(MatDenseGetColumnVecRead(B, i, &cb));
242: if (ksp->guess_zero) PetscCall(MatDenseGetColumnVecWrite(X, i, &cx));
243: else PetscCall(MatDenseGetColumnVec(X, i, &cx));
244: PetscCall(PCApplyRichardson(ksp->pc, cb, cx, richardsonP->w, ksp->rtol, ksp->abstol, ksp->divtol, maxit, ksp->guess_zero, &ksp->its, &reason));
245: if (ksp->guess_zero) PetscCall(MatDenseRestoreColumnVecWrite(X, i, &cx));
246: else PetscCall(MatDenseRestoreColumnVec(X, i, &cx));
247: PetscCall(MatDenseRestoreColumnVecRead(B, i, &cb));
248: ksp->reason = (KSPConvergedReason)reason;
249: }
250: }
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: PetscCall(PetscInfo(ksp, "Iterating on each batch of right-hand sides, by default the whole block\n"));
255: /* R and Z are cached between calls, they are rebuilt only when the shape or the type of the block of right-hand sides, the direction of the solve, or the nonzero pattern of the operator changes,
256: the local column layout is part of the shape since MatDenseGetSubMatrix() may distribute two batches of the same global width differently */
257: PetscCall(MatGetLocalSize(B, &m, &mn));
258: PetscCall(MatGetSize(B, NULL, &N));
259: PetscCall(PetscObjectGetId((PetscObject)Amat, &id));
260: PetscCall(MatGetNonzeroState(Amat, &state));
261: if (richardsonP->R) { /* a process without cached work blocks of vectors contributes PETSC_FALSE to the reduction below */
262: PetscCall(MatGetLocalSize(richardsonP->R, &n, &nn));
263: PetscCall(MatGetSize(richardsonP->R, NULL, &NN));
264: if (n == m && nn == mn && NN == N) PetscCall(PetscObjectTypeCompare((PetscObject)richardsonP->R, ((PetscObject)B)->type_name, &reuse));
265: if (richardsonP->transpose != ksp->transpose_solve || richardsonP->id != id || richardsonP->state != state) reuse = PETSC_FALSE;
266: }
267: /* the local sizes above may match on some processes only, while both branches below are collective, so the verdict must be the same on all of them */
268: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &reuse, 1, MPI_C_BOOL, MPI_LAND, PetscObjectComm((PetscObject)ksp)));
269: if (reuse) {
270: PetscCall(PetscInfo(ksp, "Reusing the cached work blocks of vectors and the symbolic phase of the product\n"));
271: /* the cached product is bound to Z between calls, so it is bound back to the block of solutions of this call */
272: PetscCall(MatProductReplaceMats(NULL, X, NULL, richardsonP->R));
273: } else {
274: PetscCall(KSPMatSolveResetBlocks_Richardson(ksp));
275: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &richardsonP->R));
276: PetscCall(MatDuplicate(B, MAT_DO_NOT_COPY_VALUES, &richardsonP->Z));
277: /* set the product A X (or A^T X) up once so that only its numeric phase is run in the loop below and in the following calls */
278: PetscCall(MatProductCreateWithMat(Amat, X, NULL, richardsonP->R));
279: PetscCall(MatProductSetType(richardsonP->R, ksp->transpose_solve ? MATPRODUCT_AtB : MATPRODUCT_AB));
280: PetscCall(MatProductSetFromOptions(richardsonP->R));
281: PetscCall(MatProductSymbolic(richardsonP->R));
282: richardsonP->transpose = ksp->transpose_solve;
283: richardsonP->id = id;
284: richardsonP->state = state;
285: }
286: R = richardsonP->R;
287: Z = richardsonP->Z;
289: if (!ksp->guess_zero) { /* R <- B - A X */
290: PetscCall(MatProductNumeric(R));
291: PetscCall(MatAYPX(R, -1.0, B, SAME_NONZERO_PATTERN));
292: } else PetscCall(MatCopy(B, R, SAME_NONZERO_PATTERN));
294: for (i = 0; i < maxit; i++) {
295: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) PetscCall(MatNorm(R, NORM_FROBENIUS, &rnorm)); /* rnorm <- ||R||_F */
296: else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
297: PetscCall(KSPMatSolvePCMatApply_Richardson(ksp, R, Z)); /* Z <- B R */
298: PetscCall(MatNorm(Z, NORM_FROBENIUS, &rnorm)); /* rnorm <- ||Z||_F */
299: } else rnorm = 0.0;
301: PetscCall(KSPMatSolveCheckNorm_Richardson(ksp, rnorm, X));
302: if (ksp->reason) break;
303: ksp->rnorm = rnorm;
304: PetscCall(KSPLogResidualHistory(ksp, rnorm));
305: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
306: if (ksp->reason) break;
307: if (ksp->normtype != KSP_NORM_PRECONDITIONED) PetscCall(KSPMatSolvePCMatApply_Richardson(ksp, R, Z)); /* Z <- B R */
309: PetscCall(MatAXPY(X, richardsonP->scale, Z, SAME_NONZERO_PATTERN)); /* X <- X + scale Z */
310: ksp->its++;
312: if (i + 1 < maxit || ksp->normtype != KSP_NORM_NONE) {
313: PetscCall(MatProductNumeric(R)); /* R <- B - A X */
314: PetscCall(MatAYPX(R, -1.0, B, SAME_NONZERO_PATTERN));
315: }
316: }
317: if (!ksp->reason) {
318: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) PetscCall(MatNorm(R, NORM_FROBENIUS, &rnorm)); /* rnorm <- ||R||_F */
319: else if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
320: PetscCall(KSPMatSolvePCMatApply_Richardson(ksp, R, Z)); /* Z <- B R */
321: PetscCall(MatNorm(Z, NORM_FROBENIUS, &rnorm)); /* rnorm <- ||Z||_F */
322: } else rnorm = 0.0;
324: PetscCall(KSPMatSolveCheckNorm_Richardson(ksp, rnorm, X));
325: if (!ksp->reason) {
326: ksp->rnorm = rnorm;
327: PetscCall(KSPLogResidualHistory(ksp, rnorm));
328: if (ksp->its >= ksp->max_it) {
329: if (ksp->normtype != KSP_NORM_NONE) {
330: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
331: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
332: } else ksp->reason = KSP_CONVERGED_ITS;
333: }
334: }
335: }
336: /* rebind the cached product to a work block owned by the KSP so that it does not keep the block of solutions of this call alive until the next one, Z has the same type and layout as X so the symbolic phase is not run again */
337: PetscCall(MatProductReplaceMats(NULL, Z, NULL, R));
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: static PetscErrorCode KSPView_Richardson(KSP ksp, PetscViewer viewer)
342: {
343: KSP_Richardson *richardsonP = (KSP_Richardson *)ksp->data;
344: PetscBool isascii;
346: PetscFunctionBegin;
347: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
348: if (isascii) {
349: if (richardsonP->selfscale) {
350: PetscCall(PetscViewerASCIIPrintf(viewer, " using self-scale best computed damping factor\n"));
351: } else {
352: PetscCall(PetscViewerASCIIPrintf(viewer, " damping factor=%g\n", (double)richardsonP->scale));
353: }
354: }
355: PetscFunctionReturn(PETSC_SUCCESS);
356: }
358: static PetscErrorCode KSPSetFromOptions_Richardson(KSP ksp, PetscOptionItems PetscOptionsObject)
359: {
360: KSP_Richardson *rich = (KSP_Richardson *)ksp->data;
361: PetscReal tmp;
362: PetscBool flg, flg2;
364: PetscFunctionBegin;
365: PetscOptionsHeadBegin(PetscOptionsObject, "KSP Richardson Options");
366: PetscCall(PetscOptionsReal("-ksp_richardson_scale", "damping factor", "KSPRichardsonSetScale", rich->scale, &tmp, &flg));
367: if (flg) PetscCall(KSPRichardsonSetScale(ksp, tmp));
368: PetscCall(PetscOptionsBool("-ksp_richardson_self_scale", "dynamically determine optimal damping factor", "KSPRichardsonSetSelfScale", rich->selfscale, &flg2, &flg));
369: if (flg) PetscCall(KSPRichardsonSetSelfScale(ksp, flg2));
370: PetscOptionsHeadEnd();
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: static PetscErrorCode KSPDestroy_Richardson(KSP ksp)
375: {
376: PetscFunctionBegin;
377: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetScale_C", NULL));
378: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetSelfScale_C", NULL));
379: PetscCall(KSPDestroyDefault(ksp));
380: PetscFunctionReturn(PETSC_SUCCESS);
381: }
383: static PetscErrorCode KSPRichardsonSetScale_Richardson(KSP ksp, PetscReal scale)
384: {
385: KSP_Richardson *richardsonP;
387: PetscFunctionBegin;
388: richardsonP = (KSP_Richardson *)ksp->data;
389: richardsonP->scale = scale;
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: static PetscErrorCode KSPRichardsonSetSelfScale_Richardson(KSP ksp, PetscBool selfscale)
394: {
395: KSP_Richardson *richardsonP;
397: PetscFunctionBegin;
398: richardsonP = (KSP_Richardson *)ksp->data;
399: if (richardsonP->selfscale != selfscale) {
400: /* KSPSetUp_Richardson() picks the number of work vectors from this flag, so setup must be run again */
401: ksp->setupstage = KSP_SETUP_NEW;
402: /* neither the block iteration nor its PCApplyRichardson() fallback is reachable once the self-scaled variant is on, so their cached work space is released */
403: if (selfscale) PetscCall(KSPReset_Richardson(ksp));
404: }
405: richardsonP->selfscale = selfscale;
406: /* the self-scaled variant has no block analog, so KSPMatSolve() falls back to solving the right-hand sides one at a time */
407: ksp->ops->matsolve = selfscale ? NULL : KSPMatSolve_Richardson;
408: PetscFunctionReturn(PETSC_SUCCESS);
409: }
411: static PetscErrorCode KSPBuildResidual_Richardson(KSP ksp, Vec t, Vec v, Vec *V)
412: {
413: PetscFunctionBegin;
414: if (ksp->normtype == KSP_NORM_NONE) {
415: PetscCall(KSPBuildResidualDefault(ksp, t, v, V));
416: } else {
417: PetscCall(VecCopy(ksp->work[0], v));
418: *V = v;
419: }
420: PetscFunctionReturn(PETSC_SUCCESS);
421: }
423: /*MC
424: KSPRICHARDSON - The preconditioned Richardson iterative method {cite}`richarson1911`
426: Options Database Key:
427: . -ksp_richardson_scale - damping factor on the correction (defaults to 1.0)
429: Level: beginner
431: Notes:
432: $ x^{n+1} = x^{n} + scale*B(b - A x^{n})$
434: Here B is the application of the preconditioner
436: This method often (usually) will not converge unless scale is very small.
438: For some preconditioners, currently `PCSOR`, the convergence test is skipped to improve speed,
439: thus it always iterates the maximum number of iterations you've selected. When -ksp_monitor
440: (or any other monitor) is turned on, the norm is computed at each iteration and so the convergence test is run unless
441: you specifically call `KSPSetNormType`(ksp,`KSP_NORM_NONE`);
443: For some preconditioners, currently `PCMG` and `PCHYPRE` with BoomerAMG if -ksp_monitor (and also
444: any other monitor) is not turned on then the convergence test is done by the preconditioner itself and
445: so the solver may run more or fewer iterations then if -ksp_monitor is selected.
447: Supports only left preconditioning
449: If using direct solvers such as `PCLU` and `PCCHOLESKY` one generally uses `KSPPREONLY` instead of this which uses exactly one iteration
451: `-ksp_type richardson -pc_type jacobi` gives one classical Jacobi preconditioning
453: `KSPMatSolve()` and `KSPMatSolveTranspose()` are supported natively, that is, the iteration is performed on a whole batch of right-hand sides at once, a batch being
454: the whole block unless `KSPSetMatSolveBatchSize()` is used, except with `KSPRichardsonSetSelfScale()`, for which the right-hand sides are solved one at a time.
455: The convergence of each batch is tested on the Frobenius norm of its block of (preconditioned) residuals and, with a nonzero initial guess, the relative tolerance
456: is by default based on the Frobenius norm of its block of (preconditioned) right-hand sides, as in `KSPSolve()`. `KSPConvergedDefaultSetUIRNorm()` can be used to
457: base it on the initial residual norm instead. Unlike `KSPSolve()`, `KSPMatSolve()` does not remove the (transpose) null space of the operator from the block of
458: right-hand sides, so the caller must make a singular system consistent by projecting `B` itself
460: As in `KSPSolve()`, the iteration is delegated to the preconditioner when it provides a fast Richardson code, the convergence test is the default one or is skipped,
461: and `Pmat` has no null space. `PCMatApplyRichardson()` is used when the preconditioner provides it, otherwise `PCApplyRichardson()` is applied
462: to one right-hand side at a time so that `KSPMatSolve()` computes the same solutions as `KSPSolve()`. Since `KSPMatSolve()` is a stripped-down version of `KSPSolve()`,
463: monitors are not called during a block iteration and so, unlike in `KSPSolve()`, they do not prevent this delegation
465: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
466: `KSPRichardsonSetScale()`, `KSPPREONLY`, `KSPRichardsonSetSelfScale()`, `KSPMatSolve()`, `PCApplyRichardson()`, `PCMatApplyRichardson()`
467: M*/
469: PETSC_EXTERN PetscErrorCode KSPCreate_Richardson(KSP ksp)
470: {
471: KSP_Richardson *richardsonP;
473: PetscFunctionBegin;
474: PetscCall(PetscNew(&richardsonP));
475: ksp->data = (void *)richardsonP;
477: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
478: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
479: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
481: ksp->ops->setup = KSPSetUp_Richardson;
482: ksp->ops->solve = KSPSolve_Richardson;
483: ksp->ops->matsolve = KSPMatSolve_Richardson;
484: ksp->ops->reset = KSPReset_Richardson;
485: ksp->ops->destroy = KSPDestroy_Richardson;
486: ksp->ops->buildsolution = KSPBuildSolutionDefault;
487: ksp->ops->buildresidual = KSPBuildResidual_Richardson;
488: ksp->ops->view = KSPView_Richardson;
489: ksp->ops->setfromoptions = KSPSetFromOptions_Richardson;
491: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetScale_C", KSPRichardsonSetScale_Richardson));
492: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPRichardsonSetSelfScale_C", KSPRichardsonSetSelfScale_Richardson));
494: richardsonP->scale = 1.0;
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }