Actual source code: pipegcr.c
petsc-3.7.5 2017-01-01
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> /*I "petscksp.h" I*/
8: #define KSPPIPEGCR_DEFAULT_MMAX 15
9: #define KSPPIPEGCR_DEFAULT_NPREALLOC 5
10: #define KSPPIPEGCR_DEFAULT_VECB 5
11: #define KSPPIPEGCR_DEFAULT_TRUNCSTRAT KSP_FCD_TRUNC_TYPE_NOTAY
12: #define KSPPIPEGCR_DEFAULT_UNROLL_W PETSC_TRUE
14: #include <petscksp.h>
18: static PetscErrorCode KSPAllocateVectors_PIPEGCR(KSP ksp, PetscInt nvecsneeded, PetscInt chunksize)
19: {
20: PetscErrorCode ierr;
21: PetscInt i;
22: KSP_PIPEGCR *pipegcr;
23: PetscInt nnewvecs, nvecsprev;
26: pipegcr = (KSP_PIPEGCR*)ksp->data;
28: /* Allocate enough new vectors to add chunksize new vectors, reach nvecsneedtotal, or to reach mmax+1, whichever is smallest */
29: if(pipegcr->nvecs < PetscMin(pipegcr->mmax+1,nvecsneeded)){
30: nvecsprev = pipegcr->nvecs;
31: nnewvecs = PetscMin(PetscMax(nvecsneeded-pipegcr->nvecs,chunksize),pipegcr->mmax+1-pipegcr->nvecs);
32: KSPCreateVecs(ksp,nnewvecs,&pipegcr->ppvecs[pipegcr->nchunks],0,NULL);
33: PetscLogObjectParents((PetscObject)ksp,nnewvecs,pipegcr->ppvecs[pipegcr->nchunks]);
34: KSPCreateVecs(ksp,nnewvecs,&pipegcr->psvecs[pipegcr->nchunks],0,NULL);
35: PetscLogObjectParents((PetscObject)ksp,nnewvecs,pipegcr->psvecs[pipegcr->nchunks]);
36: KSPCreateVecs(ksp,nnewvecs,&pipegcr->pqvecs[pipegcr->nchunks],0,NULL);
37: PetscLogObjectParents((PetscObject)ksp,nnewvecs,pipegcr->pqvecs[pipegcr->nchunks]);
38: if (pipegcr->unroll_w) {
39: KSPCreateVecs(ksp,nnewvecs,&pipegcr->ptvecs[pipegcr->nchunks],0,NULL);
40: PetscLogObjectParents((PetscObject)ksp,nnewvecs,pipegcr->ptvecs[pipegcr->nchunks]);
41: }
42: pipegcr->nvecs += nnewvecs;
43: for(i=0;i<nnewvecs;i++){
44: pipegcr->qvecs[nvecsprev+i] = pipegcr->pqvecs[pipegcr->nchunks][i];
45: pipegcr->pvecs[nvecsprev+i] = pipegcr->ppvecs[pipegcr->nchunks][i];
46: pipegcr->svecs[nvecsprev+i] = pipegcr->psvecs[pipegcr->nchunks][i];
47: if (pipegcr->unroll_w) {
48: pipegcr->tvecs[nvecsprev+i] = pipegcr->ptvecs[pipegcr->nchunks][i];
49: }
50: }
51: pipegcr->chunksizes[pipegcr->nchunks] = nnewvecs;
52: pipegcr->nchunks++;
53: }
54: return(0);
55: }
59: PetscErrorCode KSPSolve_PIPEGCR_cycle(KSP ksp)
60: {
61: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
63: Mat A, B;
64: Vec x,r,b,z,w,m,n,p,s,q,t,*redux;
65: PetscInt i,j,k,idx,kdx,mi;
66: PetscScalar alpha=0.0,gamma,*betas,*dots;
67: PetscReal rnorm=0.0, delta,*eta,*etas;
71: /* !!PS We have not checked these routines for use with complex numbers. The inner products
72: are likely not defined correctly for that case */
73: #if (defined(PETSC_USE_COMPLEX) && !defined(PETSC_SKIP_COMPLEX))
74: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"PIPEGCR has not been implemented for use with complex scalars");
75: #endif
77: KSPGetOperators(ksp, &A, &B);
78: x = ksp->vec_sol;
79: b = ksp->vec_rhs;
80: r = ksp->work[0];
81: z = ksp->work[1];
82: w = ksp->work[2]; /* w = Az = AB(r) (pipelining intermediate) */
83: m = ksp->work[3]; /* m = B(w) = B(Az) = B(AB(r)) (pipelining intermediate) */
84: n = ksp->work[4]; /* n = AB(w) = AB(Az) = AB(AB(r)) (pipelining intermediate) */
85: p = pipegcr->pvecs[0];
86: s = pipegcr->svecs[0];
87: q = pipegcr->qvecs[0];
88: t = pipegcr->unroll_w ? pipegcr->tvecs[0] : NULL;
90: redux = pipegcr->redux;
91: dots = pipegcr->dots;
92: etas = pipegcr->etas;
93: betas = dots; /* dots takes the result of all dot products of which the betas are a subset */
95: /* cycle initial residual */
96: KSP_MatMult(ksp,A,x,r);
97: VecAYPX(r,-1.0,b); /* r <- b - Ax */
98: KSP_PCApply(ksp,r,z); /* z <- B(r) */
99: KSP_MatMult(ksp,A,z,w); /* w <- Az */
101: /* initialization of other variables and pipelining intermediates */
102: VecCopy(z,p);
103: KSP_MatMult(ksp,A,p,s);
105: /* overlap initial computation of delta, gamma */
106: redux[0] = w;
107: redux[1] = r;
108: VecMDotBegin(w,2,redux,dots); /* Start split reductions for gamma = (w,r), delta = (w,w) */
109: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)s)); /* perform asynchronous reduction */
110: KSP_PCApply(ksp,s,q); /* q = B(s) */
111: if (pipegcr->unroll_w) {
112: KSP_MatMult(ksp,A,q,t); /* t = Aq */
113: }
114: VecMDotEnd(w,2,redux,dots); /* Finish split reduction */
115: delta = PetscRealPart(dots[0]);
116: etas[0] = delta;
117: gamma = dots[1];
118: alpha = gamma/delta;
120: i = 0;
121: do {
122: PetscObjectSAWsTakeAccess((PetscObject)ksp);
123: ksp->its++;
124: PetscObjectSAWsGrantAccess((PetscObject)ksp);
126: /* update solution, residuals, .. */
127: VecAXPY(x,+alpha,p);
128: VecAXPY(r,-alpha,s);
129: VecAXPY(z,-alpha,q);
130: if(pipegcr->unroll_w){
131: VecAXPY(w,-alpha,t);
132: } else {
133: KSP_MatMult(ksp,A,z,w);
134: }
136: /* Computations of current iteration done */
137: i++;
139: if (pipegcr->modifypc) {
140: (*pipegcr->modifypc)(ksp,ksp->its,ksp->rnorm,pipegcr->modifypc_ctx);
141: }
143: /* If needbe, allocate a new chunk of vectors */
144: KSPAllocateVectors_PIPEGCR(ksp,i+1,pipegcr->vecb);
146: /* Note that we wrap around and start clobbering old vectors */
147: idx = i % (pipegcr->mmax+1);
148: p = pipegcr->pvecs[idx];
149: s = pipegcr->svecs[idx];
150: q = pipegcr->qvecs[idx];
151: if (pipegcr->unroll_w) {
152: t = pipegcr->tvecs[idx];
153: }
154: eta = pipegcr->etas+idx;
156: /* number of old directions to orthogonalize against */
157: switch(pipegcr->truncstrat){
158: case KSP_FCD_TRUNC_TYPE_STANDARD:
159: mi = pipegcr->mmax;
160: break;
161: case KSP_FCD_TRUNC_TYPE_NOTAY:
162: mi = ((i-1) % pipegcr->mmax)+1;
163: break;
164: default:
165: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unrecognized Truncation Strategy");
166: }
168: /* Pick old p,s,q,zeta in a way suitable for VecMDot */
169: for(k=PetscMax(0,i-mi),j=0;k<i;j++,k++){
170: kdx = k % (pipegcr->mmax+1);
171: pipegcr->pold[j] = pipegcr->pvecs[kdx];
172: pipegcr->sold[j] = pipegcr->svecs[kdx];
173: pipegcr->qold[j] = pipegcr->qvecs[kdx];
174: if (pipegcr->unroll_w) {
175: pipegcr->told[j] = pipegcr->tvecs[kdx];
176: }
177: redux[j] = pipegcr->svecs[kdx];
178: }
179: /* If the above loop is not run redux contains only r and w => all beta_k = 0, only gamma, delta != 0 */
180: redux[j] = r;
181: redux[j+1] = w;
183: /* Dot products */
184: /* Start split reductions for beta_k = (w,s_k), gamma = (w,r), delta = (w,w) */
185: VecMDotBegin(w,j+2,redux,dots);
186: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)w));
188: /* B(w-r) + u stabilization */
189: VecWAXPY(n,-1.0,r,w); /* m = u + B(w-r): (a) ntmp = w-r */
190: KSP_PCApply(ksp,n,m); /* m = u + B(w-r): (b) mtmp = B(ntmp) = B(w-r) */
191: VecAXPY(m,1.0,z); /* m = u + B(w-r): (c) m = z + mtmp */
192: if(pipegcr->unroll_w){
193: KSP_MatMult(ksp,A,m,n); /* n = Am */
194: }
196: /* Finish split reductions for beta_k = (w,s_k), gamma = (w,r), delta = (w,w) */
197: VecMDotEnd(w,j+2,redux,dots);
198: gamma = dots[j];
199: delta = PetscRealPart(dots[j+1]);
201: /* compute new residual norm.
202: this cannot be done before this point so that the natural norm
203: is available for free and the communication involved is overlapped */
204: switch (ksp->normtype) {
205: case KSP_NORM_PRECONDITIONED:
206: VecNorm(z,NORM_2,&rnorm); /* ||r|| <- sqrt(z'*z) */
207: break;
208: case KSP_NORM_UNPRECONDITIONED:
209: VecNorm(r,NORM_2,&rnorm); /* ||r|| <- sqrt(r'*r) */
210: break;
211: case KSP_NORM_NATURAL:
212: rnorm = PetscSqrtReal(PetscAbsScalar(gamma)); /* ||r|| <- sqrt(r,w) */
213: break;
214: case KSP_NORM_NONE:
215: rnorm = 0.0;
216: break;
217: default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
218: }
220: /* Check for convergence */
221: PetscObjectSAWsTakeAccess((PetscObject)ksp);
222: ksp->rnorm = rnorm;
223: PetscObjectSAWsGrantAccess((PetscObject)ksp);
224: KSPLogResidualHistory(ksp,rnorm);
225: KSPMonitor(ksp,ksp->its,rnorm);
226: (*ksp->converged)(ksp,ksp->its+1,rnorm,&ksp->reason,ksp->cnvP);
227: if (ksp->reason) break;
229: /* compute new eta and scale beta */
230: *eta = 0.;
231: for(k=PetscMax(0,i-mi),j=0;k<i;j++,k++){
232: kdx = k % (pipegcr->mmax+1);
233: betas[j] /= -etas[kdx]; /* betak /= etak */
234: *eta -= ((PetscReal)(PetscAbsScalar(betas[j])*PetscAbsScalar(betas[j]))) * etas[kdx];
235: /* etaitmp = -betaik^2 * etak */
236: }
237: *eta += delta; /* etai = delta -betaik^2 * etak */
239: /* check breakdown of eta = (s,s) */
240: if(*eta < 0.) {
241: pipegcr->norm_breakdown = PETSC_TRUE;
242: PetscInfo1(ksp,"Restart due to square root breakdown at it = \n",ksp->its);
243: break;
244: } else {
245: alpha= gamma/(*eta); /* alpha = gamma/etai */
246: }
248: /* project out stored search directions using classical G-S */
249: VecCopy(z,p);
250: VecCopy(w,s);
251: VecCopy(m,q);
252: if(pipegcr->unroll_w){
253: VecCopy(n,t);
254: VecMAXPY(t,j,betas,pipegcr->told); /* ti <- n - sum_k beta_k t_k */
255: }
256: VecMAXPY(p,j,betas,pipegcr->pold); /* pi <- ui - sum_k beta_k p_k */
257: VecMAXPY(s,j,betas,pipegcr->sold); /* si <- wi - sum_k beta_k s_k */
258: VecMAXPY(q,j,betas,pipegcr->qold); /* qi <- m - sum_k beta_k q_k */
260: } while (ksp->its < ksp->max_it);
261: return(0);
262: }
266: PetscErrorCode KSPSolve_PIPEGCR(KSP ksp)
267: {
268: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
270: Mat A, B;
271: Vec x,b,r,z,w;
272: PetscScalar gamma;
273: PetscReal rnorm=0.0;
274: PetscBool issym;
277: KSPGetOperators(ksp, &A, &B);
278: x = ksp->vec_sol;
279: b = ksp->vec_rhs;
280: r = ksp->work[0];
281: z = ksp->work[1];
282: w = ksp->work[2]; /* w = Az = AB(r) (pipelining intermediate) */
284: /* compute initial residual */
285: if (!ksp->guess_zero) {
286: KSP_MatMult(ksp,A,x,r);
287: VecAYPX(r,-1.0,b); /* r <- b - Ax */
288: } else {
289: VecCopy(b,r); /* r <- b */
290: }
292: /* initial residual norm */
293: KSP_PCApply(ksp,r,z); /* z <- B(r) */
294: KSP_MatMult(ksp,A,z,w); /* w <- Az */
295: VecDot(r,w,&gamma); /* gamma = (r,w) */
297: switch (ksp->normtype) {
298: case KSP_NORM_PRECONDITIONED:
299: VecNorm(z,NORM_2,&rnorm); /* ||r|| <- sqrt(z'*z) */
300: break;
301: case KSP_NORM_UNPRECONDITIONED:
302: VecNorm(r,NORM_2,&rnorm); /* ||r|| <- sqrt(r'*r) */
303: break;
304: case KSP_NORM_NATURAL:
305: rnorm = PetscSqrtReal(PetscAbsScalar(gamma)); /* ||r|| <- sqrt(r,w) */
306: break;
307: case KSP_NORM_NONE:
308: rnorm = 0.0;
309: break;
310: default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
311: }
313: /* Is A symmetric? */
314: PetscObjectTypeCompareAny((PetscObject)A,&issym,"sbaij","seqsbaij","mpibaij","");
315: if (!issym) {
316: PetscInfo(A,"Matrix type is not any of MATSBAIJ,MATSEQSBAIJ,MATMPIBAIJ. Is matrix A symmetric (as required by CR methods)?");
317: }
319: /* logging */
320: PetscObjectSAWsTakeAccess((PetscObject)ksp);
321: ksp->its = 0;
322: ksp->rnorm0 = rnorm;
323: PetscObjectSAWsGrantAccess((PetscObject)ksp);
324: KSPLogResidualHistory(ksp,ksp->rnorm0);
325: KSPMonitor(ksp,ksp->its,ksp->rnorm0);
326: (*ksp->converged)(ksp,ksp->its,ksp->rnorm0,&ksp->reason,ksp->cnvP);
327: if (ksp->reason) return(0);
329: do {
330: KSPSolve_PIPEGCR_cycle(ksp);
331: if (ksp->reason) break;
332: if (pipegcr->norm_breakdown) {
333: pipegcr->n_restarts++;
334: pipegcr->norm_breakdown = PETSC_FALSE;
335: }
336: } while (ksp->its < ksp->max_it);
338: if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
339: return(0);
340: }
344: PetscErrorCode KSPView_PIPEGCR(KSP ksp, PetscViewer viewer)
345: {
346: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
348: PetscBool isascii,isstring;
349: const char *truncstr;
352: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII, &isascii);
353: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
355: if(pipegcr->truncstrat == KSP_FCD_TRUNC_TYPE_STANDARD){
356: truncstr = "Using standard truncation strategy";
357: } else if(pipegcr->truncstrat == KSP_FCD_TRUNC_TYPE_NOTAY){
358: truncstr = "Using Notay's truncation strategy";
359: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Undefined FCD truncation strategy");
360:
362: if (isascii) {
363: PetscViewerASCIIPrintf(viewer," PIPEGCR: max previous directions = %D\n",pipegcr->mmax);
364: PetscViewerASCIIPrintf(viewer," PIPEGCR: preallocated %D directions\n",PetscMin(pipegcr->nprealloc,pipegcr->mmax+1));
365: PetscViewerASCIIPrintf(viewer," PIPEGCR: %s\n",truncstr);
366: PetscViewerASCIIPrintf(viewer," PIPEGCR: w unrolling = %D \n", pipegcr->unroll_w);
367: PetscViewerASCIIPrintf(viewer," PIPEGCR: restarts performed = %D \n", pipegcr->n_restarts);
368: } else if (isstring) {
369: PetscViewerStringSPrintf(viewer, "max previous directions = %D, preallocated %D directions, %s truncation strategy", pipegcr->mmax,pipegcr->nprealloc,truncstr);
370: }
371: return(0);
372: }
377: PetscErrorCode KSPSetUp_PIPEGCR(KSP ksp)
378: {
379: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
381: Mat A;
382: PetscBool diagonalscale;
383: const PetscInt nworkstd = 5;
386: PCGetDiagonalScale(ksp->pc,&diagonalscale);
387: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
389: KSPGetOperators(ksp, &A, NULL);
391: /* Allocate "standard" work vectors */
392: KSPSetWorkVecs(ksp,nworkstd);
394: /* Allocated space for pointers to additional work vectors
395: note that mmax is the number of previous directions, so we add 1 for the current direction */
396: 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));
397: if (pipegcr->unroll_w) {
398: PetscMalloc3(pipegcr->mmax+1,&(pipegcr->tvecs),pipegcr->mmax+1,&(pipegcr->ptvecs),pipegcr->mmax+2,&(pipegcr->told));
399: }
400: PetscMalloc4(pipegcr->mmax+2,&(pipegcr->pold),pipegcr->mmax+2,&(pipegcr->sold),pipegcr->mmax+2,&(pipegcr->qold),pipegcr->mmax+2,&(pipegcr->chunksizes));
401: PetscMalloc3(pipegcr->mmax+2,&(pipegcr->dots),pipegcr->mmax+1,&(pipegcr->etas),pipegcr->mmax+2,&(pipegcr->redux));
402: /* If the requested number of preallocated vectors is greater than mmax reduce nprealloc */
403: if(pipegcr->nprealloc > pipegcr->mmax+1){
404: PetscInfo2(NULL,"Requested nprealloc=%d is greater than m_max+1=%d. Resetting nprealloc = m_max+1.\n",pipegcr->nprealloc, pipegcr->mmax+1);
405: }
407: /* Preallocate additional work vectors */
408: KSPAllocateVectors_PIPEGCR(ksp,pipegcr->nprealloc,pipegcr->nprealloc);
410: PetscLogObjectMemory(
411: (PetscObject)ksp,
412: (pipegcr->mmax + 1) * 4 * sizeof(Vec*) + /* old dirs */
413: (pipegcr->mmax + 1) * 4 * sizeof(Vec**) + /* old pdirs */
414: (pipegcr->mmax + 1) * 4 * sizeof(Vec*) + /* p/s/qold/told */
415: (pipegcr->mmax + 1) * sizeof(PetscInt) + /* chunksizes */
416: (pipegcr->mmax + 2) * sizeof(Vec*) + /* redux */
417: (pipegcr->mmax + 2) * sizeof(PetscScalar) + /* dots */
418: (pipegcr->mmax + 1) * sizeof(PetscReal) /* etas */
419: );
420: return(0);
421: }
425: PetscErrorCode KSPReset_PIPEGCR(KSP ksp)
426: {
428: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
431: if (pipegcr->modifypc_destroy) {
432: (*pipegcr->modifypc_destroy)(pipegcr->modifypc_ctx);
433: }
434: return(0);
435: }
439: PetscErrorCode KSPDestroy_PIPEGCR(KSP ksp)
440: {
442: PetscInt i;
443: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
446: VecDestroyVecs(ksp->nwork,&ksp->work); /* Destroy "standard" work vecs */
448: /* Destroy vectors for old directions and the arrays that manage pointers to them */
449: if(pipegcr->nvecs){
450: for(i=0;i<pipegcr->nchunks;i++){
451: VecDestroyVecs(pipegcr->chunksizes[i],&pipegcr->ppvecs[i]);
452: VecDestroyVecs(pipegcr->chunksizes[i],&pipegcr->psvecs[i]);
453: VecDestroyVecs(pipegcr->chunksizes[i],&pipegcr->pqvecs[i]);
454: if (pipegcr->unroll_w) {
455: VecDestroyVecs(pipegcr->chunksizes[i],&pipegcr->ptvecs[i]);
456: }
457: }
458: }
460: PetscFree6(pipegcr->pvecs,pipegcr->ppvecs,pipegcr->svecs,pipegcr->psvecs,pipegcr->qvecs,pipegcr->pqvecs);
461: PetscFree4(pipegcr->pold,pipegcr->sold,pipegcr->qold,pipegcr->chunksizes);
462: PetscFree3(pipegcr->dots,pipegcr->etas,pipegcr->redux);
463: if (pipegcr->unroll_w) {
464: PetscFree3(pipegcr->tvecs,pipegcr->ptvecs,pipegcr->told);
465: }
467: KSPReset_PIPEGCR(ksp);
468: KSPDestroyDefault(ksp);
469: return(0);
470: }
474: /*@
475: KSPPIPEGCRSetUnrollW - Set to PETSC_TRUE to use PIPEGCR with unrolling of the w vector
477: Logically Collective on KSP
479: Input Parameters:
480: + ksp - the Krylov space context
481: - unroll_w - use unrolling
483: Level: intermediate
485: Options Database:
486: . -ksp_pipegcr_unroll_w
488: .seealso: KSPPIPEGCR, KSPPIPEGCRSetTruncationType(), KSPPIPEGCRSetNprealloc(),KSPPIPEGCRGetUnrollW()
489: @*/
490: PetscErrorCode KSPPIPEGCRSetUnrollW(KSP ksp,PetscBool unroll_w)
491: {
492: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;;
497: pipegcr->unroll_w=unroll_w;
498: return(0);
499: }
503: /*@
504: KSPPIPEGCRGetUnrollW - Get information on PIPEGCR unrolling the w vector
506: Logically Collective on KSP
508: Input Parameter:
509: . ksp - the Krylov space context
511: Output Parameter:
512: . unroll_w - PIPEGCR uses unrolling (bool)
514: Level: intermediate
516: Options Database:
517: . -ksp_pipegcr_unroll_w
519: .seealso: KSPPIPEGCR, KSPPIPEGCRGetTruncationType(), KSPPIPEGCRGetNprealloc(),KSPPIPEGCRSetUnrollW()
520: @*/
521: PetscErrorCode KSPPIPEGCRGetUnrollW(KSP ksp,PetscBool *unroll_w)
522: {
523: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;;
527: *unroll_w=pipegcr->unroll_w;
528: return(0);
529: }
533: /*@
534: KSPPIPEGCRSetMmax - set the maximum number of previous directions PIPEGCR will store for orthogonalization
536: Note: mmax + 1 directions are stored (mmax previous ones along with a current one)
537: and whether all are used in each iteration also depends on the truncation strategy
538: (see KSPPIPEGCRSetTruncationType)
540: Logically Collective on KSP
542: Input Parameters:
543: + ksp - the Krylov space context
544: - mmax - the maximum number of previous directions to orthogonalize againt
546: Level: intermediate
548: Options Database:
549: . -ksp_pipegcr_mmax <N>
551: .seealso: KSPPIPEGCR, KSPPIPEGCRSetTruncationType(), KSPPIPEGCRSetNprealloc()
552: @*/
553: PetscErrorCode KSPPIPEGCRSetMmax(KSP ksp,PetscInt mmax)
554: {
555: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;;
560: pipegcr->mmax=mmax;
561: return(0);
562: }
566: /*@
567: KSPPIPEGCRGetMmax - get the maximum number of previous directions PIPEGCR will store
569: Note: PIPEGCR stores mmax+1 directions at most (mmax previous ones, and one current one)
571: Not Collective
573: Input Parameter:
574: . ksp - the Krylov space context
576: Output Parameter:
577: . mmax - the maximum number of previous directons allowed for orthogonalization
579: Options Database:
580: . -ksp_pipegcr_mmax <N>
582: Level: intermediate
584: .keywords: KSP, PIPEGCR, truncation
586: .seealso: KSPPIPEGCR, KSPPIPEGCRGetTruncationType(), KSPPIPEGCRGetNprealloc(), KSPPIPEGCRSetMmax()
587: @*/
589: PetscErrorCode KSPPIPEGCRGetMmax(KSP ksp,PetscInt *mmax)
590: {
591: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;;
595: *mmax=pipegcr->mmax;
596: return(0);
597: }
601: /*@
602: KSPPIPEGCRSetNprealloc - set the number of directions to preallocate with PIPEGCR
604: Logically Collective on KSP
606: Input Parameters:
607: + ksp - the Krylov space context
608: - nprealloc - the number of vectors to preallocate
610: Level: advanced
612: Options Database:
613: . -ksp_pipegcr_nprealloc <N>
615: .seealso: KSPPIPEGCR, KSPPIPEGCRGetTruncationType(), KSPPIPEGCRGetNprealloc()
616: @*/
617: PetscErrorCode KSPPIPEGCRSetNprealloc(KSP ksp,PetscInt nprealloc)
618: {
619: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;;
624: pipegcr->nprealloc = nprealloc;
625: return(0);
626: }
630: /*@
631: KSPPIPEGCRGetNprealloc - get the number of directions preallocate by PIPEGCR
633: Not Collective
635: Input Parameter:
636: . ksp - the Krylov space context
638: Output Parameter:
639: . nprealloc - the number of directions preallocated
641: Options Database:
642: . -ksp_pipegcr_nprealloc <N>
644: Level: advanced
646: .keywords: KSP, PIPEGCR, truncation
648: .seealso: KSPPIPEGCR, KSPPIPEGCRGetTruncationType(), KSPPIPEGCRSetNprealloc()
649: @*/
650: PetscErrorCode KSPPIPEGCRGetNprealloc(KSP ksp,PetscInt *nprealloc)
651: {
652: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;;
656: *nprealloc = pipegcr->nprealloc;
657: return(0);
658: }
662: /*@
663: KSPPIPEGCRSetTruncationType - specify how many of its stored previous directions PIPEGCR uses during orthoganalization
665: Logically Collective on KSP
667: KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions
668: KSP_FCD_TRUNC_TYPE_NOTAY uses the last max(1,mod(i,mmax)) directions at iteration i=0,1,..
670: Input Parameters:
671: + ksp - the Krylov space context
672: - truncstrat - the choice of strategy
674: Level: intermediate
676: Options Database:
677: . -ksp_pipegcr_truncation, -ksp_pipegcr_truncation_restart
679: .seealso: KSPPIPEGCR, KSPPIPEGCRSetTruncationType, KSPPIPEGCRTruncationType, KSPFCDTruncationType
680: @*/
681: PetscErrorCode KSPPIPEGCRSetTruncationType(KSP ksp,KSPFCDTruncationType truncstrat)
682: {
683: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;;
688: pipegcr->truncstrat=truncstrat;
689: return(0);
690: }
694: /*@
695: KSPPIPEGCRGetTruncationType - get the truncation strategy employed by PIPEGCR
697: Not Collective
699: KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions
700: KSP_FCD_TRUNC_TYPE_NOTAY uses the last max(1,mod(i,mmax)) directions at iteration i=0,1,..
702: Input Parameter:
703: . ksp - the Krylov space context
705: Output Parameter:
706: . truncstrat - the strategy type
708: Options Database:
709: . -ksp_pipegcr_truncation, -ksp_pipegcr_truncation_restart
711: Level: intermediate
713: .keywords: KSP, PIPEGCR, truncation
715: .seealso: KSPPIPEGCR, KSPPIPEGCRSetTruncationType, KSPPIPEGCRTruncationType, KSPFCDTruncationType
716: @*/
717: PetscErrorCode KSPPIPEGCRGetTruncationType(KSP ksp,KSPFCDTruncationType *truncstrat)
718: {
719: KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;;
723: *truncstrat=pipegcr->truncstrat;
724: return(0);
725: }
729: PetscErrorCode KSPSetFromOptions_PIPEGCR(PetscOptionItems *PetscOptionsObject,KSP ksp)
730: {
732: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
733: PetscInt mmax,nprealloc;
734: PetscBool flg;
737: PetscOptionsHead(PetscOptionsObject,"KSP PIPEGCR options");
738: PetscOptionsInt("-ksp_pipegcr_mmax","Number of search directions to storue","KSPPIPEGCRSetMmax",pipegcr->mmax,&mmax,&flg);
739: if (flg) KSPPIPEGCRSetMmax(ksp,mmax);
740: PetscOptionsInt("-ksp_pipegcr_nprealloc","Number of directions to preallocate","KSPPIPEGCRSetNprealloc",pipegcr->nprealloc,&nprealloc,&flg);
741: if (flg) { KSPPIPEGCRSetNprealloc(ksp,nprealloc); }
742: PetscOptionsEnum("-ksp_pipegcr_truncation_type","Truncation approach for directions","KSPFCGSetTruncationType",KSPFCDTruncationTypes,(PetscEnum)pipegcr->truncstrat,(PetscEnum*)&pipegcr->truncstrat,NULL);
743: PetscOptionsBool("-ksp_pipegcr_unroll_w","Use unrolling of w","KSPPIPEGCRSetUnrollW",pipegcr->unroll_w,&pipegcr->unroll_w,NULL);
744: PetscOptionsTail();
745: return(0);
746: }
749: typedef PetscErrorCode (*KSPPIPEGCRModifyPCFunction)(KSP,PetscInt,PetscReal,void*);
750: typedef PetscErrorCode (*KSPPIPEGCRDestroyFunction)(void*);
754: static PetscErrorCode KSPPIPEGCRSetModifyPC_PIPEGCR(KSP ksp,KSPPIPEGCRModifyPCFunction function,void *data,KSPPIPEGCRDestroyFunction destroy)
755: {
756: KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;
760: pipegcr->modifypc = function;
761: pipegcr->modifypc_destroy = destroy;
762: pipegcr->modifypc_ctx = data;
763: return(0);
764: }
768: /*@C
769: KSPPIPEGCRSetModifyPC - Sets the routine used by PIPEGCR to modify the preconditioner.
771: Logically Collective on KSP
773: Input Parameters:
774: + ksp - iterative context obtained from KSPCreate()
775: . function - user defined function to modify the preconditioner
776: . ctx - user provided contex for the modify preconditioner function
777: - destroy - the function to use to destroy the user provided application context.
779: Calling Sequence of function:
780: PetscErrorCode function (KSP ksp, PetscInt n, PetscReal rnorm, void *ctx)
782: ksp - iterative context
783: n - the total number of PIPEGCR iterations that have occurred
784: rnorm - 2-norm residual value
785: ctx - the user provided application context
787: Level: intermediate
789: Notes:
790: The default modifypc routine is KSPPIPEGCRModifyPCNoChange()
792: .seealso: KSPPIPEGCRModifyPCNoChange()
794: @*/
795: PetscErrorCode KSPPIPEGCRSetModifyPC(KSP ksp,PetscErrorCode (*function)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*destroy)(void*))
796: {
800: PetscUseMethod(ksp,"KSPPIPEGCRSetModifyPC_C",(KSP,PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*)(void*)),(ksp,function,data,destroy));
801: return(0);
802: }
804: /*MC
805: KSPPIPEGCR - Implements the preconditioned Generalized Conjugate Residual method with pipelining.
807: The PIPEGCR Krylov method supports non-symmetric matrices and permits the use of a preconditioner
808: which may vary from one iteration to the next. Users can can define a method to vary the
809: preconditioner between iterates via KSPPIPEGCRSetModifyPC().
810: Restarts are solves with x0 not equal to zero. When a restart occurs, the initial starting
811: solution is given by the current estimate for x which was obtained by the last restart
812: iterations of the PIPEGCR algorithm.
813: The method implemented requires at most the storage of 4 x mmax + 5 vectors, roughly twice as much as GCR.
815: Only supports left preconditioning.
817: The natural norm for this method is (u,Au). This norm is available at no computational costs. Choosing norm types preconditioned or unpreconditioned involves a blocking reduction which prevents any benefit from pipelining.
819: Options Database Keys:
820: + -ksp_pipegcr_mmax <N> - the max number of Krylov directions to orthogonalize against
821: . -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)
822: . -ksp_pipegcr_nprealloc <N> - the number of vectors to preallocated for storing Krylov directions. Once exhausted new directions are allocated blockwise (default: 5)
823: . -ksp_pipegcr_truncation - Truncate number of previous Krylov directions
824: - -ksp_pipegcr_trancation_restart - Truncation-restart strategy: Keep at most mmax Krylov directions then restart (the default)
826: Level: beginner
828: Reference:
829: Pipelined, Flexible Krylov Subspace Methods
830: Patrick Sanan, Sascha M. Schnepp, and Dave A. May
832: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
833: KSPPIPEFGMRES, KSPPIPECG, KSPPIPECR, KSPPIPEFCG
835: M*/
838: PETSC_EXTERN PetscErrorCode KSPCreate_PIPEGCR(KSP ksp)
839: {
841: KSP_PIPEGCR *pipegcr;
844: PetscNewLog(ksp,&pipegcr);
845: pipegcr->mmax = KSPPIPEGCR_DEFAULT_MMAX;
846: pipegcr->nprealloc = KSPPIPEGCR_DEFAULT_NPREALLOC;
847: pipegcr->nvecs = 0;
848: pipegcr->vecb = KSPPIPEGCR_DEFAULT_VECB;
849: pipegcr->nchunks = 0;
850: pipegcr->truncstrat = KSPPIPEGCR_DEFAULT_TRUNCSTRAT;
851: pipegcr->n_restarts = 0;
852: pipegcr->unroll_w = KSPPIPEGCR_DEFAULT_UNROLL_W;
854: ksp->data = (void*)pipegcr;
856: /* natural norm is for free, precond+unprecond norm require non-overlapped reduction */
857: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
858: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,1);
859: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);
861: ksp->ops->setup = KSPSetUp_PIPEGCR;
862: ksp->ops->solve = KSPSolve_PIPEGCR;
863: ksp->ops->reset = KSPReset_PIPEGCR;
864: ksp->ops->destroy = KSPDestroy_PIPEGCR;
865: ksp->ops->view = KSPView_PIPEGCR;
866: ksp->ops->setfromoptions = KSPSetFromOptions_PIPEGCR;
867: ksp->ops->buildsolution = KSPBuildSolutionDefault;
868: ksp->ops->buildresidual = KSPBuildResidualDefault;
870: PetscObjectComposeFunction((PetscObject)ksp,"KSPPIPEGCRSetModifyPC_C",KSPPIPEGCRSetModifyPC_PIPEGCR);
871: return(0);
872: }