Actual source code: agmres.c
1: #include <../src/ksp/ksp/impls/gmres/agmres/agmresimpl.h>
3: static PetscErrorCode KSPAGMRESBuildSoln(KSP, PetscInt);
4: static PetscErrorCode KSPAGMRESBuildBasis(KSP);
5: static PetscErrorCode KSPAGMRESBuildHessenberg(KSP);
6: extern PetscErrorCode KSPSetUp_DGMRES(KSP);
7: extern PetscErrorCode KSPBuildSolution_DGMRES(KSP, Vec, Vec *);
8: extern PetscErrorCode KSPSolve_DGMRES(KSP);
9: extern PetscErrorCode KSPDGMRESComputeDeflationData_DGMRES(KSP, PetscInt *);
10: extern PetscErrorCode KSPDGMRESComputeSchurForm_DGMRES(KSP, PetscInt *);
11: extern PetscErrorCode KSPDGMRESApplyDeflation_DGMRES(KSP, Vec, Vec);
12: extern PetscErrorCode KSPDestroy_DGMRES(KSP);
13: extern PetscErrorCode KSPSetFromOptions_DGMRES(KSP, PetscOptionItems);
14: extern PetscErrorCode KSPDGMRESSetEigen_DGMRES(KSP, PetscInt);
16: PetscLogEvent KSP_AGMRESComputeDeflationData, KSP_AGMRESBuildBasis, KSP_AGMRESComputeShifts, KSP_AGMRESRoddec;
18: /*
19: This function allocates data for the Newton basis GMRES implementation.
20: Note that most data are allocated in KSPSetUp_DGMRES and KSPSetUp_GMRES, including the space for the basis vectors,
21: the various Hessenberg matrices and the Givens rotations coefficients
22: */
23: static PetscErrorCode KSPSetUp_AGMRES(KSP ksp)
24: {
25: PetscInt hes;
26: PetscInt nloc;
27: KSP_AGMRES *agmres = (KSP_AGMRES *)ksp->data;
28: PetscInt neig = agmres->neig;
29: const PetscInt max_k = agmres->max_k;
30: PetscInt N = MAXKSPSIZE;
31: PetscInt lwork = PetscMax(8 * N + 16, 4 * neig * (N - neig));
33: PetscFunctionBegin;
34: N = MAXKSPSIZE;
35: /* Preallocate space during the call to KSPSetup_GMRES for the Krylov basis */
36: agmres->q_preallocate = PETSC_TRUE; /* No allocation on the fly */
37: /* Preallocate space to compute later the eigenvalues in GMRES */
38: ksp->calc_sings = PETSC_TRUE;
39: agmres->max_k = N; /* Set the augmented size to be allocated in KSPSetup_GMRES */
40: PetscCall(KSPSetUp_DGMRES(ksp));
41: agmres->max_k = max_k;
42: hes = (N + 1) * (N + 1);
44: /* Data for the Newton basis GMRES */
45: PetscCall(PetscCalloc4(max_k, &agmres->Rshift, max_k, &agmres->Ishift, hes, &agmres->Rloc, (N + 1) * 4, &agmres->wbufptr));
46: PetscCall(PetscMalloc3(N + 1, &agmres->tau, lwork, &agmres->work, N + 1, &agmres->nrs));
47: PetscCall(PetscCalloc4(N + 1, &agmres->Scale, N + 1, &agmres->sgn, N + 1, &agmres->tloc, N + 1, &agmres->temp));
49: /* Allocate space for the vectors in the orthogonalized basis*/
50: PetscCall(VecGetLocalSize(agmres->vecs[0], &nloc));
51: PetscCall(PetscMalloc1(nloc * (N + 1), &agmres->Qloc));
53: /* Init the ring of processors for the roddec orthogonalization */
54: PetscCall(KSPAGMRESRoddecInitNeighboor(ksp));
56: if (agmres->neig < 1) PetscFunctionReturn(PETSC_SUCCESS);
58: /* Allocate space for the deflation */
59: PetscCall(PetscMalloc1(N, &agmres->select));
60: PetscCall(VecDuplicateVecs(VEC_V(0), N, &agmres->TmpU));
61: PetscCall(PetscMalloc2(N * N, &agmres->MatEigL, N * N, &agmres->MatEigR));
62: /* PetscCall(PetscMalloc6(N*N, &agmres->Q, N*N, &agmres->Z, N, &agmres->wr, N, &agmres->wi, N, &agmres->beta, N, &agmres->modul)); */
63: PetscCall(PetscMalloc3(N * N, &agmres->Q, N * N, &agmres->Z, N, &agmres->beta));
64: PetscCall(PetscMalloc2(N + 1, &agmres->perm, 2 * neig * N, &agmres->iwork));
65: PetscFunctionReturn(PETSC_SUCCESS);
66: }
68: /*
69: Returns the current solution from the private data structure of AGMRES back to ptr.
70: */
71: static PetscErrorCode KSPBuildSolution_AGMRES(KSP ksp, Vec ptr, Vec *result)
72: {
73: KSP_AGMRES *agmres = (KSP_AGMRES *)ksp->data;
75: PetscFunctionBegin;
76: if (!ptr) {
77: if (!agmres->sol_temp) {
78: PetscCall(VecDuplicate(ksp->vec_sol, &agmres->sol_temp));
79: PetscCall(VecCopy(ksp->vec_sol, agmres->sol_temp));
80: }
81: ptr = agmres->sol_temp;
82: } else {
83: PetscCall(VecCopy(ksp->vec_sol, ptr));
84: }
85: if (result) *result = ptr;
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: /* Computes the shift values (Ritz values) needed to generate stable basis vectors
90: One cycle of DGMRES is performed to find the eigenvalues. The same data structures are used since AGMRES extends DGMRES
91: Note that when the basis is to be augmented, then this function computes the harmonic Ritz vectors from this first cycle.
92: Input :
93: - The operators (matrix, preconditioners and right-hand side) are normally required.
94: - max_k : the size of the (non augmented) basis.
95: - neig: The number of eigenvectors to augment, if deflation is needed
96: Output :
97: - The shifts as complex pair of arrays in wr and wi (size max_k).
98: - The harmonic Ritz vectors (agmres->U) if deflation is needed.
99: */
100: static PetscErrorCode KSPComputeShifts_DGMRES(KSP ksp)
101: {
102: KSP_AGMRES *agmres = (KSP_AGMRES *)ksp->data;
103: PetscInt max_k = agmres->max_k; /* size of the (non augmented) Krylov subspace */
104: PetscInt Neig = 0;
105: const PetscInt max_it = ksp->max_it;
106: PetscBool flg;
108: /* Perform one cycle of dgmres to find the eigenvalues and compute the first approximations of the eigenvectors */
110: PetscFunctionBegin;
111: PetscCall(PetscLogEventBegin(KSP_AGMRESComputeShifts, ksp, 0, 0, 0));
112: /* Send the size of the augmented basis to DGMRES */
113: ksp->max_it = max_k; /* set this to have DGMRES performing only one cycle */
114: ksp->ops->buildsolution = KSPBuildSolution_DGMRES;
115: PetscCall(KSPSolve_DGMRES(ksp));
116: ksp->guess_zero = PETSC_FALSE;
117: if (ksp->reason == KSP_CONVERGED_RTOL) {
118: PetscCall(PetscLogEventEnd(KSP_AGMRESComputeShifts, ksp, 0, 0, 0));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: } else ksp->reason = KSP_CONVERGED_ITERATING;
122: if ((agmres->r == 0) && (agmres->neig > 0)) { /* Compute the eigenvalues for the shifts and the eigenvectors (to augment the Newton basis) */
123: agmres->HasSchur = PETSC_FALSE;
124: PetscCall(KSPDGMRESComputeDeflationData_DGMRES(ksp, &Neig));
125: Neig = max_k;
126: } else { /* From DGMRES, compute only the eigenvalues needed as Shifts for the Newton Basis */
127: PetscCall(KSPDGMRESComputeSchurForm_DGMRES(ksp, &Neig));
128: }
130: /* It may happen that the Ritz values from one cycle of GMRES are not accurate enough to provide a good stability. In this case, another cycle of GMRES is performed. The two sets of values thus generated are sorted and the most accurate are kept as shifts */
131: PetscCall(PetscOptionsHasName(NULL, NULL, "-ksp_agmres_ImproveShifts", &flg));
132: if (!flg) {
133: PetscCall(KSPAGMRESLejaOrdering(agmres->wr, agmres->wi, agmres->Rshift, agmres->Ishift, max_k));
134: } else { /* Perform another cycle of DGMRES to find another set of eigenvalues */
135: PetscInt i;
136: PetscScalar *wr, *wi, *Rshift, *Ishift;
137: PetscCall(PetscMalloc4(2 * max_k, &wr, 2 * max_k, &wi, 2 * max_k, &Rshift, 2 * max_k, &Ishift));
138: for (i = 0; i < max_k; i++) {
139: wr[i] = agmres->wr[i];
140: wi[i] = agmres->wi[i];
141: }
143: PetscCall(KSPSolve_DGMRES(ksp));
145: ksp->guess_zero = PETSC_FALSE;
146: if (ksp->reason == KSP_CONVERGED_RTOL) PetscFunctionReturn(PETSC_SUCCESS);
147: else ksp->reason = KSP_CONVERGED_ITERATING;
148: if (agmres->neig > 0) { /* Compute the eigenvalues for the shifts) and the eigenvectors (to augment the Newton basis */
149: agmres->HasSchur = PETSC_FALSE;
151: PetscCall(KSPDGMRESComputeDeflationData_DGMRES(ksp, &Neig));
152: Neig = max_k;
153: } else { /* From DGMRES, compute only the eigenvalues needed as Shifts for the Newton Basis */
154: PetscCall(KSPDGMRESComputeSchurForm_DGMRES(ksp, &Neig));
155: }
156: for (i = 0; i < max_k; i++) {
157: wr[max_k + i] = agmres->wr[i];
158: wi[max_k + i] = agmres->wi[i];
159: }
160: PetscCall(KSPAGMRESLejaOrdering(wr, wi, Rshift, Ishift, 2 * max_k));
161: for (i = 0; i < max_k; i++) {
162: agmres->Rshift[i] = Rshift[i];
163: agmres->Ishift[i] = Ishift[i];
164: }
165: PetscCall(PetscFree(Rshift));
166: PetscCall(PetscFree(wr));
167: PetscCall(PetscFree(Ishift));
168: PetscCall(PetscFree(wi));
169: }
170: agmres->HasShifts = PETSC_TRUE;
171: ksp->max_it = max_it;
172: PetscCall(PetscLogEventEnd(KSP_AGMRESComputeShifts, ksp, 0, 0, 0));
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: /*
177: Generate the basis vectors from the Newton polynomials with shifts and scaling factors
178: The scaling factors are computed to obtain unit vectors. Note that this step can be avoided with the preprocessing option KSP_AGMRES_NONORM.
179: Inputs :
180: - Operators (Matrix and preconditioners and the first basis vector in VEC_V(0)
181: - Shifts values in agmres->Rshift and agmres->Ishift.
182: Output :
183: - agmres->vecs or VEC_V : basis vectors
184: - agmres->Scale : Scaling factors (equal to 1 if no scaling is done)
185: */
186: static PetscErrorCode KSPAGMRESBuildBasis(KSP ksp)
187: {
188: KSP_AGMRES *agmres = (KSP_AGMRES *)ksp->data;
189: PetscReal *Rshift = agmres->Rshift;
190: PetscReal *Ishift = agmres->Ishift;
191: PetscReal *Scale = agmres->Scale;
192: const PetscInt max_k = agmres->max_k;
193: PetscInt KspSize = KSPSIZE; /* if max_k == KspSizen then the basis should not be augmented */
194: PetscInt j = 1;
196: PetscFunctionBegin;
197: PetscCall(PetscLogEventBegin(KSP_AGMRESBuildBasis, ksp, 0, 0, 0));
198: Scale[0] = 1.0;
199: while (j <= max_k) {
200: if (Ishift[j - 1] == 0) {
201: if ((ksp->pc_side == PC_LEFT) && agmres->r && agmres->DeflPrecond) {
202: /* Apply the precond-matrix operators */
203: PetscCall(KSP_PCApplyBAorAB(ksp, VEC_V(j - 1), VEC_TMP, VEC_TMP_MATOP));
204: /* Then apply deflation as a preconditioner */
205: PetscCall(KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_TMP, VEC_V(j)));
206: } else if ((ksp->pc_side == PC_RIGHT) && agmres->r && agmres->DeflPrecond) {
207: PetscCall(KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_V(j - 1), VEC_TMP));
208: PetscCall(KSP_PCApplyBAorAB(ksp, VEC_TMP, VEC_V(j), VEC_TMP_MATOP));
209: } else {
210: PetscCall(KSP_PCApplyBAorAB(ksp, VEC_V(j - 1), VEC_V(j), VEC_TMP_MATOP));
211: }
212: PetscCall(VecAXPY(VEC_V(j), -Rshift[j - 1], VEC_V(j - 1)));
213: #if defined(KSP_AGMRES_NONORM)
214: Scale[j] = 1.0;
215: #else
216: PetscCall(VecScale(VEC_V(j), Scale[j - 1])); /* This step can be postponed until all vectors are built */
217: PetscCall(VecNorm(VEC_V(j), NORM_2, &Scale[j]));
218: Scale[j] = 1.0 / Scale[j];
219: #endif
221: agmres->matvecs += 1;
222: j++;
223: } else {
224: if ((ksp->pc_side == PC_LEFT) && agmres->r && agmres->DeflPrecond) {
225: /* Apply the precond-matrix operators */
226: PetscCall(KSP_PCApplyBAorAB(ksp, VEC_V(j - 1), VEC_TMP, VEC_TMP_MATOP));
227: /* Then apply deflation as a preconditioner */
228: PetscCall(KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_TMP, VEC_V(j)));
229: } else if ((ksp->pc_side == PC_RIGHT) && agmres->r && agmres->DeflPrecond) {
230: PetscCall(KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_V(j - 1), VEC_TMP));
231: PetscCall(KSP_PCApplyBAorAB(ksp, VEC_TMP, VEC_V(j), VEC_TMP_MATOP));
232: } else {
233: PetscCall(KSP_PCApplyBAorAB(ksp, VEC_V(j - 1), VEC_V(j), VEC_TMP_MATOP));
234: }
235: PetscCall(VecAXPY(VEC_V(j), -Rshift[j - 1], VEC_V(j - 1)));
236: #if defined(KSP_AGMRES_NONORM)
237: Scale[j] = 1.0;
238: #else
239: PetscCall(VecScale(VEC_V(j), Scale[j - 1]));
240: PetscCall(VecNorm(VEC_V(j), NORM_2, &Scale[j]));
241: Scale[j] = 1.0 / Scale[j];
242: #endif
243: agmres->matvecs += 1;
244: j++;
245: if ((ksp->pc_side == PC_LEFT) && agmres->r && agmres->DeflPrecond) {
246: /* Apply the precond-matrix operators */
247: PetscCall(KSP_PCApplyBAorAB(ksp, VEC_V(j - 1), VEC_TMP, VEC_TMP_MATOP));
248: /* Then apply deflation as a preconditioner */
249: PetscCall(KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_TMP, VEC_V(j)));
250: } else if ((ksp->pc_side == PC_RIGHT) && agmres->r && agmres->DeflPrecond) {
251: PetscCall(KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_V(j - 1), VEC_TMP));
252: PetscCall(KSP_PCApplyBAorAB(ksp, VEC_TMP, VEC_V(j), VEC_TMP_MATOP));
253: } else {
254: PetscCall(KSP_PCApplyBAorAB(ksp, VEC_V(j - 1), VEC_V(j), VEC_TMP_MATOP));
255: }
256: PetscCall(VecAXPY(VEC_V(j), -Rshift[j - 2], VEC_V(j - 1)));
257: PetscCall(VecAXPY(VEC_V(j), Scale[j - 2] * Ishift[j - 2] * Ishift[j - 2], VEC_V(j - 2)));
258: #if defined(KSP_AGMRES_NONORM)
259: Scale[j] = 1.0;
260: #else
261: PetscCall(VecNorm(VEC_V(j), NORM_2, &Scale[j]));
262: Scale[j] = 1.0 / Scale[j];
263: #endif
264: agmres->matvecs += 1;
265: j++;
266: }
267: }
268: /* Augment the subspace with the eigenvectors*/
269: while (j <= KspSize) {
270: PetscCall(KSP_PCApplyBAorAB(ksp, agmres->U[j - max_k - 1], VEC_V(j), VEC_TMP_MATOP));
271: #if defined(KSP_AGMRES_NONORM)
272: Scale[j] = 1.0;
273: #else
274: PetscCall(VecScale(VEC_V(j), Scale[j - 1]));
275: PetscCall(VecNorm(VEC_V(j), NORM_2, &Scale[j]));
276: Scale[j] = 1.0 / Scale[j];
277: #endif
278: agmres->matvecs += 1;
279: j++;
280: }
281: PetscCall(PetscLogEventEnd(KSP_AGMRESBuildBasis, ksp, 0, 0, 0));
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: /*
286: Form the Hessenberg matrix for the Arnoldi-like relation.
287: Inputs :
288: - Shifts values in agmres->Rshift and agmres->Ishift
289: - RLoc : Triangular matrix from the RODDEC orthogonalization
290: Outputs :
291: - H = agmres->hh_origin : The Hessenberg matrix.
293: NOTE: Note that the computed Hessenberg matrix is not mathematically equivalent to that in the real Arnoldi process (in KSP GMRES). If it is needed, it can be explicitly formed as H <-- H * RLoc^-1.
294: */
295: static PetscErrorCode KSPAGMRESBuildHessenberg(KSP ksp)
296: {
297: KSP_AGMRES *agmres = (KSP_AGMRES *)ksp->data;
298: PetscScalar *Rshift = agmres->Rshift;
299: PetscScalar *Ishift = agmres->Ishift;
300: PetscScalar *Scale = agmres->Scale;
301: PetscInt i = 0, j = 0;
302: const PetscInt max_k = agmres->max_k;
303: PetscInt KspSize = KSPSIZE;
304: PetscInt N = MAXKSPSIZE + 1;
306: PetscFunctionBegin;
307: PetscCall(PetscArrayzero(agmres->hh_origin, (N + 1) * N));
308: while (j < max_k) {
309: /* Real shifts */
310: if (Ishift[j] == 0) {
311: for (i = 0; i <= j; i++) *H(i, j) = *RLOC(i, j + 1) / Scale[j] + (Rshift[j] * *RLOC(i, j));
312: *H(j + 1, j) = *RLOC(j + 1, j + 1) / Scale[j];
313: j++;
314: } else if (Ishift[j] > 0) {
315: for (i = 0; i <= j; i++) *H(i, j) = *RLOC(i, j + 1) / Scale[j] + Rshift[j] * *RLOC(i, j);
316: *H(j + 1, j) = *RLOC(j + 1, j + 1) / Scale[j];
317: j++;
318: for (i = 0; i <= j; i++) *H(i, j) = (*RLOC(i, j + 1) + Rshift[j - 1] * *RLOC(i, j) - Scale[j - 1] * Ishift[j - 1] * Ishift[j - 1] * *RLOC(i, j - 1));
319: *H(j, j) = (*RLOC(j, j + 1) + Rshift[j - 1] * *RLOC(j, j));
320: *H(j + 1, j) = *RLOC(j + 1, j + 1);
321: j++;
322: } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ORDER, "BAD ORDERING OF THE SHIFTS VALUES IN THE NEWTON BASIS");
323: }
324: for (j = max_k; j < KspSize; j++) { /* take into account the norm of the augmented vectors */
325: for (i = 0; i <= j + 1; i++) *H(i, j) = *RLOC(i, j + 1) / Scale[j];
326: }
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: /*
331: Form the new approximate solution from the least-square problem
332: */
333: static PetscErrorCode KSPAGMRESBuildSoln(KSP ksp, PetscInt it)
334: {
335: KSP_AGMRES *agmres = (KSP_AGMRES *)ksp->data;
336: const PetscInt max_k = agmres->max_k; /* Size of the non-augmented Krylov basis */
337: PetscInt i, j;
338: PetscInt r = agmres->r; /* current number of augmented eigenvectors */
339: PetscBLASInt KspSize;
340: PetscBLASInt lC;
341: PetscBLASInt N;
342: PetscBLASInt ldH;
343: PetscBLASInt lwork;
344: PetscBLASInt info, nrhs = 1;
346: PetscFunctionBegin;
347: PetscCall(PetscBLASIntCast(it + 1, &ldH));
348: PetscCall(PetscBLASIntCast(KSPSIZE, &KspSize));
349: PetscCall(PetscBLASIntCast(4 * (KspSize + 1), &lwork));
350: PetscCall(PetscBLASIntCast(KspSize + 1, &lC));
351: PetscCall(PetscBLASIntCast(MAXKSPSIZE + 1, &N));
352: PetscCall(PetscBLASIntCast(N + 1, &ldH));
353: /* Save a copy of the Hessenberg matrix */
354: for (j = 0; j < N - 1; j++) {
355: for (i = 0; i < N; i++) *HS(i, j) = *H(i, j);
356: }
357: /* QR factorize the Hessenberg matrix */
358: PetscCallBLAS("LAPACKgeqrf", LAPACKgeqrf_(&lC, &KspSize, agmres->hh_origin, &ldH, agmres->tau, agmres->work, &lwork, &info));
359: PetscCheck(!info, PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB, "Error in LAPACK routine XGEQRF INFO=%" PetscBLASInt_FMT, info);
360: /* Update the right-hand side of the least square problem */
361: PetscCall(PetscArrayzero(agmres->nrs, N));
363: agmres->nrs[0] = ksp->rnorm;
364: PetscCallBLAS("LAPACKormqr", LAPACKormqr_("L", "T", &lC, &nrhs, &KspSize, agmres->hh_origin, &ldH, agmres->tau, agmres->nrs, &N, agmres->work, &lwork, &info));
365: PetscCheck(!info, PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB, "Error in LAPACK routine XORMQR INFO=%" PetscBLASInt_FMT, info);
366: ksp->rnorm = PetscAbsScalar(agmres->nrs[KspSize]);
367: /* solve the least-square problem */
368: PetscCallBLAS("LAPACKtrtrs", LAPACKtrtrs_("U", "N", "N", &KspSize, &nrhs, agmres->hh_origin, &ldH, agmres->nrs, &N, &info));
369: PetscCheck(!info, PetscObjectComm((PetscObject)ksp), PETSC_ERR_LIB, "Error in LAPACK routine XTRTRS INFO=%" PetscBLASInt_FMT, info);
370: /* Accumulate the correction to the solution of the preconditioned problem in VEC_TMP */
371: PetscCall(VecMAXPBY(VEC_TMP, max_k, agmres->nrs, 0, &VEC_V(0)));
372: if (!agmres->DeflPrecond) PetscCall(VecMAXPY(VEC_TMP, r, &agmres->nrs[max_k], agmres->U));
374: if ((ksp->pc_side == PC_RIGHT) && agmres->r && agmres->DeflPrecond) {
375: PetscCall(KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_TMP, VEC_TMP_MATOP));
376: PetscCall(VecCopy(VEC_TMP_MATOP, VEC_TMP));
377: }
378: PetscCall(KSPUnwindPreconditioner(ksp, VEC_TMP, VEC_TMP_MATOP));
379: /* add the solution to the previous one */
380: PetscCall(VecAXPY(ksp->vec_sol, 1.0, VEC_TMP));
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: /*
385: Run one cycle of the Newton-basis gmres, possibly augmented with eigenvectors.
387: Return residual history if requested.
388: Input :
389: - The vector VEC_V(0) is the initia residual
390: Output :
391: - the solution vector is in agmres->vec_sol
392: - itcount : number of inner iterations
393: - res : the new residual norm
394: .
395: NOTE: Unlike GMRES where the residual norm is available at each (inner) iteration, here it is available at the end of the cycle.
396: */
397: static PetscErrorCode KSPAGMRESCycle(PetscInt *itcount, KSP ksp)
398: {
399: KSP_AGMRES *agmres = (KSP_AGMRES *)ksp->data;
400: PetscReal res;
401: PetscInt KspSize = KSPSIZE;
403: PetscFunctionBegin;
404: /* check for the convergence */
405: res = ksp->rnorm; /* Norm of the initial residual vector */
406: if (!res) {
407: if (itcount) *itcount = 0;
408: ksp->reason = KSP_CONVERGED_ATOL;
409: PetscCall(PetscInfo(ksp, "Converged due to zero residual norm on entry\n"));
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
412: PetscCall((*ksp->converged)(ksp, ksp->its, res, &ksp->reason, ksp->cnvP));
413: /* Build the Krylov basis with Newton polynomials */
414: PetscCall(KSPAGMRESBuildBasis(ksp));
415: /* QR Factorize the basis with RODDEC */
416: PetscCall(KSPAGMRESRoddec(ksp, KspSize + 1));
418: /* Recover a (partial) Hessenberg matrix for the Arnoldi-like relation */
419: PetscCall(KSPAGMRESBuildHessenberg(ksp));
420: /* Solve the least square problem and unwind the preconditioner */
421: PetscCall(KSPAGMRESBuildSoln(ksp, KspSize));
423: res = ksp->rnorm;
424: ksp->its += KspSize;
425: agmres->it = KspSize - 1;
426: /* Test for the convergence */
427: PetscCall((*ksp->converged)(ksp, ksp->its, res, &ksp->reason, ksp->cnvP));
428: PetscCall(KSPLogResidualHistory(ksp, res));
429: PetscCall(KSPMonitor(ksp, ksp->its, res));
431: *itcount = KspSize;
432: PetscFunctionReturn(PETSC_SUCCESS);
433: }
435: static PetscErrorCode KSPSolve_AGMRES(KSP ksp)
436: {
437: PetscInt its;
438: KSP_AGMRES *agmres = (KSP_AGMRES *)ksp->data;
439: PetscBool guess_zero = ksp->guess_zero;
440: PetscReal res_old, res;
441: PetscInt test;
443: PetscFunctionBegin;
444: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
445: ksp->its = 0;
446: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
448: if (!agmres->HasShifts) { /* Compute Shifts for the Newton basis */
449: PetscCall(KSPComputeShifts_DGMRES(ksp));
450: }
451: /* NOTE: At this step, the initial guess is not equal to zero since one cycle of the classical GMRES is performed to compute the shifts */
452: PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP));
453: while (!ksp->reason) {
454: PetscCall(KSPInitialResidual(ksp, ksp->vec_sol, VEC_TMP, VEC_TMP_MATOP, VEC_V(0), ksp->vec_rhs));
455: if ((ksp->pc_side == PC_LEFT) && agmres->r && agmres->DeflPrecond) {
456: PetscCall(KSPDGMRESApplyDeflation_DGMRES(ksp, VEC_V(0), VEC_TMP));
457: PetscCall(VecCopy(VEC_TMP, VEC_V(0)));
459: agmres->matvecs += 1;
460: }
461: PetscCall(VecNormalize(VEC_V(0), &ksp->rnorm));
462: KSPCheckNorm(ksp, ksp->rnorm);
463: res_old = ksp->rnorm; /* Record the residual norm to test if deflation is needed */
465: ksp->ops->buildsolution = KSPBuildSolution_AGMRES;
467: PetscCall(KSPAGMRESCycle(&its, ksp));
468: if (ksp->its >= ksp->max_it) {
469: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
470: break;
471: }
472: /* compute the eigenvectors to augment the subspace : use an adaptive strategy */
473: res = ksp->rnorm;
474: if (!ksp->reason && agmres->neig > 0) {
475: test = agmres->max_k * PetscLogReal(ksp->rtol / res) / PetscLogReal(res / res_old); /* estimate the remaining number of steps */
476: if ((test > agmres->smv * (ksp->max_it - ksp->its)) || agmres->force) {
477: /* Augment the number of eigenvalues to deflate if the convergence is too slow */
478: if (!agmres->force && ((test > agmres->bgv * (ksp->max_it - ksp->its)) && ((agmres->r + 1) < agmres->max_neig))) agmres->neig += 1;
479: PetscCall(KSPDGMRESComputeDeflationData_DGMRES(ksp, &agmres->neig));
480: }
481: }
482: ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
483: }
484: ksp->guess_zero = guess_zero; /* restore if user has provided nonzero initial guess */
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: static PetscErrorCode KSPDestroy_AGMRES(KSP ksp)
489: {
490: KSP_AGMRES *agmres = (KSP_AGMRES *)ksp->data;
492: PetscFunctionBegin;
493: PetscCall(PetscFree(agmres->hh_origin));
495: PetscCall(PetscFree(agmres->Qloc));
496: PetscCall(PetscFree4(agmres->Rshift, agmres->Ishift, agmres->Rloc, agmres->wbufptr));
497: PetscCall(PetscFree3(agmres->tau, agmres->work, agmres->nrs));
498: PetscCall(PetscFree4(agmres->Scale, agmres->sgn, agmres->tloc, agmres->temp));
500: PetscCall(PetscFree(agmres->select));
501: PetscCall(PetscFree(agmres->wr));
502: PetscCall(PetscFree(agmres->wi));
503: if (agmres->neig) {
504: PetscCall(VecDestroyVecs(MAXKSPSIZE, &agmres->TmpU));
505: PetscCall(PetscFree(agmres->perm));
506: PetscCall(PetscFree(agmres->MatEigL));
507: PetscCall(PetscFree(agmres->MatEigR));
508: PetscCall(PetscFree(agmres->Q));
509: PetscCall(PetscFree(agmres->Z));
510: PetscCall(PetscFree(agmres->beta));
511: }
512: PetscCall(KSPDestroy_DGMRES(ksp));
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }
516: static PetscErrorCode KSPView_AGMRES(KSP ksp, PetscViewer viewer)
517: {
518: KSP_AGMRES *agmres = (KSP_AGMRES *)ksp->data;
519: const char *cstr = "RODDEC ORTHOGONOLIZATION";
520: PetscBool isascii, isstring;
521: #if defined(KSP_AGMRES_NONORM)
522: const char *Nstr = "SCALING FACTORS : NO";
523: #else
524: const char *Nstr = "SCALING FACTORS : YES";
525: #endif
527: PetscFunctionBegin;
528: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
529: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
531: if (isascii) {
532: PetscCall(PetscViewerASCIIPrintf(viewer, " restart=%" PetscInt_FMT " using %s\n", agmres->max_k, cstr));
533: PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", Nstr));
534: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of matvecs : %" PetscInt_FMT "\n", agmres->matvecs));
535: PetscCall(PetscViewerASCIIPrintf(viewer, " Adaptive strategy is used: %s\n", PetscBools[agmres->force]));
536: if (agmres->DeflPrecond) {
537: PetscCall(PetscViewerASCIIPrintf(viewer, " STRATEGY OF DEFLATION: PRECONDITIONER \n"));
538: PetscCall(PetscViewerASCIIPrintf(viewer, " Frequency of extracted eigenvalues = %" PetscInt_FMT "\n", agmres->neig));
539: PetscCall(PetscViewerASCIIPrintf(viewer, " Total number of extracted eigenvalues = %" PetscInt_FMT "\n", agmres->r));
540: PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum number of eigenvalues set to be extracted = %" PetscInt_FMT "\n", agmres->max_neig));
541: } else {
542: PetscCall(PetscViewerASCIIPrintf(viewer, " STRATEGY OF DEFLATION: AUGMENT\n"));
543: PetscCall(PetscViewerASCIIPrintf(viewer, " augmented vectors %" PetscInt_FMT " at frequency %" PetscInt_FMT " with %sRitz vectors\n", agmres->r, agmres->neig, agmres->ritz ? "" : "Harmonic "));
544: }
545: PetscCall(PetscViewerASCIIPrintf(viewer, " Minimum relaxation parameter for the adaptive strategy(smv) = %g\n", (double)agmres->smv));
546: PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum relaxation parameter for the adaptive strategy(bgv) = %g\n", (double)agmres->bgv));
547: } else if (isstring) {
548: PetscCall(PetscViewerStringSPrintf(viewer, "%s restart %" PetscInt_FMT, cstr, agmres->max_k));
549: }
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }
553: static PetscErrorCode KSPSetFromOptions_AGMRES(KSP ksp, PetscOptionItems PetscOptionsObject)
554: {
555: PetscInt neig;
556: KSP_AGMRES *agmres = (KSP_AGMRES *)ksp->data;
557: PetscBool flg;
559: PetscFunctionBegin;
560: PetscCall(KSPSetFromOptions_DGMRES(ksp, PetscOptionsObject)); /* Set common options from DGMRES and GMRES */
561: PetscOptionsHeadBegin(PetscOptionsObject, "KSP AGMRES Options");
562: PetscCall(PetscOptionsInt("-ksp_agmres_eigen", "Number of eigenvalues to deflate", "KSPDGMRESSetEigen", agmres->neig, &neig, &flg));
563: if (flg) {
564: PetscCall(KSPDGMRESSetEigen_DGMRES(ksp, neig));
565: agmres->r = 0;
566: } else agmres->neig = 0;
567: PetscCall(PetscOptionsInt("-ksp_agmres_maxeigen", "Maximum number of eigenvalues to deflate", "KSPDGMRESSetMaxEigen", agmres->max_neig, &neig, &flg));
568: if (flg) agmres->max_neig = neig + EIG_OFFSET;
569: else agmres->max_neig = agmres->neig + EIG_OFFSET;
570: PetscCall(PetscOptionsBool("-ksp_agmres_DeflPrecond", "Determine if the deflation should be applied as a preconditioner -- similar to KSP DGMRES", "KSPGMRESDeflPrecond", agmres->DeflPrecond, &agmres->DeflPrecond, NULL));
571: PetscCall(PetscOptionsBool("-ksp_agmres_ritz", "Compute the Ritz vectors instead of the Harmonic Ritz vectors ", "KSPGMRESHarmonic", agmres->ritz, &agmres->ritz, &flg));
572: PetscCall(PetscOptionsReal("-ksp_agmres_MinRatio", "Relaxation parameter in the adaptive strategy; smallest multiple of the remaining number of steps allowed", "KSPGMRESSetMinRatio", agmres->smv, &agmres->smv, NULL));
573: PetscCall(PetscOptionsReal("-ksp_agmres_MaxRatio", "Relaxation parameter in the adaptive strategy; Largest multiple of the remaining number of steps allowed", "KSPGMRESSetMaxRatio", agmres->bgv, &agmres->bgv, &flg));
574: PetscOptionsHeadEnd();
575: PetscFunctionReturn(PETSC_SUCCESS);
576: }
578: /*MC
579: KSPAGMRES - a Newton basis GMRES implementation with adaptive augmented eigenvectors.
581: Options Database Keys:
582: + -ksp_gmres_restart restart - the number of Krylov directions
583: . -ksp_gmres_krylov_monitor - plot the Krylov space generated
584: . -ksp_agmres_eigen neig - Number of eigenvalues to deflate (Number of vectors to augment)
585: . -ksp_agmres_maxeigen max_neig - Maximum number of eigenvalues to deflate
586: . -ksp_agmres_MinRatio minr - Relaxation parameter in the adaptive strategy; smallest multiple of the remaining number of steps allowed
587: . -ksp_agmres_MaxRatio maxr - Relaxation parameter in the adaptive strategy; Largest multiple of the remaining number of steps allowed
588: . -ksp_agmres_DeflPrecond - Apply deflation as a preconditioner, this is similar to `KSPDGMRES` but it rather builds a Newton basis.
589: - -ksp_dgmres_force (true|false) - Force the deflation at each restart.
591: Level: intermediate
593: Notes:
594: Left and right preconditioning are supported, but not symmetric preconditioning. Complex arithmetic is not supported
596: This code can be used as an experimental framework to combine several techniques in the particular case of GMRES.
597: For instance, the computation of the shifts can be improved with techniques described in {cite}`philippe2012generation`. The orthogonalization technique can be replaced by TSQR {cite}`demmel2012communication`.
599: The techniques used are best described in {cite}`wakam2011parallelism`. The contribution of this work is that it combines many of the previous work to reduce the amount
600: of MPI messages and improve the robustness of the global approach by using deflation techniques. Please see {cite}`wakam2011parallelism` for numerical experiments and {cite}`wakam2013memory`
601: for a description of these problems.
602: The generation of the basis can be done using s-steps approaches{cite}`mohiyuddin2009minimizing`. See also {cite}`sidje1997alternatives` and {cite}`bai1994newton`.
604: Developer Note:
605: This object is subclassed off of `KSPDGMRES`, see the source code in src/ksp/ksp/impls/gmres for comments on the structure of the code
607: Contributed by:
608: Desire NUENTSA WAKAM, INRIA <desire.nuentsa_wakam@inria.fr> with inputs from Guy Atenekeng <atenekeng@yahoo.com> and R.B. Sidje <roger.b.sidje@ua.edu>
610: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPDGMRES`, `KSPPGMRES`,
611: `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESSetOrthogonalization()`, `KSPGMRESGetOrthogonalization()`,
612: `KSPGMRESClassicalGramSchmidtOrthogonalization()`, `KSPGMRESModifiedGramSchmidtOrthogonalization()`,
613: `KSPGMRESCGSRefinementType`, `KSPGMRESSetCGSRefinementType()`, `KSPGMRESGetCGSRefinementType()`, `KSPGMRESMonitorKrylov()`, `KSPSetPCSide()`
614: M*/
616: PETSC_EXTERN PetscErrorCode KSPCreate_AGMRES(KSP ksp)
617: {
618: KSP_AGMRES *agmres;
620: PetscFunctionBegin;
621: PetscCall(PetscNew(&agmres));
622: ksp->data = (void *)agmres;
624: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
625: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
626: ksp->ops->buildsolution = KSPBuildSolution_AGMRES;
627: ksp->ops->setup = KSPSetUp_AGMRES;
628: ksp->ops->solve = KSPSolve_AGMRES;
629: ksp->ops->destroy = KSPDestroy_AGMRES;
630: ksp->ops->view = KSPView_AGMRES;
631: ksp->ops->setfromoptions = KSPSetFromOptions_AGMRES;
632: ksp->guess_zero = PETSC_TRUE;
633: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
634: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES;
636: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetPreAllocateVectors_C", KSPGMRESSetPreAllocateVectors_GMRES));
637: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetOrthogonalization_C", KSPGMRESSetOrthogonalization_GMRES));
638: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetRestart_C", KSPGMRESSetRestart_GMRES));
639: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetHapTol_C", KSPGMRESSetHapTol_GMRES));
640: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetCGSRefinementType_C", KSPGMRESSetCGSRefinementType_GMRES));
641: /* -- New functions defined in DGMRES -- */
642: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESSetEigen_C", KSPDGMRESSetEigen_DGMRES));
643: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESComputeSchurForm_C", KSPDGMRESComputeSchurForm_DGMRES));
644: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESComputeDeflationData_C", KSPDGMRESComputeDeflationData_DGMRES));
645: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPDGMRESApplyDeflation_C", KSPDGMRESApplyDeflation_DGMRES));
647: PetscCall(PetscLogEventRegister("AGMRESCompDefl", KSP_CLASSID, &KSP_AGMRESComputeDeflationData));
648: PetscCall(PetscLogEventRegister("AGMRESBuildBasis", KSP_CLASSID, &KSP_AGMRESBuildBasis));
649: PetscCall(PetscLogEventRegister("AGMRESCompShifts", KSP_CLASSID, &KSP_AGMRESComputeShifts));
650: PetscCall(PetscLogEventRegister("AGMRESOrthog", KSP_CLASSID, &KSP_AGMRESRoddec));
652: agmres->haptol = 1.0e-30;
653: agmres->q_preallocate = 0;
654: agmres->delta_allocate = AGMRES_DELTA_DIRECTIONS;
655: agmres->orthog = KSPGMRESClassicalGramSchmidtOrthogonalization;
656: agmres->nrs = NULL;
657: agmres->sol_temp = NULL;
658: agmres->max_k = AGMRES_DEFAULT_MAXK;
659: agmres->Rsvd = NULL;
660: agmres->cgstype = KSP_GMRES_CGS_REFINE_NEVER;
661: agmres->orthogwork = NULL;
663: /* Default values for the deflation */
664: agmres->r = 0;
665: agmres->neig = 0;
666: agmres->max_neig = 0;
667: agmres->lambdaN = 0.0;
668: agmres->smv = SMV;
669: agmres->bgv = 1;
670: agmres->force = PETSC_FALSE;
671: agmres->matvecs = 0;
672: agmres->improve = PETSC_FALSE;
673: agmres->HasShifts = PETSC_FALSE;
674: agmres->r = 0;
675: agmres->HasSchur = PETSC_FALSE;
676: agmres->DeflPrecond = PETSC_FALSE;
677: PetscCall(PetscObjectGetNewTag((PetscObject)ksp, &agmres->tag));
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }