Actual source code: pipegcr.c
1: /*
2: Contributed by Sascha M. Schnepp and Patrick Sanan
3: */
5: #include "petscsys.h"
6: #include <../src/ksp/ksp/impls/gcr/pipegcr/pipegcrimpl.h>
8: static PetscBool cited = PETSC_FALSE;
9: static const char citation[] =
10: "@article{SSM2016,\n"
11: " author = {P. Sanan and S.M. Schnepp and D.A. May},\n"
12: " title = {Pipelined, Flexible Krylov Subspace Methods},\n"
13: " journal = {SIAM Journal on Scientific Computing},\n"
14: " volume = {38},\n"
15: " number = {5},\n"
16: " pages = {C441-C470},\n"
17: " year = {2016},\n"
18: " doi = {10.1137/15M1049130},\n"
19: " URL = {http://dx.doi.org/10.1137/15M1049130},\n"
20: " eprint = {http://dx.doi.org/10.1137/15M1049130}\n"
21: "}\n";
23: #define KSPPIPEGCR_DEFAULT_MMAX 15
24: #define KSPPIPEGCR_DEFAULT_NPREALLOC 5
25: #define KSPPIPEGCR_DEFAULT_VECB 5
26: #define KSPPIPEGCR_DEFAULT_TRUNCSTRAT KSP_FCD_TRUNC_TYPE_NOTAY
27: #define KSPPIPEGCR_DEFAULT_UNROLL_W PETSC_TRUE
29: #include <petscksp.h>
31: static PetscErrorCode KSPAllocateVectors_PIPEGCR(KSP ksp, PetscInt nvecsneeded, PetscInt chunksize)
32: {
33: PetscInt i;
34: KSP_PIPEGCR *pipegcr;
35: PetscInt nnewvecs, nvecsprev;
37: pipegcr = (KSP_PIPEGCR*)ksp->data;
39: /* Allocate enough new vectors to add chunksize new vectors, reach nvecsneedtotal, or to reach mmax+1, whichever is smallest */
40: if (pipegcr->nvecs < PetscMin(pipegcr->mmax+1,nvecsneeded)) {
41: nvecsprev = pipegcr->nvecs;
42: nnewvecs = PetscMin(PetscMax(nvecsneeded-pipegcr->nvecs,chunksize),pipegcr->mmax+1-pipegcr->nvecs);
43: KSPCreateVecs(ksp,nnewvecs,&pipegcr->ppvecs[pipegcr->nchunks],0,NULL);
44: PetscLogObjectParents((PetscObject)ksp,nnewvecs,pipegcr->ppvecs[pipegcr->nchunks]);
45: KSPCreateVecs(ksp,nnewvecs,&pipegcr->psvecs[pipegcr->nchunks],0,NULL);
46: PetscLogObjectParents((PetscObject)ksp,nnewvecs,pipegcr->psvecs[pipegcr->nchunks]);
47: KSPCreateVecs(ksp,nnewvecs,&pipegcr->pqvecs[pipegcr->nchunks],0,NULL);
48: PetscLogObjectParents((PetscObject)ksp,nnewvecs,pipegcr->pqvecs[pipegcr->nchunks]);
49: if (pipegcr->unroll_w) {
50: KSPCreateVecs(ksp,nnewvecs,&pipegcr->ptvecs[pipegcr->nchunks],0,NULL);
51: PetscLogObjectParents((PetscObject)ksp,nnewvecs,pipegcr->ptvecs[pipegcr->nchunks]);
52: }
53: pipegcr->nvecs += nnewvecs;
54: for (i=0;i<nnewvecs;i++) {
55: pipegcr->qvecs[nvecsprev+i] = pipegcr->pqvecs[pipegcr->nchunks][i];
56: pipegcr->pvecs[nvecsprev+i] = pipegcr->ppvecs[pipegcr->nchunks][i];
57: pipegcr->svecs[nvecsprev+i] = pipegcr->psvecs[pipegcr->nchunks][i];
58: if (pipegcr->unroll_w) {
59: pipegcr->tvecs[nvecsprev+i] = pipegcr->ptvecs[pipegcr->nchunks][i];
60: }
61: }
62: pipegcr->chunksizes[pipegcr->nchunks] = nnewvecs;
63: pipegcr->nchunks++;
64: }
65: return 0;
66: }
68: static PetscErrorCode KSPSolve_PIPEGCR_cycle(KSP ksp)
69: {
70: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
71: Mat A, B;
72: Vec x,r,b,z,w,m,n,p,s,q,t,*redux;
73: PetscInt i,j,k,idx,kdx,mi;
74: PetscScalar alpha=0.0,gamma,*betas,*dots;
75: PetscReal rnorm=0.0, delta,*eta,*etas;
77: /* !!PS We have not checked these routines for use with complex numbers. The inner products
78: are likely not defined correctly for that case */
81: KSPGetOperators(ksp, &A, &B);
82: x = ksp->vec_sol;
83: b = ksp->vec_rhs;
84: r = ksp->work[0];
85: z = ksp->work[1];
86: w = ksp->work[2]; /* w = Az = AB(r) (pipelining intermediate) */
87: m = ksp->work[3]; /* m = B(w) = B(Az) = B(AB(r)) (pipelining intermediate) */
88: n = ksp->work[4]; /* n = AB(w) = AB(Az) = AB(AB(r)) (pipelining intermediate) */
89: p = pipegcr->pvecs[0];
90: s = pipegcr->svecs[0];
91: q = pipegcr->qvecs[0];
92: t = pipegcr->unroll_w ? pipegcr->tvecs[0] : NULL;
94: redux = pipegcr->redux;
95: dots = pipegcr->dots;
96: etas = pipegcr->etas;
97: betas = dots; /* dots takes the result of all dot products of which the betas are a subset */
99: /* cycle initial residual */
100: KSP_MatMult(ksp,A,x,r);
101: VecAYPX(r,-1.0,b); /* r <- b - Ax */
102: KSP_PCApply(ksp,r,z); /* z <- B(r) */
103: KSP_MatMult(ksp,A,z,w); /* w <- Az */
105: /* initialization of other variables and pipelining intermediates */
106: VecCopy(z,p);
107: KSP_MatMult(ksp,A,p,s);
109: /* overlap initial computation of delta, gamma */
110: redux[0] = w;
111: redux[1] = r;
112: VecMDotBegin(w,2,redux,dots); /* Start split reductions for gamma = (w,r), delta = (w,w) */
113: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)s)); /* perform asynchronous reduction */
114: KSP_PCApply(ksp,s,q); /* q = B(s) */
115: if (pipegcr->unroll_w) {
116: KSP_MatMult(ksp,A,q,t); /* t = Aq */
117: }
118: VecMDotEnd(w,2,redux,dots); /* Finish split reduction */
119: delta = PetscRealPart(dots[0]);
120: etas[0] = delta;
121: gamma = dots[1];
122: alpha = gamma/delta;
124: i = 0;
125: do {
126: PetscObjectSAWsTakeAccess((PetscObject)ksp);
127: ksp->its++;
128: PetscObjectSAWsGrantAccess((PetscObject)ksp);
130: /* update solution, residuals, .. */
131: VecAXPY(x,+alpha,p);
132: VecAXPY(r,-alpha,s);
133: VecAXPY(z,-alpha,q);
134: if (pipegcr->unroll_w) {
135: VecAXPY(w,-alpha,t);
136: } else {
137: KSP_MatMult(ksp,A,z,w);
138: }
140: /* Computations of current iteration done */
141: i++;
143: if (pipegcr->modifypc) {
144: (*pipegcr->modifypc)(ksp,ksp->its,ksp->rnorm,pipegcr->modifypc_ctx);
145: }
147: /* If needbe, allocate a new chunk of vectors */
148: KSPAllocateVectors_PIPEGCR(ksp,i+1,pipegcr->vecb);
150: /* Note that we wrap around and start clobbering old vectors */
151: idx = i % (pipegcr->mmax+1);
152: p = pipegcr->pvecs[idx];
153: s = pipegcr->svecs[idx];
154: q = pipegcr->qvecs[idx];
155: if (pipegcr->unroll_w) {
156: t = pipegcr->tvecs[idx];
157: }
158: eta = pipegcr->etas+idx;
160: /* number of old directions to orthogonalize against */
161: switch(pipegcr->truncstrat) {
162: case KSP_FCD_TRUNC_TYPE_STANDARD:
163: mi = pipegcr->mmax;
164: break;
165: case KSP_FCD_TRUNC_TYPE_NOTAY:
166: mi = ((i-1) % pipegcr->mmax)+1;
167: break;
168: default:
169: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unrecognized Truncation Strategy");
170: }
172: /* Pick old p,s,q,zeta in a way suitable for VecMDot */
173: for (k=PetscMax(0,i-mi),j=0;k<i;j++,k++) {
174: kdx = k % (pipegcr->mmax+1);
175: pipegcr->pold[j] = pipegcr->pvecs[kdx];
176: pipegcr->sold[j] = pipegcr->svecs[kdx];
177: pipegcr->qold[j] = pipegcr->qvecs[kdx];
178: if (pipegcr->unroll_w) {
179: pipegcr->told[j] = pipegcr->tvecs[kdx];
180: }
181: redux[j] = pipegcr->svecs[kdx];
182: }
183: /* If the above loop is not run redux contains only r and w => all beta_k = 0, only gamma, delta != 0 */
184: redux[j] = r;
185: redux[j+1] = w;
187: /* Dot products */
188: /* Start split reductions for beta_k = (w,s_k), gamma = (w,r), delta = (w,w) */
189: VecMDotBegin(w,j+2,redux,dots);
190: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)w));
192: /* B(w-r) + u stabilization */
193: VecWAXPY(n,-1.0,r,w); /* m = u + B(w-r): (a) ntmp = w-r */
194: KSP_PCApply(ksp,n,m); /* m = u + B(w-r): (b) mtmp = B(ntmp) = B(w-r) */
195: VecAXPY(m,1.0,z); /* m = u + B(w-r): (c) m = z + mtmp */
196: if (pipegcr->unroll_w) {
197: KSP_MatMult(ksp,A,m,n); /* n = Am */
198: }
200: /* Finish split reductions for beta_k = (w,s_k), gamma = (w,r), delta = (w,w) */
201: VecMDotEnd(w,j+2,redux,dots);
202: gamma = dots[j];
203: delta = PetscRealPart(dots[j+1]);
205: /* compute new residual norm.
206: this cannot be done before this point so that the natural norm
207: is available for free and the communication involved is overlapped */
208: switch (ksp->normtype) {
209: case KSP_NORM_PRECONDITIONED:
210: VecNorm(z,NORM_2,&rnorm); /* ||r|| <- sqrt(z'*z) */
211: break;
212: case KSP_NORM_UNPRECONDITIONED:
213: VecNorm(r,NORM_2,&rnorm); /* ||r|| <- sqrt(r'*r) */
214: break;
215: case KSP_NORM_NATURAL:
216: rnorm = PetscSqrtReal(PetscAbsScalar(gamma)); /* ||r|| <- sqrt(r,w) */
217: break;
218: case KSP_NORM_NONE:
219: rnorm = 0.0;
220: break;
221: default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
222: }
224: /* Check for convergence */
225: PetscObjectSAWsTakeAccess((PetscObject)ksp);
226: ksp->rnorm = rnorm;
227: PetscObjectSAWsGrantAccess((PetscObject)ksp);
228: KSPLogResidualHistory(ksp,rnorm);
229: KSPMonitor(ksp,ksp->its,rnorm);
230: (*ksp->converged)(ksp,ksp->its,rnorm,&ksp->reason,ksp->cnvP);
231: if (ksp->reason) return 0;
233: /* compute new eta and scale beta */
234: *eta = 0.;
235: for (k=PetscMax(0,i-mi),j=0;k<i;j++,k++) {
236: kdx = k % (pipegcr->mmax+1);
237: betas[j] /= -etas[kdx]; /* betak /= etak */
238: *eta -= ((PetscReal)(PetscAbsScalar(betas[j])*PetscAbsScalar(betas[j]))) * etas[kdx];
239: /* etaitmp = -betaik^2 * etak */
240: }
241: *eta += delta; /* etai = delta -betaik^2 * etak */
243: /* check breakdown of eta = (s,s) */
244: if (*eta < 0.) {
245: pipegcr->norm_breakdown = PETSC_TRUE;
246: PetscInfo(ksp,"Restart due to square root breakdown at it = \n",ksp->its);
247: break;
248: } else {
249: alpha= gamma/(*eta); /* alpha = gamma/etai */
250: }
252: /* project out stored search directions using classical G-S */
253: VecCopy(z,p);
254: VecCopy(w,s);
255: VecCopy(m,q);
256: if (pipegcr->unroll_w) {
257: VecCopy(n,t);
258: VecMAXPY(t,j,betas,pipegcr->told); /* ti <- n - sum_k beta_k t_k */
259: }
260: VecMAXPY(p,j,betas,pipegcr->pold); /* pi <- ui - sum_k beta_k p_k */
261: VecMAXPY(s,j,betas,pipegcr->sold); /* si <- wi - sum_k beta_k s_k */
262: VecMAXPY(q,j,betas,pipegcr->qold); /* qi <- m - sum_k beta_k q_k */
264: } while (ksp->its < ksp->max_it);
265: if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
266: return 0;
267: }
269: static PetscErrorCode KSPSolve_PIPEGCR(KSP ksp)
270: {
271: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
272: Mat A, B;
273: Vec x,b,r,z,w;
274: PetscScalar gamma;
275: PetscReal rnorm=0.0;
276: PetscBool issym;
278: PetscCitationsRegister(citation,&cited);
280: KSPGetOperators(ksp, &A, &B);
281: x = ksp->vec_sol;
282: b = ksp->vec_rhs;
283: r = ksp->work[0];
284: z = ksp->work[1];
285: w = ksp->work[2]; /* w = Az = AB(r) (pipelining intermediate) */
287: /* compute initial residual */
288: if (!ksp->guess_zero) {
289: KSP_MatMult(ksp,A,x,r);
290: VecAYPX(r,-1.0,b); /* r <- b - Ax */
291: } else {
292: VecCopy(b,r); /* r <- b */
293: }
295: /* initial residual norm */
296: KSP_PCApply(ksp,r,z); /* z <- B(r) */
297: KSP_MatMult(ksp,A,z,w); /* w <- Az */
298: VecDot(r,w,&gamma); /* gamma = (r,w) */
300: switch (ksp->normtype) {
301: case KSP_NORM_PRECONDITIONED:
302: VecNorm(z,NORM_2,&rnorm); /* ||r|| <- sqrt(z'*z) */
303: break;
304: case KSP_NORM_UNPRECONDITIONED:
305: VecNorm(r,NORM_2,&rnorm); /* ||r|| <- sqrt(r'*r) */
306: break;
307: case KSP_NORM_NATURAL:
308: rnorm = PetscSqrtReal(PetscAbsScalar(gamma)); /* ||r|| <- sqrt(r,w) */
309: break;
310: case KSP_NORM_NONE:
311: rnorm = 0.0;
312: break;
313: default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
314: }
316: /* Is A symmetric? */
317: PetscObjectTypeCompareAny((PetscObject)A,&issym,MATSBAIJ,MATSEQSBAIJ,MATMPISBAIJ,"");
318: if (!issym) {
319: PetscInfo(A,"Matrix type is not any of MATSBAIJ,MATSEQSBAIJ,MATMPISBAIJ. Is matrix A symmetric (as required by CR methods)?");
320: }
322: /* logging */
323: PetscObjectSAWsTakeAccess((PetscObject)ksp);
324: ksp->its = 0;
325: ksp->rnorm0 = rnorm;
326: PetscObjectSAWsGrantAccess((PetscObject)ksp);
327: KSPLogResidualHistory(ksp,ksp->rnorm0);
328: KSPMonitor(ksp,ksp->its,ksp->rnorm0);
329: (*ksp->converged)(ksp,ksp->its,ksp->rnorm0,&ksp->reason,ksp->cnvP);
330: if (ksp->reason) return 0;
332: do {
333: KSPSolve_PIPEGCR_cycle(ksp);
334: if (ksp->reason) return 0;
335: if (pipegcr->norm_breakdown) {
336: pipegcr->n_restarts++;
337: pipegcr->norm_breakdown = PETSC_FALSE;
338: }
339: } while (ksp->its < ksp->max_it);
341: if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
342: return 0;
343: }
345: static PetscErrorCode KSPView_PIPEGCR(KSP ksp, PetscViewer viewer)
346: {
347: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
348: PetscBool isascii,isstring;
349: const char *truncstr;
351: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII, &isascii);
352: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
354: if (pipegcr->truncstrat == KSP_FCD_TRUNC_TYPE_STANDARD) {
355: truncstr = "Using standard truncation strategy";
356: } else if (pipegcr->truncstrat == KSP_FCD_TRUNC_TYPE_NOTAY) {
357: truncstr = "Using Notay's truncation strategy";
358: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Undefined FCD truncation strategy");
360: if (isascii) {
361: PetscViewerASCIIPrintf(viewer," max previous directions = %D\n",pipegcr->mmax);
362: PetscViewerASCIIPrintf(viewer," preallocated %D directions\n",PetscMin(pipegcr->nprealloc,pipegcr->mmax+1));
363: PetscViewerASCIIPrintf(viewer," %s\n",truncstr);
364: PetscViewerASCIIPrintf(viewer," w unrolling = %D \n", pipegcr->unroll_w);
365: PetscViewerASCIIPrintf(viewer," restarts performed = %D \n", pipegcr->n_restarts);
366: } else if (isstring) {
367: PetscViewerStringSPrintf(viewer, "max previous directions = %D, preallocated %D directions, %s truncation strategy", pipegcr->mmax,pipegcr->nprealloc,truncstr);
368: }
369: return 0;
370: }
372: static PetscErrorCode KSPSetUp_PIPEGCR(KSP ksp)
373: {
374: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
375: Mat A;
376: PetscBool diagonalscale;
377: const PetscInt nworkstd = 5;
379: PCGetDiagonalScale(ksp->pc,&diagonalscale);
382: KSPGetOperators(ksp, &A, NULL);
384: /* Allocate "standard" work vectors */
385: KSPSetWorkVecs(ksp,nworkstd);
387: /* Allocated space for pointers to additional work vectors
388: note that mmax is the number of previous directions, so we add 1 for the current direction */
389: PetscMalloc6(pipegcr->mmax+1,&(pipegcr->pvecs),pipegcr->mmax+1,&(pipegcr->ppvecs),pipegcr->mmax+1,&(pipegcr->svecs), pipegcr->mmax+1,&(pipegcr->psvecs),pipegcr->mmax+1,&(pipegcr->qvecs),pipegcr->mmax+1,&(pipegcr->pqvecs));
390: if (pipegcr->unroll_w) {
391: PetscMalloc3(pipegcr->mmax+1,&(pipegcr->tvecs),pipegcr->mmax+1,&(pipegcr->ptvecs),pipegcr->mmax+2,&(pipegcr->told));
392: }
393: PetscMalloc4(pipegcr->mmax+2,&(pipegcr->pold),pipegcr->mmax+2,&(pipegcr->sold),pipegcr->mmax+2,&(pipegcr->qold),pipegcr->mmax+2,&(pipegcr->chunksizes));
394: PetscMalloc3(pipegcr->mmax+2,&(pipegcr->dots),pipegcr->mmax+1,&(pipegcr->etas),pipegcr->mmax+2,&(pipegcr->redux));
395: /* If the requested number of preallocated vectors is greater than mmax reduce nprealloc */
396: if (pipegcr->nprealloc > pipegcr->mmax+1) {
397: PetscInfo(NULL,"Requested nprealloc=%d is greater than m_max+1=%d. Resetting nprealloc = m_max+1.\n",pipegcr->nprealloc, pipegcr->mmax+1);
398: }
400: /* Preallocate additional work vectors */
401: KSPAllocateVectors_PIPEGCR(ksp,pipegcr->nprealloc,pipegcr->nprealloc);
403: PetscCall(PetscLogObjectMemory(
404: (PetscObject)ksp,
405: (pipegcr->mmax + 1) * 4 * sizeof(Vec*) + /* old dirs */
406: (pipegcr->mmax + 1) * 4 * sizeof(Vec**) + /* old pdirs */
407: (pipegcr->mmax + 1) * 4 * sizeof(Vec*) + /* p/s/qold/told */
408: (pipegcr->mmax + 1) * sizeof(PetscInt) + /* chunksizes */
409: (pipegcr->mmax + 2) * sizeof(Vec*) + /* redux */
410: (pipegcr->mmax + 2) * sizeof(PetscScalar) + /* dots */
411: (pipegcr->mmax + 1) * sizeof(PetscReal) /* etas */));
412: return 0;
413: }
415: static PetscErrorCode KSPReset_PIPEGCR(KSP ksp)
416: {
417: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
419: if (pipegcr->modifypc_destroy) {
420: (*pipegcr->modifypc_destroy)(pipegcr->modifypc_ctx);
421: }
422: return 0;
423: }
425: static PetscErrorCode KSPDestroy_PIPEGCR(KSP ksp)
426: {
427: PetscInt i;
428: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
430: VecDestroyVecs(ksp->nwork,&ksp->work); /* Destroy "standard" work vecs */
432: /* Destroy vectors for old directions and the arrays that manage pointers to them */
433: if (pipegcr->nvecs) {
434: for (i=0;i<pipegcr->nchunks;i++) {
435: VecDestroyVecs(pipegcr->chunksizes[i],&pipegcr->ppvecs[i]);
436: VecDestroyVecs(pipegcr->chunksizes[i],&pipegcr->psvecs[i]);
437: VecDestroyVecs(pipegcr->chunksizes[i],&pipegcr->pqvecs[i]);
438: if (pipegcr->unroll_w) {
439: VecDestroyVecs(pipegcr->chunksizes[i],&pipegcr->ptvecs[i]);
440: }
441: }
442: }
444: PetscFree6(pipegcr->pvecs,pipegcr->ppvecs,pipegcr->svecs,pipegcr->psvecs,pipegcr->qvecs,pipegcr->pqvecs);
445: PetscFree4(pipegcr->pold,pipegcr->sold,pipegcr->qold,pipegcr->chunksizes);
446: PetscFree3(pipegcr->dots,pipegcr->etas,pipegcr->redux);
447: if (pipegcr->unroll_w) {
448: PetscFree3(pipegcr->tvecs,pipegcr->ptvecs,pipegcr->told);
449: }
451: KSPReset_PIPEGCR(ksp);
452: KSPDestroyDefault(ksp);
453: return 0;
454: }
456: /*@
457: KSPPIPEGCRSetUnrollW - Set to PETSC_TRUE to use PIPEGCR with unrolling of the w vector
459: Logically Collective on ksp
461: Input Parameters:
462: + ksp - the Krylov space context
463: - unroll_w - use unrolling
465: Level: intermediate
467: Options Database:
468: . -ksp_pipegcr_unroll_w <bool> - use unrolling
470: .seealso: KSPPIPEGCR, KSPPIPEGCRSetTruncationType(), KSPPIPEGCRSetNprealloc(),KSPPIPEGCRGetUnrollW()
471: @*/
472: PetscErrorCode KSPPIPEGCRSetUnrollW(KSP ksp,PetscBool unroll_w)
473: {
474: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;
478: pipegcr->unroll_w=unroll_w;
479: return 0;
480: }
482: /*@
483: KSPPIPEGCRGetUnrollW - Get information on PIPEGCR unrolling the w vector
485: Logically Collective on ksp
487: Input Parameter:
488: . ksp - the Krylov space context
490: Output Parameter:
491: . unroll_w - PIPEGCR uses unrolling (bool)
493: Level: intermediate
495: .seealso: KSPPIPEGCR, KSPPIPEGCRGetTruncationType(), KSPPIPEGCRGetNprealloc(),KSPPIPEGCRSetUnrollW()
496: @*/
497: PetscErrorCode KSPPIPEGCRGetUnrollW(KSP ksp,PetscBool *unroll_w)
498: {
499: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;
502: *unroll_w=pipegcr->unroll_w;
503: return 0;
504: }
506: /*@
507: KSPPIPEGCRSetMmax - set the maximum number of previous directions PIPEGCR will store for orthogonalization
509: Note: mmax + 1 directions are stored (mmax previous ones along with a current one)
510: and whether all are used in each iteration also depends on the truncation strategy
511: (see KSPPIPEGCRSetTruncationType)
513: Logically Collective on ksp
515: Input Parameters:
516: + ksp - the Krylov space context
517: - mmax - the maximum number of previous directions to orthogonalize againt
519: Level: intermediate
521: Options Database:
522: . -ksp_pipegcr_mmax <N> - maximum number of previous directions
524: .seealso: KSPPIPEGCR, KSPPIPEGCRSetTruncationType(), KSPPIPEGCRSetNprealloc()
525: @*/
526: PetscErrorCode KSPPIPEGCRSetMmax(KSP ksp,PetscInt mmax)
527: {
528: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;
532: pipegcr->mmax=mmax;
533: return 0;
534: }
536: /*@
537: KSPPIPEGCRGetMmax - get the maximum number of previous directions PIPEGCR will store
539: Note: PIPEGCR stores mmax+1 directions at most (mmax previous ones, and one current one)
541: Not Collective
543: Input Parameter:
544: . ksp - the Krylov space context
546: Output Parameter:
547: . mmax - the maximum number of previous directions allowed for orthogonalization
549: Level: intermediate
551: .seealso: KSPPIPEGCR, KSPPIPEGCRGetTruncationType(), KSPPIPEGCRGetNprealloc(), KSPPIPEGCRSetMmax()
552: @*/
554: PetscErrorCode KSPPIPEGCRGetMmax(KSP ksp,PetscInt *mmax)
555: {
556: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;
559: *mmax=pipegcr->mmax;
560: return 0;
561: }
563: /*@
564: KSPPIPEGCRSetNprealloc - set the number of directions to preallocate with PIPEGCR
566: Logically Collective on ksp
568: Input Parameters:
569: + ksp - the Krylov space context
570: - nprealloc - the number of vectors to preallocate
572: Level: advanced
574: Options Database:
575: . -ksp_pipegcr_nprealloc <N> - number of vectors to preallocate
577: .seealso: KSPPIPEGCR, KSPPIPEGCRGetTruncationType(), KSPPIPEGCRGetNprealloc()
578: @*/
579: PetscErrorCode KSPPIPEGCRSetNprealloc(KSP ksp,PetscInt nprealloc)
580: {
581: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;
585: pipegcr->nprealloc = nprealloc;
586: return 0;
587: }
589: /*@
590: KSPPIPEGCRGetNprealloc - get the number of directions preallocate by PIPEGCR
592: Not Collective
594: Input Parameter:
595: . ksp - the Krylov space context
597: Output Parameter:
598: . nprealloc - the number of directions preallocated
600: Level: advanced
602: .seealso: KSPPIPEGCR, KSPPIPEGCRGetTruncationType(), KSPPIPEGCRSetNprealloc()
603: @*/
604: PetscErrorCode KSPPIPEGCRGetNprealloc(KSP ksp,PetscInt *nprealloc)
605: {
606: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;
609: *nprealloc = pipegcr->nprealloc;
610: return 0;
611: }
613: /*@
614: KSPPIPEGCRSetTruncationType - specify how many of its stored previous directions PIPEGCR uses during orthoganalization
616: Logically Collective on ksp
618: KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions
619: KSP_FCD_TRUNC_TYPE_NOTAY uses the last max(1,mod(i,mmax)) directions at iteration i=0,1,..
621: Input Parameters:
622: + ksp - the Krylov space context
623: - truncstrat - the choice of strategy
625: Level: intermediate
627: Options Database:
628: . -ksp_pipegcr_truncation_type <standard,notay> - which stored basis vectors to orthogonalize against
630: .seealso: KSPPIPEGCR, KSPPIPEGCRSetTruncationType, KSPPIPEGCRTruncationType, KSPFCDTruncationType
631: @*/
632: PetscErrorCode KSPPIPEGCRSetTruncationType(KSP ksp,KSPFCDTruncationType truncstrat)
633: {
634: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;
638: pipegcr->truncstrat=truncstrat;
639: return 0;
640: }
642: /*@
643: KSPPIPEGCRGetTruncationType - get the truncation strategy employed by PIPEGCR
645: Not Collective
647: KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions
648: KSP_FCD_TRUNC_TYPE_NOTAY uses the last max(1,mod(i,mmax)) directions at iteration i=0,1,..
650: Input Parameter:
651: . ksp - the Krylov space context
653: Output Parameter:
654: . truncstrat - the strategy type
656: Options Database:
657: . -ksp_pipegcr_truncation_type <standard,notay> - which stored basis vectors to orthogonalize against
659: Level: intermediate
661: .seealso: KSPPIPEGCR, KSPPIPEGCRSetTruncationType, KSPPIPEGCRTruncationType, KSPFCDTruncationType
662: @*/
663: PetscErrorCode KSPPIPEGCRGetTruncationType(KSP ksp,KSPFCDTruncationType *truncstrat)
664: {
665: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;
668: *truncstrat=pipegcr->truncstrat;
669: return 0;
670: }
672: static PetscErrorCode KSPSetFromOptions_PIPEGCR(PetscOptionItems *PetscOptionsObject,KSP ksp)
673: {
674: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
675: PetscInt mmax,nprealloc;
676: PetscBool flg;
678: PetscOptionsHead(PetscOptionsObject,"KSP PIPEGCR options");
679: PetscOptionsInt("-ksp_pipegcr_mmax","Number of search directions to storue","KSPPIPEGCRSetMmax",pipegcr->mmax,&mmax,&flg);
680: if (flg) KSPPIPEGCRSetMmax(ksp,mmax);
681: PetscOptionsInt("-ksp_pipegcr_nprealloc","Number of directions to preallocate","KSPPIPEGCRSetNprealloc",pipegcr->nprealloc,&nprealloc,&flg);
682: if (flg) KSPPIPEGCRSetNprealloc(ksp,nprealloc);
683: PetscOptionsEnum("-ksp_pipegcr_truncation_type","Truncation approach for directions","KSPFCGSetTruncationType",KSPFCDTruncationTypes,(PetscEnum)pipegcr->truncstrat,(PetscEnum*)&pipegcr->truncstrat,NULL);
684: PetscOptionsBool("-ksp_pipegcr_unroll_w","Use unrolling of w","KSPPIPEGCRSetUnrollW",pipegcr->unroll_w,&pipegcr->unroll_w,NULL);
685: PetscOptionsTail();
686: return 0;
687: }
690: typedef PetscErrorCode (*KSPPIPEGCRModifyPCFunction)(KSP,PetscInt,PetscReal,void*);
691: typedef PetscErrorCode (*KSPPIPEGCRDestroyFunction)(void*);
693: static PetscErrorCode KSPPIPEGCRSetModifyPC_PIPEGCR(KSP ksp,KSPPIPEGCRModifyPCFunction function,void *data,KSPPIPEGCRDestroyFunction destroy)
694: {
695: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
698: pipegcr->modifypc = function;
699: pipegcr->modifypc_destroy = destroy;
700: pipegcr->modifypc_ctx = data;
701: return 0;
702: }
704: /*@C
705: KSPPIPEGCRSetModifyPC - Sets the routine used by PIPEGCR to modify the preconditioner.
707: Logically Collective on ksp
709: Input Parameters:
710: + ksp - iterative context obtained from KSPCreate()
711: . function - user defined function to modify the preconditioner
712: . ctx - user provided context for the modify preconditioner function
713: - destroy - the function to use to destroy the user provided application context.
715: Calling Sequence of function:
716: PetscErrorCode function (KSP ksp, PetscInt n, PetscReal rnorm, void *ctx)
718: ksp - iterative context
719: n - the total number of PIPEGCR iterations that have occurred
720: rnorm - 2-norm residual value
721: ctx - the user provided application context
723: Level: intermediate
725: Notes:
726: The default modifypc routine is KSPPIPEGCRModifyPCNoChange()
728: .seealso: KSPPIPEGCRModifyPCNoChange()
730: @*/
731: PetscErrorCode KSPPIPEGCRSetModifyPC(KSP ksp,PetscErrorCode (*function)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*destroy)(void*))
732: {
733: PetscUseMethod(ksp,"KSPPIPEGCRSetModifyPC_C",(KSP,PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*)(void*)),(ksp,function,data,destroy));
734: return 0;
735: }
737: /*MC
738: KSPPIPEGCR - Implements a Pipelined Generalized Conjugate Residual method.
740: Options Database Keys:
741: + -ksp_pipegcr_mmax <N> - the max number of Krylov directions to orthogonalize against
742: . -ksp_pipegcr_unroll_w - unroll w at the storage cost of a maximum of (mmax+1) extra vectors with the benefit of better pipelining (default: PETSC_TRUE)
743: . -ksp_pipegcr_nprealloc <N> - the number of vectors to preallocated for storing Krylov directions. Once exhausted new directions are allocated blockwise (default: 5)
744: - -ksp_pipegcr_truncation_type <standard,notay> - which previous search directions to orthogonalize against
746: Notes:
747: The PIPEGCR Krylov method supports non-symmetric matrices and permits the use of a preconditioner
748: which may vary from one iteration to the next. Users can can define a method to vary the
749: preconditioner between iterates via KSPPIPEGCRSetModifyPC().
750: Restarts are solves with x0 not equal to zero. When a restart occurs, the initial starting
751: solution is given by the current estimate for x which was obtained by the last restart
752: iterations of the PIPEGCR algorithm.
753: The method implemented requires at most the storage of 4 x mmax + 5 vectors, roughly twice as much as GCR.
755: Only supports left preconditioning.
757: The natural "norm" for this method is (u,Au), where u is the preconditioned residual. This norm is available at no additional computational cost, as with standard CG. Choosing preconditioned or unpreconditioned norm types involves a blocking reduction which prevents any benefit from pipelining.
759: Reference:
760: P. Sanan, S.M. Schnepp, and D.A. May,
761: "Pipelined, Flexible Krylov Subspace Methods,"
762: SIAM Journal on Scientific Computing 2016 38:5, C441-C470,
763: DOI: 10.1137/15M1049130
765: Level: intermediate
767: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
768: KSPPIPEFGMRES, KSPPIPECG, KSPPIPECR, KSPPIPEFCG,KSPPIPEGCRSetTruncationType(),KSPPIPEGCRSetNprealloc(),KSPPIPEGCRSetUnrollW(),KSPPIPEGCRSetMmax()
770: M*/
771: PETSC_EXTERN PetscErrorCode KSPCreate_PIPEGCR(KSP ksp)
772: {
773: KSP_PIPEGCR *pipegcr;
775: PetscNewLog(ksp,&pipegcr);
776: pipegcr->mmax = KSPPIPEGCR_DEFAULT_MMAX;
777: pipegcr->nprealloc = KSPPIPEGCR_DEFAULT_NPREALLOC;
778: pipegcr->nvecs = 0;
779: pipegcr->vecb = KSPPIPEGCR_DEFAULT_VECB;
780: pipegcr->nchunks = 0;
781: pipegcr->truncstrat = KSPPIPEGCR_DEFAULT_TRUNCSTRAT;
782: pipegcr->n_restarts = 0;
783: pipegcr->unroll_w = KSPPIPEGCR_DEFAULT_UNROLL_W;
785: ksp->data = (void*)pipegcr;
787: /* natural norm is for free, precond+unprecond norm require non-overlapped reduction */
788: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
789: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,1);
790: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);
791: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
793: ksp->ops->setup = KSPSetUp_PIPEGCR;
794: ksp->ops->solve = KSPSolve_PIPEGCR;
795: ksp->ops->reset = KSPReset_PIPEGCR;
796: ksp->ops->destroy = KSPDestroy_PIPEGCR;
797: ksp->ops->view = KSPView_PIPEGCR;
798: ksp->ops->setfromoptions = KSPSetFromOptions_PIPEGCR;
799: ksp->ops->buildsolution = KSPBuildSolutionDefault;
800: ksp->ops->buildresidual = KSPBuildResidualDefault;
802: PetscObjectComposeFunction((PetscObject)ksp,"KSPPIPEGCRSetModifyPC_C",KSPPIPEGCRSetModifyPC_PIPEGCR);
803: return 0;
804: }