Actual source code: deflation.c

petsc-3.13.6 2020-09-29
Report Typos and Errors
  1:  #include <../src/ksp/pc/impls/deflation/deflation.h>

  3: const char *const PCDeflationSpaceTypes[] = {
  4:   "haar",
  5:   "db2",
  6:   "db4",
  7:   "db8",
  8:   "db16",
  9:   "biorth22",
 10:   "meyer",
 11:   "aggregation",
 12:   "user",
 13:   "PCDeflationSpaceType",
 14:   "PC_DEFLATION_SPACE_",
 15:   0
 16: };

 18: static PetscErrorCode PCDeflationSetInitOnly_Deflation(PC pc,PetscBool flg)
 19: {
 20:   PC_Deflation   *def = (PC_Deflation*)pc->data;

 23:   def->init = flg;
 24:   return(0);
 25: }

 27: /*@
 28:    PCDeflationSetInitOnly - Do only initialization step.
 29:     Sets initial guess to the solution on the deflation space but does not apply
 30:     the deflation preconditioner. The additional preconditioner is still applied.

 32:    Logically Collective

 34:    Input Parameters:
 35: +  pc  - the preconditioner context
 36: -  flg - default PETSC_FALSE

 38:    Options Database Keys:
 39: .    -pc_deflation_init_only <false> - if true computes only the special guess

 41:    Level: intermediate

 43: .seealso: PCDEFLATION
 44: @*/
 45: PetscErrorCode PCDeflationSetInitOnly(PC pc,PetscBool flg)
 46: {

 52:   PetscTryMethod(pc,"PCDeflationSetInitOnly_C",(PC,PetscBool),(pc,flg));
 53:   return(0);
 54: }


 57: static PetscErrorCode PCDeflationSetLevels_Deflation(PC pc,PetscInt current,PetscInt max)
 58: {
 59:   PC_Deflation   *def = (PC_Deflation*)pc->data;

 62:   if (current) def->lvl = current;
 63:   def->maxlvl = max;
 64:   return(0);
 65: }

 67: /*@
 68:    PCDeflationSetLevels - Set the maximum level of deflation nesting.

 70:    Logically Collective

 72:    Input Parameters:
 73: +  pc  - the preconditioner context
 74: -  max - maximum deflation level

 76:    Options Database Keys:
 77: .    -pc_deflation_max_lvl <0> - maximum number of levels for multilevel deflation

 79:    Level: intermediate

 81: .seealso: PCDeflationSetSpaceToCompute(), PCDeflationSetSpace(), PCDEFLATION
 82: @*/
 83: PetscErrorCode PCDeflationSetLevels(PC pc,PetscInt max)
 84: {

 90:   PetscTryMethod(pc,"PCDeflationSetLevels_C",(PC,PetscInt,PetscInt),(pc,0,max));
 91:   return(0);
 92: }

 94: static PetscErrorCode PCDeflationSetReductionFactor_Deflation(PC pc,PetscInt red)
 95: {
 96:   PC_Deflation   *def = (PC_Deflation*)pc->data;

 99:   def->reductionfact = red;
100:   return(0);
101: }

103: /*@
104:    PCDeflationSetReductionFactor - Set reduction factor for the bottom PCTELESCOPE coarse problem solver.

106:    Logically Collective

108:    Input Parameters:
109: +  pc  - the preconditioner context
110: -  red - reduction factor (or PETSC_DETERMINE)

112:    Options Database Keys:
113: .    -pc_deflation_reduction_factor <\-1> - reduction factor on bottom level coarse problem for PCTELESCOPE

115:    Notes:
116:      Default is computed based on the size of the coarse problem.

118:    Level: intermediate

120: .seealso: PCTELESCOPE, PCDEFLATION
121: @*/
122: PetscErrorCode PCDeflationSetReductionFactor(PC pc,PetscInt red)
123: {

129:   PetscTryMethod(pc,"PCDeflationSetReductionFactor_C",(PC,PetscInt),(pc,red));
130:   return(0);
131: }

133: static PetscErrorCode PCDeflationSetCorrectionFactor_Deflation(PC pc,PetscScalar fact)
134: {
135:   PC_Deflation   *def = (PC_Deflation*)pc->data;

138:   /* TODO PETSC_DETERMINE -> compute max eigenvalue with power method */
139:   def->correct = PETSC_TRUE;
140:   def->correctfact = fact;
141:   if (def->correct == 0.0) {
142:     def->correct = PETSC_FALSE;
143:   }
144:   return(0);
145: }

