Actual source code: dgmres.c
1: /*
2: Implements deflated GMRES.
3: */
5: #include <../src/ksp/ksp/impls/gmres/dgmres/dgmresimpl.h>
7: PetscLogEvent KSP_DGMRESComputeDeflationData, KSP_DGMRESApplyDeflation;
9: static PetscErrorCode KSPDGMRESGetNewVectors(KSP, PetscInt);
10: static PetscErrorCode KSPDGMRESUpdateHessenberg(KSP, PetscInt, PetscBool, PetscReal *);
11: static PetscErrorCode KSPDGMRESBuildSoln(PetscScalar *, Vec, Vec, KSP, PetscInt);
13: static PetscErrorCode KSPDGMRESSetEigen(KSP ksp, PetscInt nb_eig)
14: {
15: PetscFunctionBegin;
16: PetscTryMethod((ksp), "KSPDGMRESSetEigen_C", (KSP, PetscInt), (ksp, nb_eig));
17: PetscFunctionReturn(PETSC_SUCCESS);
18: }
19: static PetscErrorCode KSPDGMRESSetMaxEigen(KSP ksp, PetscInt max_neig)
20: {
21: PetscFunctionBegin;
22: PetscTryMethod((ksp), "KSPDGMRESSetMaxEigen_C", (KSP, PetscInt), (ksp, max_neig));
23: PetscFunctionReturn(PETSC_SUCCESS);
24: }
25: static PetscErrorCode KSPDGMRESComputeSchurForm(KSP ksp, PetscInt *neig)
26: {
27: PetscFunctionBegin;
28: PetscUseMethod((ksp), "KSPDGMRESComputeSchurForm_C", (KSP, PetscInt *), (ksp, neig));
29: PetscFunctionReturn(PETSC_SUCCESS);
30: }
31: PetscErrorCode KSPDGMRESComputeDeflationData(KSP ksp, PetscInt *curneigh)
32: {
33: PetscFunctionBegin;
34: PetscUseMethod((ksp), "KSPDGMRESComputeDeflationData_C", (KSP, PetscInt *), (ksp, curneigh));
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
37: static PetscErrorCode KSPDGMRESApplyDeflation(KSP ksp, Vec x, Vec y)
38: {
39: PetscFunctionBegin;
40: PetscUseMethod((ksp), "KSPDGMRESApplyDeflation_C", (KSP, Vec, Vec), (ksp, x, y));
41: PetscFunctionReturn(PETSC_SUCCESS);
42: }
44: static PetscErrorCode KSPDGMRESImproveEig(KSP ksp, PetscInt neig)
45: {
46: PetscFunctionBegin;
47: PetscUseMethod((ksp), "KSPDGMRESImproveEig_C", (KSP, PetscInt), (ksp, neig));
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: PetscErrorCode KSPSetUp_DGMRES(KSP ksp)
52: {
53: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
54: PetscInt neig = dgmres->neig + EIG_OFFSET;
55: PetscInt max_k = dgmres->max_k + 1;
57: PetscFunctionBegin;
58: PetscCall(KSPSetUp_GMRES(ksp));
59: if (!dgmres->neig) PetscFunctionReturn(PETSC_SUCCESS);
61: /* Allocate workspace for the Schur vectors*/
62: PetscCall(PetscMalloc1(neig * max_k, &SR));
63: dgmres->wr = NULL;
64: dgmres->wi = NULL;
65: dgmres->perm = NULL;
66: dgmres->modul = NULL;
67: dgmres->Q = NULL;
68: dgmres->Z = NULL;
70: UU = NULL;
71: XX = NULL;
72: MX = NULL;
73: AUU = NULL;
74: XMX = NULL;
75: XMU = NULL;
76: UMX = NULL;
77: AUAU = NULL;
78: TT = NULL;
79: TTF = NULL;
80: INVP = NULL;
81: X1 = NULL;
82: X2 = NULL;
83: MU = NULL;
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: /*
88: Run GMRES, possibly with restart. Return residual history if requested.
89: input parameters:
91: . gmres - structure containing parameters and work areas
93: output parameters:
94: . nres - residuals (from preconditioned system) at each step.
95: If restarting, consider passing nres+it. If null,
96: ignored
97: . itcount - number of iterations used. nres[0] to nres[itcount]
98: are defined. If null, ignored.
100: Notes:
101: On entry, the value in vector VEC_VV(0) should be the initial residual
102: (this allows shortcuts where the initial preconditioned residual is 0).
103: */
104: static PetscErrorCode KSPDGMRESCycle(PetscInt *itcount, KSP ksp)
105: {
106: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
107: PetscReal res_norm, res, hapbnd, tt;
108: PetscInt it = 0;
109: PetscInt max_k = dgmres->max_k;
110: PetscBool hapend = PETSC_FALSE;
111: PetscReal res_old;
112: PetscInt test = 0;
114: PetscFunctionBegin;
115: PetscCall(VecNormalize(VEC_VV(0), &res_norm));
116: KSPCheckNorm(ksp, res_norm);
117: res = res_norm;
118: *GRS(0) = res_norm;
120: /* check for the convergence */
121: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
122: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res;
123: else ksp->rnorm = 0.0;
124: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
125: dgmres->it = (it - 1);
126: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
127: PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
128: if (!res) {
129: if (itcount) *itcount = 0;
130: ksp->reason = KSP_CONVERGED_ATOL;
131: PetscCall(PetscInfo(ksp, "Converged due to zero residual norm on entry\n"));
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
134: /* record the residual norm to test if deflation is needed */
135: res_old = res;
137: PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP));
138: while (!ksp->reason && it < max_k && ksp->its < ksp->max_it) {
139: if (it) {
140: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
141: PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
142: }
143: dgmres->it = (it - 1);
144: if (dgmres->vv_allocated <= it + VEC_OFFSET + 1) PetscCall(KSPDGMRESGetNewVectors(ksp, it + 1));
145: if (dgmres->r > 0) {
146: if (ksp->pc_side == PC_LEFT) {
147: /* Apply the first preconditioner */
148: PetscCall(KSP_PCApplyBAorAB(ksp, VEC_VV(it), VEC_TEMP, VEC_TEMP_MATOP));
149: /* Then apply Deflation as a preconditioner */
150: PetscCall(KSPDGMRESApplyDeflation(ksp, VEC_TEMP, VEC_VV(1 + it)));
151: } else if (ksp->pc_side == PC_RIGHT) {
152: PetscCall(KSPDGMRESApplyDeflation(ksp, VEC_VV(it), VEC_TEMP));
153: PetscCall(KSP_PCApplyBAorAB(ksp, VEC_TEMP, VEC_VV(1 + it), VEC_TEMP_MATOP));
154: }
155: } else {
156: PetscCall(KSP_PCApplyBAorAB(ksp, VEC_VV(it), VEC_VV(1 + it), VEC_TEMP_MATOP));
157: }
158: dgmres->matvecs += 1;
159: /* update Hessenberg matrix and do Gram-Schmidt */
160: PetscCall((*dgmres->orthog)(ksp, it));
162: /* vv(i+1) . vv(i+1) */
163: PetscCall(VecNormalize(VEC_VV(it + 1), &tt));
164: /* save the magnitude */
165: *HH(it + 1, it) = tt;
166: *HES(it + 1, it) = tt;
168: /* check for the happy breakdown */
169: hapbnd = PetscAbsScalar(tt / *GRS(it));
170: if (hapbnd > dgmres->haptol) hapbnd = dgmres->haptol;
171: if (tt < hapbnd) {
172: PetscCall(PetscInfo(ksp, "Detected happy breakdown, current hapbnd = %g tt = %g\n", (double)hapbnd, (double)tt));
173: hapend = PETSC_TRUE;
174: }
175: PetscCall(KSPDGMRESUpdateHessenberg(ksp, it, hapend, &res));
177: it++;
178: dgmres->it = (it - 1); /* For converged */
179: ksp->its++;
180: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res;
181: else ksp->rnorm = 0.0;
182: if (ksp->reason) break;
184: PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP));
186: /* Catch error in happy breakdown and signal convergence and break from loop */
187: if (hapend) {
188: if (!ksp->reason) {
189: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Reached happy break down, but convergence was not indicated. Residual norm = %g", (double)res);
190: ksp->reason = KSP_DIVERGED_BREAKDOWN;
191: break;
192: }
193: }
194: }
196: if (itcount) *itcount = it;
198: /*
199: Down here we have to solve for the "best" coefficients of the Krylov
200: columns, add the solution values together, and possibly unwind the
201: preconditioning from the solution
202: */
203: /* Form the solution (or the solution so far) */
204: PetscCall(KSPDGMRESBuildSoln(GRS(0), ksp->vec_sol, ksp->vec_sol, ksp, it - 1));
206: /* Monitor if we know that we will not return for a restart */
207: if (ksp->reason == KSP_CONVERGED_ITERATING && ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
208: if (it && ksp->reason) {
209: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
210: PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
211: }
213: /* Compute data for the deflation to be used during the next restart */
214: if (!ksp->reason && ksp->its < ksp->max_it) {
215: test = max_k * PetscLogReal(ksp->rtol / res) / PetscLogReal(res / res_old);
216: /* Compute data for the deflation if the residual rtol will not be reached in the remaining number of steps allowed */
217: if ((test > dgmres->smv * (ksp->max_it - ksp->its)) || dgmres->force) PetscCall(KSPDGMRESComputeDeflationData(ksp, NULL));
218: }
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: PetscErrorCode KSPSolve_DGMRES(KSP ksp)
223: {
224: PetscInt i, its = 0, itcount;
225: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
226: PetscBool guess_zero = ksp->guess_zero;
228: PetscFunctionBegin;
229: PetscCheck(!ksp->calc_sings || dgmres->Rsvd, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ORDER, "Must call KSPSetComputeSingularValues() before KSPSetUp() is called");
231: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
232: ksp->its = 0;
233: dgmres->matvecs = 0;
234: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
236: itcount = 0;
237: while (!ksp->reason) {
238: PetscCall(KSPInitialResidual(ksp, ksp->vec_sol, VEC_TEMP, VEC_TEMP_MATOP, VEC_VV(0), ksp->vec_rhs));
239: if (ksp->pc_side == PC_LEFT) {
240: dgmres->matvecs += 1;
241: if (dgmres->r > 0) {
242: PetscCall(KSPDGMRESApplyDeflation(ksp, VEC_VV(0), VEC_TEMP));
243: PetscCall(VecCopy(VEC_TEMP, VEC_VV(0)));
244: }
245: }
247: PetscCall(KSPDGMRESCycle(&its, ksp));
248: itcount += its;
249: if (itcount >= ksp->max_it) {
250: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
251: break;
252: }
253: ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
254: }
255: ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */
257: for (i = 0; i < dgmres->r; i++) PetscCall(VecViewFromOptions(UU[i], (PetscObject)ksp, "-ksp_dgmres_view_deflation_vecs"));
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: PetscErrorCode KSPDestroy_DGMRES(KSP ksp)
262: {
263: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
264: PetscInt neig1 = dgmres->neig + EIG_OFFSET;
265: PetscInt max_neig = dgmres->max_neig;
267: PetscFunctionBegin;
268: if (dgmres->r) {
269: PetscCall(VecDestroyVecs(max_neig, &UU));
270: PetscCall(VecDestroyVecs(max_neig, &MU));
271: if (XX) {
272: PetscCall(VecDestroyVecs(neig1, &XX));
273: PetscCall(VecDestroyVecs(neig1, &MX));
274: }
275: PetscCall(PetscFree(TT));
276: PetscCall(PetscFree(TTF));
277: PetscCall(PetscFree(INVP));
278: PetscCall(PetscFree(XMX));
279: PetscCall(PetscFree(UMX));
280: PetscCall(PetscFree(XMU));
281: PetscCall(PetscFree(X1));
282: PetscCall(PetscFree(X2));
283: PetscCall(PetscFree(dgmres->work));
284: PetscCall(PetscFree(dgmres->iwork));
285: PetscCall(PetscFree(dgmres->wr));
286: PetscCall(PetscFree(dgmres->wi));
287: PetscCall(PetscFree(dgmres->modul));
288: PetscCall(PetscFree(dgmres->Q));
289: PetscCall(PetscFree(ORTH));
290: PetscCall(PetscFree(AUAU));
291: PetscCall(PetscFree(AUU));
292: PetscCall(PetscFree(SR2));
293: }
294: PetscCall(PetscFree(SR));
295: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESSetEigen_C", NULL));
296: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESSetMaxEigen_C", NULL));
297: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESSetRatio_C", NULL));
298: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESForce_C", NULL));
299: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESComputeSchurForm_C", NULL));
300: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESComputeDeflationData_C", NULL));
301: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESApplyDeflation_C", NULL));
302: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESImproveEig_C", NULL));
303: PetscCall(KSPDestroy_GMRES(ksp));
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }
307: /*
308: KSPDGMRESBuildSoln - create the solution from the starting vector and the
309: current iterates.
311: Input parameters:
312: nrs - work area of size it + 1.
313: vs - index of initial guess
314: vdest - index of result. Note that vs may == vdest (replace
315: guess with the solution).
317: This is an internal routine that knows about the GMRES internals.
318: */
319: static PetscErrorCode KSPDGMRESBuildSoln(PetscScalar *nrs, Vec vs, Vec vdest, KSP ksp, PetscInt it)
320: {
321: PetscScalar tt;
322: PetscInt ii, k, j;
323: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
325: /* Solve for solution vector that minimizes the residual */
327: PetscFunctionBegin;
328: /* If it is < 0, no gmres steps have been performed */
329: if (it < 0) {
330: PetscCall(VecCopy(vs, vdest)); /* VecCopy() is smart, exists immediately if vguess == vdest */
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
333: PetscCheck(*HH(it, it) != 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Likely your matrix is the zero operator. HH(it,it) is identically zero; it = %" PetscInt_FMT " GRS(it) = %g", it, (double)PetscAbsScalar(*GRS(it)));
334: if (*HH(it, it) != 0.0) nrs[it] = *GRS(it) / *HH(it, it);
335: else nrs[it] = 0.0;
337: for (ii = 1; ii <= it; ii++) {
338: k = it - ii;
339: tt = *GRS(k);
340: for (j = k + 1; j <= it; j++) tt = tt - *HH(k, j) * nrs[j];
341: PetscCheck(*HH(k, k) != 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Likely your matrix is singular. HH(k,k) is identically zero; it = %" PetscInt_FMT " k = %" PetscInt_FMT, it, k);
342: nrs[k] = tt / *HH(k, k);
343: }
345: /* Accumulate the correction to the solution of the preconditioned problem in TEMP */
346: PetscCall(VecMAXPBY(VEC_TEMP, it + 1, nrs, 0, &VEC_VV(0)));
348: /* Apply deflation */
349: if (ksp->pc_side == PC_RIGHT && dgmres->r > 0) {
350: PetscCall(KSPDGMRESApplyDeflation(ksp, VEC_TEMP, VEC_TEMP_MATOP));
351: PetscCall(VecCopy(VEC_TEMP_MATOP, VEC_TEMP));
352: }
353: PetscCall(KSPUnwindPreconditioner(ksp, VEC_TEMP, VEC_TEMP_MATOP));
355: /* add solution to previous solution */
356: if (vdest != vs) PetscCall(VecCopy(vs, vdest));
357: PetscCall(VecAXPY(vdest, 1.0, VEC_TEMP));
358: PetscFunctionReturn(PETSC_SUCCESS);
359: }
361: /*
362: Do the scalar work for the orthogonalization. Return new residual norm.
363: */
364: static PetscErrorCode KSPDGMRESUpdateHessenberg(KSP ksp, PetscInt it, PetscBool hapend, PetscReal *res)
365: {
366: PetscScalar *hh, *cc, *ss, tt;
367: PetscInt j;
368: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
370: PetscFunctionBegin;
371: hh = HH(0, it);
372: cc = CC(0);
373: ss = SS(0);
375: /* Apply all the previously computed plane rotations to the new column
376: of the Hessenberg matrix */
377: for (j = 1; j <= it; j++) {
378: tt = *hh;
379: *hh = PetscConj(*cc) * tt + *ss * *(hh + 1);
380: hh++;
381: *hh = *cc++ * *hh - (*ss++ * tt);
382: }
384: /*
385: compute the new plane rotation, and apply it to:
386: 1) the right-hand side of the Hessenberg system
387: 2) the new column of the Hessenberg matrix
388: thus obtaining the updated value of the residual
389: */
390: if (!hapend) {
391: tt = PetscSqrtScalar(PetscConj(*hh) * *hh + PetscConj(*(hh + 1)) * *(hh + 1));
392: if (tt == 0.0) {
393: ksp->reason = KSP_DIVERGED_NULL;
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
396: *cc = *hh / tt;
397: *ss = *(hh + 1) / tt;
398: *GRS(it + 1) = -(*ss * *GRS(it));
399: *GRS(it) = PetscConj(*cc) * *GRS(it);
400: *hh = PetscConj(*cc) * *hh + *ss * *(hh + 1);
401: *res = PetscAbsScalar(*GRS(it + 1));
402: } else {
403: /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
404: another rotation matrix (so RH doesn't change). The new residual is
405: always the new sine term times the residual from last time (GRS(it)),
406: but now the new sine rotation would be zero...so the residual should
407: be zero...so we will multiply "zero" by the last residual. This might
408: not be exactly what we want to do here -could just return "zero". */
409: *res = 0.0;
410: }
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: /*
415: Allocates more work vectors, starting from VEC_VV(it).
416: */
417: static PetscErrorCode KSPDGMRESGetNewVectors(KSP ksp, PetscInt it)
418: {
419: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
420: PetscInt nwork = dgmres->nwork_alloc, k, nalloc;
422: PetscFunctionBegin;
423: nalloc = PetscMin(ksp->max_it, dgmres->delta_allocate);
424: /* Adjust the number to allocate to make sure that we don't exceed the
425: number of available slots */
426: if (it + VEC_OFFSET + nalloc >= dgmres->vecs_allocated) nalloc = dgmres->vecs_allocated - it - VEC_OFFSET;
427: if (!nalloc) PetscFunctionReturn(PETSC_SUCCESS);
429: dgmres->vv_allocated += nalloc;
431: PetscCall(KSPCreateVecs(ksp, nalloc, &dgmres->user_work[nwork], 0, NULL));
433: dgmres->mwork_alloc[nwork] = nalloc;
434: for (k = 0; k < nalloc; k++) dgmres->vecs[it + VEC_OFFSET + k] = dgmres->user_work[nwork][k];
435: dgmres->nwork_alloc++;
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }
439: PetscErrorCode KSPBuildSolution_DGMRES(KSP ksp, Vec ptr, Vec *result)
440: {
441: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
443: PetscFunctionBegin;
444: if (!ptr) {
445: if (!dgmres->sol_temp) PetscCall(VecDuplicate(ksp->vec_sol, &dgmres->sol_temp));
446: ptr = dgmres->sol_temp;
447: }
448: if (!dgmres->nrs) {
449: /* allocate the work area */
450: PetscCall(PetscMalloc1(dgmres->max_k, &dgmres->nrs));
451: }
452: PetscCall(KSPDGMRESBuildSoln(dgmres->nrs, ksp->vec_sol, ptr, ksp, dgmres->it));
453: if (result) *result = ptr;
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: static PetscErrorCode KSPView_DGMRES(KSP ksp, PetscViewer viewer)
458: {
459: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
460: PetscBool iascii, isharmonic;
462: PetscFunctionBegin;
463: PetscCall(KSPView_GMRES(ksp, viewer));
464: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
465: if (iascii) {
466: PetscCall(PetscViewerASCIIPrintf(viewer, " Adaptive strategy is used: %s\n", PetscBools[dgmres->force]));
467: PetscCall(PetscOptionsHasName(((PetscObject)ksp)->options, ((PetscObject)ksp)->prefix, "-ksp_dgmres_harmonic_ritz", &isharmonic));
468: if (isharmonic) {
469: PetscCall(PetscViewerASCIIPrintf(viewer, " Frequency of extracted eigenvalues = %" PetscInt_FMT " using Harmonic Ritz values \n", dgmres->neig));
470: } else {
471: PetscCall(PetscViewerASCIIPrintf(viewer, " Frequency of extracted eigenvalues = %" PetscInt_FMT " using Ritz values \n", dgmres->neig));
472: }
473: PetscCall(PetscViewerASCIIPrintf(viewer, " Total number of extracted eigenvalues = %" PetscInt_FMT "\n", dgmres->r));
474: PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum number of eigenvalues set to be extracted = %" PetscInt_FMT "\n", dgmres->max_neig));
475: PetscCall(PetscViewerASCIIPrintf(viewer, " relaxation parameter for the adaptive strategy(smv) = %g\n", (double)dgmres->smv));
476: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of matvecs : %" PetscInt_FMT "\n", dgmres->matvecs));
477: }
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: PetscErrorCode KSPDGMRESSetEigen_DGMRES(KSP ksp, PetscInt neig)
482: {
483: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
485: PetscFunctionBegin;
486: PetscCheck(neig >= 0 && neig <= dgmres->max_k, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "The value of neig must be positive and less than the restart value ");
487: dgmres->neig = neig;
488: PetscFunctionReturn(PETSC_SUCCESS);
489: }
491: static PetscErrorCode KSPDGMRESSetMaxEigen_DGMRES(KSP ksp, PetscInt max_neig)
492: {
493: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
495: PetscFunctionBegin;
496: PetscCheck(max_neig >= 0 && max_neig <= dgmres->max_k, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "The value of max_neig must be positive and less than the restart value ");
497: dgmres->max_neig = max_neig;
498: PetscFunctionReturn(PETSC_SUCCESS);
499: }
501: static PetscErrorCode KSPDGMRESSetRatio_DGMRES(KSP ksp, PetscReal ratio)
502: {
503: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
505: PetscFunctionBegin;
506: PetscCheck(ratio > 0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "The relaxation parameter value must be positive");
507: dgmres->smv = ratio;
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: static PetscErrorCode KSPDGMRESForce_DGMRES(KSP ksp, PetscBool force)
512: {
513: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
515: PetscFunctionBegin;
516: dgmres->force = force;
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: PetscErrorCode KSPSetFromOptions_DGMRES(KSP ksp, PetscOptionItems *PetscOptionsObject)
521: {
522: PetscInt neig;
523: PetscInt max_neig;
524: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
525: PetscBool flg;
527: PetscFunctionBegin;
528: PetscCall(KSPSetFromOptions_GMRES(ksp, PetscOptionsObject));
529: PetscOptionsHeadBegin(PetscOptionsObject, "KSP DGMRES Options");
530: PetscCall(PetscOptionsInt("-ksp_dgmres_eigen", "Number of smallest eigenvalues to extract at each restart", "KSPDGMRESSetEigen", dgmres->neig, &neig, &flg));
531: if (flg) PetscCall(KSPDGMRESSetEigen(ksp, neig));
532: PetscCall(PetscOptionsInt("-ksp_dgmres_max_eigen", "Maximum Number of smallest eigenvalues to extract ", "KSPDGMRESSetMaxEigen", dgmres->max_neig, &max_neig, &flg));
533: if (flg) PetscCall(KSPDGMRESSetMaxEigen(ksp, max_neig));
534: PetscCall(PetscOptionsReal("-ksp_dgmres_ratio", "Relaxation parameter for the smaller number of matrix-vectors product allowed", "KSPDGMRESSetRatio", dgmres->smv, &dgmres->smv, NULL));
535: PetscCall(PetscOptionsBool("-ksp_dgmres_improve", "Improve the computation of eigenvalues by solving a new generalized eigenvalue problem (experimental - not stable at this time)", NULL, dgmres->improve, &dgmres->improve, NULL));
536: PetscCall(PetscOptionsBool("-ksp_dgmres_force", "Sets DGMRES always at restart active, i.e do not use the adaptive strategy", "KSPDGMRESForce", dgmres->force, &dgmres->force, NULL));
537: PetscOptionsHeadEnd();
538: PetscFunctionReturn(PETSC_SUCCESS);
539: }
541: PetscErrorCode KSPDGMRESComputeDeflationData_DGMRES(KSP ksp, PetscInt *ExtrNeig)
542: {
543: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
544: PetscInt i, j, k;
545: PetscBLASInt nr, bmax;
546: PetscInt r = dgmres->r;
547: PetscInt neig; /* number of eigenvalues to extract at each restart */
548: PetscInt neig1 = dgmres->neig + EIG_OFFSET; /* max number of eig that can be extracted at each restart */
549: PetscInt max_neig = dgmres->max_neig; /* Max number of eigenvalues to extract during the iterative process */
550: PetscInt N = dgmres->max_k + 1;
551: PetscInt n = dgmres->it + 1;
552: PetscReal alpha;
554: PetscFunctionBegin;
555: PetscCall(PetscLogEventBegin(KSP_DGMRESComputeDeflationData, ksp, 0, 0, 0));
556: if (dgmres->neig == 0 || (max_neig < (r + neig1) && !dgmres->improve)) {
557: PetscCall(PetscLogEventEnd(KSP_DGMRESComputeDeflationData, ksp, 0, 0, 0));
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: PetscCall(KSPDGMRESComputeSchurForm(ksp, &neig));
562: /* Form the extended Schur vectors X=VV*Sr */
563: if (!XX) PetscCall(VecDuplicateVecs(VEC_VV(0), neig1, &XX));
564: for (j = 0; j < neig; j++) PetscCall(VecMAXPBY(XX[j], n, &SR[j * N], 0, &VEC_VV(0)));
566: /* Orthogonalize X against U */
567: if (!ORTH) PetscCall(PetscMalloc1(max_neig, &ORTH));
568: if (r > 0) {
569: /* modified Gram-Schmidt */
570: for (j = 0; j < neig; j++) {
571: for (i = 0; i < r; i++) {
572: /* First, compute U'*X[j] */
573: PetscCall(VecDot(XX[j], UU[i], &alpha));
574: /* Then, compute X(j)=X(j)-U*U'*X(j) */
575: PetscCall(VecAXPY(XX[j], -alpha, UU[i]));
576: }
577: }
578: }
579: /* Compute MX = M^{-1}*A*X */
580: if (!MX) PetscCall(VecDuplicateVecs(VEC_VV(0), neig1, &MX));
581: for (j = 0; j < neig; j++) PetscCall(KSP_PCApplyBAorAB(ksp, XX[j], MX[j], VEC_TEMP_MATOP));
582: dgmres->matvecs += neig;
584: if ((r + neig1) > max_neig && dgmres->improve) { /* Improve the approximate eigenvectors in X by solving a new generalized eigenvalue -- expensive to do this */
585: PetscCall(KSPDGMRESImproveEig(ksp, neig));
586: PetscCall(PetscLogEventEnd(KSP_DGMRESComputeDeflationData, ksp, 0, 0, 0));
587: PetscFunctionReturn(PETSC_SUCCESS); /* We return here since data for M have been improved in KSPDGMRESImproveEig()*/
588: }
590: /* Compute XMX = X'*M^{-1}*A*X -- size (neig, neig) */
591: if (!XMX) PetscCall(PetscMalloc1(neig1 * neig1, &XMX));
592: for (j = 0; j < neig; j++) PetscCall(VecMDot(MX[j], neig, XX, &XMX[j * neig1]));
594: if (r > 0) {
595: /* Compute UMX = U'*M^{-1}*A*X -- size (r, neig) */
596: if (!UMX) PetscCall(PetscMalloc1(max_neig * neig1, &UMX));
597: for (j = 0; j < neig; j++) PetscCall(VecMDot(MX[j], r, UU, &UMX[j * max_neig]));
598: /* Compute XMU = X'*M^{-1}*A*U -- size(neig, r) */
599: if (!XMU) PetscCall(PetscMalloc1(max_neig * neig1, &XMU));
600: for (j = 0; j < r; j++) PetscCall(VecMDot(MU[j], neig, XX, &XMU[j * neig1]));
601: }
603: /* Form the new matrix T = [T UMX; XMU XMX]; */
604: if (!TT) PetscCall(PetscMalloc1(max_neig * max_neig, &TT));
605: if (r > 0) {
606: /* Add XMU to T */
607: for (j = 0; j < r; j++) PetscCall(PetscArraycpy(&TT[max_neig * j + r], &XMU[neig1 * j], neig));
608: /* Add [UMX; XMX] to T */
609: for (j = 0; j < neig; j++) {
610: k = r + j;
611: PetscCall(PetscArraycpy(&TT[max_neig * k], &UMX[max_neig * j], r));
612: PetscCall(PetscArraycpy(&TT[max_neig * k + r], &XMX[neig1 * j], neig));
613: }
614: } else { /* Add XMX to T */
615: for (j = 0; j < neig; j++) PetscCall(PetscArraycpy(&TT[max_neig * j], &XMX[neig1 * j], neig));
616: }
618: dgmres->r += neig;
619: r = dgmres->r;
620: PetscCall(PetscBLASIntCast(r, &nr));
621: /*LU Factorize T with Lapack xgetrf routine */
623: PetscCall(PetscBLASIntCast(max_neig, &bmax));
624: if (!TTF) PetscCall(PetscMalloc1(bmax * bmax, &TTF));
625: PetscCall(PetscArraycpy(TTF, TT, bmax * r));
626: if (!INVP) PetscCall(PetscMalloc1(bmax, &INVP));
627: {
628: PetscBLASInt info;
629: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&nr, &nr, TTF, &bmax, INVP, &info));
630: PetscCheck(!info, PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB, "Error in LAPACK routine XGETRF INFO=%d", (int)info);
631: }
633: /* Save X in U and MX in MU for the next cycles and increase the size of the invariant subspace */
634: if (!UU) {
635: PetscCall(VecDuplicateVecs(VEC_VV(0), max_neig, &UU));
636: PetscCall(VecDuplicateVecs(VEC_VV(0), max_neig, &MU));
637: }
638: for (j = 0; j < neig; j++) {
639: PetscCall(VecCopy(XX[j], UU[r - neig + j]));
640: PetscCall(VecCopy(MX[j], MU[r - neig + j]));
641: }
642: PetscCall(PetscLogEventEnd(KSP_DGMRESComputeDeflationData, ksp, 0, 0, 0));
643: PetscFunctionReturn(PETSC_SUCCESS);
644: }
646: PetscErrorCode KSPDGMRESComputeSchurForm_DGMRES(KSP ksp, PetscInt *neig)
647: {
648: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
649: PetscInt N = dgmres->max_k + 1, n = dgmres->it + 1;
650: PetscBLASInt bn;
651: PetscReal *A;
652: PetscBLASInt ihi;
653: PetscBLASInt ldA = 0; /* leading dimension of A */
654: PetscBLASInt ldQ; /* leading dimension of Q */
655: PetscReal *Q; /* orthogonal matrix of (left) Schur vectors */
656: PetscReal *work; /* working vector */
657: PetscBLASInt lwork; /* size of the working vector */
658: PetscInt *perm; /* Permutation vector to sort eigenvalues */
659: PetscInt i, j;
660: PetscBLASInt NbrEig; /* Number of eigenvalues really extracted */
661: PetscReal *wr, *wi, *modul; /* Real and imaginary part and modulus of the eigenvalues of A */
662: PetscBLASInt *select;
663: PetscBLASInt *iwork;
664: PetscBLASInt liwork;
665: PetscScalar *Ht; /* Transpose of the Hessenberg matrix */
666: PetscScalar *t; /* Store the result of the solution of H^T*t=h_{m+1,m}e_m */
667: PetscBLASInt *ipiv; /* Permutation vector to be used in LAPACK */
668: PetscBool flag; /* determine whether to use Ritz vectors or harmonic Ritz vectors */
670: PetscFunctionBegin;
671: PetscCall(PetscBLASIntCast(n, &bn));
672: PetscCall(PetscBLASIntCast(N, &ldA));
673: ihi = ldQ = bn;
674: PetscCall(PetscBLASIntCast(5 * N, &lwork));
676: #if defined(PETSC_USE_COMPLEX)
677: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "No support for complex numbers.");
678: #endif
680: PetscCall(PetscMalloc1(ldA * ldA, &A));
681: PetscCall(PetscMalloc1(ldQ * n, &Q));
682: PetscCall(PetscMalloc1(lwork, &work));
683: if (!dgmres->wr) {
684: PetscCall(PetscMalloc1(n, &dgmres->wr));
685: PetscCall(PetscMalloc1(n, &dgmres->wi));
686: }
687: wr = dgmres->wr;
688: wi = dgmres->wi;
689: PetscCall(PetscMalloc1(n, &modul));
690: PetscCall(PetscMalloc1(n, &perm));
691: /* copy the Hessenberg matrix to work space */
692: PetscCall(PetscArraycpy(A, dgmres->hes_origin, ldA * ldA));
693: PetscCall(PetscOptionsHasName(((PetscObject)ksp)->options, ((PetscObject)ksp)->prefix, "-ksp_dgmres_harmonic_ritz", &flag));
694: if (flag) {
695: /* Compute the matrix H + H^{-T}*h^2_{m+1,m}e_m*e_m^T */
696: /* Transpose the Hessenberg matrix */
697: PetscCall(PetscMalloc1(bn * bn, &Ht));
698: for (i = 0; i < bn; i++) {
699: for (j = 0; j < bn; j++) Ht[i * bn + j] = dgmres->hes_origin[j * ldA + i];
700: }
702: /* Solve the system H^T*t = h_{m+1,m}e_m */
703: PetscCall(PetscCalloc1(bn, &t));
704: t[bn - 1] = dgmres->hes_origin[(bn - 1) * ldA + bn]; /* Pick the last element H(m+1,m) */
705: PetscCall(PetscMalloc1(bn, &ipiv));
706: /* Call the LAPACK routine dgesv to solve the system Ht^-1 * t */
707: {
708: PetscBLASInt info;
709: PetscBLASInt nrhs = 1;
710: PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&bn, &nrhs, Ht, &bn, ipiv, t, &bn, &info));
711: PetscCheck(!info, PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB, "Error while calling the Lapack routine DGESV");
712: }
713: /* Now form H + H^{-T}*h^2_{m+1,m}e_m*e_m^T */
714: for (i = 0; i < bn; i++) A[(bn - 1) * bn + i] += t[i];
715: PetscCall(PetscFree(t));
716: PetscCall(PetscFree(Ht));
717: }
718: /* Compute eigenvalues with the Schur form */
719: {
720: PetscBLASInt info = 0;
721: PetscBLASInt ilo = 1;
722: PetscCallBLAS("LAPACKhseqr", LAPACKhseqr_("S", "I", &bn, &ilo, &ihi, A, &ldA, wr, wi, Q, &ldQ, work, &lwork, &info));
723: PetscCheck(!info, PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB, "Error in LAPACK routine XHSEQR %d", (int)info);
724: }
725: PetscCall(PetscFree(work));
727: /* sort the eigenvalues */
728: for (i = 0; i < n; i++) modul[i] = PetscSqrtReal(wr[i] * wr[i] + wi[i] * wi[i]);
729: for (i = 0; i < n; i++) perm[i] = i;
731: PetscCall(PetscSortRealWithPermutation(n, modul, perm));
732: /* save the complex modulus of the largest eigenvalue in magnitude */
733: if (dgmres->lambdaN < modul[perm[n - 1]]) dgmres->lambdaN = modul[perm[n - 1]];
734: /* count the number of extracted eigenvalues (with complex conjugates) */
735: NbrEig = 0;
736: while (NbrEig < dgmres->neig) {
737: if (wi[perm[NbrEig]] != 0) NbrEig += 2;
738: else NbrEig += 1;
739: }
740: /* Reorder the Schur decomposition so that the cluster of smallest eigenvalues appears in the leading diagonal blocks of A */
742: PetscCall(PetscCalloc1(n, &select));
744: if (!dgmres->GreatestEig) {
745: for (j = 0; j < NbrEig; j++) select[perm[j]] = 1;
746: } else {
747: for (j = 0; j < NbrEig; j++) select[perm[n - j - 1]] = 1;
748: }
749: /* call Lapack dtrsen */
750: lwork = PetscMax(1, 4 * NbrEig * (bn - NbrEig));
751: liwork = PetscMax(1, 2 * NbrEig * (bn - NbrEig));
752: PetscCall(PetscMalloc1(lwork, &work));
753: PetscCall(PetscMalloc1(liwork, &iwork));
754: {
755: PetscBLASInt info = 0;
756: PetscReal CondEig; /* lower bound on the reciprocal condition number for the selected cluster of eigenvalues */
757: PetscReal CondSub; /* estimated reciprocal condition number of the specified invariant subspace. */
758: PetscCallBLAS("LAPACKtrsen", LAPACKtrsen_("B", "V", select, &bn, A, &ldA, Q, &ldQ, wr, wi, &NbrEig, &CondEig, &CondSub, work, &lwork, iwork, &liwork, &info));
759: PetscCheck(info != 1, PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB, "Unable to reorder the eigenvalues with the LAPACK routine: ILL-CONDITIONED PROBLEM");
760: }
761: PetscCall(PetscFree(select));
763: /* Extract the Schur vectors */
764: for (j = 0; j < NbrEig; j++) PetscCall(PetscArraycpy(&SR[j * N], &Q[j * ldQ], n));
765: *neig = NbrEig;
766: PetscCall(PetscFree(A));
767: PetscCall(PetscFree(work));
768: PetscCall(PetscFree(perm));
769: PetscCall(PetscFree(work));
770: PetscCall(PetscFree(iwork));
771: PetscCall(PetscFree(modul));
772: PetscCall(PetscFree(Q));
773: PetscFunctionReturn(PETSC_SUCCESS);
774: }
776: PetscErrorCode KSPDGMRESApplyDeflation_DGMRES(KSP ksp, Vec x, Vec y)
777: {
778: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
779: PetscInt i, r = dgmres->r;
780: PetscReal alpha = 1.0;
781: PetscInt max_neig = dgmres->max_neig;
782: PetscBLASInt br, bmax;
783: PetscReal lambda = dgmres->lambdaN;
785: PetscFunctionBegin;
786: PetscCall(PetscBLASIntCast(r, &br));
787: PetscCall(PetscBLASIntCast(max_neig, &bmax));
788: PetscCall(PetscLogEventBegin(KSP_DGMRESApplyDeflation, ksp, 0, 0, 0));
789: if (!r) {
790: PetscCall(VecCopy(x, y));
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
793: /* Compute U'*x */
794: if (!X1) {
795: PetscCall(PetscMalloc1(bmax, &X1));
796: PetscCall(PetscMalloc1(bmax, &X2));
797: }
798: PetscCall(VecMDot(x, r, UU, X1));
800: /* Solve T*X1=X2 for X1*/
801: PetscCall(PetscArraycpy(X2, X1, br));
802: {
803: PetscBLASInt info;
804: PetscBLASInt nrhs = 1;
805: PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_("N", &br, &nrhs, TTF, &bmax, INVP, X1, &bmax, &info));
806: PetscCheck(!info, PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB, "Error in LAPACK routine XGETRS %d", (int)info);
807: }
808: /* Iterative refinement -- is it really necessary ?? */
809: if (!WORK) {
810: PetscCall(PetscMalloc1(3 * bmax, &WORK));
811: PetscCall(PetscMalloc1(bmax, &IWORK));
812: }
813: {
814: PetscBLASInt info;
815: PetscReal berr, ferr;
816: PetscBLASInt nrhs = 1;
817: PetscCallBLAS("LAPACKgerfs", LAPACKgerfs_("N", &br, &nrhs, TT, &bmax, TTF, &bmax, INVP, X2, &bmax, X1, &bmax, &ferr, &berr, WORK, IWORK, &info));
818: PetscCheck(!info, PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB, "Error in LAPACK routine XGERFS %d", (int)info);
819: }
821: for (i = 0; i < r; i++) X2[i] = X1[i] / lambda - X2[i];
823: /* Compute X2=U*X2 */
824: PetscCall(VecMAXPBY(y, r, X2, 0, UU));
825: PetscCall(VecAXPY(y, alpha, x));
827: PetscCall(PetscLogEventEnd(KSP_DGMRESApplyDeflation, ksp, 0, 0, 0));
828: PetscFunctionReturn(PETSC_SUCCESS);
829: }
831: static PetscErrorCode KSPDGMRESImproveEig_DGMRES(KSP ksp, PetscInt neig)
832: {
833: KSP_DGMRES *dgmres = (KSP_DGMRES *)ksp->data;
834: PetscInt j, r_old, r = dgmres->r;
835: PetscBLASInt i = 0;
836: PetscInt neig1 = dgmres->neig + EIG_OFFSET;
837: PetscInt bmax = dgmres->max_neig;
838: PetscInt aug = r + neig; /* actual size of the augmented invariant basis */
839: PetscInt aug1 = bmax + neig1; /* maximum size of the augmented invariant basis */
840: PetscBLASInt ldA; /* leading dimension of AUAU and AUU*/
841: PetscBLASInt N; /* size of AUAU */
842: PetscReal *Q; /* orthogonal matrix of (left) schur vectors */
843: PetscReal *Z; /* orthogonal matrix of (right) schur vectors */
844: PetscReal *work; /* working vector */
845: PetscBLASInt lwork; /* size of the working vector */
846: PetscInt *perm; /* Permutation vector to sort eigenvalues */
847: PetscReal *wr, *wi, *beta, *modul; /* Real and imaginary part and modulus of the eigenvalues of A*/
848: PetscBLASInt NbrEig = 0, nr, bm;
849: PetscBLASInt *select;
850: PetscBLASInt liwork, *iwork;
852: PetscFunctionBegin;
853: /* Block construction of the matrices AUU=(AU)'*U and (AU)'*AU*/
854: if (!AUU) {
855: PetscCall(PetscMalloc1(aug1 * aug1, &AUU));
856: PetscCall(PetscMalloc1(aug1 * aug1, &AUAU));
857: }
858: /* AUU = (AU)'*U = [(MU)'*U (MU)'*X; (MX)'*U (MX)'*X]
859: * Note that MU and MX have been computed previously either in ComputeDataDeflation() or down here in a previous call to this function */
860: /* (MU)'*U size (r x r) -- store in the <r> first columns of AUU*/
861: for (j = 0; j < r; j++) PetscCall(VecMDot(UU[j], r, MU, &AUU[j * aug1]));
862: /* (MU)'*X size (r x neig) -- store in AUU from the column <r>*/
863: for (j = 0; j < neig; j++) PetscCall(VecMDot(XX[j], r, MU, &AUU[(r + j) * aug1]));
864: /* (MX)'*U size (neig x r) -- store in the <r> first columns of AUU from the row <r>*/
865: for (j = 0; j < r; j++) PetscCall(VecMDot(UU[j], neig, MX, &AUU[j * aug1 + r]));
866: /* (MX)'*X size (neig neig) -- store in AUU from the column <r> and the row <r>*/
867: for (j = 0; j < neig; j++) PetscCall(VecMDot(XX[j], neig, MX, &AUU[(r + j) * aug1 + r]));
869: /* AUAU = (AU)'*AU = [(MU)'*MU (MU)'*MX; (MX)'*MU (MX)'*MX] */
870: /* (MU)'*MU size (r x r) -- store in the <r> first columns of AUAU*/
871: for (j = 0; j < r; j++) PetscCall(VecMDot(MU[j], r, MU, &AUAU[j * aug1]));
872: /* (MU)'*MX size (r x neig) -- store in AUAU from the column <r>*/
873: for (j = 0; j < neig; j++) PetscCall(VecMDot(MX[j], r, MU, &AUAU[(r + j) * aug1]));
874: /* (MX)'*MU size (neig x r) -- store in the <r> first columns of AUAU from the row <r>*/
875: for (j = 0; j < r; j++) PetscCall(VecMDot(MU[j], neig, MX, &AUAU[j * aug1 + r]));
876: /* (MX)'*MX size (neig neig) -- store in AUAU from the column <r> and the row <r>*/
877: for (j = 0; j < neig; j++) PetscCall(VecMDot(MX[j], neig, MX, &AUAU[(r + j) * aug1 + r]));
879: /* Computation of the eigenvectors */
880: PetscCall(PetscBLASIntCast(aug1, &ldA));
881: PetscCall(PetscBLASIntCast(aug, &N));
882: lwork = 8 * N + 20; /* sizeof the working space */
883: PetscCall(PetscMalloc1(N, &wr));
884: PetscCall(PetscMalloc1(N, &wi));
885: PetscCall(PetscMalloc1(N, &beta));
886: PetscCall(PetscMalloc1(N, &modul));
887: PetscCall(PetscMalloc1(N, &perm));
888: PetscCall(PetscMalloc1(N * N, &Q));
889: PetscCall(PetscMalloc1(N * N, &Z));
890: PetscCall(PetscMalloc1(lwork, &work));
891: {
892: PetscBLASInt info = 0;
893: PetscCallBLAS("LAPACKgges", LAPACKgges_("V", "V", "N", NULL, &N, AUAU, &ldA, AUU, &ldA, &i, wr, wi, beta, Q, &N, Z, &N, work, &lwork, NULL, &info));
894: PetscCheck(!info, PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB, "Error in LAPACK routine XGGES %d", (int)info);
895: }
896: for (i = 0; i < N; i++) {
897: if (beta[i] != 0.0) {
898: wr[i] /= beta[i];
899: wi[i] /= beta[i];
900: }
901: }
902: /* sort the eigenvalues */
903: for (i = 0; i < N; i++) modul[i] = PetscSqrtReal(wr[i] * wr[i] + wi[i] * wi[i]);
904: for (i = 0; i < N; i++) perm[i] = i;
905: PetscCall(PetscSortRealWithPermutation(N, modul, perm));
906: /* Save the norm of the largest eigenvalue */
907: if (dgmres->lambdaN < modul[perm[N - 1]]) dgmres->lambdaN = modul[perm[N - 1]];
908: /* Allocate space to extract the first r schur vectors */
909: if (!SR2) PetscCall(PetscMalloc1(aug1 * bmax, &SR2));
910: /* count the number of extracted eigenvalues (complex conjugates count as 2) */
911: while (NbrEig < bmax) {
912: if (wi[perm[NbrEig]] == 0) NbrEig += 1;
913: else NbrEig += 2;
914: }
915: if (NbrEig > bmax) NbrEig = bmax - 1;
916: r_old = r; /* previous size of r */
917: dgmres->r = r = NbrEig;
919: /* Select the eigenvalues to reorder */
920: PetscCall(PetscCalloc1(N, &select));
921: if (!dgmres->GreatestEig) {
922: for (j = 0; j < NbrEig; j++) select[perm[j]] = 1;
923: } else {
924: for (j = 0; j < NbrEig; j++) select[perm[N - j - 1]] = 1;
925: }
926: /* Reorder and extract the new <r> schur vectors */
927: lwork = PetscMax(4 * N + 16, 2 * NbrEig * (N - NbrEig));
928: liwork = PetscMax(N + 6, 2 * NbrEig * (N - NbrEig));
929: PetscCall(PetscFree(work));
930: PetscCall(PetscMalloc1(lwork, &work));
931: PetscCall(PetscMalloc1(liwork, &iwork));
932: {
933: PetscBLASInt info = 0;
934: PetscReal Dif[2];
935: PetscBLASInt ijob = 2;
936: PetscBLASInt wantQ = 1, wantZ = 1;
937: PetscCallBLAS("LAPACKtgsen", LAPACKtgsen_(&ijob, &wantQ, &wantZ, select, &N, AUAU, &ldA, AUU, &ldA, wr, wi, beta, Q, &N, Z, &N, &NbrEig, NULL, NULL, &Dif[0], work, &lwork, iwork, &liwork, &info));
938: PetscCheck(info != 1, PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB, "Unable to reorder the eigenvalues with the LAPACK routine: ill-conditioned problem.");
939: }
940: PetscCall(PetscFree(select));
942: for (j = 0; j < r; j++) PetscCall(PetscArraycpy(&SR2[j * aug1], &Z[j * N], N));
944: /* Multiply the Schur vectors SR2 by U (and X) to get a new U
945: -- save it temporarily in MU */
946: for (j = 0; j < r; j++) {
947: PetscCall(VecMAXPBY(MU[j], r_old, &SR2[j * aug1], 0, UU));
948: PetscCall(VecMAXPY(MU[j], neig, &SR2[j * aug1 + r_old], XX));
949: }
950: /* Form T = U'*MU*U */
951: for (j = 0; j < r; j++) {
952: PetscCall(VecCopy(MU[j], UU[j]));
953: PetscCall(KSP_PCApplyBAorAB(ksp, UU[j], MU[j], VEC_TEMP_MATOP));
954: }
955: dgmres->matvecs += r;
956: for (j = 0; j < r; j++) PetscCall(VecMDot(MU[j], r, UU, &TT[j * bmax]));
957: /* Factorize T */
958: PetscCall(PetscArraycpy(TTF, TT, bmax * r));
959: PetscCall(PetscBLASIntCast(r, &nr));
960: PetscCall(PetscBLASIntCast(bmax, &bm));
961: {
962: PetscBLASInt info;
963: PetscCallBLAS("LAPACKgetrf", LAPACKgetrf_(&nr, &nr, TTF, &bm, INVP, &info));
964: PetscCheck(!info, PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB, "Error in LAPACK routine XGETRF INFO=%d", (int)info);
965: }
966: /* Free Memory */
967: PetscCall(PetscFree(wr));
968: PetscCall(PetscFree(wi));
969: PetscCall(PetscFree(beta));
970: PetscCall(PetscFree(modul));
971: PetscCall(PetscFree(perm));
972: PetscCall(PetscFree(Q));
973: PetscCall(PetscFree(Z));
974: PetscCall(PetscFree(work));
975: PetscCall(PetscFree(iwork));
976: PetscFunctionReturn(PETSC_SUCCESS);
977: }
979: /*MC
980: KSPDGMRES - Implements the deflated GMRES as defined in {cite}`erhel1996restarted` and {cite}`wakam2013memory`
982: Options Database Keys:
983: GMRES Options (inherited):
984: + -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against
985: . -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
986: . -ksp_gmres_preallocate - preallocate all the Krylov search directions initially
987: (otherwise groups of vectors are allocated as needed)
988: . -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against
989: the Krylov space (fast) (the default)
990: . -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
991: . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the
992: stability of the classical Gram-Schmidt orthogonalization.
993: - -ksp_gmres_krylov_monitor - plot the Krylov space generated
995: DGMRES Options Database Keys:
996: + -ksp_dgmres_eigen <neig> - number of smallest eigenvalues to extract at each restart
997: . -ksp_dgmres_max_eigen <max_neig> - maximum number of eigenvalues that can be extracted during the iterative process
998: . -ksp_dgmres_force - use the deflation at each restart; switch off the adaptive strategy.
999: - -ksp_dgmres_view_deflation_vecs <viewerspec> - View the deflation vectors, where viewerspec is a key that can be
1000: parsed by `PetscOptionsCreateViewer()`. If neig > 1, viewerspec should
1001: end with ":append". No vectors will be viewed if the adaptive
1002: strategy chooses not to deflate, so -ksp_dgmres_force should also
1003: be given.
1004: The deflation vectors span a subspace that may be a good
1005: approximation of the subspace of smallest eigenvectors of the
1006: preconditioned operator, so this option can aid in understanding
1007: the performance of a preconditioner.
1009: Level: beginner
1011: Notes:
1012: Left and right preconditioning are supported, but not symmetric preconditioning. Complex arithmetic is not supported
1014: In this implementation, the adaptive strategy allows switching to deflated GMRES when the stagnation occurs.
1016: Contributed by:
1017: Desire NUENTSA WAKAM, INRIA
1019: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPFGMRES`, `KSPLGMRES`,
1020: `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`,
1021: `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`,
1022: `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()`, `KSPSetPCSide()`
1023: M*/
1025: PETSC_EXTERN PetscErrorCode KSPCreate_DGMRES(KSP ksp)
1026: {
1027: KSP_DGMRES *dgmres;
1029: PetscFunctionBegin;
1030: PetscCall(PetscNew(&dgmres));
1031: ksp->data = (void *)dgmres;
1033: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
1034: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
1035: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
1037: ksp->ops->buildsolution = KSPBuildSolution_DGMRES;
1038: ksp->ops->setup = KSPSetUp_DGMRES;
1039: ksp->ops->solve = KSPSolve_DGMRES;
1040: ksp->ops->destroy = KSPDestroy_DGMRES;
1041: ksp->ops->view = KSPView_DGMRES;
1042: ksp->ops->setfromoptions = KSPSetFromOptions_DGMRES;
1043: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
1044: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES;
1046: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetPreAllocateVectors_C", KSPGMRESSetPreAllocateVectors_GMRES));
1047: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetOrthogonalization_C", KSPGMRESSetOrthogonalization_GMRES));
1048: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetRestart_C", KSPGMRESSetRestart_GMRES));
1049: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetHapTol_C", KSPGMRESSetHapTol_GMRES));
1050: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetCGSRefinementType_C", KSPGMRESSetCGSRefinementType_GMRES));
1051: /* -- New functions defined in DGMRES -- */
1052: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESSetEigen_C", KSPDGMRESSetEigen_DGMRES));
1053: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESSetMaxEigen_C", KSPDGMRESSetMaxEigen_DGMRES));
1054: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESSetRatio_C", KSPDGMRESSetRatio_DGMRES));
1055: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESForce_C", KSPDGMRESForce_DGMRES));
1056: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESComputeSchurForm_C", KSPDGMRESComputeSchurForm_DGMRES));
1057: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESComputeDeflationData_C", KSPDGMRESComputeDeflationData_DGMRES));
1058: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESApplyDeflation_C", KSPDGMRESApplyDeflation_DGMRES));
1059: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESImproveEig_C", KSPDGMRESImproveEig_DGMRES));
1061: PetscCall(PetscLogEventRegister("DGMRESCompDefl", KSP_CLASSID, &KSP_DGMRESComputeDeflationData));
1062: PetscCall(PetscLogEventRegister("DGMRESApplyDefl", KSP_CLASSID, &KSP_DGMRESApplyDeflation));
1064: dgmres->haptol = 1.0e-30;
1065: dgmres->q_preallocate = 0;
1066: dgmres->delta_allocate = GMRES_DELTA_DIRECTIONS;
1067: dgmres->orthog = KSPGMRESClassicalGramSchmidtOrthogonalization;
1068: dgmres->nrs = NULL;
1069: dgmres->sol_temp = NULL;
1070: dgmres->max_k = GMRES_DEFAULT_MAXK;
1071: dgmres->Rsvd = NULL;
1072: dgmres->cgstype = KSP_GMRES_CGS_REFINE_NEVER;
1073: dgmres->orthogwork = NULL;
1075: /* Default values for the deflation */
1076: dgmres->r = 0;
1077: dgmres->neig = DGMRES_DEFAULT_EIG;
1078: dgmres->max_neig = DGMRES_DEFAULT_MAXEIG - 1;
1079: dgmres->lambdaN = 0.0;
1080: dgmres->smv = SMV;
1081: dgmres->matvecs = 0;
1082: dgmres->GreatestEig = PETSC_FALSE; /* experimental */
1083: dgmres->HasSchur = PETSC_FALSE;
1084: PetscFunctionReturn(PETSC_SUCCESS);
1085: }