Actual source code: eksm.c
1: /*
2: This file implements the Extended Krylov Subspace Method (EKSM)
3: for multiple shifted linear systems.
4: */
6: #include <../src/ksp/ksp/impls/eksm/eksmimpl.h>
7: #include <petscblaslapack.h>
9: #define EKSM_DELTA_DIRECTIONS 10
11: static PetscErrorCode KSPReset_EKSM(KSP);
13: /*
14: Compute the default shift varsigma as the centroid of all provided shifts sigma_i
15: */
16: static PetscErrorCode KSPEKSMSetDefaultShift(KSP ksp, PetscInt nshift, PetscScalar *sigma, PetscBool *cmplx)
17: {
18: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
20: PetscFunctionBegin;
21: eksm->shift = 0.0;
22: for (PetscInt i = 0; i < nshift; i++) {
23: eksm->shift += sigma[i];
24: if (SHIFT_IS_COMPLEX(cmplx, i)) { // Add complex conjugate pair
25: eksm->shift += sigma[i];
26: i++;
27: }
28: }
29: eksm->shift /= nshift;
30: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: static PetscErrorCode KSPSetUp_EKSM(KSP ksp)
34: {
35: PetscInt hh, hes, tt, sol, rs, max_k;
36: PetscBool ispcnone;
37: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
38: Mat Anest, Pnest, A, P;
39: Vec rhs, b;
40: Mat_MultiShift actx, pctx;
42: PetscFunctionBegin;
43: PetscCall(PetscObjectTypeCompare((PetscObject)ksp->pc, PCNONE, &ispcnone));
44: PetscCheck(ispcnone, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "In KSPEKSM the PC type must be PCNONE; use KSPEKSMGetKSP() to set an appropriate preconditioner");
46: PetscCall(KSPReset_EKSM(ksp));
47: PetscCall(PCGetOperators(ksp->pc, &Anest, &Pnest));
48: MatCheckMultiShift(Anest, &actx);
49: if (Pnest != Anest) MatCheckMultiShift(Pnest, &pctx);
51: PetscCall(MatCreateVecs(Anest, NULL, &rhs));
52: PetscCall(VecNestGetSubVec(rhs, 0, &b));
53: PetscCall(KSPEKSMGetKSP(ksp, &eksm->ksps, actx->M ? &eksm->kspm : NULL));
54: if (!eksm->shift_set) PetscCall(KSPEKSMSetDefaultShift(ksp, actx->nshift, actx->sigma, actx->cmplx));
56: PetscCall(MatMultiShiftBuildShiftedMatrix_Internal(actx->K, eksm->shift, actx->M, actx->str, PETSC_TRUE, &A));
57: P = A;
58: if (Pnest != Anest) PetscCall(MatMultiShiftBuildShiftedMatrix_Internal(pctx->K, eksm->shift, pctx->M, pctx->str, PETSC_TRUE, &P));
60: PetscCall(KSPSetOperators(eksm->ksps, A, P));
61: if (ksp->setfromoptionscalled) PetscCall(KSPSetFromOptions(eksm->ksps));
62: PetscCall(KSPSetUp(eksm->ksps));
63: PetscCall(MatDestroy(&A));
64: if (Pnest != Anest) PetscCall(MatDestroy(&P));
65: if (actx->M) {
66: if (Pnest != Anest) PetscCall(KSPSetOperators(eksm->kspm, actx->M, pctx->M));
67: else PetscCall(KSPSetOperators(eksm->kspm, actx->M, actx->M));
68: if (ksp->setfromoptionscalled) PetscCall(KSPSetFromOptions(eksm->kspm));
69: PetscCall(KSPSetUp(eksm->kspm));
70: }
72: eksm->factor = 1;
73: #if !PetscDefined(USE_COMPLEX)
74: // check for complex conjugate pairs of shifts
75: for (PetscInt k = 0; k < actx->nshift; k++) {
76: if (SHIFT_IS_COMPLEX(actx->cmplx, k)) {
77: eksm->factor = 2; // the projected matrix is twice as large
78: break;
79: }
80: }
81: #endif
83: max_k = 2 * ksp->max_it; /* two vectors per iteration */
84: eksm->ldh = max_k + 1;
85: eksm->ldk = max_k;
86: eksm->ldt = eksm->factor * (max_k + 1);
88: hh = eksm->ldh * max_k;
89: hes = eksm->ldk * max_k;
90: tt = eksm->ldt * eksm->factor * max_k;
91: sol = eksm->ldk * actx->nshift;
92: rs = eksm->ldt;
94: eksm->lwork = eksm->ldt + PetscMax(eksm->ldt, actx->nshift);
96: PetscCall(PetscCalloc6(hh, &eksm->hh_origin, hes, &eksm->kk_origin, tt, &eksm->tt_origin, sol, &eksm->yy_origin, rs, &eksm->rs_origin, eksm->lwork, &eksm->work));
98: /* Allocate array to hold pointers to basis vectors */
99: eksm->vecs_allocated = VEC_OFFSET + 2 + max_k;
100: PetscCall(PetscMalloc1(eksm->vecs_allocated, &eksm->vecs));
101: PetscCall(PetscMalloc1(VEC_OFFSET + 2 + max_k, &eksm->user_work));
102: PetscCall(PetscMalloc1(VEC_OFFSET + 2 + max_k, &eksm->mwork_alloc));
103: eksm->vv_allocated = VEC_OFFSET + 2 + PetscMin(2 * eksm->delta_allocate, max_k);
104: PetscCall(VecDuplicateVecs(b, eksm->vv_allocated, &eksm->user_work[0]));
105: eksm->mwork_alloc[0] = eksm->vv_allocated;
106: eksm->nwork_alloc = 1;
107: for (PetscInt k = 0; k < eksm->vv_allocated; k++) eksm->vecs[k] = eksm->user_work[0][k];
108: PetscCall(VecDestroy(&rhs));
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
112: /*
113: Allocate more work vectors, starting from VEC_VV(2*it)
114: */
115: static PetscErrorCode KSPEKSMGetNewVectors(KSP ksp, PetscInt it)
116: {
117: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
118: PetscInt nwork = eksm->nwork_alloc, nalloc;
120: PetscFunctionBegin;
121: nalloc = 2 * eksm->delta_allocate;
122: /* Adjust the number to allocate to make sure that we don't exceed the
123: number of available slots */
124: if (2 * it + VEC_OFFSET + nalloc >= eksm->vecs_allocated) nalloc = eksm->vecs_allocated - 2 * it - VEC_OFFSET;
125: if (!nalloc) PetscFunctionReturn(PETSC_SUCCESS);
127: eksm->vv_allocated += nalloc;
128: PetscCall(VecDuplicateVecs(VEC_VV(0), nalloc, &eksm->user_work[nwork]));
129: eksm->mwork_alloc[nwork] = nalloc;
130: for (PetscInt k = 0; k < nalloc; k++) eksm->vecs[2 * it + VEC_OFFSET + k] = eksm->user_work[nwork][k];
131: eksm->nwork_alloc++;
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
135: /*
136: Solve the projected problem for each shift
137: */
138: static PetscErrorCode KSPEKSMProjectedProblem(KSP ksp, PetscInt m, PetscReal *res, PetscInt nshift, PetscScalar *sigma, PetscBool *cmplx)
139: {
140: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
141: PetscScalar sone = 1.0, szero = 0.0;
142: PetscScalar *e;
143: PetscBLASInt n, n1, n2, n22, one = 1, ldk, ldt, lwork;
145: PetscFunctionBegin;
146: e = eksm->rs_origin;
147: PetscCall(PetscBLASIntCast(eksm->ldk, &ldk));
148: PetscCall(PetscBLASIntCast(eksm->ldt, &ldt));
149: PetscCall(PetscBLASIntCast(eksm->lwork, &lwork));
150: PetscCall(PetscBLASIntCast(m, &n));
151: PetscCall(PetscBLASIntCast(m + 1, &n1));
152: PetscCall(PetscBLASIntCast(2 * m, &n2));
153: PetscCall(PetscBLASIntCast(2 * m + 2, &n22));
154: *res = 0.0;
155: for (PetscInt k = 0; k < nshift; k++) {
156: PetscCall(PetscArrayzero(eksm->tt_origin, eksm->ldt * eksm->factor * m));
157: if (!SHIFT_IS_COMPLEX(cmplx, k)) {
158: /*
159: Usual case: real shift or complex arithmetic, process just one shift.
160: z = argmin || (H + sigma*K) z - e ||, y = K*z
161: */
162: PetscCall(PetscArrayzero(e, m + 1));
163: e[0] = eksm->v0norm;
164: /* T = H + sigma K */
165: for (PetscInt i = 0; i < m + 1; i++)
166: for (PetscInt j = 0; j < m; j++) *TT(i, j) = *HH(i, j);
167: for (PetscInt i = 0; i < m; i++)
168: for (PetscInt j = 0; j < m; j++) *TT(i, j) += sigma[k] * *KK(i, j);
169: /* z = T\e (least-squares) */
170: PetscCallLAPACKInfo("LAPACKgels", LAPACKgels_("N", &n1, &n, &one, TT(0, 0), &ldt, e, &ldt, eksm->work, &lwork, &info));
171: *res = PetscMax(*res, PetscAbsScalar(e[m]));
172: /* y(:,k) = K z */
173: PetscCallBLAS("BLASgemv", BLASgemv_("N", &n, &n, &sone, KK(0, 0), &ldk, e, &one, &szero, YY(0, k), &one));
174: } else {
175: /*
176: Complex conjugate pair of shifts in real arithmetic, process two iterations at once.
177: Same least squares problem as before, but coefficient matrix is
178: [ Re(H + sigma*K) -Im(H + sigma*K ] = [ H + Re(sigma)*K -Im(sigma)*K ]
179: [ Im(H + sigma*K) Re(H + sigma*K ] [ Im(sigma)*K H + Re(sigma)*K ]
180: */
181: PetscCall(PetscArrayzero(e, 2 * m + 2));
182: e[0] = eksm->v0norm;
183: /* form T */
184: for (PetscInt i = 0; i < m + 1; i++)
185: for (PetscInt j = 0; j < m; j++) {
186: *TT(i, j) = *HH(i, j);
187: *TT(i + m + 1, j + m) = *HH(i, j);
188: }
189: for (PetscInt i = 0; i < m; i++)
190: for (PetscInt j = 0; j < m; j++) {
191: *TT(i, j) += sigma[k] * *KK(i, j);
192: *TT(i, j + m) = -sigma[k + 1] * *KK(i, j);
193: *TT(i + m + 1, j) = sigma[k + 1] * *KK(i, j);
194: *TT(i + m + 1, j + m) += sigma[k] * *KK(i, j);
195: }
196: /* z = T\e (least-squares) */
197: PetscCallLAPACKInfo("LAPACKgels", LAPACKgels_("N", &n22, &n2, &one, TT(0, 0), &ldt, e, &ldt, eksm->work, &lwork, &info));
198: *res = PetscMax(*res, PetscHypotReal(PetscAbsScalar(e[2 * m]), PetscAbsScalar(e[2 * m + 1])));
199: /* y(:,k) = K z */
200: PetscCallBLAS("BLASgemv", BLASgemv_("N", &n, &n, &sone, KK(0, 0), &ldk, e, &one, &szero, YY(0, k), &one));
201: PetscCallBLAS("BLASgemv", BLASgemv_("N", &n, &n, &sone, KK(0, 0), &ldk, e + m, &one, &szero, YY(0, k + 1), &one));
202: k++; // skip next shift, result for conj(sigma) is conj(y)
203: }
204: }
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: /*
209: Build the solution vectors for each shift, i.e., subvectors of vec_sol
210: */
211: static PetscErrorCode KSPEKSMBuildSoln(KSP ksp, Vec vec_sol, PetscInt m, PetscBool *cmplx)
212: {
213: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
214: Vec x, xi;
215: PetscInt nshift;
217: PetscFunctionBegin;
218: PetscCall(VecSet(vec_sol, 0.0));
219: PetscCall(VecNestGetSize(vec_sol, &nshift));
220: for (PetscInt k = 0; k < nshift; k++) {
221: // x = V*y
222: PetscCall(VecNestGetSubVec(vec_sol, k, &x));
223: PetscCall(VecMAXPY(x, m, YY(0, k), &VEC_VV(0)));
224: if (SHIFT_IS_COMPLEX(cmplx, k)) { // imaginary part
225: PetscCall(VecNestGetSubVec(vec_sol, k + 1, &xi));
226: PetscCall(VecMAXPY(xi, m, YY(0, k + 1), &VEC_VV(0)));
227: k++; // skip next shift, result for conj(sigma) is conj(x)
228: }
229: }
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: /*
234: Stop the outer iteration if one of the internal linear solves failed
235: */
236: static PetscErrorCode KSPEKSMCheckInnerSolve(KSP ksp, KSP inner)
237: {
238: KSPConvergedReason reason;
240: PetscFunctionBegin;
241: PetscCall(KSPGetConvergedReason(inner, &reason));
242: if (reason < 0) {
243: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Internal KSPEKSM solve diverged with %s", KSPConvergedReasons[reason]);
244: PetscCall(PetscInfo(ksp, "Internal KSPEKSM solve diverged with %s\n", KSPConvergedReasons[reason]));
245: ksp->reason = KSP_DIVERGED_INNER_SOLVE_FAILED;
246: }
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: static PetscErrorCode KSPSolve_EKSM(KSP ksp)
251: {
252: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
253: PetscBool hapend = PETSC_FALSE;
254: PetscReal res, nrm;
255: PetscInt max_k, m = 0;
256: Mat Anest;
257: Vec b;
258: Mat_MultiShift actx;
260: PetscFunctionBegin;
261: PetscCheck(ksp->transpose_solve == PETSC_FALSE, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "This solver does not support transpose solve");
262: PetscCheck(ksp->guess_zero, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "This solver does not support nonzero initial guess");
263: PetscCall(VecNestGetSubVec(ksp->vec_rhs, 0, &b));
264: PetscCall(VecCopy(b, VEC_VV(0)));
265: PetscCall(PCGetOperators(ksp->pc, &Anest, NULL));
266: MatCheckMultiShift(Anest, &actx);
267: if (actx->M) {
268: PetscCall(VecCopy(VEC_VV(0), VEC_TEMP));
269: PetscCall(KSPSolve(eksm->kspm, VEC_TEMP, VEC_VV(0)));
270: PetscCall(KSPEKSMCheckInnerSolve(ksp, eksm->kspm));
271: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
272: }
273: PetscCall(VecNormalize(VEC_VV(0), &eksm->v0norm));
275: /*
276: Initial residual norm. The initial guess is zero, so this is the norm of the right-hand side of the
277: projected least-squares problem: ||b|| when M is the identity and ||M^{-1} b|| otherwise. Using it
278: keeps iteration 0 in the same norm as the projected residual reported at later iterations, so the
279: relative convergence test compares like with like.
280: */
281: res = eksm->v0norm;
282: KSPCheckNorm(ksp, res);
283: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
284: ksp->its = 0;
285: ksp->rnorm = res;
286: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
288: max_k = eksm->ldh - 1;
289: PetscCheck(max_k == 2 * ksp->max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_WRONGSTATE, "KSPSetTolerances() must be called before KSPSetUp() for KSPEKSM, and not between successive calls to KSPSolve()");
291: PetscCall(KSPLogResidualHistory(ksp, res));
292: PetscCall(KSPLogErrorHistory(ksp));
293: PetscCall(KSPMonitor(ksp, ksp->its, res));
294: if (!res) {
295: ksp->reason = KSP_CONVERGED_ATOL;
296: PetscCall(PetscInfo(ksp, "Converged due to zero residual norm on entry\n"));
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: /* check for the convergence */
301: PetscCall((*ksp->converged)(ksp, ksp->its, res, &ksp->reason, ksp->cnvP));
303: PetscCall(PetscArrayzero(eksm->hh_origin, eksm->ldh * max_k));
305: while (!ksp->reason && ksp->its < ksp->max_it) {
306: if (eksm->vv_allocated <= 2 * ksp->its + VEC_OFFSET + 2) PetscCall(KSPEKSMGetNewVectors(ksp, ksp->its + 1));
308: if (ksp->its == 0) { /* first iteration */
309: if (actx->M) PetscCall(MatMult(actx->M, VEC_VV(0), VEC_TEMP));
310: else PetscCall(VecCopy(VEC_VV(0), VEC_TEMP));
311: *HH(0, 0) = 1.0;
312: } else {
313: if (actx->M) PetscCall(MatMult(actx->M, VEC_VV(2 * ksp->its - 1), VEC_TEMP));
314: else PetscCall(VecCopy(VEC_VV(2 * ksp->its - 1), VEC_TEMP));
315: *HH(2 * ksp->its - 1, 2 * ksp->its) = 1.0;
316: }
317: PetscCall(KSPSolve(eksm->ksps, VEC_TEMP, VEC_VV(2 * ksp->its + 1)));
318: PetscCall(KSPEKSMCheckInnerSolve(ksp, eksm->ksps));
319: if (ksp->reason) break;
320: /* update Hessenberg matrix and do Gram-Schmidt */
321: PetscCall((*ksp->orthog)(ksp, &VEC_VV(0), 2 * ksp->its + 1, NULL, KK(0, 2 * ksp->its)));
322: if (ksp->reason) break;
324: /* vv(i+1) . vv(i+1) */
325: PetscCall(VecNormalize(VEC_VV(2 * ksp->its + 1), &nrm));
326: KSPCheckNorm(ksp, nrm);
328: /* save the magnitude */
329: *KK(2 * ksp->its + 1, 2 * ksp->its) = nrm;
330: for (PetscInt i = 0; i < 2 * ksp->its + 2; i++) *HH(i, 2 * ksp->its) -= eksm->shift * *KK(i, 2 * ksp->its);
332: /* check for the happy breakdown */
333: if (nrm <= eksm->haptol) {
334: PetscCall(PetscInfo(ksp, "Detected happy breakdown, nrm = %14.12e\n", (double)nrm));
335: hapend = PETSC_TRUE;
337: m = 2 * ksp->its + 1;
338: PetscCall(KSPEKSMProjectedProblem(ksp, m, &res, actx->nshift, actx->sigma, actx->cmplx));
339: } else {
340: PetscCall(MatMult(actx->K, VEC_VV(2 * ksp->its), VEC_TEMP));
341: if (actx->M) {
342: PetscCall(KSPSolve(eksm->kspm, VEC_TEMP, VEC_VV(2 * ksp->its + 2)));
343: PetscCall(KSPEKSMCheckInnerSolve(ksp, eksm->kspm));
344: if (ksp->reason) break;
345: } else PetscCall(VecCopy(VEC_TEMP, VEC_VV(2 * ksp->its + 2)));
346: *KK(2 * ksp->its, 2 * ksp->its + 1) = 1.0;
348: /* update Hessenberg matrix and do Gram-Schmidt */
349: PetscCall((*ksp->orthog)(ksp, &VEC_VV(0), 2 * ksp->its + 2, NULL, HH(0, 2 * ksp->its + 1)));
350: if (ksp->reason) break;
352: /* vv(i+1) . vv(i+1) */
353: PetscCall(VecNormalize(VEC_VV(2 * ksp->its + 2), &nrm));
354: KSPCheckNorm(ksp, nrm);
356: /* save the magnitude */
357: *HH(2 * ksp->its + 2, 2 * ksp->its + 1) = nrm;
359: /* check for the happy breakdown */
360: if (nrm <= eksm->haptol) {
361: PetscCall(PetscInfo(ksp, "Detected happy breakdown, nrm = %14.12e\n", (double)nrm));
362: hapend = PETSC_TRUE;
363: }
365: m = 2 * ksp->its + 2;
366: PetscCall(KSPEKSMProjectedProblem(ksp, m, &res, actx->nshift, actx->sigma, actx->cmplx));
367: }
369: ksp->its++;
370: ksp->rnorm = res;
371: if (ksp->reason) break;
373: PetscCall((*ksp->converged)(ksp, ksp->its, res, &ksp->reason, ksp->cnvP));
375: /* Catch error in happy breakdown and signal convergence and break from loop */
376: if (hapend) {
377: if (ksp->normtype == KSP_NORM_NONE) ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN;
378: else if (!ksp->reason) {
379: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Reached happy breakdown, but convergence was not indicated. Residual norm = %g", (double)res);
380: ksp->reason = KSP_DIVERGED_BREAKDOWN;
381: break;
382: }
383: }
384: PetscCall(KSPLogResidualHistory(ksp, res));
385: PetscCall(KSPLogErrorHistory(ksp));
386: PetscCall(KSPMonitor(ksp, ksp->its, res));
387: }
389: /* Form the solution */
390: PetscCall(KSPEKSMBuildSoln(ksp, ksp->vec_sol, m, actx->cmplx));
392: if (ksp->reason == KSP_CONVERGED_ITERATING && ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: static PetscErrorCode KSPReset_EKSM(KSP ksp)
397: {
398: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
400: PetscFunctionBegin;
401: /* Free the Hessenberg matrices */
402: PetscCall(PetscFree6(eksm->hh_origin, eksm->kk_origin, eksm->tt_origin, eksm->yy_origin, eksm->rs_origin, eksm->work));
404: /* free work vectors */
405: PetscCall(PetscFree(eksm->vecs));
406: for (PetscInt i = 0; i < eksm->nwork_alloc; i++) PetscCall(VecDestroyVecs(eksm->mwork_alloc[i], &eksm->user_work[i]));
407: eksm->nwork_alloc = 0;
409: PetscCall(PetscFree(eksm->user_work));
410: PetscCall(PetscFree(eksm->mwork_alloc));
412: if (eksm->ksps) PetscCall(KSPReset(eksm->ksps));
413: if (eksm->kspm) PetscCall(KSPReset(eksm->kspm));
415: eksm->vv_allocated = 0;
416: eksm->vecs_allocated = 0;
417: PetscFunctionReturn(PETSC_SUCCESS);
418: }
420: static PetscErrorCode KSPDestroy_EKSM(KSP ksp)
421: {
422: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
424: PetscFunctionBegin;
425: PetscCall(KSPReset_EKSM(ksp));
426: PetscCall(KSPDestroy(&eksm->ksps));
427: PetscCall(KSPDestroy(&eksm->kspm));
428: PetscCall(PetscFree(ksp->data));
429: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPEKSMSetHapTol_C", NULL));
430: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPEKSMGetHapTol_C", NULL));
431: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPEKSMSetKSP_C", NULL));
432: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPEKSMGetKSP_C", NULL));
433: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPEKSMSetShift_C", NULL));
434: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPEKSMGetShift_C", NULL));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: static PetscErrorCode KSPView_EKSM(KSP ksp, PetscViewer viewer)
439: {
440: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
441: const char *cstr;
442: PetscBool isascii, isstring;
444: PetscFunctionBegin;
445: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
446: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
447: if (ksp->orthog == KSPOrthogonalizationClassicalGramSchmidt) {
448: switch (ksp->cgstype) {
449: case KSP_ORTHOGONALIZATION_CGS_REFINE_NEVER:
450: cstr = "classical (unmodified) Gram-Schmidt orthogonalization with no iterative refinement";
451: break;
452: case KSP_ORTHOGONALIZATION_CGS_REFINE_ALWAYS:
453: cstr = "classical (unmodified) Gram-Schmidt orthogonalization with one step of iterative refinement";
454: break;
455: case KSP_ORTHOGONALIZATION_CGS_REFINE_IFNEEDED:
456: cstr = "classical (unmodified) Gram-Schmidt orthogonalization with one step of iterative refinement when needed";
457: break;
458: default:
459: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Unknown orthogonalization");
460: }
461: } else if (ksp->orthog == KSPOrthogonalizationModifiedGramSchmidt) cstr = "modified Gram-Schmidt orthogonalization";
462: else cstr = "unknown orthogonalization";
463: if (isascii) {
464: PetscCall(PetscViewerASCIIPrintf(viewer, " using %s\n", cstr));
465: #if !PetscDefined(USE_COMPLEX)
466: PetscCall(PetscViewerASCIIPrintf(viewer, " shift=%g\n", (double)eksm->shift));
467: #else
468: PetscCall(PetscViewerASCIIPrintf(viewer, " shift=%g%+gi\n", (double)PetscRealPart(eksm->shift), (double)PetscImaginaryPart(eksm->shift)));
469: #endif
470: PetscCall(PetscViewerASCIIPrintf(viewer, " happy breakdown tolerance=%g\n", (double)eksm->haptol));
471: if (eksm->ksps) {
472: PetscCall(PetscViewerASCIIPushTab(viewer));
473: PetscCall(KSPView(eksm->ksps, viewer));
474: PetscCall(PetscViewerASCIIPopTab(viewer));
475: }
476: if (eksm->kspm) {
477: PetscCall(PetscViewerASCIIPushTab(viewer));
478: PetscCall(KSPView(eksm->kspm, viewer));
479: PetscCall(PetscViewerASCIIPopTab(viewer));
480: }
481: } else if (isstring) {
482: #if !PetscDefined(USE_COMPLEX)
483: PetscCall(PetscViewerStringSPrintf(viewer, "%s shift %g", cstr, (double)eksm->shift));
484: #else
485: PetscCall(PetscViewerStringSPrintf(viewer, "%s shift %g%+gi", cstr, (double)PetscRealPart(eksm->shift), (double)PetscImaginaryPart(eksm->shift)));
486: #endif
487: }
488: PetscFunctionReturn(PETSC_SUCCESS);
489: }
491: static PetscErrorCode KSPSetFromOptions_EKSM(KSP ksp, PetscOptionItems PetscOptionsObject)
492: {
493: PetscScalar shift;
494: PetscReal haptol;
495: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
496: PetscBool flg;
498: PetscFunctionBegin;
499: PetscOptionsHeadBegin(PetscOptionsObject, "KSP EKSM Options");
500: PetscCall(PetscOptionsScalar("-ksp_eksm_shift", "Scalar shift for the internal linear solves", "KSPEKSMSetShift", eksm->shift, &shift, &flg));
501: if (flg) PetscCall(KSPEKSMSetShift(ksp, shift));
502: PetscCall(PetscOptionsReal("-ksp_eksm_haptol", "Tolerance for exact convergence (happy breakdown)", "KSPEKSMSetHapTol", eksm->haptol, &haptol, &flg));
503: if (flg) PetscCall(KSPEKSMSetHapTol(ksp, haptol));
504: PetscOptionsHeadEnd();
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: static PetscErrorCode KSPEKSMSetHapTol_EKSM(KSP ksp, PetscReal tol)
509: {
510: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
512: PetscFunctionBegin;
513: PetscCheck(tol >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Tolerance must be non-negative");
514: eksm->haptol = tol;
515: PetscFunctionReturn(PETSC_SUCCESS);
516: }
518: /*@
519: KSPEKSMSetHapTol - Sets the tolerance for detecting a happy breakdown in EKSM.
521: Logically Collective
523: Input Parameters:
524: + ksp - the Krylov space solver context
525: - tol - the tolerance for detecting a happy breakdown
527: Options Database Key:
528: . -ksp_eksm_haptol tol - set tolerance for determining happy breakdown
530: Level: intermediate
532: Notes:
533: Happy breakdown is the rare case in `KSPEKSM` where a very near zero matrix entry is generated in the upper Hessenberg matrix indicating
534: an 'exact' solution has been obtained.
536: The default tolerance value for detecting a happy breakdown with EKSM in PETSc is 1.0e-30.
538: .seealso: [](ch_ksp), `KSPEKSM`, `KSPSetTolerances()`, `KSPGMRESSetHapTol()`
539: @*/
540: PetscErrorCode KSPEKSMSetHapTol(KSP ksp, PetscReal tol)
541: {
542: PetscFunctionBegin;
545: PetscTryMethod(ksp, "KSPEKSMSetHapTol_C", (KSP, PetscReal), (ksp, tol));
546: PetscFunctionReturn(PETSC_SUCCESS);
547: }
549: static PetscErrorCode KSPEKSMGetHapTol_EKSM(KSP ksp, PetscReal *tol)
550: {
551: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
553: PetscFunctionBegin;
554: *tol = eksm->haptol;
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: /*@
559: KSPEKSMGetHapTol - Gets the tolerance for detecting a happy breakdown in EKSM.
561: Not Collective
563: Input Parameter:
564: . ksp - the Krylov space solver context
566: Output Parameter:
567: . tol - the tolerance for detecting a happy breakdown
569: Level: intermediate
571: .seealso: [](ch_ksp), `KSPEKSM`, `KSPEKSMSetHapTol()`
572: @*/
573: PetscErrorCode KSPEKSMGetHapTol(KSP ksp, PetscReal *tol)
574: {
575: PetscFunctionBegin;
577: PetscAssertPointer(tol, 2);
578: PetscUseMethod(ksp, "KSPEKSMGetHapTol_C", (KSP, PetscReal *), (ksp, tol));
579: PetscFunctionReturn(PETSC_SUCCESS);
580: }
582: static PetscErrorCode KSPEKSMSetKSP_EKSM(KSP ksp, KSP ksps, KSP kspm)
583: {
584: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
586: PetscFunctionBegin;
587: if (ksps) {
588: PetscCall(PetscObjectReference((PetscObject)ksps));
589: PetscCall(PetscObjectDereference((PetscObject)eksm->ksps));
590: eksm->ksps = ksps;
591: }
592: if (kspm) {
593: PetscCall(PetscObjectReference((PetscObject)kspm));
594: PetscCall(PetscObjectDereference((PetscObject)eksm->kspm));
595: eksm->kspm = kspm;
596: }
597: if (ksp->setupstage) ksp->setupstage = KSP_SETUP_NEWMATRIX;
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: /*@
602: KSPEKSMSetKSP - Sets the internal `KSP` objects to be used by EKSM.
604: Not Collective, but all `KSP` objects must live on the same `MPI_Comm`
606: Input Parameters:
607: + ksp - the Krylov space solver context
608: . ksps - the `KSP` context used internally for linear systems with coefficient matrix $K + \varsigma M$, or `NULL` to keep the current one
609: - kspm - the `KSP` context used internally for linear systems with coefficient matrix $M$, or `NULL` to keep the current one
611: Level: intermediate
613: Notes:
614: This routine is rarely needed. The most common usage is to call `KSPEKSMGetKSP()` to extract the internal
615: objects and then configure them.
617: Each `KSP` object actually passed in has its reference count increased by one and the internal object it
618: replaces has its reference count decreased by one. An argument passed as `NULL` leaves the corresponding
619: internal object, and its reference count, unchanged.
621: .seealso: [](ch_ksp), `KSPEKSM`, `KSPEKSMGetKSP()`
622: @*/
623: PetscErrorCode KSPEKSMSetKSP(KSP ksp, KSP ksps, KSP kspm)
624: {
625: PetscFunctionBegin;
627: if (ksps) {
629: PetscCheckSameComm(ksp, 1, ksps, 2);
630: }
631: if (kspm) {
633: PetscCheckSameComm(ksp, 1, kspm, 3);
634: }
635: PetscTryMethod(ksp, "KSPEKSMSetKSP_C", (KSP, KSP, KSP), (ksp, ksps, kspm));
636: PetscFunctionReturn(PETSC_SUCCESS);
637: }
639: static PetscErrorCode KSPEKSMGetKSP_EKSM(KSP ksp, KSP *ksps, KSP *kspm)
640: {
641: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
643: PetscFunctionBegin;
644: if (ksps) {
645: if (!eksm->ksps) {
646: PetscCall(KSPCreate(PetscObjectComm((PetscObject)ksp), &eksm->ksps));
647: PetscCall(PetscObjectIncrementTabLevel((PetscObject)eksm->ksps, (PetscObject)ksp, 1));
648: PetscCall(PetscObjectSetOptions((PetscObject)eksm->ksps, ((PetscObject)ksp)->options));
649: PetscCall(KSPSetOptionsPrefix(eksm->ksps, ((PetscObject)ksp)->prefix));
650: PetscCall(KSPAppendOptionsPrefix(eksm->ksps, "eksm_s_"));
651: PetscCall(KSPSetTolerances(eksm->ksps, ksp->rtol, ksp->abstol, ksp->divtol, PETSC_DETERMINE));
652: }
653: *ksps = eksm->ksps;
654: }
655: if (kspm) {
656: if (!eksm->kspm) {
657: PetscCall(KSPCreate(PetscObjectComm((PetscObject)ksp), &eksm->kspm));
658: PetscCall(PetscObjectIncrementTabLevel((PetscObject)eksm->kspm, (PetscObject)ksp, 1));
659: PetscCall(PetscObjectSetOptions((PetscObject)eksm->kspm, ((PetscObject)ksp)->options));
660: PetscCall(KSPSetOptionsPrefix(eksm->kspm, ((PetscObject)ksp)->prefix));
661: PetscCall(KSPAppendOptionsPrefix(eksm->kspm, "eksm_m_"));
662: PetscCall(KSPSetTolerances(eksm->kspm, ksp->rtol, ksp->abstol, ksp->divtol, PETSC_DETERMINE));
663: }
664: *kspm = eksm->kspm;
665: }
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: /*@
670: KSPEKSMGetKSP - Returns the internal `KSP` objects used by EKSM.
672: Not Collective, but if `ksp` is parallel, then `ksps` and `kspm` are parallel
674: Input Parameter:
675: . ksp - the Krylov space solver context
677: Output Parameters:
678: + ksps - the `KSP` context used internally for linear systems with coefficient matrix $K + \varsigma M$
679: - kspm - the `KSP` context used internally for linear systems with coefficient matrix $M$
681: Level: intermediate
683: Notes:
684: The extended Krylov subspace method needs to solve linear systems with matrices $K + \varsigma M$ as well as $M$ (in case $M$ is not the
685: identity). Here, $K$ and $M$ are the matrices provided in `MatCreateNestFromMultipleShifts()` and $\varsigma$ is the shift
686: set with `KSPEKSMSetShift()`.
688: This function is provided to allow the user to configure the internal `KSP` objects that will do most of the work.
689: Pass `NULL` to any of the arguments if not needed. In particular, if $M=I$ the `kspm` object is not needed.
691: To configure the linear solvers from the command line, use the `-eksm_s_` and `-eksm_m_` prefixes, respectively.
692: For instance `-ksp_type eksm -eksm_s_ksp_type bcgs -eksm_s_pc_type bjacobi -eksm_m_pc_type none`.
694: .seealso: [](ch_ksp), `KSPEKSM`, `MatCreateNestFromMultipleShifts()`, `KSPEKSMSetShift()`, `KSPEKSMSetKSP()`
695: @*/
696: PetscErrorCode KSPEKSMGetKSP(KSP ksp, KSP *ksps, KSP *kspm)
697: {
698: PetscFunctionBegin;
700: PetscUseMethod(ksp, "KSPEKSMGetKSP_C", (KSP, KSP *, KSP *), (ksp, ksps, kspm));
701: PetscFunctionReturn(PETSC_SUCCESS);
702: }
704: static PetscErrorCode KSPEKSMSetShift_EKSM(KSP ksp, PetscScalar shift)
705: {
706: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
708: PetscFunctionBegin;
709: eksm->shift_set = PETSC_TRUE;
710: if (!ksp->setupstage) eksm->shift = shift;
711: else if (eksm->shift != shift) {
712: eksm->shift = shift;
713: ksp->setupstage = KSP_SETUP_NEW;
714: /* free the data structures, then create them again */
715: PetscCall(KSPReset_EKSM(ksp));
716: }
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: /*@
721: KSPEKSMSetShift - Sets the value of the shift $\varsigma$ that will be used for the internal linear solves in the EKSM solver.
723: Logically Collective
725: Input Parameters:
726: + ksp - the Krylov space solver context
727: - shift - the value of the shift, $\varsigma$
729: Options Database Key:
730: . -ksp_eksm_shift shift - the value of the shift
732: Level: intermediate
734: Notes:
735: The extended Krylov subspace method needs to solve linear systems with the matrix $K + \varsigma M$ using an auxiliary `KSP` object,
736: see `KSPEKSMGetKSP()`. The scalar value $\varsigma$ is provided with this function. If not set, the default
737: is to use a value of $\varsigma$ equal to the centroid (or center of mass) of all shifts $\sigma_i$ provided
738: by the user in `MatCreateNestFromMultipleShifts()`.
740: .seealso: [](ch_ksp), `KSPEKSM`, `KSPEKSMGetKSP()`, `KSPEKSMGetShift()`, `MatCreateNestFromMultipleShifts()`
741: @*/
742: PetscErrorCode KSPEKSMSetShift(KSP ksp, PetscScalar shift)
743: {
744: PetscFunctionBegin;
747: PetscTryMethod(ksp, "KSPEKSMSetShift_C", (KSP, PetscScalar), (ksp, shift));
748: PetscFunctionReturn(PETSC_SUCCESS);
749: }
751: static PetscErrorCode KSPEKSMGetShift_EKSM(KSP ksp, PetscScalar *shift)
752: {
753: KSP_EKSM *eksm = (KSP_EKSM *)ksp->data;
755: PetscFunctionBegin;
756: *shift = eksm->shift;
757: PetscFunctionReturn(PETSC_SUCCESS);
758: }
760: /*@
761: KSPEKSMGetShift - Gets the value of the shift $\varsigma$ used for the internal linear solves.
763: Not Collective
765: Input Parameter:
766: . ksp - the Krylov space solver context
768: Output Parameter:
769: . shift - the value of the shift, $\varsigma$
771: Level: intermediate
773: Note:
774: The returned value is the one set with `KSPEKSMSetShift()`. If that routine was not called, the shift is
775: computed internally during `KSPSetUp()` as the centroid of the shifts passed to `MatCreateNestFromMultipleShifts()`,
776: so this routine returns 0 when called before the solver has been set up.
778: .seealso: [](ch_ksp), `KSPEKSM`, `KSPEKSMSetShift()`
779: @*/
780: PetscErrorCode KSPEKSMGetShift(KSP ksp, PetscScalar *shift)
781: {
782: PetscFunctionBegin;
784: PetscAssertPointer(shift, 2);
785: PetscUseMethod(ksp, "KSPEKSMGetShift_C", (KSP, PetscScalar *), (ksp, shift));
786: PetscFunctionReturn(PETSC_SUCCESS);
787: }
789: /*MC
790: KSPEKSM - Implements the Extended Krylov Subspace Method (EKSM) for solving multiple shifted linear systems simultaneously using `KSP`.
792: Options Database Keys:
793: + -ksp_eksm_shift shift - the scalar shift for the internal linear solves
794: - -ksp_eksm_haptol tol - the tolerance for happy breakdown (exact convergence) of `KSPEKSM`
796: Level: beginner
798: Note:
799: This solver should be used with no preconditioning, `PCNONE`. It is possible to use a preconditioner in the internal
800: linear solvers, see `KSPEKSMGetKSP()`. To configure the internal linear solvers from the command line, use the
801: `-eksm_s_` and `-eksm_m_` prefixes, respectively.
802: For instance `-ksp_type eksm -eksm_s_ksp_type bcgs -eksm_s_pc_type bjacobi -eksm_m_pc_type none`.
804: Unlike restarted methods such as `KSPGMRES`, this solver keeps the whole extended Krylov basis, so the maximum
805: iteration count set with `KSPSetTolerances()` also fixes the size of the workspace allocated in `KSPSetUp()`\:
806: two basis vectors per iteration and $O(\textrm{max\_it}^2)$ scalars of dense storage. Choose it accordingly, and
807: set it before `KSPSetUp()`; changing it afterwards makes the next `KSPSolve()` fail.
809: When a mass matrix $M$ is given to `MatCreateNestFromMultipleShifts()`, the residual norm used by the convergence
810: test and reported by `-ksp_monitor` is $\|M^{-1}(b-(K+\sigma_i M)x_i)\|$ maximized over the shifts, not the
811: unpreconditioned residual norm.
813: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPEKSM`, `KSPEKSMSetHapTol()`, `KSPEKSMSetKSP()`, `KSPEKSMSetShift()`,
814: `MatCreateNestFromMultipleShifts()`, `MatCreateVecNestFromMultipleShifts()`
815: M*/
817: PETSC_EXTERN PetscErrorCode KSPCreate_EKSM(KSP ksp)
818: {
819: KSP_EKSM *eksm;
820: PC pc;
822: PetscFunctionBegin;
823: PetscCall(PetscNew(&eksm));
824: ksp->data = (void *)eksm;
826: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 4));
827: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 3));
828: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_SYMMETRIC, 2));
829: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
830: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
832: PetscObjectParameterSetDefault(ksp, max_it, 50);
834: ksp->setupnewmatrix = PETSC_TRUE;
836: ksp->ops->setup = KSPSetUp_EKSM;
837: ksp->ops->solve = KSPSolve_EKSM;
838: ksp->ops->reset = KSPReset_EKSM;
839: ksp->ops->destroy = KSPDestroy_EKSM;
840: ksp->ops->view = KSPView_EKSM;
841: ksp->ops->setfromoptions = KSPSetFromOptions_EKSM;
843: PetscCall(KSPGetPC(ksp, &pc));
844: PetscCall(PCSetType(pc, PCNONE));
846: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPEKSMSetHapTol_C", KSPEKSMSetHapTol_EKSM));
847: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPEKSMGetHapTol_C", KSPEKSMGetHapTol_EKSM));
848: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPEKSMSetKSP_C", KSPEKSMSetKSP_EKSM));
849: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPEKSMGetKSP_C", KSPEKSMGetKSP_EKSM));
850: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPEKSMSetShift_C", KSPEKSMSetShift_EKSM));
851: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPEKSMGetShift_C", KSPEKSMGetShift_EKSM));
853: eksm->haptol = 1.0e-30;
854: eksm->delta_allocate = EKSM_DELTA_DIRECTIONS;
855: PetscFunctionReturn(PETSC_SUCCESS);
856: }