147: /*@
148:    PCDeflationSetCorrectionFactor - Set coarse problem correction factor.
149:     The Preconditioner becomes P*M^{-1} + fact*Q.

151:    Logically Collective

153:    Input Parameters:
154: +  pc   - the preconditioner context
155: -  fact - correction factor

157:    Options Database Keys:
158: +    -pc_deflation_correction        <false> - if true apply coarse problem correction
159: -    -pc_deflation_correction_factor <1.0>   - sets coarse problem correction factor

161:    Notes:
162:     Any non-zero fact enables the coarse problem correction.

164:    Level: intermediate

166: .seealso: PCDEFLATION
167: @*/
168: PetscErrorCode PCDeflationSetCorrectionFactor(PC pc,PetscScalar fact)
169: {

175:   PetscTryMethod(pc,"PCDeflationSetCorrectionFactor_C",(PC,PetscScalar),(pc,fact));
176:   return(0);
177: }

179: static PetscErrorCode PCDeflationSetSpaceToCompute_Deflation(PC pc,PCDeflationSpaceType type,PetscInt size)
180: {
181:   PC_Deflation   *def = (PC_Deflation*)pc->data;

184:   if (type) def->spacetype = type;
185:   if (size > 0) def->spacesize = size;
186:   return(0);
187: }

189: /*@
190:    PCDeflationSetSpaceToCompute - Set deflation space type and size to compute.

192:    Logically Collective

194:    Input Parameters:
195: +  pc   - the preconditioner context
196: .  type - deflation space type to compute (or PETSC_IGNORE)
197: -  size - size of the space to compute (or PETSC_DEFAULT)

199:    Options Database Keys:
200: +    -pc_deflation_compute_space      <haar> - compute PCDeflationSpaceType deflation space
201: -    -pc_deflation_compute_space_size <1>    - size of the deflation space

203:    Notes:
204:     For wavelet-based deflation, size represents number of levels.

206:     The deflation space is computed in PCSetUp().

208:    Level: intermediate

210: .seealso: PCDeflationSetLevels(), PCDEFLATION
211: @*/
212: PetscErrorCode PCDeflationSetSpaceToCompute(PC pc,PCDeflationSpaceType type,PetscInt size)
213: {

220:   PetscTryMethod(pc,"PCDeflationSetSpaceToCompute_C",(PC,PCDeflationSpaceType,PetscInt),(pc,type,size));
221:   return(0);
222: }

224: static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc,Mat W,PetscBool transpose)
225: {
226:   PC_Deflation   *def = (PC_Deflation*)pc->data;

230:   /* possibly allows W' = Wt (which is valid but not tested) */
231:   PetscObjectReference((PetscObject)W);
232:   if (transpose) {
233:     MatDestroy(&def->Wt);
234:     def->Wt = W;
235:   } else {
236:     MatDestroy(&def->W);
237:     def->W = W;
238:   }
239:   return(0);
240: }

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

245:    Logically Collective

247:    Input Parameters:
248: +  pc        - the preconditioner context
249: .  W         - deflation matrix
250: -  transpose - indicates that W is an explicit transpose of the deflation matrix

