Actual source code: deflation.c
petsc-3.12.5 2020-03-29
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 application of the
837: preconditioner consists only of application of the additional preconditioner M^{-1}. Otherwise, the preconditioner
838: application 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 valueproblems�, 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: }