Actual source code: deflation.c

  1: #include <../src/ksp/pc/impls/deflation/deflation.h>

  3: const char *const PCDeflationSpaceTypes[] = {"haar", "db2", "db4", "db8", "db16", "biorth22", "meyer", "aggregation", "user", "PCDeflationSpaceType", "PC_DEFLATION_SPACE_", NULL};

  5: static PetscErrorCode PCDeflationSetInitOnly_Deflation(PC pc, PetscBool flg)
  6: {
  7:   PC_Deflation *def = (PC_Deflation *)pc->data;

  9:   def->init = flg;
 10:   return 0;
 11: }

 13: /*@
 14:    PCDeflationSetInitOnly - Do only initialization step.
 15:     Sets initial guess to the solution on the deflation space but does not apply
 16:     the deflation preconditioner. The additional preconditioner is still applied.

 18:    Logically Collective

 20:    Input Parameters:
 21: +  pc  - the preconditioner context
 22: -  flg - default `PETSC_FALSE`

 24:    Options Database Key:
 25: .    -pc_deflation_init_only <false> - if true computes only the special guess

 27:    Level: intermediate

 29: .seealso: `PCDEFLATION`
 30: @*/
 31: PetscErrorCode PCDeflationSetInitOnly(PC pc, PetscBool flg)
 32: {
 35:   PetscTryMethod(pc, "PCDeflationSetInitOnly_C", (PC, PetscBool), (pc, flg));
 36:   return 0;
 37: }

 39: static PetscErrorCode PCDeflationSetLevels_Deflation(PC pc, PetscInt current, PetscInt max)
 40: {
 41:   PC_Deflation *def = (PC_Deflation *)pc->data;

 43:   if (current) def->lvl = current;
 44:   def->maxlvl = max;
 45:   return 0;
 46: }

 48: /*@
 49:    PCDeflationSetLevels - Set the maximum level of deflation nesting.

 51:    Logically Collective

 53:    Input Parameters:
 54: +  pc  - the preconditioner context
 55: -  max - maximum deflation level

 57:    Options Database Key:
 58: .    -pc_deflation_max_lvl <0> - maximum number of levels for multilevel deflation

 60:    Level: intermediate

 62: .seealso: `PCDeflationSetSpaceToCompute()`, `PCDeflationSetSpace()`, `PCDEFLATION`
 63: @*/
 64: PetscErrorCode PCDeflationSetLevels(PC pc, PetscInt max)
 65: {
 68:   PetscTryMethod(pc, "PCDeflationSetLevels_C", (PC, PetscInt, PetscInt), (pc, 0, max));
 69:   return 0;
 70: }

 72: static PetscErrorCode PCDeflationSetReductionFactor_Deflation(PC pc, PetscInt red)
 73: {
 74:   PC_Deflation *def = (PC_Deflation *)pc->data;

 76:   def->reductionfact = red;
 77:   return 0;
 78: }

 80: /*@
 81:    PCDeflationSetReductionFactor - Set reduction factor for the `PCDEFLATION`

 83:    Logically Collective

 85:    Input Parameters:
 86: +  pc  - the preconditioner context
 87: -  red - reduction factor (or `PETSC_DETERMINE`)

 89:    Options Database Key:
 90: .    -pc_deflation_reduction_factor <\-1> - reduction factor on bottom level coarse problem for `PCDEFLATION`

 92:    Note:
 93:    Default is computed based on the size of the coarse problem.

 95:    Level: intermediate

 97: .seealso: `PCTELESCOPE`, `PCDEFLATION`, `PCDeflationSetLevels()`
 98: @*/
 99: PetscErrorCode PCDeflationSetReductionFactor(PC pc, PetscInt red)
100: {
103:   PetscTryMethod(pc, "PCDeflationSetReductionFactor_C", (PC, PetscInt), (pc, red));
104:   return 0;
105: }

107: static PetscErrorCode PCDeflationSetCorrectionFactor_Deflation(PC pc, PetscScalar fact)
108: {
109:   PC_Deflation *def = (PC_Deflation *)pc->data;

111:   /* TODO PETSC_DETERMINE -> compute max eigenvalue with power method */
112:   def->correct     = PETSC_TRUE;
113:   def->correctfact = fact;
114:   if (def->correct == 0.0) def->correct = PETSC_FALSE;
115:   return 0;
116: }

118: /*@
119:    PCDeflationSetCorrectionFactor - Set coarse problem correction factor.
120:     The preconditioner becomes P*M^{-1} + fact*Q.

122:    Logically Collective

124:    Input Parameters:
125: +  pc   - the preconditioner context
126: -  fact - correction factor

128:    Options Database Keys:
129: +    -pc_deflation_correction        <false> - if true apply coarse problem correction
130: -    -pc_deflation_correction_factor <1.0>   - sets coarse problem correction factor

132:    Note:
133:     Any non-zero fact enables the coarse problem correction.

135:    Level: intermediate

137: .seealso: `PCDEFLATION`, `PCDeflationSetLevels()`, `PCDeflationSetReductionFactor()`
138: @*/
139: PetscErrorCode PCDeflationSetCorrectionFactor(PC pc, PetscScalar fact)
140: {
143:   PetscTryMethod(pc, "PCDeflationSetCorrectionFactor_C", (PC, PetscScalar), (pc, fact));
144:   return 0;
145: }

147: static PetscErrorCode PCDeflationSetSpaceToCompute_Deflation(PC pc, PCDeflationSpaceType type, PetscInt size)
148: {
149:   PC_Deflation *def = (PC_Deflation *)pc->data;

151:   if (type) def->spacetype = type;
152:   if (size > 0) def->spacesize = size;
153:   return 0;
154: }

156: /*@
157:    PCDeflationSetSpaceToCompute - Set deflation space type and size to compute.

159:    Logically Collective

161:    Input Parameters:
162: +  pc   - the preconditioner context
163: .  type - deflation space type to compute (or `PETSC_IGNORE`)
164: -  size - size of the space to compute (or `PETSC_DEFAULT`)

166:    Options Database Keys:
167: +    -pc_deflation_compute_space      <haar> - compute `PCDeflationSpaceType` deflation space
168: -    -pc_deflation_compute_space_size <1>    - size of the deflation space

170:    Notes:
171:     For wavelet-based deflation, size represents number of levels.

173:     The deflation space is computed in `PCSetUp()`.

175:    Level: intermediate

177: .seealso: `PCDeflationSetLevels()`, `PCDEFLATION`
178: @*/
179: PetscErrorCode PCDeflationSetSpaceToCompute(PC pc, PCDeflationSpaceType type, PetscInt size)
180: {
184:   PetscTryMethod(pc, "PCDeflationSetSpaceToCompute_C", (PC, PCDeflationSpaceType, PetscInt), (pc, type, size));
185:   return 0;
186: }

188: static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc, Mat W, PetscBool transpose)
189: {
190:   PC_Deflation *def = (PC_Deflation *)pc->data;

192:   /* possibly allows W' = Wt (which is valid but not tested) */
193:   PetscObjectReference((PetscObject)W);
194:   if (transpose) {
195:     MatDestroy(&def->Wt);
196:     def->Wt = W;
197:   } else {
198:     MatDestroy(&def->W);
199:     def->W = W;
200:   }
201:   return 0;
202: }

204: /*@
205:    PCDeflationSetSpace - Set the deflation space matrix (or its (Hermitian) transpose).

207:    Logically Collective

209:    Input Parameters:
210: +  pc        - the preconditioner context
211: .  W         - deflation matrix
212: -  transpose - indicates that W is an explicit transpose of the deflation matrix

214:    Notes:
215:     Setting W as a multipliplicative `MATCOMPOSITE` enables use of the multilevel
216:     deflation. If W = W0*W1*W2*...*Wn, W0 is taken as the first deflation space and
217:     the coarse problem (W0'*A*W0)^{-1} is again preconditioned by deflation with
218:     W1 as the deflation matrix. This repeats until the maximum level set by
219:     PCDeflationSetLevels() is reached or there are no more matrices available.
220:     If there are matrices left after reaching the maximum level,
221:     they are merged into a deflation matrix ...*W{n-1}*Wn.

223:    Level: intermediate

225: .seealso: `PCDeflationSetLevels()`, `PCDEFLATION`, `PCDeflationSetProjectionNullSpaceMat()`
226: @*/
227: PetscErrorCode PCDeflationSetSpace(PC pc, Mat W, PetscBool transpose)
228: {
232:   PetscTryMethod(pc, "PCDeflationSetSpace_C", (PC, Mat, PetscBool), (pc, W, transpose));
233:   return 0;
234: }

236: static PetscErrorCode PCDeflationSetProjectionNullSpaceMat_Deflation(PC pc, Mat mat)
237: {
238:   PC_Deflation *def = (PC_Deflation *)pc->data;

240:   PetscObjectReference((PetscObject)mat);
241:   MatDestroy(&def->WtA);
242:   def->WtA = mat;
243:   return 0;
244: }

246: /*@
247:    PCDeflationSetProjectionNullSpaceMat - Set the projection null space matrix (W'*A).

249:    Collective

251:    Input Parameters:
252: +  pc  - preconditioner context
253: -  mat - projection null space matrix

255:    Level: developer

257: .seealso: `PCDEFLATION`, `PCDeflationSetSpace()`
258: @*/
259: PetscErrorCode PCDeflationSetProjectionNullSpaceMat(PC pc, Mat mat)
260: {
263:   PetscTryMethod(pc, "PCDeflationSetProjectionNullSpaceMat_C", (PC, Mat), (pc, mat));
264:   return 0;
265: }

267: static PetscErrorCode PCDeflationSetCoarseMat_Deflation(PC pc, Mat mat)
268: {
269:   PC_Deflation *def = (PC_Deflation *)pc->data;

271:   PetscObjectReference((PetscObject)mat);
272:   MatDestroy(&def->WtAW);
273:   def->WtAW = mat;
274:   return 0;
275: }

277: /*@
278:    PCDeflationSetCoarseMat - Set the coarse problem `Mat`.

280:    Collective

282:    Input Parameters:
283: +  pc  - preconditioner context
284: -  mat - coarse problem mat

286:    Level: developer

288: .seealso: `PCDEFLATION`, `PCDeflationGetCoarseKSP()`
289: @*/
290: PetscErrorCode PCDeflationSetCoarseMat(PC pc, Mat mat)
291: {
294:   PetscTryMethod(pc, "PCDeflationSetCoarseMat_C", (PC, Mat), (pc, mat));
295:   return 0;
296: }

298: static PetscErrorCode PCDeflationGetCoarseKSP_Deflation(PC pc, KSP *ksp)
299: {
300:   PC_Deflation *def = (PC_Deflation *)pc->data;

302:   *ksp = def->WtAWinv;
303:   return 0;
304: }

306: /*@
307:    PCDeflationGetCoarseKSP - Returns the coarse problem `KSP`.

309:    Not Collective

311:    Input Parameter:
312: .  pc - preconditioner context

314:    Output Parameter:
315: .  ksp - coarse problem `KSP` context

317:    Level: advanced

319: .seealso: `PCDEFLATION`, `PCDeflationSetCoarseMat()`
320: @*/
321: PetscErrorCode PCDeflationGetCoarseKSP(PC pc, KSP *ksp)
322: {
325:   PetscTryMethod(pc, "PCDeflationGetCoarseKSP_C", (PC, KSP *), (pc, ksp));
326:   return 0;
327: }

329: static PetscErrorCode PCDeflationGetPC_Deflation(PC pc, PC *apc)
330: {
331:   PC_Deflation *def = (PC_Deflation *)pc->data;

333:   *apc = def->pc;
334:   return 0;
335: }

337: /*@
338:    PCDeflationGetPC - Returns the additional preconditioner M^{-1}.

340:    Not Collective

342:    Input Parameter:
343: .  pc  - the preconditioner context

345:    Output Parameter:
346: .  apc - additional preconditioner

348:    Level: advanced

350: .seealso: `PCDEFLATION`, `PCDeflationGetCoarseKSP()`
351: @*/
352: PetscErrorCode PCDeflationGetPC(PC pc, PC *apc)
353: {
356:   PetscTryMethod(pc, "PCDeflationGetPC_C", (PC, PC *), (pc, apc));
357:   return 0;
358: }

360: /*
361:   x <- x + W*(W'*A*W)^{-1}*W'*r  = x + Q*r
362: */
363: static PetscErrorCode PCPreSolve_Deflation(PC pc, KSP ksp, Vec b, Vec x)
364: {
365:   PC_Deflation *def = (PC_Deflation *)pc->data;
366:   Mat           A;
367:   Vec           r, w1, w2;
368:   PetscBool     nonzero;

370:   w1 = def->workcoarse[0];
371:   w2 = def->workcoarse[1];
372:   r  = def->work;
373:   PCGetOperators(pc, NULL, &A);

375:   KSPGetInitialGuessNonzero(ksp, &nonzero);
376:   KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
377:   if (nonzero) {
378:     MatMult(A, x, r); /*    r  <- b - Ax              */
379:     VecAYPX(r, -1.0, b);
380:   } else {
381:     VecCopy(b, r); /*    r  <- b (x is 0)          */
382:   }

384:   if (def->Wt) {
385:     MatMult(def->Wt, r, w1); /*    w1 <- W'*r                */
386:   } else {
387:     MatMultHermitianTranspose(def->W, r, w1); /*    w1 <- W'*r                */
388:   }
389:   KSPSolve(def->WtAWinv, w1, w2); /*    w2 <- (W'*A*W)^{-1}*w1    */
390:   MatMult(def->W, w2, r);         /*    r  <- W*w2                */
391:   VecAYPX(x, 1.0, r);
392:   return 0;
393: }

395: /*
396:   if (def->correct) {
397:     z <- M^{-1}r - W*(W'*A*W)^{-1}*(W'*A*M^{-1}r - l*W'*r) = (P*M^{-1} + l*Q)*r
398:   } else {
399:     z <- M^{-1}*r - W*(W'*A*W)^{-1}*W'*A*M{-1}*r = P*M^{-1}*r
400:   }
401: */
402: static PetscErrorCode PCApply_Deflation(PC pc, Vec r, Vec z)
403: {
404:   PC_Deflation *def = (PC_Deflation *)pc->data;
405:   Mat           A;
406:   Vec           u, w1, w2;

408:   w1 = def->workcoarse[0];
409:   w2 = def->workcoarse[1];
410:   u  = def->work;
411:   PCGetOperators(pc, NULL, &A);

413:   PCApply(def->pc, r, z); /*    z <- M^{-1}*r             */
414:   if (!def->init) {
415:     MatMult(def->WtA, z, w1); /*    w1 <- W'*A*z              */
416:     if (def->correct) {
417:       if (def->Wt) {
418:         MatMult(def->Wt, r, w2); /*    w2 <- W'*r                */
419:       } else {
420:         MatMultHermitianTranspose(def->W, r, w2); /*    w2 <- W'*r                */
421:       }
422:       VecAXPY(w1, -1.0 * def->correctfact, w2); /*    w1 <- w1 - l*w2           */
423:     }
424:     KSPSolve(def->WtAWinv, w1, w2); /*    w2 <- (W'*A*W)^{-1}*w1    */
425:     MatMult(def->W, w2, u);         /*    u  <- W*w2                */
426:     VecAXPY(z, -1.0, u);            /*    z  <- z - u               */
427:   }
428:   return 0;
429: }