252:    Notes:
253:     Setting W as a multipliplicative MATCOMPOSITE enables use of the multilevel
254:     deflation. If W = W0*W1*W2*...*Wn, W0 is taken as the first deflation space and
255:     the coarse problem (W0'*A*W0)^{-1} is again preconditioned by deflation with
256:     W1 as the deflation matrix. This repeats until the maximum level set by
257:     PCDeflationSetLevels() is reached or there are no more matrices available.
258:     If there are matrices left after reaching the maximum level,
259:     they are merged into a deflation matrix ...*W{n-1}*Wn.

261:    Level: intermediate

263: .seealso: PCDeflationSetLevels(), PCDEFLATION
264: @*/
265: PetscErrorCode PCDeflationSetSpace(PC pc,Mat W,PetscBool transpose)
266: {

273:   PetscTryMethod(pc,"PCDeflationSetSpace_C",(PC,Mat,PetscBool),(pc,W,transpose));
274:   return(0);
275: }

277: static PetscErrorCode PCDeflationSetProjectionNullSpaceMat_Deflation(PC pc,Mat mat)
278: {
279:   PC_Deflation     *def = (PC_Deflation*)pc->data;
280:   PetscErrorCode   ierr;

283:   PetscObjectReference((PetscObject)mat);
284:   MatDestroy(&def->WtA);
285:   def->WtA = mat;
286:   PetscLogObjectParent((PetscObject)pc,(PetscObject)def->WtA);
287:   return(0);
288: }

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

293:    Collective

295:    Input Parameters:
296: +  pc  - preconditioner context
297: -  mat - projection null space matrix

299:    Level: developer

301: .seealso: PCDEFLATION
302: @*/
303: PetscErrorCode  PCDeflationSetProjectionNullSpaceMat(PC pc,Mat mat)
304: {

310:   PetscTryMethod(pc,"PCDeflationSetProjectionNullSpaceMat_C",(PC,Mat),(pc,mat));
311:   return(0);
312: }

314: static PetscErrorCode PCDeflationSetCoarseMat_Deflation(PC pc,Mat mat)
315: {
316:   PC_Deflation     *def = (PC_Deflation*)pc->data;
317:   PetscErrorCode   ierr;

320:   PetscObjectReference((PetscObject)mat);
321:   MatDestroy(&def->WtAW);
322:   def->WtAW = mat;
323:   PetscLogObjectParent((PetscObject)pc,(PetscObject)def->WtAW);
324:   return(0);
325: }

327: /*@
328:    PCDeflationSetCoarseMat - Set the coarse problem Mat.

330:    Collective

332:    Input Parameters:
333: +  pc  - preconditioner context
334: -  mat - coarse problem mat

336:    Level: developer

338: .seealso: PCDEFLATION
339: @*/
340: PetscErrorCode  PCDeflationSetCoarseMat(PC pc,Mat mat)
341: {

347:   PetscTryMethod(pc,"PCDeflationSetCoarseMat_C",(PC,Mat),(pc,mat));
348:   return(0);
349: }

351: static PetscErrorCode PCDeflationGetCoarseKSP_Deflation(PC pc,KSP *ksp)
352: {
353:   PC_Deflation     *def = (PC_Deflation*)pc->data;

356:   *ksp = def->WtAWinv;
357:   return(0);
358: }

360: /*@
361:    PCDeflationGetCoarseKSP - Returns the coarse problem KSP.

363:    Not Collective

365:    Input Parameters:
366: .  pc - preconditioner context

368:    Output Parameters:
369: .  ksp - coarse problem KSP context

371:    Level: advanced

373: .seealso: PCDEFLATION
374: @*/
375: PetscErrorCode  PCDeflationGetCoarseKSP(PC pc,KSP *ksp)
376: {

382:   PetscTryMethod(pc,"PCDeflationGetCoarseKSP_C",(PC,KSP*),(pc,ksp));
383:   return(0);
384: }

386: static PetscErrorCode PCDeflationGetPC_Deflation(PC pc,PC *apc)
387: {
388:   PC_Deflation   *def = (PC_Deflation*)pc->data;

391:   *apc = def->pc;
392:   return(0);
393: }

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

398:    Not Collective

400:    Input Parameters:
401: .  pc  - the preconditioner context

403:    Output Parameters:
404: .  apc - additional preconditioner

406:    Level: advanced

408: .seealso: PCDEFLATION
409: @*/
410: PetscErrorCode PCDeflationGetPC(PC pc,PC *apc)
411: {

417:   PetscTryMethod(pc,"PCDeflationGetPC_C",(PC,PC*),(pc,apc));
418:   return(0);
419: }

421: /*
422:   x <- x + W*(W'*A*W)^{-1}*W'*r  = x + Q*r
423: */
424: static PetscErrorCode PCPreSolve_Deflation(PC pc,KSP ksp,Vec b, Vec x)
425: {
426:   PC_Deflation     *def = (PC_Deflation*)pc->data;
427:   Mat              A;
428:   Vec              r,w1,w2;
429:   PetscBool        nonzero;
430:   PetscErrorCode   ierr;

433:   w1 = def->workcoarse[0];
434:   w2 = def->workcoarse[1];
435:   r  = def->work;
436:   PCGetOperators(pc,NULL,&A);

438:   KSPGetInitialGuessNonzero(ksp,&nonzero);
439:   KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
440:   if (nonzero) {
441:     MatMult(A,x,r);                          /*    r  <- b - Ax              */
442:     VecAYPX(r,-1.0,b);
443:   } else {
444:     VecCopy(b,r);                            /*    r  <- b (x is 0)          */
445:   }

447:   if (def->Wt) {
448:     MatMult(def->Wt,r,w1);                   /*    w1 <- W'*r                */
449:   } else {
450:     MatMultHermitianTranspose(def->W,r,w1);  /*    w1 <- W'*r                */
451:   }
452:   KSPSolve(def->WtAWinv,w1,w2);              /*    w2 <- (W'*A*W)^{-1}*w1    */
453:   MatMult(def->W,w2,r);                      /*    r  <- W*w2                */
454:   VecAYPX(x,1.0,r);
455:   return(0);
456: }

458: /*
459:   if (def->correct) {
460:     z <- M^{-1}r - W*(W'*A*W)^{-1}*(W'*A*M^{-1}r - l*W'*r) = (P*M^{-1} + l*Q)*r
461:   } else {
462:     z <- M^{-1}*r - W*(W'*A*W)^{-1}*W'*A*M{-1}*r = P*M^{-1}*r
463:   }
464: */
465: static PetscErrorCode PCApply_Deflation(PC pc,Vec r,Vec z)
466: {
467:   PC_Deflation     *def = (PC_Deflation*)pc->data;
468:   Mat              A;
469:   Vec              u,w1,w2;
470:   PetscErrorCode   ierr;

473:   w1 = def->workcoarse[0];
474:   w2 = def->workcoarse[1];
475:   u  = def->work;
476:   PCGetOperators(pc,NULL,&A);

478:   PCApply(def->pc,r,z);                          /*    z <- M^{-1}*r             */
479:   if (!def->init) {
480:     MatMult(def->WtA,z,w1);                      /*    w1 <- W'*A*z              */
481:     if (def->correct) {
482:       if (def->Wt) {
483:         MatMult(def->Wt,r,w2);                   /*    w2 <- W'*r                */
484:       } else {
485:         MatMultHermitianTranspose(def->W,r,w2);  /*    w2 <- W'*r                */
486:       }
487:       VecAXPY(w1,-1.0*def->correctfact,w2);      /*    w1 <- w1 - l*w2           */
488:     }
489:     KSPSolve(def->WtAWinv,w1,w2);                /*    w2 <- (W'*A*W)^{-1}*w1    */
490:     MatMult(def->W,w2,u);                        /*    u  <- W*w2                */
491:     VecAXPY(z,-1.0,u);                           /*    z  <- z - u               */
492:   }
493:   return(0);
494: }

496: static PetscErrorCode PCSetUp_Deflation(PC pc)
497: {
498:   PC_Deflation     *def = (PC_Deflation*)pc->data;
499:   KSP              innerksp;
500:   PC               pcinner;
501:   Mat              Amat,nextDef=NULL,*mats;
502:   PetscInt         i,m,red,size;
503:   PetscMPIInt      commsize;
504:   PetscBool        match,flgspd,transp=PETSC_FALSE;
505:   MatCompositeType ctype;
506:   MPI_Comm         comm;
507:   char             prefix[128]="";
508:   PetscErrorCode   ierr;

511:   if (pc->setupcalled) return(0);
512:   PetscObjectGetComm((PetscObject)pc,&comm);
513:   PCGetOperators(pc,NULL,&Amat);
514:   if (!def->lvl && !def->prefix) {
515:     PCGetOptionsPrefix(pc,&def->prefix);
516:   }
517:   if (def->lvl) {
518:     PetscSNPrintf(prefix,sizeof(prefix),"%d_",(int)def->lvl);
519:   }

521:   /* compute a deflation space */
522:   if (def->W || def->Wt) {
523:     def->spacetype = PC_DEFLATION_SPACE_USER;
524:   } else {
525:     PCDeflationComputeSpace(pc);
526:   }

528:   /* nested deflation */
529:   if (def->W) {
530:     PetscObjectTypeCompare((PetscObject)def->W,MATCOMPOSITE,&match);
531:     if (match) {
532:       MatCompositeGetType(def->W,&ctype);
533:       MatCompositeGetNumberMat(def->W,&size);
534:     }
535:   } else {
536:     MatCreateHermitianTranspose(def->Wt,&def->W);
537:     PetscObjectTypeCompare((PetscObject)def->Wt,MATCOMPOSITE,&match);
538:     if (match) {
539:       MatCompositeGetType(def->Wt,&ctype);
540:       MatCompositeGetNumberMat(def->Wt,&size);
541:     }
542:     transp = PETSC_TRUE;
543:   }
544:   if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) {
545:     if (!transp) {
546:       if (def->lvl < def->maxlvl) {
547:         PetscMalloc1(size,&mats);
548:         for (i=0; i<size; i++) {
549:           MatCompositeGetMat(def->W,i,&mats[i]);
550:         }
551:         size -= 1;
552:         MatDestroy(&def->W);
553:         def->W = mats[size];
554:         PetscObjectReference((PetscObject)mats[size]);
555:         if (size > 1) {
556:           MatCreateComposite(comm,size,mats,&nextDef);
557:           MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);
558:         } else {
559:           nextDef = mats[0];
560:           PetscObjectReference((PetscObject)mats[0]);
561:         }
562:         PetscFree(mats);
563:       } else {
564:         /* TODO test merge side performance */
565:         /* MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT); */
566:         MatCompositeMerge(def->W);
567:       }
568:     } else {
569:       if (def->lvl < def->maxlvl) {
570:         PetscMalloc1(size,&mats);
571:         for (i=0; i<size; i++) {
572:           MatCompositeGetMat(def->Wt,i,&mats[i]);
573:         }
574:         size -= 1;
575:         MatDestroy(&def->Wt);
576:         def->Wt = mats[0];
577:         PetscObjectReference((PetscObject)mats[0]);
578:         if (size > 1) {
579:           MatCreateComposite(comm,size,&mats[1],&nextDef);
580:           MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);
581:         } else {
582:           nextDef = mats[1];
583:           PetscObjectReference((PetscObject)mats[1]);
584:         }
585:         PetscFree(mats);
586:       } else {
587:         /* MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT); */
588:         MatCompositeMerge(def->Wt);
589:       }
590:     }
591:   }

