Actual source code: deflation.c
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: NULL
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: }
56: static PetscErrorCode PCDeflationSetLevels_Deflation(PC pc,PetscInt current,PetscInt max)
57: {
58: PC_Deflation *def = (PC_Deflation*)pc->data;
61: if (current) def->lvl = current;
62: def->maxlvl = max;
63: return(0);
64: }
66: /*@
67: PCDeflationSetLevels - Set the maximum level of deflation nesting.
69: Logically Collective
71: Input Parameters:
72: + pc - the preconditioner context
73: - max - maximum deflation level
75: Options Database Keys:
76: . -pc_deflation_max_lvl <0> - maximum number of levels for multilevel deflation
78: Level: intermediate
80: .seealso: PCDeflationSetSpaceToCompute(), PCDeflationSetSpace(), PCDEFLATION
81: @*/
82: PetscErrorCode PCDeflationSetLevels(PC pc,PetscInt max)
83: {
89: PetscTryMethod(pc,"PCDeflationSetLevels_C",(PC,PetscInt,PetscInt),(pc,0,max));
90: return(0);
91: }
93: static PetscErrorCode PCDeflationSetReductionFactor_Deflation(PC pc,PetscInt red)
94: {
95: PC_Deflation *def = (PC_Deflation*)pc->data;
98: def->reductionfact = red;
99: return(0);
100: }
102: /*@
103: PCDeflationSetReductionFactor - Set reduction factor for the bottom PCTELESCOPE coarse problem solver.
105: Logically Collective
107: Input Parameters:
108: + pc - the preconditioner context
109: - red - reduction factor (or PETSC_DETERMINE)
111: Options Database Keys:
112: . -pc_deflation_reduction_factor <\-1> - reduction factor on bottom level coarse problem for PCTELESCOPE
114: Notes:
115: Default is computed based on the size of the coarse problem.
117: Level: intermediate
119: .seealso: PCTELESCOPE, PCDEFLATION
120: @*/
121: PetscErrorCode PCDeflationSetReductionFactor(PC pc,PetscInt red)
122: {
128: PetscTryMethod(pc,"PCDeflationSetReductionFactor_C",(PC,PetscInt),(pc,red));
129: return(0);
130: }
132: static PetscErrorCode PCDeflationSetCorrectionFactor_Deflation(PC pc,PetscScalar fact)
133: {
134: PC_Deflation *def = (PC_Deflation*)pc->data;
137: /* TODO PETSC_DETERMINE -> compute max eigenvalue with power method */
138: def->correct = PETSC_TRUE;
139: def->correctfact = fact;
140: if (def->correct == 0.0) {
141: def->correct = PETSC_FALSE;
142: }
143: return(0);
144: }
146: /*@
147: PCDeflationSetCorrectionFactor - Set coarse problem correction factor.
148: The Preconditioner becomes P*M^{-1} + fact*Q.
150: Logically Collective
152: Input Parameters:
153: + pc - the preconditioner context
154: - fact - correction factor
156: Options Database Keys:
157: + -pc_deflation_correction <false> - if true apply coarse problem correction
158: - -pc_deflation_correction_factor <1.0> - sets coarse problem correction factor
160: Notes:
161: Any non-zero fact enables the coarse problem correction.
163: Level: intermediate
165: .seealso: PCDEFLATION
166: @*/
167: PetscErrorCode PCDeflationSetCorrectionFactor(PC pc,PetscScalar fact)
168: {
174: PetscTryMethod(pc,"PCDeflationSetCorrectionFactor_C",(PC,PetscScalar),(pc,fact));
175: return(0);
176: }
178: static PetscErrorCode PCDeflationSetSpaceToCompute_Deflation(PC pc,PCDeflationSpaceType type,PetscInt size)
179: {
180: PC_Deflation *def = (PC_Deflation*)pc->data;
183: if (type) def->spacetype = type;
184: if (size > 0) def->spacesize = size;
185: return(0);
186: }
188: /*@
189: PCDeflationSetSpaceToCompute - Set deflation space type and size to compute.
191: Logically Collective
193: Input Parameters:
194: + pc - the preconditioner context
195: . type - deflation space type to compute (or PETSC_IGNORE)
196: - size - size of the space to compute (or PETSC_DEFAULT)
198: Options Database Keys:
199: + -pc_deflation_compute_space <haar> - compute PCDeflationSpaceType deflation space
200: - -pc_deflation_compute_space_size <1> - size of the deflation space
202: Notes:
203: For wavelet-based deflation, size represents number of levels.
205: The deflation space is computed in PCSetUp().
207: Level: intermediate
209: .seealso: PCDeflationSetLevels(), PCDEFLATION
210: @*/
211: PetscErrorCode PCDeflationSetSpaceToCompute(PC pc,PCDeflationSpaceType type,PetscInt size)
212: {
219: PetscTryMethod(pc,"PCDeflationSetSpaceToCompute_C",(PC,PCDeflationSpaceType,PetscInt),(pc,type,size));
220: return(0);
221: }
223: static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc,Mat W,PetscBool transpose)
224: {
225: PC_Deflation *def = (PC_Deflation*)pc->data;
229: /* possibly allows W' = Wt (which is valid but not tested) */
230: PetscObjectReference((PetscObject)W);
231: if (transpose) {
232: MatDestroy(&def->Wt);
233: def->Wt = W;
234: } else {
235: MatDestroy(&def->W);
236: def->W = W;
237: }
238: return(0);
239: }
241: /*@
242: PCDeflationSetSpace - Set the deflation space matrix (or its (Hermitian) transpose).
244: Logically Collective
246: Input Parameters:
247: + pc - the preconditioner context
248: . W - deflation matrix
249: - transpose - indicates that W is an explicit transpose of the deflation matrix
251: Notes:
252: Setting W as a multipliplicative MATCOMPOSITE enables use of the multilevel
253: deflation. If W = W0*W1*W2*...*Wn, W0 is taken as the first deflation space and
254: the coarse problem (W0'*A*W0)^{-1} is again preconditioned by deflation with
255: W1 as the deflation matrix. This repeats until the maximum level set by
256: PCDeflationSetLevels() is reached or there are no more matrices available.
257: If there are matrices left after reaching the maximum level,
258: they are merged into a deflation matrix ...*W{n-1}*Wn.
260: Level: intermediate
262: .seealso: PCDeflationSetLevels(), PCDEFLATION
263: @*/
264: PetscErrorCode PCDeflationSetSpace(PC pc,Mat W,PetscBool transpose)
265: {
272: PetscTryMethod(pc,"PCDeflationSetSpace_C",(PC,Mat,PetscBool),(pc,W,transpose));
273: return(0);
274: }
276: static PetscErrorCode PCDeflationSetProjectionNullSpaceMat_Deflation(PC pc,Mat mat)
277: {
278: PC_Deflation *def = (PC_Deflation*)pc->data;
279: PetscErrorCode ierr;
282: PetscObjectReference((PetscObject)mat);
283: MatDestroy(&def->WtA);
284: def->WtA = mat;
285: PetscLogObjectParent((PetscObject)pc,(PetscObject)def->WtA);
286: return(0);
287: }
289: /*@
290: PCDeflationSetProjectionNullSpaceMat - Set the projection null space matrix (W'*A).
292: Collective
294: Input Parameters:
295: + pc - preconditioner context
296: - mat - projection null space matrix
298: Level: developer
300: .seealso: PCDEFLATION
301: @*/
302: PetscErrorCode PCDeflationSetProjectionNullSpaceMat(PC pc,Mat mat)
303: {
309: PetscTryMethod(pc,"PCDeflationSetProjectionNullSpaceMat_C",(PC,Mat),(pc,mat));
310: return(0);
311: }
313: static PetscErrorCode PCDeflationSetCoarseMat_Deflation(PC pc,Mat mat)
314: {
315: PC_Deflation *def = (PC_Deflation*)pc->data;
316: PetscErrorCode ierr;
319: PetscObjectReference((PetscObject)mat);
320: MatDestroy(&def->WtAW);
321: def->WtAW = mat;
322: PetscLogObjectParent((PetscObject)pc,(PetscObject)def->WtAW);
323: return(0);
324: }
326: /*@
327: PCDeflationSetCoarseMat - Set the coarse problem Mat.
329: Collective
331: Input Parameters:
332: + pc - preconditioner context
333: - mat - coarse problem mat
335: Level: developer
337: .seealso: PCDEFLATION
338: @*/
339: PetscErrorCode PCDeflationSetCoarseMat(PC pc,Mat mat)
340: {
346: PetscTryMethod(pc,"PCDeflationSetCoarseMat_C",(PC,Mat),(pc,mat));
347: return(0);
348: }
350: static PetscErrorCode PCDeflationGetCoarseKSP_Deflation(PC pc,KSP *ksp)
351: {
352: PC_Deflation *def = (PC_Deflation*)pc->data;
355: *ksp = def->WtAWinv;
356: return(0);
357: }
359: /*@
360: PCDeflationGetCoarseKSP - Returns the coarse problem KSP.
362: Not Collective
364: Input Parameters:
365: . pc - preconditioner context
367: Output Parameters:
368: . ksp - coarse problem KSP context
370: Level: advanced
372: .seealso: PCDEFLATION
373: @*/
374: PetscErrorCode PCDeflationGetCoarseKSP(PC pc,KSP *ksp)
375: {
381: PetscTryMethod(pc,"PCDeflationGetCoarseKSP_C",(PC,KSP*),(pc,ksp));
382: return(0);
383: }
385: static PetscErrorCode PCDeflationGetPC_Deflation(PC pc,PC *apc)
386: {
387: PC_Deflation *def = (PC_Deflation*)pc->data;
390: *apc = def->pc;
391: return(0);
392: }
394: /*@
395: PCDeflationGetPC - Returns the additional preconditioner M^{-1}.
397: Not Collective
399: Input Parameters:
400: . pc - the preconditioner context
402: Output Parameters:
403: . apc - additional preconditioner
405: Level: advanced
407: .seealso: PCDEFLATION
408: @*/
409: PetscErrorCode PCDeflationGetPC(PC pc,PC *apc)
410: {
416: PetscTryMethod(pc,"PCDeflationGetPC_C",(PC,PC*),(pc,apc));
417: return(0);
418: }
420: /*
421: x <- x + W*(W'*A*W)^{-1}*W'*r = x + Q*r
422: */
423: static PetscErrorCode PCPreSolve_Deflation(PC pc,KSP ksp,Vec b, Vec x)
424: {
425: PC_Deflation *def = (PC_Deflation*)pc->data;
426: Mat A;
427: Vec r,w1,w2;
428: PetscBool nonzero;
429: PetscErrorCode ierr;
432: w1 = def->workcoarse[0];
433: w2 = def->workcoarse[1];
434: r = def->work;
435: PCGetOperators(pc,NULL,&A);
437: KSPGetInitialGuessNonzero(ksp,&nonzero);
438: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
439: if (nonzero) {
440: MatMult(A,x,r); /* r <- b - Ax */
441: VecAYPX(r,-1.0,b);
442: } else {
443: VecCopy(b,r); /* r <- b (x is 0) */
444: }
446: if (def->Wt) {
447: MatMult(def->Wt,r,w1); /* w1 <- W'*r */
448: } else {
449: MatMultHermitianTranspose(def->W,r,w1); /* w1 <- W'*r */
450: }
451: KSPSolve(def->WtAWinv,w1,w2); /* w2 <- (W'*A*W)^{-1}*w1 */
452: MatMult(def->W,w2,r); /* r <- W*w2 */
453: VecAYPX(x,1.0,r);
454: return(0);
455: }
457: /*
458: if (def->correct) {
459: z <- M^{-1}r - W*(W'*A*W)^{-1}*(W'*A*M^{-1}r - l*W'*r) = (P*M^{-1} + l*Q)*r
460: } else {
461: z <- M^{-1}*r - W*(W'*A*W)^{-1}*W'*A*M{-1}*r = P*M^{-1}*r
462: }
463: */
464: static PetscErrorCode PCApply_Deflation(PC pc,Vec r,Vec z)
465: {
466: PC_Deflation *def = (PC_Deflation*)pc->data;
467: Mat A;
468: Vec u,w1,w2;
469: PetscErrorCode ierr;
472: w1 = def->workcoarse[0];
473: w2 = def->workcoarse[1];
474: u = def->work;
475: PCGetOperators(pc,NULL,&A);
477: PCApply(def->pc,r,z); /* z <- M^{-1}*r */
478: if (!def->init) {
479: MatMult(def->WtA,z,w1); /* w1 <- W'*A*z */
480: if (def->correct) {
481: if (def->Wt) {
482: MatMult(def->Wt,r,w2); /* w2 <- W'*r */
483: } else {
484: MatMultHermitianTranspose(def->W,r,w2); /* w2 <- W'*r */
485: }
486: VecAXPY(w1,-1.0*def->correctfact,w2); /* w1 <- w1 - l*w2 */
487: }
488: KSPSolve(def->WtAWinv,w1,w2); /* w2 <- (W'*A*W)^{-1}*w1 */
489: MatMult(def->W,w2,u); /* u <- W*w2 */
490: VecAXPY(z,-1.0,u); /* z <- z - u */
491: }
492: return(0);
493: }
495: static PetscErrorCode PCSetUp_Deflation(PC pc)
496: {
497: PC_Deflation *def = (PC_Deflation*)pc->data;
498: KSP innerksp;
499: PC pcinner;
500: Mat Amat,nextDef=NULL,*mats;
501: PetscInt i,m,red,size;
502: PetscMPIInt commsize;
503: PetscBool match,flgspd,transp=PETSC_FALSE;
504: MatCompositeType ctype;
505: MPI_Comm comm;
506: char prefix[128]="";
507: PetscErrorCode ierr;
510: if (pc->setupcalled) return(0);
511: PetscObjectGetComm((PetscObject)pc,&comm);
512: PCGetOperators(pc,NULL,&Amat);
513: if (!def->lvl && !def->prefix) {
514: PCGetOptionsPrefix(pc,&def->prefix);
515: }
516: if (def->lvl) {
517: PetscSNPrintf(prefix,sizeof(prefix),"%d_",(int)def->lvl);
518: }
520: /* compute a deflation space */
521: if (def->W || def->Wt) {
522: def->spacetype = PC_DEFLATION_SPACE_USER;
523: } else {
524: PCDeflationComputeSpace(pc);
525: }
527: /* nested deflation */
528: if (def->W) {
529: PetscObjectTypeCompare((PetscObject)def->W,MATCOMPOSITE,&match);
530: if (match) {
531: MatCompositeGetType(def->W,&ctype);
532: MatCompositeGetNumberMat(def->W,&size);
533: }
534: } else {
535: MatCreateHermitianTranspose(def->Wt,&def->W);
536: PetscObjectTypeCompare((PetscObject)def->Wt,MATCOMPOSITE,&match);
537: if (match) {
538: MatCompositeGetType(def->Wt,&ctype);
539: MatCompositeGetNumberMat(def->Wt,&size);
540: }
541: transp = PETSC_TRUE;
542: }
543: if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) {
544: if (!transp) {
545: if (def->lvl < def->maxlvl) {
546: PetscMalloc1(size,&mats);
547: for (i=0; i<size; i++) {
548: MatCompositeGetMat(def->W,i,&mats[i]);
549: }
550: size -= 1;
551: MatDestroy(&def->W);
552: def->W = mats[size];
553: PetscObjectReference((PetscObject)mats[size]);
554: if (size > 1) {
555: MatCreateComposite(comm,size,mats,&nextDef);
556: MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);
557: } else {
558: nextDef = mats[0];
559: PetscObjectReference((PetscObject)mats[0]);
560: }
561: PetscFree(mats);
562: } else {
563: /* TODO test merge side performance */
564: /* MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT); */
565: MatCompositeMerge(def->W);
566: }
567: } else {
568: if (def->lvl < def->maxlvl) {
569: PetscMalloc1(size,&mats);
570: for (i=0; i<size; i++) {
571: MatCompositeGetMat(def->Wt,i,&mats[i]);
572: }
573: size -= 1;
574: MatDestroy(&def->Wt);
575: def->Wt = mats[0];
576: PetscObjectReference((PetscObject)mats[0]);
577: if (size > 1) {
578: MatCreateComposite(comm,size,&mats[1],&nextDef);
579: MatCompositeSetType(nextDef,MAT_COMPOSITE_MULTIPLICATIVE);
580: } else {
581: nextDef = mats[1];
582: PetscObjectReference((PetscObject)mats[1]);
583: }
584: PetscFree(mats);
585: } else {
586: /* MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT); */
587: MatCompositeMerge(def->Wt);
588: }
589: }
590: }
592: if (transp) {
593: MatDestroy(&def->W);
594: MatHermitianTranspose(def->Wt,MAT_INITIAL_MATRIX,&def->W);
595: }
597: /* assemble WtA */
598: if (!def->WtA) {
599: if (def->Wt) {
600: MatMatMult(def->Wt,Amat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtA);
601: } else {
602: #if defined(PETSC_USE_COMPLEX)
603: MatHermitianTranspose(def->W,MAT_INITIAL_MATRIX,&def->Wt);
604: MatMatMult(def->Wt,Amat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtA);
605: #else
606: MatTransposeMatMult(def->W,Amat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtA);
607: #endif
608: }
609: }
610: /* setup coarse problem */
611: if (!def->WtAWinv) {
612: MatGetSize(def->W,NULL,&m);
613: if (!def->WtAW) {
614: MatMatMult(def->WtA,def->W,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&def->WtAW);
615: /* TODO create MatInheritOption(Mat,MatOption) */
616: MatGetOption(Amat,MAT_SPD,&flgspd);
617: MatSetOption(def->WtAW,MAT_SPD,flgspd);
618: if (PetscDefined(USE_DEBUG)) {
619: /* Check columns of W are not in kernel of A */
620: PetscReal *norms;
621: PetscMalloc1(m,&norms);
622: MatGetColumnNorms(def->WtAW,NORM_INFINITY,norms);
623: for (i=0; i<m; i++) {
624: if (norms[i] < 100*PETSC_MACHINE_EPSILON) {
625: SETERRQ1(comm,PETSC_ERR_SUP,"Column %D of W is in kernel of A.",i);
626: }
627: }
628: PetscFree(norms);
629: }
630: } else {
631: MatGetOption(def->WtAW,MAT_SPD,&flgspd);
632: }
633: /* TODO use MATINV ? */
634: KSPCreate(comm,&def->WtAWinv);
635: KSPSetOperators(def->WtAWinv,def->WtAW,def->WtAW);
636: KSPGetPC(def->WtAWinv,&pcinner);
637: /* Setup KSP and PC */
638: if (nextDef) { /* next level for multilevel deflation */
639: innerksp = def->WtAWinv;
640: /* set default KSPtype */
641: if (!def->ksptype) {
642: def->ksptype = KSPFGMRES;
643: if (flgspd) { /* SPD system */
644: def->ksptype = KSPFCG;
645: }
646: }
647: KSPSetType(innerksp,def->ksptype); /* TODO iherit from KSP + tolerances */
648: PCSetType(pcinner,PCDEFLATION); /* TODO create coarse preconditinoner M_c = WtMW ? */
649: PCDeflationSetSpace(pcinner,nextDef,transp);
650: PCDeflationSetLevels_Deflation(pcinner,def->lvl+1,def->maxlvl);
651: /* inherit options */
652: if (def->prefix) ((PC_Deflation*)(pcinner->data))->prefix = def->prefix;
653: ((PC_Deflation*)(pcinner->data))->init = def->init;
654: ((PC_Deflation*)(pcinner->data))->ksptype = def->ksptype;
655: ((PC_Deflation*)(pcinner->data))->correct = def->correct;
656: ((PC_Deflation*)(pcinner->data))->correctfact = def->correctfact;
657: ((PC_Deflation*)(pcinner->data))->reductionfact = def->reductionfact;
658: MatDestroy(&nextDef);
659: } else { /* the last level */
660: KSPSetType(def->WtAWinv,KSPPREONLY);
661: PCSetType(pcinner,PCTELESCOPE);
662: /* do not overwrite PCTELESCOPE */
663: if (def->prefix) {
664: KSPSetOptionsPrefix(def->WtAWinv,def->prefix);
665: }
666: KSPAppendOptionsPrefix(def->WtAWinv,"deflation_tel_");
667: PCSetFromOptions(pcinner);
668: PetscObjectTypeCompare((PetscObject)pcinner,PCTELESCOPE,&match);
669: if (!match) SETERRQ(comm,PETSC_ERR_SUP,"User can not owerwrite PCTELESCOPE on bottom level, use reduction factor = 1 instead.");
670: /* Reduction factor choice */
671: red = def->reductionfact;
672: if (red < 0) {
673: MPI_Comm_size(comm,&commsize);
674: red = ceil((float)commsize/ceil((float)m/commsize));
675: PetscObjectTypeCompareAny((PetscObject)(def->WtAW),&match,MATSEQDENSE,MATMPIDENSE,MATDENSE,"");
676: if (match) red = commsize;
677: PetscInfo1(pc,"Auto choosing reduction factor %D\n",red);
678: }
679: PCTelescopeSetReductionFactor(pcinner,red);
680: PCSetUp(pcinner);
681: PCTelescopeGetKSP(pcinner,&innerksp);
682: if (innerksp) {
683: KSPGetPC(innerksp,&pcinner);
684: PCSetType(pcinner,PCLU);
685: #if defined(PETSC_HAVE_SUPERLU)
686: MatGetFactorAvailable(def->WtAW,MATSOLVERSUPERLU,MAT_FACTOR_LU,&match);
687: if (match) {
688: PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU);
689: }
690: #endif
691: #if defined(PETSC_HAVE_SUPERLU_DIST)
692: MatGetFactorAvailable(def->WtAW,MATSOLVERSUPERLU_DIST,MAT_FACTOR_LU,&match);
693: if (match) {
694: PCFactorSetMatSolverType(pcinner,MATSOLVERSUPERLU_DIST);
695: }
696: #endif
697: }
698: }
700: if (innerksp) {
701: if (def->prefix) {
702: KSPSetOptionsPrefix(innerksp,def->prefix);
703: KSPAppendOptionsPrefix(innerksp,"deflation_");
704: } else {
705: KSPSetOptionsPrefix(innerksp,"deflation_");
706: }
707: KSPAppendOptionsPrefix(innerksp,prefix);
708: KSPSetFromOptions(innerksp);
709: KSPSetUp(innerksp);
710: }
711: }
712: KSPSetFromOptions(def->WtAWinv);
713: KSPSetUp(def->WtAWinv);
715: /* create preconditioner */
716: if (!def->pc) {
717: PCCreate(comm,&def->pc);
718: PCSetOperators(def->pc,Amat,Amat);
719: PCSetType(def->pc,PCNONE);
720: if (def->prefix) {
721: PCSetOptionsPrefix(def->pc,def->prefix);
722: }
723: PCAppendOptionsPrefix(def->pc,"deflation_");
724: PCAppendOptionsPrefix(def->pc,prefix);
725: PCAppendOptionsPrefix(def->pc,"pc_");
726: PCSetFromOptions(def->pc);
727: PCSetUp(def->pc);
728: }
730: /* create work vecs */
731: MatCreateVecs(Amat,NULL,&def->work);
732: KSPCreateVecs(def->WtAWinv,2,&def->workcoarse,0,NULL);
733: return(0);
734: }
736: static PetscErrorCode PCReset_Deflation(PC pc)
737: {
738: PC_Deflation *def = (PC_Deflation*)pc->data;
739: PetscErrorCode ierr;
742: VecDestroy(&def->work);
743: VecDestroyVecs(2,&def->workcoarse);
744: MatDestroy(&def->W);
745: MatDestroy(&def->Wt);
746: MatDestroy(&def->WtA);
747: MatDestroy(&def->WtAW);
748: KSPDestroy(&def->WtAWinv);
749: PCDestroy(&def->pc);
750: return(0);
751: }
753: static PetscErrorCode PCDestroy_Deflation(PC pc)
754: {
758: PCReset_Deflation(pc);
759: PetscFree(pc->data);
760: return(0);
761: }
763: static PetscErrorCode PCView_Deflation(PC pc,PetscViewer viewer)
764: {
765: PC_Deflation *def = (PC_Deflation*)pc->data;
766: PetscInt its;
767: PetscBool iascii;
768: PetscErrorCode ierr;
771: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
772: if (iascii) {
773: if (def->correct) {
774: PetscViewerASCIIPrintf(viewer,"using CP correction, factor = %g+%gi\n",
775: (double)PetscRealPart(def->correctfact),
776: (double)PetscImaginaryPart(def->correctfact));
777: }
778: if (!def->lvl) {
779: PetscViewerASCIIPrintf(viewer,"deflation space type: %s\n",PCDeflationSpaceTypes[def->spacetype]);
780: }
782: PetscViewerASCIIPrintf(viewer,"--- Additional PC:\n");
783: PetscViewerASCIIPushTab(viewer);
784: PCView(def->pc,viewer);
785: PetscViewerASCIIPopTab(viewer);
787: PetscViewerASCIIPrintf(viewer,"--- Coarse problem solver:\n");
788: PetscViewerASCIIPushTab(viewer);
789: KSPGetTotalIterations(def->WtAWinv,&its);
790: PetscViewerASCIIPrintf(viewer,"total number of iterations: %D\n",its);
791: KSPView(def->WtAWinv,viewer);
792: PetscViewerASCIIPopTab(viewer);
793: }
794: return(0);
795: }
797: static PetscErrorCode PCSetFromOptions_Deflation(PetscOptionItems *PetscOptionsObject,PC pc)
798: {
799: PC_Deflation *def = (PC_Deflation*)pc->data;
800: PetscErrorCode ierr;
803: PetscOptionsHead(PetscOptionsObject,"Deflation options");
804: PetscOptionsBool("-pc_deflation_init_only","Use only initialization step - Initdef","PCDeflationSetInitOnly",def->init,&def->init,NULL);
805: PetscOptionsInt("-pc_deflation_levels","Maximum of deflation levels","PCDeflationSetLevels",def->maxlvl,&def->maxlvl,NULL);
806: PetscOptionsInt("-pc_deflation_reduction_factor","Reduction factor for coarse problem solution using PCTELESCOPE","PCDeflationSetReductionFactor",def->reductionfact,&def->reductionfact,NULL);
807: PetscOptionsBool("-pc_deflation_correction","Add coarse problem correction Q to P","PCDeflationSetCorrectionFactor",def->correct,&def->correct,NULL);
808: PetscOptionsScalar("-pc_deflation_correction_factor","Set multiple of Q to use as coarse problem correction","PCDeflationSetCorrectionFactor",def->correctfact,&def->correctfact,NULL);
809: PetscOptionsEnum("-pc_deflation_compute_space","Compute deflation space","PCDeflationSetSpace",PCDeflationSpaceTypes,(PetscEnum)def->spacetype,(PetscEnum*)&def->spacetype,NULL);
810: PetscOptionsInt("-pc_deflation_compute_space_size","Set size of the deflation space to compute","PCDeflationSetSpace",def->spacesize,&def->spacesize,NULL);
811: PetscOptionsBool("-pc_deflation_space_extend","Extend deflation space instead of truncating (wavelets)","PCDeflation",def->extendsp,&def->extendsp,NULL);
812: PetscOptionsTail();
813: return(0);
814: }
816: /*MC
817: PCDEFLATION - Deflation preconditioner shifts (deflates) part of the spectrum to zero or to a predefined value.
819: Options Database Keys:
820: + -pc_deflation_init_only <false> - if true computes only the special guess
821: . -pc_deflation_max_lvl <0> - maximum number of levels for multilevel deflation
822: . -pc_deflation_reduction_factor <\-1> - reduction factor on bottom level coarse problem for PCTELESCOPE (default based on the size of the coarse problem)
823: . -pc_deflation_correction <false> - if true apply coarse problem correction
824: . -pc_deflation_correction_factor <1.0> - sets coarse problem correction factor
825: . -pc_deflation_compute_space <haar> - compute PCDeflationSpaceType deflation space
826: - -pc_deflation_compute_space_size <1> - size of the deflation space (corresponds to number of levels for wavelet-based deflation)
828: Notes:
829: Given a (complex - transpose is always Hermitian) full rank deflation matrix W, the deflation (introduced in [1,2])
830: preconditioner uses projections Q = W*(W'*A*W)^{-1}*W' and P = I - Q*A, where A is pmat.
832: The deflation computes initial guess x0 = x_{-1} - Q*r_{-1}, which is the solution on the deflation space.
833: If PCDeflationSetInitOnly() or -pc_deflation_init_only is set to PETSC_TRUE (InitDef scheme), the application of the
834: preconditioner consists only of application of the additional preconditioner M^{-1}. Otherwise, the preconditioner
835: application consists of P*M^{-1} + factor*Q. The first part of the preconditioner (PM^{-1}) shifts some eigenvalues
836: to zero while the addition of the coarse problem correction (factor*Q) makes the preconditioner to shift some
837: eigenvalues to the given factor. The InitDef scheme is recommended for deflation using high accuracy estimates
838: of eigenvectors of A when it exhibits similar convergence to the full deflation but is cheaper.
840: The deflation matrix is by default automatically computed. The type of deflation matrix and its size to compute can
841: be controlled by PCDeflationSetSpaceToCompute() or -pc_deflation_compute_space and -pc_deflation_compute_space_size.
842: User can set an arbitrary deflation space matrix with PCDeflationSetSpace(). If the deflation matrix
843: is a multiplicative MATCOMPOSITE, a multilevel deflation [3] is used. The first matrix in the composite is used as the
844: deflation matrix, and the coarse problem (W'*A*W)^{-1} is solved by KSPFCG (if A is MAT_SPD) or KSPFGMRES preconditioned
845: by deflation with deflation matrix being the next matrix in the MATCOMPOSITE. This scheme repeats until the maximum
846: level is reached or there are no more matrices. If the maximum level is reached, the remaining matrices are merged
847: (multiplied) to create the last deflation matrix. The maximum level defaults to 0 and can be set by
848: PCDeflationSetLevels() or by -pc_deflation_levels.
850: The coarse problem KSP can be controlled from the command line with prefix -deflation_ for the first level and -deflation_[lvl-1]
851: from the second level onward. You can also use
852: PCDeflationGetCoarseKSP() to control it from code. The bottom level KSP defaults to
853: KSPPREONLY with PCLU direct solver (MATSOLVERSUPERLU/MATSOLVERSUPERLU_DIST if available) wrapped into PCTELESCOPE.
854: For convenience, the reduction factor can be set by PCDeflationSetReductionFactor()
855: or -pc_deflation_recduction_factor. The default is chosen heuristically based on the coarse problem size.
857: The additional preconditioner can be controlled from command line with prefix -deflation_[lvl]_pc (same rules used for
858: coarse problem KSP apply for [lvl]_ part of prefix), e.g., -deflation_1_pc_pc_type bjacobi. You can also use
859: PCDeflationGetPC() to control the additional preconditioner from code. It defaults to PCNONE.
861: The coarse problem correction term (factor*Q) can be turned on by -pc_deflation_correction and the factor value can
862: be set by pc_deflation_correction_factor or by PCDeflationSetCorrectionFactor(). The coarse problem can
863: significantly improve convergence when the deflation coarse problem is not solved with high enough accuracy. We
864: recommend setting factor to some eigenvalue, e.g., the largest eigenvalue so that the preconditioner does not create
865: an isolated eigenvalue.
867: The options are automatically inherited from the previous deflation level.
869: The preconditioner supports KSPMonitorDynamicTolerance(). This is useful for the multilevel scheme for which we also
870: recommend limiting the number of iterations for the coarse problems.
872: See section 3 of [4] for additional references and decription of the algorithm when used for conjugate gradients.
873: Section 4 describes some possible choices for the deflation space.
875: Developer Notes:
876: Contributed by Jakub Kruzik (PERMON), Institute of Geonics of the Czech
877: Academy of Sciences and VSB - TU Ostrava.
879: Developed from PERMON code used in [4] while on a research stay with
880: Prof. Reinhard Nabben at the Institute of Mathematics, TU Berlin.
882: References:
883: + [1] - A. Nicolaides. "Deflation of conjugate gradients with applications to boundary value problems", SIAM J. Numer. Anal. 24.2, 1987.
884: . [2] - Z. Dostal. "Conjugate gradient method with preconditioning by projector", Int J. Comput. Math. 23.3-4, 1988.
885: . [3] - Y. A. Erlangga and R. Nabben. "Multilevel Projection-Based Nested Krylov Iteration for Boundary Value Problems", SIAM J. Sci. Comput. 30.3, 2008.
886: - [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
888: Level: intermediate
890: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
891: PCDeflationSetInitOnly(), PCDeflationSetLevels(), PCDeflationSetReductionFactor(),
892: PCDeflationSetCorrectionFactor(), PCDeflationSetSpaceToCompute(),
893: PCDeflationSetSpace(), PCDeflationSpaceType, PCDeflationSetProjectionNullSpaceMat(),
894: PCDeflationSetCoarseMat(), PCDeflationGetCoarseKSP(), PCDeflationGetPC()
895: M*/
897: PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc)
898: {
899: PC_Deflation *def;
903: PetscNewLog(pc,&def);
904: pc->data = (void*)def;
906: def->init = PETSC_FALSE;
907: def->correct = PETSC_FALSE;
908: def->correctfact = 1.0;
909: def->reductionfact = -1;
910: def->spacetype = PC_DEFLATION_SPACE_HAAR;
911: def->spacesize = 1;
912: def->extendsp = PETSC_FALSE;
913: def->lvl = 0;
914: def->maxlvl = 0;
915: def->W = NULL;
916: def->Wt = NULL;
918: pc->ops->apply = PCApply_Deflation;
919: pc->ops->presolve = PCPreSolve_Deflation;
920: pc->ops->setup = PCSetUp_Deflation;
921: pc->ops->reset = PCReset_Deflation;
922: pc->ops->destroy = PCDestroy_Deflation;
923: pc->ops->setfromoptions = PCSetFromOptions_Deflation;
924: pc->ops->view = PCView_Deflation;
926: PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetInitOnly_C",PCDeflationSetInitOnly_Deflation);
927: PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetLevels_C",PCDeflationSetLevels_Deflation);
928: PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetReductionFactor_C",PCDeflationSetReductionFactor_Deflation);
929: PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetCorrectionFactor_C",PCDeflationSetCorrectionFactor_Deflation);
930: PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpaceToCompute_C",PCDeflationSetSpaceToCompute_Deflation);
931: PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetSpace_C",PCDeflationSetSpace_Deflation);
932: PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetProjectionNullSpaceMat_C",PCDeflationSetProjectionNullSpaceMat_Deflation);
933: PetscObjectComposeFunction((PetscObject)pc,"PCDeflationSetCoarseMat_C",PCDeflationSetCoarseMat_Deflation);
934: PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetCoarseKSP_C",PCDeflationGetCoarseKSP_Deflation);
935: PetscObjectComposeFunction((PetscObject)pc,"PCDeflationGetPC_C",PCDeflationGetPC_Deflation);
936: return(0);
937: }