431: static PetscErrorCode PCSetUp_Deflation(PC pc)
432: {
433:   PC_Deflation    *def = (PC_Deflation *)pc->data;
434:   KSP              innerksp;
435:   PC               pcinner;
436:   Mat              Amat, nextDef = NULL, *mats;
437:   PetscInt         i, m, red, size;
438:   PetscMPIInt      commsize;
439:   PetscBool        match, flgspd, isset, transp = PETSC_FALSE;
440:   MatCompositeType ctype;
441:   MPI_Comm         comm;
442:   char             prefix[128] = "";

444:   if (pc->setupcalled) return 0;
445:   PetscObjectGetComm((PetscObject)pc, &comm);
446:   PCGetOperators(pc, NULL, &Amat);
447:   if (!def->lvl && !def->prefix) PCGetOptionsPrefix(pc, &def->prefix);
448:   if (def->lvl) PetscSNPrintf(prefix, sizeof(prefix), "%d_", (int)def->lvl);

450:   /* compute a deflation space */
451:   if (def->W || def->Wt) {
452:     def->spacetype = PC_DEFLATION_SPACE_USER;
453:   } else {
454:     PCDeflationComputeSpace(pc);
455:   }

457:   /* nested deflation */
458:   if (def->W) {
459:     PetscObjectTypeCompare((PetscObject)def->W, MATCOMPOSITE, &match);
460:     if (match) {
461:       MatCompositeGetType(def->W, &ctype);
462:       MatCompositeGetNumberMat(def->W, &size);
463:     }
464:   } else {
465:     MatCreateHermitianTranspose(def->Wt, &def->W);
466:     PetscObjectTypeCompare((PetscObject)def->Wt, MATCOMPOSITE, &match);
467:     if (match) {
468:       MatCompositeGetType(def->Wt, &ctype);
469:       MatCompositeGetNumberMat(def->Wt, &size);
470:     }
471:     transp = PETSC_TRUE;
472:   }
473:   if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) {
474:     if (!transp) {
475:       if (def->lvl < def->maxlvl) {
476:         PetscMalloc1(size, &mats);
477:         for (i = 0; i < size; i++) MatCompositeGetMat(def->W, i, &mats[i]);
478:         size -= 1;
479:         MatDestroy(&def->W);
480:         def->W = mats[size];
481:         PetscObjectReference((PetscObject)mats[size]);
482:         if (size > 1) {
483:           MatCreateComposite(comm, size, mats, &nextDef);
484:           MatCompositeSetType(nextDef, MAT_COMPOSITE_MULTIPLICATIVE);
485:         } else {
486:           nextDef = mats[0];
487:           PetscObjectReference((PetscObject)mats[0]);
488:         }
489:         PetscFree(mats);
490:       } else {
491:         /* TODO test merge side performance */
492:         /* MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT); */
493:         MatCompositeMerge(def->W);
494:       }
495:     } else {
496:       if (def->lvl < def->maxlvl) {
497:         PetscMalloc1(size, &mats);
498:         for (i = 0; i < size; i++) MatCompositeGetMat(def->Wt, i, &mats[i]);
499:         size -= 1;
500:         MatDestroy(&def->Wt);
501:         def->Wt = mats[0];
502:         PetscObjectReference((PetscObject)mats[0]);
503:         if (size > 1) {
504:           MatCreateComposite(comm, size, &mats[1], &nextDef);
505:           MatCompositeSetType(nextDef, MAT_COMPOSITE_MULTIPLICATIVE);
506:         } else {
507:           nextDef = mats[1];
508:           PetscObjectReference((PetscObject)mats[1]);
509:         }
510:         PetscFree(mats);
511:       } else {
512:         /* MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT); */
513:         MatCompositeMerge(def->Wt);
514:       }
515:     }
516:   }

518:   if (transp) {
519:     MatDestroy(&def->W);
520:     MatHermitianTranspose(def->Wt, MAT_INITIAL_MATRIX, &def->W);
521:   }

523:   /* assemble WtA */
524:   if (!def->WtA) {
525:     if (def->Wt) {
526:       MatMatMult(def->Wt, Amat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &def->WtA);
527:     } else {
528: #if defined(PETSC_USE_COMPLEX)
529:       MatHermitianTranspose(def->W, MAT_INITIAL_MATRIX, &def->Wt);
530:       MatMatMult(def->Wt, Amat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &def->WtA);
531: #else
532:       MatTransposeMatMult(def->W, Amat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &def->WtA);
533: #endif
534:     }
535:   }
536:   /* setup coarse problem */
537:   if (!def->WtAWinv) {
538:     MatGetSize(def->W, NULL, &m);
539:     if (!def->WtAW) {
540:       MatMatMult(def->WtA, def->W, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &def->WtAW);
541:       /* TODO create MatInheritOption(Mat,MatOption) */
542:       MatIsSPDKnown(Amat, &isset, &flgspd);
543:       if (isset) MatSetOption(def->WtAW, MAT_SPD, flgspd);
544:       if (PetscDefined(USE_DEBUG)) {
545:         /* Check columns of W are not in kernel of A */
546:         PetscReal *norms;

548:         PetscMalloc1(m, &norms);
549:         MatGetColumnNorms(def->WtAW, NORM_INFINITY, norms);
551:         PetscFree(norms);
552:       }
553:     } else MatIsSPDKnown(def->WtAW, &isset, &flgspd);

555:     /* TODO use MATINV ? */
556:     KSPCreate(comm, &def->WtAWinv);
557:     KSPSetOperators(def->WtAWinv, def->WtAW, def->WtAW);
558:     KSPGetPC(def->WtAWinv, &pcinner);
559:     /* Setup KSP and PC */
560:     if (nextDef) { /* next level for multilevel deflation */
561:       innerksp = def->WtAWinv;
562:       /* set default KSPtype */
563:       if (!def->ksptype) {
564:         def->ksptype = KSPFGMRES;
565:         if (isset && flgspd) { /* SPD system */
566:           def->ksptype = KSPFCG;
567:         }
568:       }
569:       KSPSetType(innerksp, def->ksptype); /* TODO iherit from KSP + tolerances */
570:       PCSetType(pcinner, PCDEFLATION);    /* TODO create coarse preconditinoner M_c = WtMW ? */
571:       PCDeflationSetSpace(pcinner, nextDef, transp);
572:       PCDeflationSetLevels_Deflation(pcinner, def->lvl + 1, def->maxlvl);
573:       /* inherit options */
574:       if (def->prefix) ((PC_Deflation *)(pcinner->data))->prefix = def->prefix;
575:       ((PC_Deflation *)(pcinner->data))->init          = def->init;
576:       ((PC_Deflation *)(pcinner->data))->ksptype       = def->ksptype;
577:       ((PC_Deflation *)(pcinner->data))->correct       = def->correct;
578:       ((PC_Deflation *)(pcinner->data))->correctfact   = def->correctfact;
579:       ((PC_Deflation *)(pcinner->data))->reductionfact = def->reductionfact;
580:       MatDestroy(&nextDef);
581:     } else { /* the last level */
582:       KSPSetType(def->WtAWinv, KSPPREONLY);
583:       PCSetType(pcinner, PCTELESCOPE);
584:       /* do not overwrite PCTELESCOPE */
585:       if (def->prefix) KSPSetOptionsPrefix(def->WtAWinv, def->prefix);
586:       KSPAppendOptionsPrefix(def->WtAWinv, "deflation_tel_");
587:       PCSetFromOptions(pcinner);
588:       PetscObjectTypeCompare((PetscObject)pcinner, PCTELESCOPE, &match);
590:       /* Reduction factor choice */
591:       red = def->reductionfact;
592:       if (red < 0) {
593:         MPI_Comm_size(comm, &commsize);
594:         red = PetscCeilInt(commsize, PetscCeilInt(m, commsize));
595:         PetscObjectTypeCompareAny((PetscObject)(def->WtAW), &match, MATSEQDENSE, MATMPIDENSE, MATDENSE, "");
596:         if (match) red = commsize;
597:         PetscInfo(pc, "Auto choosing reduction factor %" PetscInt_FMT "\n", red);
598:       }
599:       PCTelescopeSetReductionFactor(pcinner, red);
600:       PCSetUp(pcinner);
601:       PCTelescopeGetKSP(pcinner, &innerksp);
602:       if (innerksp) {
603:         KSPGetPC(innerksp, &pcinner);
604:         PCSetType(pcinner, PCLU);
605: #if defined(PETSC_HAVE_SUPERLU)
606:         MatGetFactorAvailable(def->WtAW, MATSOLVERSUPERLU, MAT_FACTOR_LU, &match);
607:         if (match) PCFactorSetMatSolverType(pcinner, MATSOLVERSUPERLU);
608: #endif
609: #if defined(PETSC_HAVE_SUPERLU_DIST)
610:         MatGetFactorAvailable(def->WtAW, MATSOLVERSUPERLU_DIST, MAT_FACTOR_LU, &match);
611:         if (match) PCFactorSetMatSolverType(pcinner, MATSOLVERSUPERLU_DIST);
612: #endif
613:       }
614:     }

616:     if (innerksp) {
617:       if (def->prefix) {
618:         KSPSetOptionsPrefix(innerksp, def->prefix);
619:         KSPAppendOptionsPrefix(innerksp, "deflation_");
620:       } else {
621:         KSPSetOptionsPrefix(innerksp, "deflation_");
622:       }
623:       KSPAppendOptionsPrefix(innerksp, prefix);
624:       KSPSetFromOptions(innerksp);
625:       KSPSetUp(innerksp);
626:     }
627:   }
628:   KSPSetFromOptions(def->WtAWinv);
629:   KSPSetUp(def->WtAWinv);