593:   if (transp) {
594:     MatDestroy(&def->W);
595:     MatHermitianTranspose(def->Wt,MAT_INITIAL_MATRIX,&def->W);
596:   }

598:   /* assemble WtA */
599:   if (!def->WtA) {
600:     if (def->Wt) {
601:       MatMatMult(def->Wt,Amat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtA);
602:     } else {
603: #if defined(PETSC_USE_COMPLEX)
604:       MatHermitianTranspose(def->W,MAT_INITIAL_MATRIX,&def->Wt);
605:       MatMatMult(def->Wt,Amat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtA);
606: #else
607:       MatTransposeMatMult(def->W,Amat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtA);
608: #endif
609:     }
610:   }
611:   /* setup coarse problem */
612:   if (!def->WtAWinv) {
613:     MatGetSize(def->W,NULL,&m);
614:     if (!def->WtAW) {
615:       MatMatMult(def->WtA,def->W,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);
616:       /* TODO create MatInheritOption(Mat,MatOption) */
617:       MatGetOption(Amat,MAT_SPD,&flgspd);
618:       MatSetOption(def->WtAW,MAT_SPD,flgspd);
619: #if defined(PETSC_USE_DEBUG)
620:       /* Check columns of W are not in kernel of A */
621:       {
622:         PetscReal *norms;
623:         PetscMalloc1(m,&norms);
624:         MatGetColumnNorms(def->WtAW,NORM_INFINITY,norms);
625:         for (i=0; i<m; i++) {
626:           if (norms[i] < 100*PETSC_MACHINE_EPSILON) {
627:             SETERRQ1(comm,PETSC_ERR_SUP,"Column %D of W is in kernel of A.",i);
628:           }
629:         }
630:         PetscFree(norms);
631:       }
632: #endif
633:     } else {
634:       MatGetOption(def->WtAW,MAT_SPD,&flgspd);
635:     }
636:     /* TODO use MATINV ? */
637:     KSPCreate(comm,&def->WtAWinv);
638:     KSPSetOperators(def->WtAWinv,def->WtAW,def->WtAW);
639:     KSPGetPC(def->WtAWinv,&pcinner);
640:     /* Setup KSP and PC */
641:     if (nextDef) { /* next level for multilevel deflation */
642:       innerksp = def->WtAWinv;
643:       /* set default KSPtype */
644:       if (!def->ksptype) {
645:         def->ksptype = KSPFGMRES;
646:         if (flgspd) { /* SPD system */
647:           def->ksptype = KSPFCG;
648:         }
649:       }
650:       KSPSetType(innerksp,def->ksptype); /* TODO iherit from KSP + tolerances */
651:       PCSetType(pcinner,PCDEFLATION); /* TODO create coarse preconditinoner M_c = WtMW ? */
652:       PCDeflationSetSpace(pcinner,nextDef,transp);
653:       PCDeflationSetLevels_Deflation(pcinner,def->lvl+1,def->maxlvl);
654:       /* inherit options */
655:       if (def->prefix) ((PC_Deflation*)(pcinner->data))->prefix = def->prefix;
656:       ((PC_Deflation*)(pcinner->data))->init          = def->init;
657:       ((PC_Deflation*)(pcinner->data))->ksptype       = def->ksptype;
658:       ((PC_Deflation*)(pcinner->data))->correct       = def->correct;
659:       ((PC_Deflation*)(pcinner->data))->correctfact   = def->correctfact;
660:       ((PC_Deflation*)(pcinner->data))->reductionfact = def->reductionfact;
661:       MatDestroy(&nextDef);
662:     } else { /* the last level */
663:       KSPSetType(def->WtAWinv,KSPPREONLY);
664:       PCSetType(pcinner,PCTELESCOPE);
665:       /* do not overwrite PCTELESCOPE */
666:       if (def->prefix) {
667:         KSPSetOptionsPrefix(def->WtAWinv,def->prefix);
668:       }
669:       KSPAppendOptionsPrefix(def->WtAWinv,"deflation_tel_");
670:       PCSetFromOptions(pcinner);
671:       PetscObjectTypeCompare((PetscObject)pcinner,PCTELESCOPE,&match);
672:       if (!match) SETERRQ(comm,PETSC_ERR_SUP,"User can not owerwrite PCTELESCOPE on bottom level, use reduction factor = 1 instead.");
673:       /* Reduction factor choice */
674:       red = def->reductionfact;
675:       if (red < 0) {
676:         MPI_Comm_size(comm,&commsize);
677:         red  = ceil((float)commsize/ceil((float)m/commsize));
678:         PetscObjectTypeCompareAny((PetscObject)(def->WtAW),&match,MATSEQDENSE,MATMPIDENSE,MATDENSE,"");
679:         if (match) red = commsize;
680:         PetscInfo1(pc,"Auto choosing reduction factor %D\n",red);
681:       }
682:       PCTelescopeSetReductionFactor(pcinner,red);
683:       PCSetUp(pcinner);
684:       PCTelescopeGetKSP(pcinner,&innerksp);
685:       if (innerksp) {
686:         KSPGetPC(innerksp,&pcinner);
687:         PCSetType(pcinner,PCLU);
688: #if defined(PETSC_HAVE_SUPERLU)
689:         MatGetFactorAvailable(def->WtAW,MATSOLVERSUPERLU,MAT_FACTOR_LU,&match);
690:         if (match) {
691:           PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU);
692:         }
693: #endif
694: #if defined(PETSC_HAVE_SUPERLU_DIST)
695:         MatGetFactorAvailable(def->WtAW,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&match);
696:         if (match) {
697:           PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU_DIST);
698:         }
699: #endif
700:       }
701:     }

