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