631:   /* create preconditioner */
632:   if (!def->pc) {
633:     PCCreate(comm, &def->pc);
634:     PCSetOperators(def->pc, Amat, Amat);
635:     PCSetType(def->pc, PCNONE);
636:     if (def->prefix) PCSetOptionsPrefix(def->pc, def->prefix);
637:     PCAppendOptionsPrefix(def->pc, "deflation_");
638:     PCAppendOptionsPrefix(def->pc, prefix);
639:     PCAppendOptionsPrefix(def->pc, "pc_");
640:     PCSetFromOptions(def->pc);
641:     PCSetUp(def->pc);
642:   }

644:   /* create work vecs */
645:   MatCreateVecs(Amat, NULL, &def->work);
646:   KSPCreateVecs(def->WtAWinv, 2, &def->workcoarse, 0, NULL);
647:   return 0;
648: }

650: static PetscErrorCode PCReset_Deflation(PC pc)
651: {
652:   PC_Deflation *def = (PC_Deflation *)pc->data;

654:   VecDestroy(&def->work);
655:   VecDestroyVecs(2, &def->workcoarse);
656:   MatDestroy(&def->W);
657:   MatDestroy(&def->Wt);
658:   MatDestroy(&def->WtA);
659:   MatDestroy(&def->WtAW);
660:   KSPDestroy(&def->WtAWinv);
661:   PCDestroy(&def->pc);
662:   return 0;
663: }

665: static PetscErrorCode PCDestroy_Deflation(PC pc)
666: {
667:   PCReset_Deflation(pc);
668:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetInitOnly_C", NULL);
669:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetLevels_C", NULL);
670:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetReductionFactor_C", NULL);
671:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetCorrectionFactor_C", NULL);
672:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetSpaceToCompute_C", NULL);
673:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetSpace_C", NULL);
674:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetProjectionNullSpaceMat_C", NULL);
675:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetCoarseMat_C", NULL);
676:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationGetCoarseKSP_C", NULL);
677:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationGetPC_C", NULL);
678:   PetscFree(pc->data);
679:   return 0;
680: }

682: static PetscErrorCode PCView_Deflation(PC pc, PetscViewer viewer)
683: {
684:   PC_Deflation *def = (PC_Deflation *)pc->data;
685:   PetscInt      its;
686:   PetscBool     iascii;

688:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
689:   if (iascii) {
690:     if (def->correct) PetscViewerASCIIPrintf(viewer, "using CP correction, factor = %g+%gi\n", (double)PetscRealPart(def->correctfact), (double)PetscImaginaryPart(def->correctfact));
691:     if (!def->lvl) PetscViewerASCIIPrintf(viewer, "deflation space type: %s\n", PCDeflationSpaceTypes[def->spacetype]);

693:     PetscViewerASCIIPrintf(viewer, "--- Additional PC:\n");
694:     PetscViewerASCIIPushTab(viewer);
695:     PCView(def->pc, viewer);
696:     PetscViewerASCIIPopTab(viewer);

698:     PetscViewerASCIIPrintf(viewer, "--- Coarse problem solver:\n");
699:     PetscViewerASCIIPushTab(viewer);
700:     KSPGetTotalIterations(def->WtAWinv, &its);
701:     PetscViewerASCIIPrintf(viewer, "total number of iterations: %" PetscInt_FMT "\n", its);
702:     KSPView(def->WtAWinv, viewer);
703:     PetscViewerASCIIPopTab(viewer);
704:   }
705:   return 0;
706: }

708: static PetscErrorCode PCSetFromOptions_Deflation(PC pc, PetscOptionItems *PetscOptionsObject)
709: {
710:   PC_Deflation *def = (PC_Deflation *)pc->data;

712:   PetscOptionsHeadBegin(PetscOptionsObject, "Deflation options");
713:   PetscOptionsBool("-pc_deflation_init_only", "Use only initialization step - Initdef", "PCDeflationSetInitOnly", def->init, &def->init, NULL);
714:   PetscOptionsInt("-pc_deflation_levels", "Maximum of deflation levels", "PCDeflationSetLevels", def->maxlvl, &def->maxlvl, NULL);
715:   PetscOptionsInt("-pc_deflation_reduction_factor", "Reduction factor for coarse problem solution using PCTELESCOPE", "PCDeflationSetReductionFactor", def->reductionfact, &def->reductionfact, NULL);
716:   PetscOptionsBool("-pc_deflation_correction", "Add coarse problem correction Q to P", "PCDeflationSetCorrectionFactor", def->correct, &def->correct, NULL);
717:   PetscOptionsScalar("-pc_deflation_correction_factor", "Set multiple of Q to use as coarse problem correction", "PCDeflationSetCorrectionFactor", def->correctfact, &def->correctfact, NULL);
718:   PetscOptionsEnum("-pc_deflation_compute_space", "Compute deflation space", "PCDeflationSetSpace", PCDeflationSpaceTypes, (PetscEnum)def->spacetype, (PetscEnum *)&def->spacetype, NULL);
719:   PetscOptionsInt("-pc_deflation_compute_space_size", "Set size of the deflation space to compute", "PCDeflationSetSpace", def->spacesize, &def->spacesize, NULL);
720:   PetscOptionsBool("-pc_deflation_space_extend", "Extend deflation space instead of truncating (wavelets)", "PCDeflation", def->extendsp, &def->extendsp, NULL);
721:   PetscOptionsHeadEnd();
722:   return 0;
723: }

725: /*MC
726:    PCDEFLATION - Deflation preconditioner shifts (deflates) part of the spectrum to zero or to a predefined value.

728:    Options Database Keys:
729: +    -pc_deflation_init_only          <false> - if true computes only the special guess
730: .    -pc_deflation_max_lvl            <0>     - maximum number of levels for multilevel deflation
731: .    -pc_deflation_reduction_factor <\-1>     - reduction factor on bottom level coarse problem for PCTELESCOPE (default based on the size of the coarse problem)
732: .    -pc_deflation_correction         <false> - if true apply coarse problem correction
733: .    -pc_deflation_correction_factor  <1.0>   - sets coarse problem correction factor
734: .    -pc_deflation_compute_space      <haar>  - compute PCDeflationSpaceType deflation space
735: -    -pc_deflation_compute_space_size <1>     - size of the deflation space (corresponds to number of levels for wavelet-based deflation)

737:    Notes:
738:     Given a (complex - transpose is always Hermitian) full rank deflation matrix W, the deflation (introduced in [1,2])
739:     preconditioner uses projections Q = W*(W'*A*W)^{-1}*W' and P = I - Q*A, where A is pmat.

741:     The deflation computes initial guess x0 = x_{-1} - Q*r_{-1}, which is the solution on the deflation space.
742:     If `PCDeflationSetInitOnly()` or -pc_deflation_init_only is set to `PETSC_TRUE` (InitDef scheme), the application of the
743:     preconditioner consists only of application of the additional preconditioner M^{-1}. Otherwise, the preconditioner
744:     application consists of P*M^{-1} + factor*Q. The first part of the preconditioner (PM^{-1}) shifts some eigenvalues
745:     to zero while the addition of the coarse problem correction (factor*Q) makes the preconditioner to shift some
746:     eigenvalues to the given factor. The InitDef scheme is recommended for deflation using high accuracy estimates
747:     of eigenvectors of A when it exhibits similar convergence to the full deflation but is cheaper.

749:     The deflation matrix is by default automatically computed. The type of deflation matrix and its size to compute can
750:     be controlled by `PCDeflationSetSpaceToCompute()` or -pc_deflation_compute_space and -pc_deflation_compute_space_size.
751:     User can set an arbitrary deflation space matrix with `PCDeflationSetSpace()`. If the deflation matrix
752:     is a multiplicative `MATCOMPOSITE`, a multilevel deflation [3] is used. The first matrix in the composite is used as the
753:     deflation matrix, and the coarse problem (W'*A*W)^{-1} is solved by `KSPFCG` (if A is `MAT_SPD`) or `KSPFGMRES` preconditioned
754:     by deflation with deflation matrix being the next matrix in the `MATCOMPOSITE`. This scheme repeats until the maximum
755:     level is reached or there are no more matrices. If the maximum level is reached, the remaining matrices are merged
756:     (multiplied) to create the last deflation matrix. The maximum level defaults to 0 and can be set by
757:     `PCDeflationSetLevels()` or by -pc_deflation_levels.

759:     The coarse problem `KSP` can be controlled from the command line with prefix -deflation_ for the first level and -deflation_[lvl-1]
760:     from the second level onward. You can also use
761:     `PCDeflationGetCoarseKSP()` to control it from code. The bottom level KSP defaults to
762:     `KSPPREONLY` with `PCLU` direct solver (`MATSOLVERSUPERLU`/`MATSOLVERSUPERLU_DIST` if available) wrapped into `PCTELESCOPE`.
763:     For convenience, the reduction factor can be set by `PCDeflationSetReductionFactor()`
764:     or -pc_deflation_recduction_factor. The default is chosen heuristically based on the coarse problem size.

766:     The additional preconditioner can be controlled from command line with prefix -deflation_[lvl]_pc (same rules used for
767:     coarse problem `KSP` apply for [lvl]_ part of prefix), e.g., -deflation_1_pc_pc_type bjacobi. You can also use
768:     `PCDeflationGetPC()` to control the additional preconditioner from code. It defaults to `PCNONE`.

770:     The coarse problem correction term (factor*Q) can be turned on by -pc_deflation_correction and the factor value can
771:     be set by pc_deflation_correction_factor or by `PCDeflationSetCorrectionFactor()`. The coarse problem can
772:     significantly improve convergence when the deflation coarse problem is not solved with high enough accuracy. We
773:     recommend setting factor to some eigenvalue, e.g., the largest eigenvalue so that the preconditioner does not create
774:     an isolated eigenvalue.

776:     The options are automatically inherited from the previous deflation level.

778:     The preconditioner supports `KSPMonitorDynamicTolerance()`. This is useful for the multilevel scheme for which we also
779:     recommend limiting the number of iterations for the coarse problems.

781:     See section 3 of [4] for additional references and decription of the algorithm when used for conjugate gradients.
782:     Section 4 describes some possible choices for the deflation space.

784:      Contributed by Jakub Kruzik (PERMON), Institute of Geonics of the Czech
785:      Academy of Sciences and VSB - TU Ostrava.

787:      Developed from PERMON code used in [4] while on a research stay with
788:      Prof. Reinhard Nabben at the Institute of Mathematics, TU Berlin.

790:    References:
791: +  * - A. Nicolaides. "Deflation of conjugate gradients with applications to boundary value problems", SIAM J. Numer. Anal. 24.2, 1987.
792: .  * - Z. Dostal. "Conjugate gradient method with preconditioning by projector", Int J. Comput. Math. 23.3-4, 1988.
793: .  * - Y. A. Erlangga and R. Nabben. "Multilevel Projection-Based Nested Krylov Iteration for Boundary Value Problems", SIAM J. Sci. Comput. 30.3, 2008.
794: -  * - J. Kruzik "Implementation of the Deflated Variants of the Conjugate Gradient Method", Master's thesis, VSB-TUO, 2018 - http://dspace5.vsb.cz/bitstream/handle/10084/130303/KRU0097_USP_N2658_2612T078_2018.pdf

796:    Level: intermediate

798: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
799:           `PCDeflationSetInitOnly()`, `PCDeflationSetLevels()`, `PCDeflationSetReductionFactor()`,
800:           `PCDeflationSetCorrectionFactor()`, `PCDeflationSetSpaceToCompute()`,
801:           `PCDeflationSetSpace()`, `PCDeflationSpaceType`, `PCDeflationSetProjectionNullSpaceMat()`,
802:           `PCDeflationSetCoarseMat()`, `PCDeflationGetCoarseKSP()`, `PCDeflationGetPC()`
803: M*/

805: PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc)
806: {
807:   PC_Deflation *def;

809:   PetscNew(&def);
810:   pc->data = (void *)def;

812:   def->init          = PETSC_FALSE;
813:   def->correct       = PETSC_FALSE;
814:   def->correctfact   = 1.0;
815:   def->reductionfact = -1;
816:   def->spacetype     = PC_DEFLATION_SPACE_HAAR;
817:   def->spacesize     = 1;
818:   def->extendsp      = PETSC_FALSE;
819:   def->lvl           = 0;
820:   def->maxlvl        = 0;
821:   def->W             = NULL;
822:   def->Wt            = NULL;

824:   pc->ops->apply          = PCApply_Deflation;
825:   pc->ops->presolve       = PCPreSolve_Deflation;
826:   pc->ops->setup          = PCSetUp_Deflation;
827:   pc->ops->reset          = PCReset_Deflation;
828:   pc->ops->destroy        = PCDestroy_Deflation;
829:   pc->ops->setfromoptions = PCSetFromOptions_Deflation;
830:   pc->ops->view           = PCView_Deflation;

832:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetInitOnly_C", PCDeflationSetInitOnly_Deflation);
833:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetLevels_C", PCDeflationSetLevels_Deflation);
834:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetReductionFactor_C", PCDeflationSetReductionFactor_Deflation);
835:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetCorrectionFactor_C", PCDeflationSetCorrectionFactor_Deflation);
836:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetSpaceToCompute_C", PCDeflationSetSpaceToCompute_Deflation);
837:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetSpace_C", PCDeflationSetSpace_Deflation);
838:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetProjectionNullSpaceMat_C", PCDeflationSetProjectionNullSpaceMat_Deflation);
839:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetCoarseMat_C", PCDeflationSetCoarseMat_Deflation);
840:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationGetCoarseKSP_C", PCDeflationGetCoarseKSP_Deflation);
841:   PetscObjectComposeFunction((PetscObject)pc, "PCDeflationGetPC_C", PCDeflationGetPC_Deflation);
842:   return 0;
843: }