703:     if (innerksp) {
704:       if (def->prefix) {
705:         KSPSetOptionsPrefix(innerksp,def->prefix);
706:         KSPAppendOptionsPrefix(innerksp,"deflation_");
707:       } else {
708:         KSPSetOptionsPrefix(innerksp,"deflation_");
709:       }
710:       KSPAppendOptionsPrefix(innerksp,prefix);
711:       KSPSetFromOptions(innerksp);
712:       KSPSetUp(innerksp);
713:     }
714:   }
715:   KSPSetFromOptions(def->WtAWinv);
716:   KSPSetUp(def->WtAWinv);

718:   /* create preconditioner */
719:   if (!def->pc) {
720:     PCCreate(comm,&def->pc);
721:     PCSetOperators(def->pc,Amat,Amat);
722:     PCSetType(def->pc,PCNONE);
723:     if (def->prefix) {
724:       PCSetOptionsPrefix(def->pc,def->prefix);
725:     }
726:     PCAppendOptionsPrefix(def->pc,"deflation_");
727:     PCAppendOptionsPrefix(def->pc,prefix);
728:     PCAppendOptionsPrefix(def->pc,"pc_");
729:     PCSetFromOptions(def->pc);
730:     PCSetUp(def->pc);
731:   }

733:   /* create work vecs */
734:   MatCreateVecs(Amat,NULL,&def->work);
735:   KSPCreateVecs(def->WtAWinv,2,&def->workcoarse,0,NULL);
736:   return(0);
737: }

739: static PetscErrorCode PCReset_Deflation(PC pc)
740: {
741:   PC_Deflation      *def = (PC_Deflation*)pc->data;
742:   PetscErrorCode    ierr;

745:   VecDestroy(&def->work);
746:   VecDestroyVecs(2,&def->workcoarse);
747:   MatDestroy(&def->W);
748:   MatDestroy(&def->Wt);
749:   MatDestroy(&def->WtA);
750:   MatDestroy(&def->WtAW);
751:   KSPDestroy(&def->WtAWinv);
752:   PCDestroy(&def->pc);
753:   return(0);
754: }

756: static PetscErrorCode PCDestroy_Deflation(PC pc)
757: {

761:   PCReset_Deflation(pc);
762:   PetscFree(pc->data);
763:   return(0);
764: }

766: static PetscErrorCode PCView_Deflation(PC pc,PetscViewer viewer)
767: {
768:   PC_Deflation      *def = (PC_Deflation*)pc->data;
769:   PetscInt          its;
770:   PetscBool         iascii;
771:   PetscErrorCode    ierr;

774:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
775:   if (iascii) {
776:     if (def->correct) {
777:       PetscViewerASCIIPrintf(viewer,"using CP correction, factor = %g+%gi\n",
778:                                     (double)PetscRealPart(def->correctfact),
779:                                     (double)PetscImaginaryPart(def->correctfact));
780:     }
781:     if (!def->lvl) {
782:       PetscViewerASCIIPrintf(viewer,"deflation space type: %s\n",PCDeflationSpaceTypes[def->spacetype]);
783:     }

785:     PetscViewerASCIIPrintf(viewer,"--- Additional PC:\n");
786:     PetscViewerASCIIPushTab(viewer);
787:     PCView(def->pc,viewer);
788:     PetscViewerASCIIPopTab(viewer);

790:     PetscViewerASCIIPrintf(viewer,"--- Coarse problem solver:\n");
791:     PetscViewerASCIIPushTab(viewer);
792:     KSPGetTotalIterations(def->WtAWinv,&its);
793:     PetscViewerASCIIPrintf(viewer,"total number of iterations: %D\n",its);
794:     KSPView(def->WtAWinv,viewer);
795:     PetscViewerASCIIPopTab(viewer);
796:   }
797:   return(0);
798: }

800: static PetscErrorCode PCSetFromOptions_Deflation(PetscOptionItems *PetscOptionsObject,PC pc)
801: {
802:   PC_Deflation      *def = (PC_Deflation*)pc->data;
803:   PetscErrorCode    ierr;

806:   PetscOptionsHead(PetscOptionsObject,"Deflation options");
807:   PetscOptionsBool("-pc_deflation_init_only","Use only initialization step - Initdef","PCDeflationSetInitOnly",def->init,&def->init,NULL);
808:   PetscOptionsInt("-pc_deflation_levels","Maximum of deflation levels","PCDeflationSetLevels",def->maxlvl,&def->maxlvl,NULL);
809:   PetscOptionsInt("-pc_deflation_reduction_factor","Reduction factor for coarse problem solution using PCTELESCOPE","PCDeflationSetReductionFactor",def->reductionfact,&def->reductionfact,NULL);
810:   PetscOptionsBool("-pc_deflation_correction","Add coarse problem correction Q to P","PCDeflationSetCorrectionFactor",def->correct,&def->correct,NULL);
811:   PetscOptionsScalar("-pc_deflation_correction_factor","Set multiple of Q to use as coarse problem correction","PCDeflationSetCorrectionFactor",def->correctfact,&def->correctfact,NULL);
812:   PetscOptionsEnum("-pc_deflation_compute_space","Compute deflation space","PCDeflationSetSpace",PCDeflationSpaceTypes,(PetscEnum)def->spacetype,(PetscEnum*)&def->spacetype,NULL);
813:   PetscOptionsInt("-pc_deflation_compute_space_size","Set size of the deflation space to compute","PCDeflationSetSpace",def->spacesize,&def->spacesize,NULL);
814:   PetscOptionsBool("-pc_deflation_space_extend","Extend deflation space instead of truncating (wavelets)","PCDeflation",def->extendsp,&def->extendsp,NULL);
815:   PetscOptionsTail();
816:   return(0);
817: }

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

822:    Options Database Keys:
823: +    -pc_deflation_init_only          <false> - if true computes only the special guess
824: .    -pc_deflation_max_lvl            <0>     - maximum number of levels for multilevel deflation
825: .    -pc_deflation_reduction_factor <\-1>     - reduction factor on bottom level coarse problem for PCTELESCOPE (default based on the size of the coarse problem)
826: .    -pc_deflation_correction         <false> - if true apply coarse problem correction
827: .    -pc_deflation_correction_factor  <1.0>   - sets coarse problem correction factor
828: .    -pc_deflation_compute_space      <haar>  - compute PCDeflationSpaceType deflation space
829: -    -pc_deflation_compute_space_size <1>     - size of the deflation space (corresponds to number of levels for wavelet-based deflation)

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

835:     The deflation computes initial guess x0 = x_{-1} - Q*r_{-1}, which is the solution on the deflation space.
836:     If PCDeflationSetInitOnly() or -pc_deflation_init_only is set to PETSC_TRUE (InitDef scheme), the Section 1.5 Writing Application Codes with PETSc of the
837:     preconditioner consists only of Section 1.5 Writing Application Codes with PETSc of the additional preconditioner M^{-1}. Otherwise, the preconditioner
838:     Section 1.5 Writing Application Codes with PETSc consists of P*M^{-1} + factor*Q. The first part of the preconditioner (PM^{-1}) shifts some eigenvalues
839:     to zero while the addition of the coarse problem correction (factor*Q) makes the preconditioner to shift some
840:     eigenvalues to the given factor. The InitDef scheme is recommended for deflation using high accuracy estimates
841:     of eigenvectors of A when it exhibits similar convergence to the full deflation but is cheaper.

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

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

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

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

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

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

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

878:    Developer Notes:
879:      Contributed by Jakub Kruzik (PERMON), Institute of Geonics of the Czech
880:      Academy of Sciences and VSB - TU Ostrava.

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

885:    References:
886: +    [1] - A. Nicolaides. "Deflation of conjugate gradients with applications to boundary value problems", SIAM J. Numer. Anal. 24.2, 1987.
887: .    [2] - Z. Dostal. "Conjugate gradient method with preconditioning by projector", Int J. Comput. Math. 23.3-4, 1988.
888: .    [3] - Y. A. Erlangga and R. Nabben. "Multilevel Projection-Based Nested Krylov Iteration for Boundary Value Problems", SIAM J. Sci. Comput. 30.3, 2008.
889: -    [4] - 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

891:    Level: intermediate

893: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
894:            PCDeflationSetInitOnly(), PCDeflationSetLevels(), PCDeflationSetReductionFactor(),
895:            PCDeflationSetCorrectionFactor(), PCDeflationSetSpaceToCompute(),
896:            PCDeflationSetSpace(), PCDeflationSpaceType, PCDeflationSetProjectionNullSpaceMat(),
897:            PCDeflationSetCoarseMat(), PCDeflationGetCoarseKSP(), PCDeflationGetPC()
898: M*/

900: PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc)
901: {
902:   PC_Deflation   *def;

906:   PetscNewLog(pc,&def);
907:   pc->data = (void*)def;

909:   def->init          = PETSC_FALSE;
910:   def->correct       = PETSC_FALSE;
911:   def->correctfact   = 1.0;
912:   def->reductionfact = -1;
913:   def->spacetype     = PC_DEFLATION_SPACE_HAAR;
914:   def->spacesize     = 1;
915:   def->extendsp      = PETSC_FALSE;
916:   def->lvl           = 0;
917:   def->maxlvl        = 0;
918:   def->W             = NULL;
919:   def->Wt            = NULL;

921:   pc->ops->apply          = PCApply_Deflation;
922:   pc->ops->presolve       = PCPreSolve_Deflation;
923:   pc->ops->setup          = PCSetUp_Deflation;
924:   pc->ops->reset          = PCReset_Deflation;
925:   pc->ops->destroy        = PCDestroy_Deflation;
926:   pc->ops->setfromoptions = PCSetFromOptions_Deflation;
927:   pc->ops->view           = PCView_Deflation;

929:   PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetInitOnly_C",PCDeflationSetInitOnly_Deflation);
930:   PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetLevels_C",PCDeflationSetLevels_Deflation);
931:   PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetReductionFactor_C",PCDeflationSetReductionFactor_Deflation);
932:   PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetCorrectionFactor_C",PCDeflationSetCorrectionFactor_Deflation);
933:   PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpaceToCompute_C",PCDeflationSetSpaceToCompute_Deflation);
934:   PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpace_C",PCDeflationSetSpace_Deflation);
935:   PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetProjectionNullSpaceMat_C",PCDeflationSetProjectionNullSpaceMat_Deflation);
936:   PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetCoarseMat_C",PCDeflationSetCoarseMat_Deflation);
937:   PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetCoarseKSP_C",PCDeflationGetCoarseKSP_Deflation);
938:   PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetPC_C",PCDeflationGetPC_Deflation);
939:   return(0);
940: }