Actual source code: idr.c
1: #include <petsc/private/kspimpl.h>
3: typedef struct {
4: PetscInt s; /* shadow space dimension; default 4 */
5: PetscObjectId workid; /* id of the work vectors used to initialize the shadow space */
6: Vec *GG; /* s direction vectors G[0..s-1] */
7: Vec *UU; /* s update vectors U[0..s-1] */
8: Vec *PP; /* s shadow vectors P[0..s-1] (fixed, random orthonormal) */
9: Vec r; /* current residual */
10: Vec v; /* work vector */
11: Vec t; /* work vector (preconditioned operator result) */
12: Vec guess; /* saved initial guess, used with right preconditioning */
13: PetscScalar *M; /* s*s matrix M[j,k] = <G[k],P[j]>, column-major */
14: PetscScalar *f; /* length s: P^T r */
15: PetscScalar *c; /* length s: solution of M c = f */
16: PetscReal cth; /* omega stabilization threshold (0 = off, default 0.7) */
17: PetscRandom rand; /* random context to initialize shadow vectors */
18: } KSP_IDR;
20: /*
21: KSPIDRInitShadowSpace_IDR - Fill shadow space P[0..s-1] with random
22: orthonormal vectors (Gram-Schmidt). Called when the work vectors are
23: created or recreated.
24: */
25: static PetscErrorCode KSPIDRInitShadowSpace_IDR(KSP ksp)
26: {
27: KSP_IDR *idr = (KSP_IDR *)ksp->data;
28: PetscScalar *dot;
29: PetscInt k;
31: PetscFunctionBegin;
32: PetscCall(KSPIDRGetRandom(ksp, &idr->rand));
33: PetscCall(PetscMalloc1(idr->s, &dot));
34: for (k = 0; k < idr->s; k++) PetscCall(VecSetRandom(idr->PP[k], idr->rand));
35: /* Gram-Schmidt orthonormalization */
36: for (k = 0; k < idr->s; k++) {
37: PetscCall((*ksp->orthog)(ksp, idr->PP, k, NULL, dot));
38: PetscCall(VecNormalize(idr->PP[k], NULL));
39: }
40: PetscCall(PetscFree(dot));
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
44: /*
45: KSPSetUp_IDR - Allocate the (3s+4) work vectors and the s*s+2s scalar
46: arrays, then initialize the shadow space P.
47: */
48: static PetscErrorCode KSPSetUp_IDR(KSP ksp)
49: {
50: KSP_IDR *idr = (KSP_IDR *)ksp->data;
52: PetscFunctionBegin;
53: PetscCall(KSPSetWorkVecs(ksp, 4 + 3 * idr->s));
54: idr->r = ksp->work[0];
55: idr->v = ksp->work[1];
56: idr->t = ksp->work[2];
57: idr->GG = ksp->work + 3;
58: idr->UU = ksp->work + 3 + idr->s;
59: idr->PP = ksp->work + 3 + 2 * idr->s;
60: idr->guess = ksp->work[3 + 3 * idr->s];
61: if (idr->M) PetscCall(PetscFree3(idr->M, idr->f, idr->c));
62: PetscCall(PetscMalloc3(idr->s * idr->s, &idr->M, idr->s, &idr->f, idr->s, &idr->c));
63: PetscCall(KSPIDRInitShadowSpace_IDR(ksp));
64: PetscCall(PetscObjectGetId((PetscObject)ksp->work[0], &idr->workid));
65: PetscFunctionReturn(PETSC_SUCCESS);
66: }
68: /*
69: KSPSolve_IDR - IDR(s) biorthogonal solve kernel.
71: This implements the biorthogonal IDR(s) recurrence (Algorithm 2 of
72: {cite}`gijzen:2011`) applied to the
73: preconditioned operator A' (= K^{-1}A for left, AK^{-1} for right
74: preconditioning, evaluated by KSP_PCApplyBAorAB()). Working on A'
75: keeps x and r consistent through left/right preconditioning and makes
76: IDR(1) reduce to BiCGSTAB.
78: The s-by-s matrix M (column-major, M[i + j*s] = p_i^H g_j) is kept
79: lower triangular: it is initialized to the identity and only its lower
80: part is updated, so the small system M[k:s-1,k:s-1] c = f[k:s-1] is a
81: forward substitution.
82: */
83: static PetscErrorCode KSPSolve_IDR(KSP ksp)
84: {
85: KSP_IDR *idr = (KSP_IDR *)ksp->data;
86: PetscInt s = idr->s, i, j, k;
87: PetscScalar *M = idr->M, *f = idr->f, *c = idr->c;
88: PetscScalar alpha, beta, om, tr, sum;
89: PetscReal dp = 0.0, nr, nt, rho;
90: PetscObjectId workid;
91: Vec X, B, R, V, T;
92: Vec *G, *U, *P;
94: PetscFunctionBegin;
95: PetscCheck(ksp->nwork == 4 + 3 * s, PetscObjectComm((PetscObject)ksp), PETSC_ERR_COR, "Unexpected number of work vectors %" PetscInt_FMT " != %" PetscInt_FMT, ksp->nwork, 4 + 3 * s);
96: idr->r = ksp->work[0];
97: idr->v = ksp->work[1];
98: idr->t = ksp->work[2];
99: idr->GG = ksp->work + 3;
100: idr->UU = ksp->work + 3 + s;
101: idr->PP = ksp->work + 3 + 2 * s;
102: idr->guess = ksp->work[3 + 3 * s];
103: PetscCall(PetscObjectGetId((PetscObject)ksp->work[0], &workid));
104: if (workid != idr->workid) {
105: PetscCall(KSPIDRInitShadowSpace_IDR(ksp));
106: idr->workid = workid;
107: }
108: X = ksp->vec_sol;
109: B = ksp->vec_rhs;
110: R = idr->r;
111: V = idr->v;
112: T = idr->t;
113: G = idr->GG;
114: U = idr->UU;
115: P = idr->PP;
116: nr = 1.0;
118: /* Compute initial (preconditioned for left PC) residual R */
119: PetscCall(KSPInitialResidual(ksp, X, V, T, R, B));
121: if (ksp->pc_side == PC_RIGHT && !ksp->guess_zero) {
122: PetscCall(VecCopy(X, idr->guess));
123: PetscCall(VecSet(X, 0.0));
124: }
126: if (ksp->normtype != KSP_NORM_NONE) {
127: PetscCall(VecNorm(R, NORM_2, &dp));
128: KSPCheckNorm(ksp, dp);
129: nr = dp;
130: }
131: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
132: ksp->its = 0;
133: ksp->rnorm = dp;
134: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
135: PetscCall(KSPLogResidualHistory(ksp, dp));
136: PetscCall(KSPMonitor(ksp, 0, dp));
137: PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP));
139: if (!ksp->reason) {
140: /* Initialize the IDR data: G = U = 0, M = I, omega = 1 */
141: for (k = 0; k < s; k++) {
142: PetscCall(VecSet(G[k], 0.0));
143: PetscCall(VecSet(U[k], 0.0));
144: }
145: PetscCall(PetscArrayzero(M, s * s));
146: for (k = 0; k < s; k++) M[k + k * s] = 1.0;
147: om = 1.0;
149: while (ksp->its < ksp->max_it && !ksp->reason) {
150: /* f = P^H r */
151: PetscCall(VecMDot(R, s, P, f));
153: for (k = 0; k < s; k++) {
154: /* Forward substitution: solve lower-triangular M[k:s-1,k:s-1] c[k:s-1] = f[k:s-1] */
155: for (i = k; i < s; i++) {
156: sum = f[i];
157: for (j = k; j < i; j++) sum -= M[i + j * s] * c[j];
158: c[i] = sum / M[i + i * s];
159: }
161: /* v = r - sum_{j=k}^{s-1} c[j] G[j] */
162: PetscCall(VecCopy(R, V));
163: for (j = k; j < s; j++) c[j] = -c[j];
164: PetscCall(VecMAXPY(V, s - k, c + k, G + k));
165: for (j = k; j < s; j++) c[j] = -c[j];
167: /* U[k] = omega*v + sum_{j=k}^{s-1} c[j] U[j] (scale U[k] by c[k] first to avoid aliasing) */
168: PetscCall(VecAXPBY(U[k], om, c[k], V));
169: if (s - k - 1 > 0) PetscCall(VecMAXPY(U[k], s - k - 1, c + k + 1, U + k + 1));
171: /* G[k] = A' U[k] */
172: PetscCall(KSP_PCApplyBAorAB(ksp, U[k], G[k], T));
174: /* Bi-orthogonalize G[k], U[k] against p_0,...,p_{k-1} */
175: for (i = 0; i < k; i++) {
176: PetscCall(VecDot(G[k], P[i], &alpha));
177: alpha /= M[i + i * s];
178: PetscCall(VecAXPY(G[k], -alpha, G[i]));
179: PetscCall(VecAXPY(U[k], -alpha, U[i]));
180: }
182: /* Update column k of M: M[k:s-1][k] = P[k:s-1]^H G[k] */
183: PetscCall(VecMDot(G[k], s - k, &P[k], &M[k + k * s]));
185: if (PetscAbsScalar(M[k + k * s]) < 10 * PETSC_MACHINE_EPSILON * nr) {
186: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve breakdown due to zero M[k,k] in IDR(s)");
187: ksp->reason = KSP_DIVERGED_BREAKDOWN;
188: PetscCall(PetscInfo(ksp, "Breakdown in IDR(s) half-step: M[k,k] = 0\n"));
189: break;
190: }
192: /* Make r orthogonal to p_0,...,p_k: r -= beta g_k, x += beta u_k */
193: beta = f[k] / M[k + k * s];
194: PetscCall(VecAXPY(R, -beta, G[k]));
195: PetscCall(VecAXPY(X, beta, U[k]));
197: /* With right preconditioning: R doubles as both the residual for x_0 and the
198: RHS for the shifted system A K^{-1} y = R iterated from y = 0 */
199: if (ksp->normtype != KSP_NORM_NONE) {
200: PetscCall(VecNorm(R, NORM_2, &dp));
201: KSPCheckNorm(ksp, dp);
202: }
203: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
204: ksp->its++;
205: ksp->rnorm = dp;
206: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
207: PetscCall(KSPLogResidualHistory(ksp, dp));
208: PetscCall(KSPMonitor(ksp, ksp->its, dp));
209: PetscCall((*ksp->converged)(ksp, ksp->its, dp, &ksp->reason, ksp->cnvP));
210: if (ksp->reason || ksp->its >= ksp->max_it) break;
212: /* Update remaining shadow projections: f[j] -= beta M[j][k], j > k */
213: for (i = k + 1; i < s; i++) f[i] -= beta * M[i + k * s];
214: }
215: if (ksp->reason || ksp->its >= ksp->max_it) break;
217: /* Minimal-residual (omega) step with angle stabilization.
218: Batch ||r||, ||t||, (r,t) into a single MPI collective via Begin/End. */
219: PetscCall(KSP_PCApplyBAorAB(ksp, R, T, V));
220: PetscCall(VecNormBegin(R, NORM_2, &nr));
221: PetscCall(VecNormBegin(T, NORM_2, &nt));
222: PetscCall(VecDotBegin(R, T, &tr));
223: PetscCall(VecNormEnd(R, NORM_2, &nr));
224: PetscCall(VecNormEnd(T, NORM_2, &nt));
225: PetscCall(VecDotEnd(R, T, &tr));
226: if (nt == 0.0) {
227: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve breakdown: zero ||A'r|| in IDR(s) omega step");
228: ksp->reason = KSP_DIVERGED_BREAKDOWN;
229: PetscCall(PetscInfo(ksp, "Breakdown in IDR(s) omega step: ||A'r|| = 0\n"));
230: break;
231: }
232: om = tr / (nt * nt);
233: if (idr->cth > 0.0) { /* 0 is a flag */
234: rho = PetscAbsScalar(tr) / (nt * nr);
235: if (rho < idr->cth) om *= idr->cth / rho;
236: }
237: PetscCall(VecAXPY(X, om, R));
238: PetscCall(VecAXPY(R, -om, T));
240: if (ksp->normtype != KSP_NORM_NONE) {
241: PetscCall(VecNorm(R, NORM_2, &dp));
242: KSPCheckNorm(ksp, dp);
243: }
244: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
245: ksp->its++;
246: ksp->rnorm = dp;
247: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
248: PetscCall(KSPLogResidualHistory(ksp, dp));
249: PetscCall(KSPMonitor(ksp, ksp->its, dp));
250: PetscCall((*ksp->converged)(ksp, ksp->its, dp, &ksp->reason, ksp->cnvP));
251: }
252: if (!ksp->reason && ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
253: }
255: /* Recover the true solution: unwind right preconditioning and add back guess */
256: PetscCall(KSPUnwindPreconditioner(ksp, X, T));
257: if (ksp->pc_side == PC_RIGHT && !ksp->guess_zero) PetscCall(VecAXPY(X, 1.0, idr->guess));
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: static PetscErrorCode KSPDestroy_IDR(KSP ksp)
262: {
263: KSP_IDR *idr = (KSP_IDR *)ksp->data;
265: PetscFunctionBegin;
266: PetscCall(PetscRandomDestroy(&idr->rand));
267: PetscCall(PetscFree3(idr->M, idr->f, idr->c));
268: PetscCall(KSPDestroyDefault(ksp));
269: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPIDRSetS_C", NULL));
270: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPIDRGetS_C", NULL));
271: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPIDRSetCosine_C", NULL));
272: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPIDRGetCosine_C", NULL));
273: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPIDRSetRandom_C", NULL));
274: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPIDRGetRandom_C", NULL));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: static PetscErrorCode KSPView_IDR(KSP ksp, PetscViewer viewer)
279: {
280: KSP_IDR *idr = (KSP_IDR *)ksp->data;
281: PetscBool isascii;
283: PetscFunctionBegin;
284: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
285: if (isascii) {
286: PetscCall(PetscViewerASCIIPrintf(viewer, " s (shadow space dimension) = %" PetscInt_FMT "\n", idr->s));
287: PetscCall(PetscViewerASCIIPrintf(viewer, " omega stabilization cosine threshold = %g\n", (double)idr->cth));
288: if (idr->rand) {
289: PetscCall(PetscViewerASCIIPushTab(viewer));
290: PetscCall(PetscRandomView(idr->rand, viewer));
291: PetscCall(PetscViewerASCIIPopTab(viewer));
292: }
293: }
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: static PetscErrorCode KSPSetFromOptions_IDR(KSP ksp, PetscOptionItems PetscOptionsObject)
298: {
299: KSP_IDR *idr = (KSP_IDR *)ksp->data;
300: PetscReal cth;
301: PetscInt s;
302: PetscBool flg;
304: PetscFunctionBegin;
305: PetscOptionsHeadBegin(PetscOptionsObject, "KSP IDR(s) options");
306: PetscCall(PetscOptionsBoundedInt("-ksp_idr_s", "Shadow space dimension", "KSPIDRSetS", idr->s, &s, &flg, 1));
307: if (flg) PetscCall(KSPIDRSetS(ksp, s));
308: PetscCall(PetscOptionsRangeReal("-ksp_idr_cosine", "Omega stabilization cosine threshold (0 = off)", "KSPIDRSetCosine", idr->cth, &cth, &flg, 0.0, 1.0));
309: if (flg) PetscCall(KSPIDRSetCosine(ksp, cth));
310: PetscOptionsHeadEnd();
311: PetscCall(KSPIDRGetRandom(ksp, &idr->rand));
312: PetscCall(PetscRandomSetFromOptions(idr->rand));
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: static PetscErrorCode KSPIDRSetS_IDR(KSP ksp, PetscInt s)
317: {
318: KSP_IDR *idr = (KSP_IDR *)ksp->data;
320: PetscFunctionBegin;
321: PetscCheck(s >= 1, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Shadow space dimension s must be >= 1, got %" PetscInt_FMT, s);
322: if (idr->s != s) {
323: idr->s = s;
324: ksp->setupstage = KSP_SETUP_NEW;
325: }
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: static PetscErrorCode KSPIDRGetS_IDR(KSP ksp, PetscInt *s)
330: {
331: KSP_IDR *idr = (KSP_IDR *)ksp->data;
333: PetscFunctionBegin;
334: *s = idr->s;
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: static PetscErrorCode KSPIDRSetCosine_IDR(KSP ksp, PetscReal cth)
339: {
340: KSP_IDR *idr = (KSP_IDR *)ksp->data;
342: PetscFunctionBegin;
343: PetscCheck(cth >= 0.0 && cth < 1.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Omega stabilization cosine threshold must be in [0,1), got %g", (double)cth);
344: idr->cth = cth;
345: PetscFunctionReturn(PETSC_SUCCESS);
346: }
348: static PetscErrorCode KSPIDRGetCosine_IDR(KSP ksp, PetscReal *cth)
349: {
350: KSP_IDR *idr = (KSP_IDR *)ksp->data;
352: PetscFunctionBegin;
353: *cth = idr->cth;
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: static PetscErrorCode KSPIDRSetRandom_IDR(KSP ksp, PetscRandom rand)
358: {
359: KSP_IDR *idr = (KSP_IDR *)ksp->data;
361: PetscFunctionBegin;
362: PetscCall(PetscObjectReference((PetscObject)rand));
363: PetscCall(PetscRandomDestroy(&idr->rand));
364: idr->rand = rand;
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: static PetscErrorCode KSPIDRGetRandom_IDR(KSP ksp, PetscRandom *rand)
369: {
370: KSP_IDR *idr = (KSP_IDR *)ksp->data;
372: PetscFunctionBegin;
373: if (!idr->rand) {
374: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)ksp), &idr->rand));
375: PetscCall(PetscObjectIncrementTabLevel((PetscObject)idr->rand, (PetscObject)ksp, 1));
376: PetscCall(PetscRandomSetOptionsPrefix(idr->rand, ((PetscObject)ksp)->prefix));
377: PetscCall(PetscRandomAppendOptionsPrefix(idr->rand, "ksp_idr_"));
378: PetscCall(PetscObjectSetOptions((PetscObject)idr->rand, ((PetscObject)ksp)->options));
379: }
380: *rand = idr->rand;
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: /*@
385: KSPIDRSetS - Sets the shadow space dimension s for the `KSPIDR` solver.
387: Logically Collective
389: Input Parameters:
390: + ksp - the Krylov solver context
391: - s - shadow space dimension (default 4); must be >= 1
393: Options Database Key:
394: . -ksp_idr_s s - shadow space dimension
396: Level: intermediate
398: Notes:
399: Increasing `s` generally improves convergence but requires `s` additional
400: vectors. If `s` is changed after `KSPSetUp()` has been called, the solver
401: is reset automatically.
403: .seealso: [](ch_ksp), `KSPIDR`, `KSPIDRGetS()`, `KSPIDRSetRandom()`
404: @*/
405: PetscErrorCode KSPIDRSetS(KSP ksp, PetscInt s)
406: {
407: PetscFunctionBegin;
410: PetscTryMethod(ksp, "KSPIDRSetS_C", (KSP, PetscInt), (ksp, s));
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: /*@
415: KSPIDRGetS - Gets the shadow space dimension s used by the `KSPIDR` solver.
417: Not Collective
419: Input Parameter:
420: . ksp - the Krylov solver context
422: Output Parameter:
423: . s - the shadow space dimension
425: Level: intermediate
427: .seealso: [](ch_ksp), `KSPIDR`, `KSPIDRSetS()`
428: @*/
429: PetscErrorCode KSPIDRGetS(KSP ksp, PetscInt *s)
430: {
431: PetscFunctionBegin;
433: PetscAssertPointer(s, 2);
434: PetscUseMethod(ksp, "KSPIDRGetS_C", (KSP, PetscInt *), (ksp, s));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: /*@
439: KSPIDRSetCosine - Sets the omega stabilization cosine threshold for the `KSPIDR` solver.
441: Logically Collective
443: Input Parameters:
444: + ksp - the Krylov solver context
445: - cth - stabilization cosine threshold in [0,1) (default 0.7, 0 = off)
447: Options Database Key:
448: . -ksp_idr_cosine cth - omega stabilization cosine threshold
450: Level: intermediate
452: Notes:
453: When the cosine of the angle between the residual and the preconditioned
454: residual drops below this threshold, omega is scaled to prevent the
455: near-orthogonality stalling described in {cite}`sleijpen:1993,sleijpen:1995`.
456: Setting `cth` to 0 disables stabilization.
458: .seealso: [](ch_ksp), `KSPIDR`, `KSPIDRGetCosine()`, `KSPIDRSetS()`, `KSPIDRSetRandom()`
459: @*/
460: PetscErrorCode KSPIDRSetCosine(KSP ksp, PetscReal cth)
461: {
462: PetscFunctionBegin;
465: PetscTryMethod(ksp, "KSPIDRSetCosine_C", (KSP, PetscReal), (ksp, cth));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: /*@
470: KSPIDRGetCosine - Gets the omega stabilization cosine threshold used by the `KSPIDR` solver.
472: Not Collective
474: Input Parameter:
475: . ksp - the Krylov solver context
477: Output Parameter:
478: . cth - the stabilization cosine threshold
480: Level: intermediate
482: .seealso: [](ch_ksp), `KSPIDR`, `KSPIDRSetCosine()`, `KSPIDRGetS()`
483: @*/
484: PetscErrorCode KSPIDRGetCosine(KSP ksp, PetscReal *cth)
485: {
486: PetscFunctionBegin;
488: PetscAssertPointer(cth, 2);
489: PetscUseMethod(ksp, "KSPIDRGetCosine_C", (KSP, PetscReal *), (ksp, cth));
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }
493: /*@
494: KSPIDRSetRandom - Sets the `PetscRandom` object used by the `KSPIDR` solver
495: to initialize the shadow vectors.
497: Collective
499: Input Parameters:
500: + ksp - the Krylov solver context
501: - rand - the random number generator context
503: Level: advanced
505: Note:
506: `KSPIDR` creates its own random number generator internally that can be accessed
507: with `KSPIDRGetRandom()` and controlled from the options database with the options
508: prefix of the `KSP` object.
510: .seealso: [](ch_ksp), `KSPIDR`, `KSPIDRGetRandom()`, `PetscRandomCreate()`, `KSPIDRSetS()`, `KSPIDRSetCosine()`
511: @*/
512: PetscErrorCode KSPIDRSetRandom(KSP ksp, PetscRandom rand)
513: {
514: PetscFunctionBegin;
517: PetscCheckSameComm(ksp, 1, rand, 2);
518: PetscTryMethod(ksp, "KSPIDRSetRandom_C", (KSP, PetscRandom), (ksp, rand));
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }
522: /*@
523: KSPIDRGetRandom - Gets the `PetscRandom` object used by the `KSPIDR` solver.
525: Collective
527: Input Parameter:
528: . ksp - the Krylov solver context
530: Output Parameter:
531: . rand - the random number generator context
533: Level: advanced
535: .seealso: [](ch_ksp), `KSPIDR`, `KSPIDRSetRandom()`, `KSPIDRSetCosine()`, `KSPIDRGetS()`
536: @*/
537: PetscErrorCode KSPIDRGetRandom(KSP ksp, PetscRandom *rand)
538: {
539: PetscFunctionBegin;
541: PetscAssertPointer(rand, 2);
542: PetscUseMethod(ksp, "KSPIDRGetRandom_C", (KSP, PetscRandom *), (ksp, rand));
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: /*MC
547: KSPIDR - IDR(s): Induced Dimension Reduction method for general nonsymmetric
548: linear systems {cite}`gijzen:2011`.
550: Options Database Keys:
551: + -ksp_idr_s s - shadow space dimension (default 4); larger `s` improves convergence
552: at the cost of `s` additional vectors and `s` extra inner products per
553: step, see `KSPIDRSetS()`
554: . -ksp_idr_cosine cth - omega stabilization cosine threshold (default 0.7, 0 = off); prevents
555: near-orthogonality stalling in the minimal-residual omega step
556: - -ksp_idr_random_type (rander48|rand|rand48|sprng|random123|curand) - set the random number generator
558: Level: intermediate
560: Notes:
561: IDR(s) is a short-recurrence, non-restarting Krylov method for general
562: nonsymmetric linear systems. It requires no growing subspace and avoids
563: the restart stagnation of `KSPGMRES`. The parameter `s` controls the
564: trade-off between memory and convergence speed\: s=1 is mathematically
565: equivalent to `KSPBCGS`; s=4 typically converges as fast as
566: GMRES(50); s=8 often outperforms GMRES(100).
567: Memory usage is (3s+3) vectors plus an s-by-s dense matrix.
568: This implements the biorthogonal variant described in {cite}`gijzen:2011`.
570: `KSPIDR` uses a `PetscRandom` which may be obtained with `KSPIDRGetRandom()`
571: (see also `KSPIDRSetRandom()`). The `PetscRandom` may be controlled from the
572: options database with the options prefix of the `KSP` object.
574: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
575: `KSPBCGS`, `KSPBCGSL`, `KSPGMRES`, `KSPIDRSetS()`, `KSPIDRGetS()`,
576: `KSPIDRSetCosine()`, `KSPIDRGetCosine()`, `KSPIDRSetRandom()`, `KSPIDRGetRandom()`
577: M*/
578: PETSC_EXTERN PetscErrorCode KSPCreate_IDR(KSP ksp)
579: {
580: KSP_IDR *idr;
582: PetscFunctionBegin;
583: PetscCall(PetscNew(&idr));
584: idr->s = 4;
585: idr->cth = 0.7;
586: ksp->data = (void *)idr;
588: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
589: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
590: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
591: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
593: ksp->ops->setup = KSPSetUp_IDR;
594: ksp->ops->solve = KSPSolve_IDR;
595: ksp->ops->destroy = KSPDestroy_IDR;
596: ksp->ops->view = KSPView_IDR;
597: ksp->ops->setfromoptions = KSPSetFromOptions_IDR;
598: ksp->ops->buildsolution = KSPBuildSolutionDefault;
599: ksp->ops->buildresidual = KSPBuildResidualDefault;
601: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPIDRSetS_C", KSPIDRSetS_IDR));
602: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPIDRGetS_C", KSPIDRGetS_IDR));
603: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPIDRSetCosine_C", KSPIDRSetCosine_IDR));
604: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPIDRGetCosine_C", KSPIDRGetCosine_IDR));
605: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPIDRSetRandom_C", KSPIDRSetRandom_IDR));
606: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPIDRGetRandom_C", KSPIDRGetRandom_IDR));
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }