Actual source code: pipegcr.c

petsc-master 2020-09-24
Report Typos and Errors
  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:   PetscErrorCode  ierr;
 34:   PetscInt        i;
 35:   KSP_PIPEGCR     *pipegcr;
 36:   PetscInt        nnewvecs, nvecsprev;

 39:   pipegcr = (KSP_PIPEGCR*)ksp->data;

 41:   /* Allocate enough new vectors to add chunksize new vectors, reach nvecsneedtotal, or to reach mmax+1, whichever is smallest */
 42:   if (pipegcr->nvecs < PetscMin(pipegcr->mmax+1,nvecsneeded)){
 43:     nvecsprev = pipegcr->nvecs;
 44:     nnewvecs = PetscMin(PetscMax(nvecsneeded-pipegcr->nvecs,chunksize),pipegcr->mmax+1-pipegcr->nvecs);
 45:     KSPCreateVecs(ksp,nnewvecs,&pipegcr->ppvecs[pipegcr->nchunks],0,NULL);
 46:     PetscLogObjectParents((PetscObject)ksp,nnewvecs,pipegcr->ppvecs[pipegcr->nchunks]);
 47:     KSPCreateVecs(ksp,nnewvecs,&pipegcr->psvecs[pipegcr->nchunks],0,NULL);
 48:     PetscLogObjectParents((PetscObject)ksp,nnewvecs,pipegcr->psvecs[pipegcr->nchunks]);
 49:     KSPCreateVecs(ksp,nnewvecs,&pipegcr->pqvecs[pipegcr->nchunks],0,NULL);
 50:     PetscLogObjectParents((PetscObject)ksp,nnewvecs,pipegcr->pqvecs[pipegcr->nchunks]);
 51:     if (pipegcr->unroll_w) {
 52:       KSPCreateVecs(ksp,nnewvecs,&pipegcr->ptvecs[pipegcr->nchunks],0,NULL);
 53:       PetscLogObjectParents((PetscObject)ksp,nnewvecs,pipegcr->ptvecs[pipegcr->nchunks]);
 54:     }
 55:     pipegcr->nvecs += nnewvecs;
 56:     for (i=0;i<nnewvecs;i++){
 57:       pipegcr->qvecs[nvecsprev+i] = pipegcr->pqvecs[pipegcr->nchunks][i];
 58:       pipegcr->pvecs[nvecsprev+i] = pipegcr->ppvecs[pipegcr->nchunks][i];
 59:       pipegcr->svecs[nvecsprev+i] = pipegcr->psvecs[pipegcr->nchunks][i];
 60:       if (pipegcr->unroll_w) {
 61:         pipegcr->tvecs[nvecsprev+i] = pipegcr->ptvecs[pipegcr->nchunks][i];
 62:       }
 63:     }
 64:     pipegcr->chunksizes[pipegcr->nchunks] = nnewvecs;
 65:     pipegcr->nchunks++;
 66:   }
 67:   return(0);
 68: }

 70: static PetscErrorCode KSPSolve_PIPEGCR_cycle(KSP ksp)
 71: {
 72:   KSP_PIPEGCR    *pipegcr = (KSP_PIPEGCR*)ksp->data;
 74:   Mat            A, B;
 75:   Vec            x,r,b,z,w,m,n,p,s,q,t,*redux;
 76:   PetscInt       i,j,k,idx,kdx,mi;
 77:   PetscScalar    alpha=0.0,gamma,*betas,*dots;
 78:   PetscReal      rnorm=0.0, delta,*eta,*etas;


 82:   /* !!PS We have not checked these routines for use with complex numbers. The inner products
 83:      are likely not defined correctly for that case */
 84: #if (defined(PETSC_USE_COMPLEX) && !defined(PETSC_SKIP_COMPLEX))
 85:   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"PIPEGCR has not been implemented for use with complex scalars");
 86: #endif

 88:   KSPGetOperators(ksp, &A, &B);
 89:   x = ksp->vec_sol;
 90:   b = ksp->vec_rhs;
 91:   r = ksp->work[0];
 92:   z = ksp->work[1];
 93:   w = ksp->work[2]; /* w = Az = AB(r)                 (pipelining intermediate) */
 94:   m = ksp->work[3]; /* m = B(w) = B(Az) = B(AB(r))    (pipelining intermediate) */
 95:   n = ksp->work[4]; /* n = AB(w) = AB(Az) = AB(AB(r)) (pipelining intermediate) */
 96:   p = pipegcr->pvecs[0];
 97:   s = pipegcr->svecs[0];
 98:   q = pipegcr->qvecs[0];
 99:   t = pipegcr->unroll_w ? pipegcr->tvecs[0] : NULL;

101:   redux = pipegcr->redux;
102:   dots  = pipegcr->dots;
103:   etas  = pipegcr->etas;
104:   betas = dots;        /* dots takes the result of all dot products of which the betas are a subset */

106:   /* cycle initial residual */
107:   KSP_MatMult(ksp,A,x,r);
108:   VecAYPX(r,-1.0,b);                   /* r <- b - Ax */
109:   KSP_PCApply(ksp,r,z);                /* z <- B(r)   */
110:   KSP_MatMult(ksp,A,z,w);              /* w <- Az     */

112:   /* initialization of other variables and pipelining intermediates */
113:   VecCopy(z,p);
114:   KSP_MatMult(ksp,A,p,s);

116:   /* overlap initial computation of delta, gamma */
117:   redux[0] = w;
118:   redux[1] = r;
119:   VecMDotBegin(w,2,redux,dots);    /* Start split reductions for gamma = (w,r), delta = (w,w) */
120:   PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)s)); /* perform asynchronous reduction */
121:   KSP_PCApply(ksp,s,q);            /* q = B(s) */
122:   if (pipegcr->unroll_w) {
123:     KSP_MatMult(ksp,A,q,t);        /* t = Aq   */
124:   }
125:   VecMDotEnd(w,2,redux,dots);      /* Finish split reduction */
126:   delta    = PetscRealPart(dots[0]);
127:   etas[0]  = delta;
128:   gamma    = dots[1];
129:   alpha    = gamma/delta;

131:   i = 0;
132:   do {
133:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
134:     ksp->its++;
135:     PetscObjectSAWsGrantAccess((PetscObject)ksp);

137:     /* update solution, residuals, .. */
138:     VecAXPY(x,+alpha,p);
139:     VecAXPY(r,-alpha,s);
140:     VecAXPY(z,-alpha,q);
141:     if (pipegcr->unroll_w){
142:       VecAXPY(w,-alpha,t);
143:     } else {
144:       KSP_MatMult(ksp,A,z,w);
145:     }

147:     /* Computations of current iteration done */
148:     i++;

150:     if (pipegcr->modifypc) {
151:       (*pipegcr->modifypc)(ksp,ksp->its,ksp->rnorm,pipegcr->modifypc_ctx);
152:     }

154:     /* If needbe, allocate a new chunk of vectors */
155:     KSPAllocateVectors_PIPEGCR(ksp,i+1,pipegcr->vecb);

157:     /* Note that we wrap around and start clobbering old vectors */
158:     idx = i % (pipegcr->mmax+1);
159:     p   = pipegcr->pvecs[idx];
160:     s   = pipegcr->svecs[idx];
161:     q   = pipegcr->qvecs[idx];
162:     if (pipegcr->unroll_w) {
163:       t   = pipegcr->tvecs[idx];
164:     }
165:     eta = pipegcr->etas+idx;

167:     /* number of old directions to orthogonalize against */
168:     switch(pipegcr->truncstrat){
169:       case KSP_FCD_TRUNC_TYPE_STANDARD:
170:         mi = pipegcr->mmax;
171:         break;
172:       case KSP_FCD_TRUNC_TYPE_NOTAY:
173:         mi = ((i-1) % pipegcr->mmax)+1;
174:         break;
175:       default:
176:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unrecognized Truncation Strategy");
177:     }

179:     /* Pick old p,s,q,zeta in a way suitable for VecMDot */
180:     for (k=PetscMax(0,i-mi),j=0;k<i;j++,k++){
181:       kdx = k % (pipegcr->mmax+1);
182:       pipegcr->pold[j] = pipegcr->pvecs[kdx];
183:       pipegcr->sold[j] = pipegcr->svecs[kdx];
184:       pipegcr->qold[j] = pipegcr->qvecs[kdx];
185:       if (pipegcr->unroll_w) {
186:         pipegcr->told[j] = pipegcr->tvecs[kdx];
187:       }
188:       redux[j]         = pipegcr->svecs[kdx];
189:     }
190:     /* If the above loop is not run redux contains only r and w => all beta_k = 0, only gamma, delta != 0 */
191:     redux[j]   = r;
192:     redux[j+1] = w;

194:     /* Dot products */
195:     /* Start split reductions for beta_k = (w,s_k), gamma = (w,r), delta = (w,w) */
196:     VecMDotBegin(w,j+2,redux,dots);
197:     PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)w));

199:     /* B(w-r) + u stabilization */
200:     VecWAXPY(n,-1.0,r,w);              /* m = u + B(w-r): (a) ntmp = w-r              */
201:     KSP_PCApply(ksp,n,m);              /* m = u + B(w-r): (b) mtmp = B(ntmp) = B(w-r) */
202:     VecAXPY(m,1.0,z);                  /* m = u + B(w-r): (c) m = z + mtmp            */
203:     if (pipegcr->unroll_w){
204:       KSP_MatMult(ksp,A,m,n);          /* n = Am                                      */
205:     }

207:     /* Finish split reductions for beta_k = (w,s_k), gamma = (w,r), delta = (w,w) */
208:     VecMDotEnd(w,j+2,redux,dots);
209:     gamma = dots[j];
210:     delta = PetscRealPart(dots[j+1]);

212:     /* compute new residual norm.
213:        this cannot be done before this point so that the natural norm
214:        is available for free and the communication involved is overlapped */
215:     switch (ksp->normtype) {
216:     case KSP_NORM_PRECONDITIONED:
217:       VecNorm(z,NORM_2,&rnorm);        /* ||r|| <- sqrt(z'*z) */
218:       break;
219:     case KSP_NORM_UNPRECONDITIONED:
220:       VecNorm(r,NORM_2,&rnorm);        /* ||r|| <- sqrt(r'*r) */
221:       break;
222:     case KSP_NORM_NATURAL:
223:       rnorm = PetscSqrtReal(PetscAbsScalar(gamma));         /* ||r|| <- sqrt(r,w)  */
224:       break;
225:     case KSP_NORM_NONE:
226:       rnorm = 0.0;
227:       break;
228:     default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
229:     }

231:     /* Check for convergence */
232:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
233:     ksp->rnorm = rnorm;
234:     PetscObjectSAWsGrantAccess((PetscObject)ksp);
235:     KSPLogResidualHistory(ksp,rnorm);
236:     KSPMonitor(ksp,ksp->its,rnorm);
237:     (*ksp->converged)(ksp,ksp->its,rnorm,&ksp->reason,ksp->cnvP);
238:     if (ksp->reason) return(0);

240:     /* compute new eta and scale beta */
241:     *eta = 0.;
242:     for (k=PetscMax(0,i-mi),j=0;k<i;j++,k++){
243:       kdx = k % (pipegcr->mmax+1);
244:       betas[j] /= -etas[kdx];                               /* betak  /= etak */
245:       *eta -= ((PetscReal)(PetscAbsScalar(betas[j])*PetscAbsScalar(betas[j]))) * etas[kdx];
246:                                                             /* etaitmp = -betaik^2 * etak */
247:     }
248:     *eta += delta;                                          /* etai    = delta -betaik^2 * etak */

250:     /* check breakdown of eta = (s,s) */
251:     if (*eta < 0.) {
252:       pipegcr->norm_breakdown = PETSC_TRUE;
253:       PetscInfo1(ksp,"Restart due to square root breakdown at it = \n",ksp->its);
254:       break;
255:     } else {
256:       alpha= gamma/(*eta);                                  /* alpha = gamma/etai */
257:     }

259:     /* project out stored search directions using classical G-S */
260:     VecCopy(z,p);
261:     VecCopy(w,s);
262:     VecCopy(m,q);
263:     if (pipegcr->unroll_w) {
264:       VecCopy(n,t);
265:       VecMAXPY(t,j,betas,pipegcr->told); /* ti <- n  - sum_k beta_k t_k */
266:     }
267:     VecMAXPY(p,j,betas,pipegcr->pold); /* pi <- ui - sum_k beta_k p_k */
268:     VecMAXPY(s,j,betas,pipegcr->sold); /* si <- wi - sum_k beta_k s_k */
269:     VecMAXPY(q,j,betas,pipegcr->qold); /* qi <- m  - sum_k beta_k q_k */

271:   } while (ksp->its < ksp->max_it);
272:   if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
273:   return(0);
274: }

276: static PetscErrorCode KSPSolve_PIPEGCR(KSP ksp)
277: {
278:   KSP_PIPEGCR    *pipegcr = (KSP_PIPEGCR*)ksp->data;
280:   Mat            A, B;
281:   Vec            x,b,r,z,w;
282:   PetscScalar    gamma;
283:   PetscReal      rnorm=0.0;
284:   PetscBool      issym;

287:   PetscCitationsRegister(citation,&cited);

289:   KSPGetOperators(ksp, &A, &B);
290:   x = ksp->vec_sol;
291:   b = ksp->vec_rhs;
292:   r = ksp->work[0];
293:   z = ksp->work[1];
294:   w = ksp->work[2]; /* w = Az = AB(r)                 (pipelining intermediate) */

296:   /* compute initial residual */
297:   if (!ksp->guess_zero) {
298:     KSP_MatMult(ksp,A,x,r);
299:     VecAYPX(r,-1.0,b);                 /* r <- b - Ax       */
300:   } else {
301:     VecCopy(b,r);                      /* r <- b            */
302:   }

304:   /* initial residual norm */
305:   KSP_PCApply(ksp,r,z);                /* z <- B(r)         */
306:   KSP_MatMult(ksp,A,z,w);              /* w <- Az           */
307:   VecDot(r,w,&gamma);                  /* gamma = (r,w)     */

309:   switch (ksp->normtype) {
310:     case KSP_NORM_PRECONDITIONED:
311:       VecNorm(z,NORM_2,&rnorm);        /* ||r|| <- sqrt(z'*z) */
312:       break;
313:     case KSP_NORM_UNPRECONDITIONED:
314:       VecNorm(r,NORM_2,&rnorm);        /* ||r|| <- sqrt(r'*r) */
315:       break;
316:     case KSP_NORM_NATURAL:
317:       rnorm = PetscSqrtReal(PetscAbsScalar(gamma));         /* ||r|| <- sqrt(r,w)  */
318:       break;
319:     case KSP_NORM_NONE:
320:       rnorm = 0.0;
321:       break;
322:     default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
323:   }

325:   /* Is A symmetric? */
326:   PetscObjectTypeCompareAny((PetscObject)A,&issym,MATSBAIJ,MATSEQSBAIJ,MATMPISBAIJ,"");
327:   if (!issym) {
328:     PetscInfo(A,"Matrix type is not any of MATSBAIJ,MATSEQSBAIJ,MATMPISBAIJ. Is matrix A symmetric (as required by CR methods)?");
329:   }

331:   /* logging */
332:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
333:   ksp->its    = 0;
334:   ksp->rnorm0 = rnorm;
335:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
336:   KSPLogResidualHistory(ksp,ksp->rnorm0);
337:   KSPMonitor(ksp,ksp->its,ksp->rnorm0);
338:   (*ksp->converged)(ksp,ksp->its,ksp->rnorm0,&ksp->reason,ksp->cnvP);
339:   if (ksp->reason) return(0);

341:   do {
342:     KSPSolve_PIPEGCR_cycle(ksp);
343:     if (ksp->reason) return(0);
344:     if (pipegcr->norm_breakdown) {
345:       pipegcr->n_restarts++;
346:       pipegcr->norm_breakdown = PETSC_FALSE;
347:     }
348:   } while (ksp->its < ksp->max_it);

350:   if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
351:   return(0);
352: }

354: static PetscErrorCode KSPView_PIPEGCR(KSP ksp, PetscViewer viewer)
355: {
356:   KSP_PIPEGCR    *pipegcr = (KSP_PIPEGCR*)ksp->data;
358:   PetscBool      isascii,isstring;
359:   const char     *truncstr;

362:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII, &isascii);
363:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);

365:   if (pipegcr->truncstrat == KSP_FCD_TRUNC_TYPE_STANDARD){
366:     truncstr = "Using standard truncation strategy";
367:   } else if (pipegcr->truncstrat == KSP_FCD_TRUNC_TYPE_NOTAY){
368:     truncstr = "Using Notay's truncation strategy";
369:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Undefined FCD truncation strategy");


372:   if (isascii) {
373:     PetscViewerASCIIPrintf(viewer,"  max previous directions = %D\n",pipegcr->mmax);
374:     PetscViewerASCIIPrintf(viewer,"  preallocated %D directions\n",PetscMin(pipegcr->nprealloc,pipegcr->mmax+1));
375:     PetscViewerASCIIPrintf(viewer,"  %s\n",truncstr);
376:     PetscViewerASCIIPrintf(viewer,"  w unrolling = %D \n", pipegcr->unroll_w);
377:     PetscViewerASCIIPrintf(viewer,"  restarts performed = %D \n", pipegcr->n_restarts);
378:   } else if (isstring) {
379:     PetscViewerStringSPrintf(viewer, "max previous directions = %D, preallocated %D directions, %s truncation strategy", pipegcr->mmax,pipegcr->nprealloc,truncstr);
380:   }
381:   return(0);
382: }


385: static PetscErrorCode KSPSetUp_PIPEGCR(KSP ksp)
386: {
387:   KSP_PIPEGCR   *pipegcr = (KSP_PIPEGCR*)ksp->data;
389:   Mat            A;
390:   PetscBool      diagonalscale;
391:   const PetscInt nworkstd = 5;

394:   PCGetDiagonalScale(ksp->pc,&diagonalscale);
395:   if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);

397:   KSPGetOperators(ksp, &A, NULL);

399:   /* Allocate "standard" work vectors */
400:   KSPSetWorkVecs(ksp,nworkstd);

402:   /* Allocated space for pointers to additional work vectors
403:     note that mmax is the number of previous directions, so we add 1 for the current direction */
404:   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));
405:   if (pipegcr->unroll_w) {
406:     PetscMalloc3(pipegcr->mmax+1,&(pipegcr->tvecs),pipegcr->mmax+1,&(pipegcr->ptvecs),pipegcr->mmax+2,&(pipegcr->told));
407:   }
408:   PetscMalloc4(pipegcr->mmax+2,&(pipegcr->pold),pipegcr->mmax+2,&(pipegcr->sold),pipegcr->mmax+2,&(pipegcr->qold),pipegcr->mmax+2,&(pipegcr->chunksizes));
409:   PetscMalloc3(pipegcr->mmax+2,&(pipegcr->dots),pipegcr->mmax+1,&(pipegcr->etas),pipegcr->mmax+2,&(pipegcr->redux));
410:   /* If the requested number of preallocated vectors is greater than mmax reduce nprealloc */
411:   if (pipegcr->nprealloc > pipegcr->mmax+1){
412:     PetscInfo2(NULL,"Requested nprealloc=%d is greater than m_max+1=%d. Resetting nprealloc = m_max+1.\n",pipegcr->nprealloc, pipegcr->mmax+1);
413:   }

415:   /* Preallocate additional work vectors */
416:   KSPAllocateVectors_PIPEGCR(ksp,pipegcr->nprealloc,pipegcr->nprealloc);

418:   PetscLogObjectMemory(
419:     (PetscObject)ksp,
420:     (pipegcr->mmax + 1) * 4 * sizeof(Vec*) +        /* old dirs  */
421:     (pipegcr->mmax + 1) * 4 * sizeof(Vec**) +       /* old pdirs */
422:     (pipegcr->mmax + 1) * 4 * sizeof(Vec*) +        /* p/s/qold/told */
423:     (pipegcr->mmax + 1) *     sizeof(PetscInt) +    /* chunksizes */
424:     (pipegcr->mmax + 2) *     sizeof(Vec*) +        /* redux */
425:     (pipegcr->mmax + 2) *     sizeof(PetscScalar) + /* dots */
426:     (pipegcr->mmax + 1) *     sizeof(PetscReal)     /* etas */
427: );
428:   return(0);
429: }

431: static PetscErrorCode KSPReset_PIPEGCR(KSP ksp)
432: {
434:   KSP_PIPEGCR    *pipegcr = (KSP_PIPEGCR*)ksp->data;

437:   if (pipegcr->modifypc_destroy) {
438:     (*pipegcr->modifypc_destroy)(pipegcr->modifypc_ctx);
439:   }
440:   return(0);
441: }

443: static PetscErrorCode KSPDestroy_PIPEGCR(KSP ksp)
444: {
446:   PetscInt       i;
447:   KSP_PIPEGCR    *pipegcr = (KSP_PIPEGCR*)ksp->data;

450:   VecDestroyVecs(ksp->nwork,&ksp->work); /* Destroy "standard" work vecs */

452:   /* Destroy vectors for old directions and the arrays that manage pointers to them */
453:   if (pipegcr->nvecs){
454:     for (i=0;i<pipegcr->nchunks;i++){
455:       VecDestroyVecs(pipegcr->chunksizes[i],&pipegcr->ppvecs[i]);
456:       VecDestroyVecs(pipegcr->chunksizes[i],&pipegcr->psvecs[i]);
457:       VecDestroyVecs(pipegcr->chunksizes[i],&pipegcr->pqvecs[i]);
458:       if (pipegcr->unroll_w) {
459:         VecDestroyVecs(pipegcr->chunksizes[i],&pipegcr->ptvecs[i]);
460:       }
461:     }
462:   }

464:   PetscFree6(pipegcr->pvecs,pipegcr->ppvecs,pipegcr->svecs,pipegcr->psvecs,pipegcr->qvecs,pipegcr->pqvecs);
465:   PetscFree4(pipegcr->pold,pipegcr->sold,pipegcr->qold,pipegcr->chunksizes);
466:   PetscFree3(pipegcr->dots,pipegcr->etas,pipegcr->redux);
467:   if (pipegcr->unroll_w) {
468:     PetscFree3(pipegcr->tvecs,pipegcr->ptvecs,pipegcr->told);
469:   }

471:   KSPReset_PIPEGCR(ksp);
472:   KSPDestroyDefault(ksp);
473:   return(0);
474: }

476: /*@
477:   KSPPIPEGCRSetUnrollW - Set to PETSC_TRUE to use PIPEGCR with unrolling of the w vector

479:   Logically Collective on ksp

481:   Input Parameters:
482: +  ksp - the Krylov space context
483: -  unroll_w - use unrolling

485:   Level: intermediate

487:   Options Database:
488: . -ksp_pipegcr_unroll_w

490: .seealso: KSPPIPEGCR, KSPPIPEGCRSetTruncationType(), KSPPIPEGCRSetNprealloc(),KSPPIPEGCRGetUnrollW()
491: @*/
492: PetscErrorCode KSPPIPEGCRSetUnrollW(KSP ksp,PetscBool unroll_w)
493: {
494:   KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;

499:   pipegcr->unroll_w=unroll_w;
500:   return(0);
501: }

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: }

531: /*@
532:   KSPPIPEGCRSetMmax - set the maximum number of previous directions PIPEGCR will store for orthogonalization

534:   Note: mmax + 1 directions are stored (mmax previous ones along with a current one)
535:   and whether all are used in each iteration also depends on the truncation strategy
536:   (see KSPPIPEGCRSetTruncationType)

538:   Logically Collective on ksp

540:   Input Parameters:
541: +  ksp - the Krylov space context
542: -  mmax - the maximum number of previous directions to orthogonalize againt

544:   Level: intermediate

546:   Options Database:
547: . -ksp_pipegcr_mmax <N>

549: .seealso: KSPPIPEGCR, KSPPIPEGCRSetTruncationType(), KSPPIPEGCRSetNprealloc()
550: @*/
551: PetscErrorCode KSPPIPEGCRSetMmax(KSP ksp,PetscInt mmax)
552: {
553:   KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;

558:   pipegcr->mmax=mmax;
559:   return(0);
560: }

562: /*@
563:   KSPPIPEGCRGetMmax - get the maximum number of previous directions PIPEGCR will store

565:   Note: PIPEGCR stores mmax+1 directions at most (mmax previous ones, and one current one)

567:    Not Collective

569:    Input Parameter:
570: .  ksp - the Krylov space context

572:    Output Parameter:
573: .  mmax - the maximum number of previous directons allowed for orthogonalization

575:   Options Database:
576: . -ksp_pipegcr_mmax <N>

578:    Level: intermediate

580: .seealso: KSPPIPEGCR, KSPPIPEGCRGetTruncationType(), KSPPIPEGCRGetNprealloc(), KSPPIPEGCRSetMmax()
581: @*/

583: PetscErrorCode KSPPIPEGCRGetMmax(KSP ksp,PetscInt *mmax)
584: {
585:   KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;

589:   *mmax=pipegcr->mmax;
590:   return(0);
591: }

593: /*@
594:   KSPPIPEGCRSetNprealloc - set the number of directions to preallocate with PIPEGCR

596:   Logically Collective on ksp

598:   Input Parameters:
599: +  ksp - the Krylov space context
600: -  nprealloc - the number of vectors to preallocate

602:   Level: advanced

604:   Options Database:
605: . -ksp_pipegcr_nprealloc <N>

607: .seealso: KSPPIPEGCR, KSPPIPEGCRGetTruncationType(), KSPPIPEGCRGetNprealloc()
608: @*/
609: PetscErrorCode KSPPIPEGCRSetNprealloc(KSP ksp,PetscInt nprealloc)
610: {
611:   KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;

616:   pipegcr->nprealloc = nprealloc;
617:   return(0);
618: }

620: /*@
621:   KSPPIPEGCRGetNprealloc - get the number of directions preallocate by PIPEGCR

623:    Not Collective

625:    Input Parameter:
626: .  ksp - the Krylov space context

628:    Output Parameter:
629: .  nprealloc - the number of directions preallocated

631:   Options Database:
632: . -ksp_pipegcr_nprealloc <N>

634:    Level: advanced

636: .seealso: KSPPIPEGCR, KSPPIPEGCRGetTruncationType(), KSPPIPEGCRSetNprealloc()
637: @*/
638: PetscErrorCode KSPPIPEGCRGetNprealloc(KSP ksp,PetscInt *nprealloc)
639: {
640:   KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;

644:   *nprealloc = pipegcr->nprealloc;
645:   return(0);
646: }

648: /*@
649:   KSPPIPEGCRSetTruncationType - specify how many of its stored previous directions PIPEGCR uses during orthoganalization

651:   Logically Collective on ksp

653:   KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions
654:   KSP_FCD_TRUNC_TYPE_NOTAY uses the last max(1,mod(i,mmax)) directions at iteration i=0,1,..

656:   Input Parameters:
657: +  ksp - the Krylov space context
658: -  truncstrat - the choice of strategy

660:   Level: intermediate

662:   Options Database:
663: . -ksp_pipegcr_truncation_type <standard,notay> - which stored basis vectors to orthogonalize against

665: .seealso: KSPPIPEGCR, KSPPIPEGCRSetTruncationType, KSPPIPEGCRTruncationType, KSPFCDTruncationType
666: @*/
667: PetscErrorCode KSPPIPEGCRSetTruncationType(KSP ksp,KSPFCDTruncationType truncstrat)
668: {
669:   KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;

674:   pipegcr->truncstrat=truncstrat;
675:   return(0);
676: }

678: /*@
679:   KSPPIPEGCRGetTruncationType - get the truncation strategy employed by PIPEGCR

681:   Not Collective

683:   KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions
684:   KSP_FCD_TRUNC_TYPE_NOTAY uses the last max(1,mod(i,mmax)) directions at iteration i=0,1,..

686:    Input Parameter:
687: .  ksp - the Krylov space context

689:    Output Parameter:
690: .  truncstrat - the strategy type

692:   Options Database:
693: . -ksp_pipegcr_truncation_type <standard,notay> - which stored basis vectors to orthogonalize against

695:    Level: intermediate

697: .seealso: KSPPIPEGCR, KSPPIPEGCRSetTruncationType, KSPPIPEGCRTruncationType, KSPFCDTruncationType
698: @*/
699: PetscErrorCode KSPPIPEGCRGetTruncationType(KSP ksp,KSPFCDTruncationType *truncstrat)
700: {
701:   KSP_PIPEGCR *pipegcr=(KSP_PIPEGCR*)ksp->data;

705:   *truncstrat=pipegcr->truncstrat;
706:   return(0);
707: }

709: static PetscErrorCode KSPSetFromOptions_PIPEGCR(PetscOptionItems *PetscOptionsObject,KSP ksp)
710: {
712:   KSP_PIPEGCR    *pipegcr = (KSP_PIPEGCR*)ksp->data;
713:   PetscInt       mmax,nprealloc;
714:   PetscBool      flg;

717:   PetscOptionsHead(PetscOptionsObject,"KSP PIPEGCR options");
718:   PetscOptionsInt("-ksp_pipegcr_mmax","Number of search directions to storue","KSPPIPEGCRSetMmax",pipegcr->mmax,&mmax,&flg);
719:   if (flg) KSPPIPEGCRSetMmax(ksp,mmax);
720:   PetscOptionsInt("-ksp_pipegcr_nprealloc","Number of directions to preallocate","KSPPIPEGCRSetNprealloc",pipegcr->nprealloc,&nprealloc,&flg);
721:   if (flg) { KSPPIPEGCRSetNprealloc(ksp,nprealloc); }
722:   PetscOptionsEnum("-ksp_pipegcr_truncation_type","Truncation approach for directions","KSPFCGSetTruncationType",KSPFCDTruncationTypes,(PetscEnum)pipegcr->truncstrat,(PetscEnum*)&pipegcr->truncstrat,NULL);
723:   PetscOptionsBool("-ksp_pipegcr_unroll_w","Use unrolling of w","KSPPIPEGCRSetUnrollW",pipegcr->unroll_w,&pipegcr->unroll_w,NULL);
724:   PetscOptionsTail();
725:   return(0);
726: }

729: typedef PetscErrorCode (*KSPPIPEGCRModifyPCFunction)(KSP,PetscInt,PetscReal,void*);
730: typedef PetscErrorCode (*KSPPIPEGCRDestroyFunction)(void*);

732: static PetscErrorCode  KSPPIPEGCRSetModifyPC_PIPEGCR(KSP ksp,KSPPIPEGCRModifyPCFunction function,void *data,KSPPIPEGCRDestroyFunction destroy)
733: {
734:   KSP_PIPEGCR *pipegcr = (KSP_PIPEGCR*)ksp->data;

738:   pipegcr->modifypc         = function;
739:   pipegcr->modifypc_destroy = destroy;
740:   pipegcr->modifypc_ctx     = data;
741:   return(0);
742: }

744: /*@C
745:  KSPPIPEGCRSetModifyPC - Sets the routine used by PIPEGCR to modify the preconditioner.

747:  Logically Collective on ksp

749:  Input Parameters:
750:  +  ksp      - iterative context obtained from KSPCreate()
751:  .  function - user defined function to modify the preconditioner
752:  .  ctx      - user provided contex for the modify preconditioner function
753:  -  destroy  - the function to use to destroy the user provided application context.

755:  Calling Sequence of function:
756:   PetscErrorCode function (KSP ksp, PetscInt n, PetscReal rnorm, void *ctx)

758:  ksp   - iterative context
759:  n     - the total number of PIPEGCR iterations that have occurred
760:  rnorm - 2-norm residual value
761:  ctx   - the user provided application context

763:  Level: intermediate

765:  Notes:
766:  The default modifypc routine is KSPPIPEGCRModifyPCNoChange()

768:  .seealso: KSPPIPEGCRModifyPCNoChange()

770:  @*/
771: PetscErrorCode  KSPPIPEGCRSetModifyPC(KSP ksp,PetscErrorCode (*function)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*destroy)(void*))
772: {

776:   PetscUseMethod(ksp,"KSPPIPEGCRSetModifyPC_C",(KSP,PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*),void *data,PetscErrorCode (*)(void*)),(ksp,function,data,destroy));
777:   return(0);
778: }

780: /*MC
781:      KSPPIPEGCR - Implements a Pipelined Generalized Conjugate Residual method.

783:   Options Database Keys:
784: +   -ksp_pipegcr_mmax <N>  - the max number of Krylov directions to orthogonalize against
785: .   -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)
786: .   -ksp_pipegcr_nprealloc <N> - the number of vectors to preallocated for storing Krylov directions. Once exhausted new directions are allocated blockwise (default: 5)
787: -   -ksp_pipegcr_truncation_type <standard,notay> - which previous search directions to orthogonalize against


790:   Notes:
791:     The PIPEGCR Krylov method supports non-symmetric matrices and permits the use of a preconditioner
792:     which may vary from one iteration to the next. Users can can define a method to vary the
793:     preconditioner between iterates via KSPPIPEGCRSetModifyPC().
794:     Restarts are solves with x0 not equal to zero. When a restart occurs, the initial starting
795:     solution is given by the current estimate for x which was obtained by the last restart
796:     iterations of the PIPEGCR algorithm.
797:     The method implemented requires at most the storage of 4 x mmax + 5 vectors, roughly twice as much as GCR.

799:     Only supports left preconditioning.

801:     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.

803:   Reference:
804:     P. Sanan, S.M. Schnepp, and D.A. May,
805:     "Pipelined, Flexible Krylov Subspace Methods,"
806:     SIAM Journal on Scientific Computing 2016 38:5, C441-C470,
807:     DOI: 10.1137/15M1049130

809:    Level: intermediate

811: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
812:            KSPPIPEFGMRES, KSPPIPECG, KSPPIPECR, KSPPIPEFCG,KSPPIPEGCRSetTruncationType(),KSPPIPEGCRSetNprealloc(),KSPPIPEGCRSetUnrollW(),KSPPIPEGCRSetMmax()


815: M*/
816: PETSC_EXTERN PetscErrorCode KSPCreate_PIPEGCR(KSP ksp)
817: {
819:   KSP_PIPEGCR    *pipegcr;

822:   PetscNewLog(ksp,&pipegcr);
823:   pipegcr->mmax       = KSPPIPEGCR_DEFAULT_MMAX;
824:   pipegcr->nprealloc  = KSPPIPEGCR_DEFAULT_NPREALLOC;
825:   pipegcr->nvecs      = 0;
826:   pipegcr->vecb       = KSPPIPEGCR_DEFAULT_VECB;
827:   pipegcr->nchunks    = 0;
828:   pipegcr->truncstrat = KSPPIPEGCR_DEFAULT_TRUNCSTRAT;
829:   pipegcr->n_restarts = 0;
830:   pipegcr->unroll_w   = KSPPIPEGCR_DEFAULT_UNROLL_W;

832:   ksp->data       = (void*)pipegcr;

834:   /* natural norm is for free, precond+unprecond norm require non-overlapped reduction */
835:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
836:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,1);
837:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);
838:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);

840:   ksp->ops->setup          = KSPSetUp_PIPEGCR;
841:   ksp->ops->solve          = KSPSolve_PIPEGCR;
842:   ksp->ops->reset          = KSPReset_PIPEGCR;
843:   ksp->ops->destroy        = KSPDestroy_PIPEGCR;
844:   ksp->ops->view           = KSPView_PIPEGCR;
845:   ksp->ops->setfromoptions = KSPSetFromOptions_PIPEGCR;
846:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
847:   ksp->ops->buildresidual  = KSPBuildResidualDefault;

849:   PetscObjectComposeFunction((PetscObject)ksp,"KSPPIPEGCRSetModifyPC_C",KSPPIPEGCRSetModifyPC_PIPEGCR);
850:   return(0);
851: }