Actual source code: pgmres.c
2: /*
3: This file implements PGMRES (a Pipelined Generalized Minimal Residual method)
4: */
6: #include <../src/ksp/ksp/impls/gmres/pgmres/pgmresimpl.h>
7: #define PGMRES_DELTA_DIRECTIONS 10
8: #define PGMRES_DEFAULT_MAXK 30
10: static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP,PetscInt,PetscBool*,PetscReal*);
11: static PetscErrorCode KSPPGMRESBuildSoln(PetscScalar*,Vec,Vec,KSP,PetscInt);
13: /*
15: KSPSetUp_PGMRES - Sets up the workspace needed by pgmres.
17: This is called once, usually automatically by KSPSolve() or KSPSetUp(),
18: but can be called directly by KSPSetUp().
20: */
21: static PetscErrorCode KSPSetUp_PGMRES(KSP ksp)
22: {
26: KSPSetUp_GMRES(ksp);
27: return(0);
28: }
30: /*
32: KSPPGMRESCycle - Run pgmres, possibly with restart. Return residual
33: history if requested.
35: input parameters:
36: . pgmres - structure containing parameters and work areas
38: output parameters:
39: . itcount - number of iterations used. If null, ignored.
40: . converged - 0 if not converged
42: Notes:
43: On entry, the value in vector VEC_VV(0) should be
44: the initial residual.
46: */
47: static PetscErrorCode KSPPGMRESCycle(PetscInt *itcount,KSP ksp)
48: {
49: KSP_PGMRES *pgmres = (KSP_PGMRES*)(ksp->data);
50: PetscReal res_norm,res,newnorm;
52: PetscInt it = 0,j,k;
53: PetscBool hapend = PETSC_FALSE;
56: if (itcount) *itcount = 0;
57: VecNormalize(VEC_VV(0),&res_norm);
58: KSPCheckNorm(ksp,res_norm);
59: res = res_norm;
60: *RS(0) = res_norm;
62: /* check for the convergence */
63: PetscObjectSAWsTakeAccess((PetscObject)ksp);
64: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res;
65: else ksp->rnorm = 0;
66: PetscObjectSAWsGrantAccess((PetscObject)ksp);
67: pgmres->it = it-2;
68: KSPLogResidualHistory(ksp,ksp->rnorm);
69: KSPMonitor(ksp,ksp->its,ksp->rnorm);
70: if (!res) {
71: ksp->reason = KSP_CONVERGED_ATOL;
72: PetscInfo(ksp,"Converged due to zero residual norm on entry\n");
73: return(0);
74: }
76: (*ksp->converged)(ksp,ksp->its,ksp->rnorm,&ksp->reason,ksp->cnvP);
77: for (; !ksp->reason; it++) {
78: Vec Zcur,Znext;
79: if (pgmres->vv_allocated <= it + VEC_OFFSET + 1) {
80: KSPGMRESGetNewVectors(ksp,it+1);
81: }
82: /* VEC_VV(it-1) is orthogonal, it will be normalized once the VecNorm arrives. */
83: Zcur = VEC_VV(it); /* Zcur is not yet orthogonal, but the VecMDot to orthogonalize it has been started. */
84: Znext = VEC_VV(it+1); /* This iteration will compute Znext, update with a deferred correction once we know how
85: * Zcur relates to the previous vectors, and start the reduction to orthogonalize it. */
87: if (it < pgmres->max_k+1 && ksp->its+1 < PetscMax(2,ksp->max_it)) { /* We don't know whether what we have computed is enough, so apply the matrix. */
88: KSP_PCApplyBAorAB(ksp,Zcur,Znext,VEC_TEMP_MATOP);
89: }
91: if (it > 1) { /* Complete the pending reduction */
92: VecNormEnd(VEC_VV(it-1),NORM_2,&newnorm);
93: *HH(it-1,it-2) = newnorm;
94: }
95: if (it > 0) { /* Finish the reduction computing the latest column of H */
96: VecMDotEnd(Zcur,it,&(VEC_VV(0)),HH(0,it-1));
97: }
99: if (it > 1) {
100: /* normalize the base vector from two iterations ago, basis is complete up to here */
101: VecScale(VEC_VV(it-1),1./ *HH(it-1,it-2));
103: KSPPGMRESUpdateHessenberg(ksp,it-2,&hapend,&res);
104: pgmres->it = it-2;
105: ksp->its++;
106: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res;
107: else ksp->rnorm = 0;
109: (*ksp->converged)(ksp,ksp->its,ksp->rnorm,&ksp->reason,ksp->cnvP);
110: if (it < pgmres->max_k+1 || ksp->reason || ksp->its == ksp->max_it) { /* Monitor if we are done or still iterating, but not before a restart. */
111: KSPLogResidualHistory(ksp,ksp->rnorm);
112: KSPMonitor(ksp,ksp->its,ksp->rnorm);
113: }
114: if (ksp->reason) break;
115: /* Catch error in happy breakdown and signal convergence and break from loop */
116: if (hapend) {
117: if (ksp->errorifnotconverged) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"You reached the happy break down, but convergence was not indicated. Residual norm = %g",(double)res);
118: else {
119: ksp->reason = KSP_DIVERGED_BREAKDOWN;
120: break;
121: }
122: }
124: if (!(it < pgmres->max_k+1 && ksp->its < ksp->max_it)) break;
126: /* The it-2 column of H was not scaled when we computed Zcur, apply correction */
127: VecScale(Zcur,1./ *HH(it-1,it-2));
128: /* And Znext computed in this iteration was computed using the under-scaled Zcur */
129: VecScale(Znext,1./ *HH(it-1,it-2));
131: /* In the previous iteration, we projected an unnormalized Zcur against the Krylov basis, so we need to fix the column of H resulting from that projection. */
132: for (k=0; k<it; k++) *HH(k,it-1) /= *HH(it-1,it-2);
133: /* When Zcur was projected against the Krylov basis, VV(it-1) was still not normalized, so fix that too. This
134: * column is complete except for HH(it,it-1) which we won't know until the next iteration. */
135: *HH(it-1,it-1) /= *HH(it-1,it-2);
136: }
138: if (it > 0) {
139: PetscScalar *work;
140: if (!pgmres->orthogwork) {PetscMalloc1(pgmres->max_k + 2,&pgmres->orthogwork);}
141: work = pgmres->orthogwork;
142: /* Apply correction computed by the VecMDot in the last iteration to Znext. The original form is
143: *
144: * Znext -= sum_{j=0}^{i-1} Z[j+1] * H[j,i-1]
145: *
146: * where
147: *
148: * Z[j] = sum_{k=0}^j V[k] * H[k,j-1]
149: *
150: * substituting
151: *
152: * Znext -= sum_{j=0}^{i-1} sum_{k=0}^{j+1} V[k] * H[k,j] * H[j,i-1]
153: *
154: * rearranging the iteration space from row-column to column-row
155: *
156: * Znext -= sum_{k=0}^i sum_{j=k-1}^{i-1} V[k] * H[k,j] * H[j,i-1]
157: *
158: * Note that column it-1 of HH is correct. For all previous columns, we must look at HES because HH has already
159: * been transformed to upper triangular form.
160: */
161: for (k=0; k<it+1; k++) {
162: work[k] = 0;
163: for (j=PetscMax(0,k-1); j<it-1; j++) work[k] -= *HES(k,j) * *HH(j,it-1);
164: }
165: VecMAXPY(Znext,it+1,work,&VEC_VV(0));
166: VecAXPY(Znext,-*HH(it-1,it-1),Zcur);
168: /* Orthogonalize Zcur against existing basis vectors. */
169: for (k=0; k<it; k++) work[k] = -*HH(k,it-1);
170: VecMAXPY(Zcur,it,work,&VEC_VV(0));
171: /* Zcur is now orthogonal, and will be referred to as VEC_VV(it) again, though it is still not normalized. */
172: /* Begin computing the norm of the new vector, will be normalized after the MatMult in the next iteration. */
173: VecNormBegin(VEC_VV(it),NORM_2,&newnorm);
174: }
176: /* Compute column of H (to the diagonal, but not the subdiagonal) to be able to orthogonalize the newest vector. */
177: VecMDotBegin(Znext,it+1,&VEC_VV(0),HH(0,it));
179: /* Start an asynchronous split-mode reduction, the result of the MDot and Norm will be collected on the next iteration. */
180: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)Znext));
181: }
183: if (itcount) *itcount = it-1; /* Number of iterations actually completed. */
185: /*
186: Down here we have to solve for the "best" coefficients of the Krylov
187: columns, add the solution values together, and possibly unwind the
188: preconditioning from the solution
189: */
190: /* Form the solution (or the solution so far) */
191: KSPPGMRESBuildSoln(RS(0),ksp->vec_sol,ksp->vec_sol,ksp,it-2);
192: return(0);
193: }
195: /*
196: KSPSolve_PGMRES - This routine applies the PGMRES method.
198: Input Parameter:
199: . ksp - the Krylov space object that was set to use pgmres
201: Output Parameter:
202: . outits - number of iterations used
204: */
205: static PetscErrorCode KSPSolve_PGMRES(KSP ksp)
206: {
208: PetscInt its,itcount;
209: KSP_PGMRES *pgmres = (KSP_PGMRES*)ksp->data;
210: PetscBool guess_zero = ksp->guess_zero;
213: if (ksp->calc_sings && !pgmres->Rsvd) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ORDER,"Must call KSPSetComputeSingularValues() before KSPSetUp() is called");
214: PetscObjectSAWsTakeAccess((PetscObject)ksp);
215: ksp->its = 0;
216: PetscObjectSAWsGrantAccess((PetscObject)ksp);
218: itcount = 0;
219: ksp->reason = KSP_CONVERGED_ITERATING;
220: while (!ksp->reason) {
221: KSPInitialResidual(ksp,ksp->vec_sol,VEC_TEMP,VEC_TEMP_MATOP,VEC_VV(0),ksp->vec_rhs);
222: KSPPGMRESCycle(&its,ksp);
223: itcount += its;
224: if (itcount >= ksp->max_it) {
225: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
226: break;
227: }
228: ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
229: }
230: ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */
231: return(0);
232: }
234: static PetscErrorCode KSPDestroy_PGMRES(KSP ksp)
235: {
239: KSPDestroy_GMRES(ksp);
240: return(0);
241: }
243: /*
244: KSPPGMRESBuildSoln - create the solution from the starting vector and the
245: current iterates.
247: Input parameters:
248: nrs - work area of size it + 1.
249: vguess - index of initial guess
250: vdest - index of result. Note that vguess may == vdest (replace
251: guess with the solution).
252: it - HH upper triangular part is a block of size (it+1) x (it+1)
254: This is an internal routine that knows about the PGMRES internals.
255: */
256: static PetscErrorCode KSPPGMRESBuildSoln(PetscScalar *nrs,Vec vguess,Vec vdest,KSP ksp,PetscInt it)
257: {
258: PetscScalar tt;
260: PetscInt k,j;
261: KSP_PGMRES *pgmres = (KSP_PGMRES*)(ksp->data);
264: /* Solve for solution vector that minimizes the residual */
266: if (it < 0) { /* no pgmres steps have been performed */
267: VecCopy(vguess,vdest); /* VecCopy() is smart, exits immediately if vguess == vdest */
268: return(0);
269: }
271: /* solve the upper triangular system - RS is the right side and HH is
272: the upper triangular matrix - put soln in nrs */
273: if (*HH(it,it) != 0.0) nrs[it] = *RS(it) / *HH(it,it);
274: else nrs[it] = 0.0;
276: for (k=it-1; k>=0; k--) {
277: tt = *RS(k);
278: for (j=k+1; j<=it; j++) tt -= *HH(k,j) * nrs[j];
279: nrs[k] = tt / *HH(k,k);
280: }
282: /* Accumulate the correction to the solution of the preconditioned problem in TEMP */
283: VecZeroEntries(VEC_TEMP);
284: VecMAXPY(VEC_TEMP,it+1,nrs,&VEC_VV(0));
285: KSPUnwindPreconditioner(ksp,VEC_TEMP,VEC_TEMP_MATOP);
286: /* add solution to previous solution */
287: if (vdest == vguess) {
288: VecAXPY(vdest,1.0,VEC_TEMP);
289: } else {
290: VecWAXPY(vdest,1.0,VEC_TEMP,vguess);
291: }
292: return(0);
293: }
295: /*
297: KSPPGMRESUpdateHessenberg - Do the scalar work for the orthogonalization.
298: Return new residual.
300: input parameters:
302: . ksp - Krylov space object
303: . it - plane rotations are applied to the (it+1)th column of the
304: modified hessenberg (i.e. HH(:,it))
305: . hapend - PETSC_FALSE not happy breakdown ending.
307: output parameters:
308: . res - the new residual
310: */
311: /*
312: . it - column of the Hessenberg that is complete, PGMRES is actually computing two columns ahead of this
313: */
314: static PetscErrorCode KSPPGMRESUpdateHessenberg(KSP ksp,PetscInt it,PetscBool *hapend,PetscReal *res)
315: {
316: PetscScalar *hh,*cc,*ss,*rs;
317: PetscInt j;
318: PetscReal hapbnd;
319: KSP_PGMRES *pgmres = (KSP_PGMRES*)(ksp->data);
323: hh = HH(0,it); /* pointer to beginning of column to update */
324: cc = CC(0); /* beginning of cosine rotations */
325: ss = SS(0); /* beginning of sine rotations */
326: rs = RS(0); /* right hand side of least squares system */
328: /* The Hessenberg matrix is now correct through column it, save that form for possible spectral analysis */
329: for (j=0; j<=it+1; j++) *HES(j,it) = hh[j];
331: /* check for the happy breakdown */
332: hapbnd = PetscMin(PetscAbsScalar(hh[it+1] / rs[it]),pgmres->haptol);
333: if (PetscAbsScalar(hh[it+1]) < hapbnd) {
334: PetscInfo4(ksp,"Detected happy breakdown, current hapbnd = %14.12e H(%D,%D) = %14.12e\n",(double)hapbnd,it+1,it,(double)PetscAbsScalar(*HH(it+1,it)));
335: *hapend = PETSC_TRUE;
336: }
338: /* Apply all the previously computed plane rotations to the new column
339: of the Hessenberg matrix */
340: /* Note: this uses the rotation [conj(c) s ; -s c], c= cos(theta), s= sin(theta),
341: and some refs have [c s ; -conj(s) c] (don't be confused!) */
343: for (j=0; j<it; j++) {
344: PetscScalar hhj = hh[j];
345: hh[j] = PetscConj(cc[j])*hhj + ss[j]*hh[j+1];
346: hh[j+1] = -ss[j] *hhj + cc[j]*hh[j+1];
347: }
349: /*
350: compute the new plane rotation, and apply it to:
351: 1) the right-hand-side of the Hessenberg system (RS)
352: note: it affects RS(it) and RS(it+1)
353: 2) the new column of the Hessenberg matrix
354: note: it affects HH(it,it) which is currently pointed to
355: by hh and HH(it+1, it) (*(hh+1))
356: thus obtaining the updated value of the residual...
357: */
359: /* compute new plane rotation */
361: if (!*hapend) {
362: PetscReal delta = PetscSqrtReal(PetscSqr(PetscAbsScalar(hh[it])) + PetscSqr(PetscAbsScalar(hh[it+1])));
363: if (delta == 0.0) {
364: ksp->reason = KSP_DIVERGED_NULL;
365: return(0);
366: }
368: cc[it] = hh[it] / delta; /* new cosine value */
369: ss[it] = hh[it+1] / delta; /* new sine value */
371: hh[it] = PetscConj(cc[it])*hh[it] + ss[it]*hh[it+1];
372: rs[it+1] = -ss[it]*rs[it];
373: rs[it] = PetscConj(cc[it])*rs[it];
374: *res = PetscAbsScalar(rs[it+1]);
375: } else { /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
376: another rotation matrix (so RH doesn't change). The new residual is
377: always the new sine term times the residual from last time (RS(it)),
378: but now the new sine rotation would be zero...so the residual should
379: be zero...so we will multiply "zero" by the last residual. This might
380: not be exactly what we want to do here -could just return "zero". */
382: *res = 0.0;
383: }
384: return(0);
385: }
387: /*
388: KSPBuildSolution_PGMRES
390: Input Parameter:
391: . ksp - the Krylov space object
392: . ptr-
394: Output Parameter:
395: . result - the solution
397: Note: this calls KSPPGMRESBuildSoln - the same function that KSPPGMRESCycle
398: calls directly.
400: */
401: PetscErrorCode KSPBuildSolution_PGMRES(KSP ksp,Vec ptr,Vec *result)
402: {
403: KSP_PGMRES *pgmres = (KSP_PGMRES*)ksp->data;
407: if (!ptr) {
408: if (!pgmres->sol_temp) {
409: VecDuplicate(ksp->vec_sol,&pgmres->sol_temp);
410: PetscLogObjectParent((PetscObject)ksp,(PetscObject)pgmres->sol_temp);
411: }
412: ptr = pgmres->sol_temp;
413: }
414: if (!pgmres->nrs) {
415: /* allocate the work area */
416: PetscMalloc1(pgmres->max_k,&pgmres->nrs);
417: PetscLogObjectMemory((PetscObject)ksp,pgmres->max_k*sizeof(PetscScalar));
418: }
420: KSPPGMRESBuildSoln(pgmres->nrs,ksp->vec_sol,ptr,ksp,pgmres->it);
421: if (result) *result = ptr;
422: return(0);
423: }
425: PetscErrorCode KSPSetFromOptions_PGMRES(PetscOptionItems *PetscOptionsObject,KSP ksp)
426: {
430: KSPSetFromOptions_GMRES(PetscOptionsObject,ksp);
431: PetscOptionsHead(PetscOptionsObject,"KSP pipelined GMRES Options");
432: PetscOptionsTail();
433: return(0);
434: }
436: PetscErrorCode KSPReset_PGMRES(KSP ksp)
437: {
441: KSPReset_GMRES(ksp);
442: return(0);
443: }
445: /*MC
446: KSPPGMRES - Implements the Pipelined Generalized Minimal Residual method.
448: Options Database Keys:
449: + -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against
450: . -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
451: . -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of
452: vectors are allocated as needed)
453: . -ksp_gmres_classicalgramschmidt - use classical (unmodified) Gram-Schmidt to orthogonalize against the Krylov space (fast) (the default)
454: . -ksp_gmres_modifiedgramschmidt - use modified Gram-Schmidt in the orthogonalization (more stable, but slower)
455: . -ksp_gmres_cgs_refinement_type <refine_never,refine_ifneeded,refine_always> - determine if iterative refinement is used to increase the
456: stability of the classical Gram-Schmidt orthogonalization.
457: - -ksp_gmres_krylov_monitor - plot the Krylov space generated
459: Level: beginner
461: Notes:
462: MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods.
463: See the FAQ on the PETSc website for details.
465: Reference:
466: Ghysels, Ashby, Meerbergen, Vanroose, Hiding global communication latencies in the GMRES algorithm on massively parallel machines, 2012.
468: Developer Notes:
469: This object is subclassed off of KSPGMRES
471: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGMRES, KSPLGMRES, KSPPIPECG, KSPPIPECR,
472: KSPGMRESSetRestart(), KSPGMRESSetHapTol(), KSPGMRESSetPreAllocateVectors(), KSPGMRESSetOrthogonalization(), KSPGMRESGetOrthogonalization(),
473: KSPGMRESClassicalGramSchmidtOrthogonalization(), KSPGMRESModifiedGramSchmidtOrthogonalization(),
474: KSPGMRESCGSRefinementType, KSPGMRESSetCGSRefinementType(), KSPGMRESGetCGSRefinementType(), KSPGMRESMonitorKrylov()
475: M*/
477: PETSC_EXTERN PetscErrorCode KSPCreate_PGMRES(KSP ksp)
478: {
479: KSP_PGMRES *pgmres;
483: PetscNewLog(ksp,&pgmres);
485: ksp->data = (void*)pgmres;
486: ksp->ops->buildsolution = KSPBuildSolution_PGMRES;
487: ksp->ops->setup = KSPSetUp_PGMRES;
488: ksp->ops->solve = KSPSolve_PGMRES;
489: ksp->ops->reset = KSPReset_PGMRES;
490: ksp->ops->destroy = KSPDestroy_PGMRES;
491: ksp->ops->view = KSPView_GMRES;
492: ksp->ops->setfromoptions = KSPSetFromOptions_PGMRES;
493: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
494: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES;
496: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
497: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
498: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
500: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetPreAllocateVectors_C",KSPGMRESSetPreAllocateVectors_GMRES);
501: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetOrthogonalization_C",KSPGMRESSetOrthogonalization_GMRES);
502: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetOrthogonalization_C",KSPGMRESGetOrthogonalization_GMRES);
503: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetRestart_C",KSPGMRESSetRestart_GMRES);
504: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetRestart_C",KSPGMRESGetRestart_GMRES);
505: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESSetCGSRefinementType_C",KSPGMRESSetCGSRefinementType_GMRES);
506: PetscObjectComposeFunction((PetscObject)ksp,"KSPGMRESGetCGSRefinementType_C",KSPGMRESGetCGSRefinementType_GMRES);
508: pgmres->nextra_vecs = 1;
509: pgmres->haptol = 1.0e-30;
510: pgmres->q_preallocate = 0;
511: pgmres->delta_allocate = PGMRES_DELTA_DIRECTIONS;
512: pgmres->orthog = KSPGMRESClassicalGramSchmidtOrthogonalization;
513: pgmres->nrs = NULL;
514: pgmres->sol_temp = NULL;
515: pgmres->max_k = PGMRES_DEFAULT_MAXK;
516: pgmres->Rsvd = NULL;
517: pgmres->orthogwork = NULL;
518: pgmres->cgstype = KSP_GMRES_CGS_REFINE_NEVER;
519: return(0);
520: }