Actual source code: deflation.c
1: #include <../src/ksp/pc/impls/deflation/deflation.h>
3: const char *const PCDeflationSpaceTypes[] = {"haar", "db2", "db4", "db8", "db16", "biorth22", "meyer", "aggregation", "user", "PCDeflationSpaceType", "PC_DEFLATION_SPACE_", NULL};
5: static PetscErrorCode PCDeflationSetInitOnly_Deflation(PC pc, PetscBool flg)
6: {
7: PC_Deflation *def = (PC_Deflation *)pc->data;
9: def->init = flg;
10: return 0;
11: }
13: /*@
14: PCDeflationSetInitOnly - Do only initialization step.
15: Sets initial guess to the solution on the deflation space but does not apply
16: the deflation preconditioner. The additional preconditioner is still applied.
18: Logically Collective
20: Input Parameters:
21: + pc - the preconditioner context
22: - flg - default `PETSC_FALSE`
24: Options Database Key:
25: . -pc_deflation_init_only <false> - if true computes only the special guess
27: Level: intermediate
29: .seealso: `PCDEFLATION`
30: @*/
31: PetscErrorCode PCDeflationSetInitOnly(PC pc, PetscBool flg)
32: {
35: PetscTryMethod(pc, "PCDeflationSetInitOnly_C", (PC, PetscBool), (pc, flg));
36: return 0;
37: }
39: static PetscErrorCode PCDeflationSetLevels_Deflation(PC pc, PetscInt current, PetscInt max)
40: {
41: PC_Deflation *def = (PC_Deflation *)pc->data;
43: if (current) def->lvl = current;
44: def->maxlvl = max;
45: return 0;
46: }
48: /*@
49: PCDeflationSetLevels - Set the maximum level of deflation nesting.
51: Logically Collective
53: Input Parameters:
54: + pc - the preconditioner context
55: - max - maximum deflation level
57: Options Database Key:
58: . -pc_deflation_max_lvl <0> - maximum number of levels for multilevel deflation
60: Level: intermediate
62: .seealso: `PCDeflationSetSpaceToCompute()`, `PCDeflationSetSpace()`, `PCDEFLATION`
63: @*/
64: PetscErrorCode PCDeflationSetLevels(PC pc, PetscInt max)
65: {
68: PetscTryMethod(pc, "PCDeflationSetLevels_C", (PC, PetscInt, PetscInt), (pc, 0, max));
69: return 0;
70: }
72: static PetscErrorCode PCDeflationSetReductionFactor_Deflation(PC pc, PetscInt red)
73: {
74: PC_Deflation *def = (PC_Deflation *)pc->data;
76: def->reductionfact = red;
77: return 0;
78: }
80: /*@
81: PCDeflationSetReductionFactor - Set reduction factor for the `PCDEFLATION`
83: Logically Collective
85: Input Parameters:
86: + pc - the preconditioner context
87: - red - reduction factor (or `PETSC_DETERMINE`)
89: Options Database Key:
90: . -pc_deflation_reduction_factor <\-1> - reduction factor on bottom level coarse problem for `PCDEFLATION`
92: Note:
93: Default is computed based on the size of the coarse problem.
95: Level: intermediate
97: .seealso: `PCTELESCOPE`, `PCDEFLATION`, `PCDeflationSetLevels()`
98: @*/
99: PetscErrorCode PCDeflationSetReductionFactor(PC pc, PetscInt red)
100: {
103: PetscTryMethod(pc, "PCDeflationSetReductionFactor_C", (PC, PetscInt), (pc, red));
104: return 0;
105: }
107: static PetscErrorCode PCDeflationSetCorrectionFactor_Deflation(PC pc, PetscScalar fact)
108: {
109: PC_Deflation *def = (PC_Deflation *)pc->data;
111: /* TODO PETSC_DETERMINE -> compute max eigenvalue with power method */
112: def->correct = PETSC_TRUE;
113: def->correctfact = fact;
114: if (def->correct == 0.0) def->correct = PETSC_FALSE;
115: return 0;
116: }
118: /*@
119: PCDeflationSetCorrectionFactor - Set coarse problem correction factor.
120: The preconditioner becomes P*M^{-1} + fact*Q.
122: Logically Collective
124: Input Parameters:
125: + pc - the preconditioner context
126: - fact - correction factor
128: Options Database Keys:
129: + -pc_deflation_correction <false> - if true apply coarse problem correction
130: - -pc_deflation_correction_factor <1.0> - sets coarse problem correction factor
132: Note:
133: Any non-zero fact enables the coarse problem correction.
135: Level: intermediate
137: .seealso: `PCDEFLATION`, `PCDeflationSetLevels()`, `PCDeflationSetReductionFactor()`
138: @*/
139: PetscErrorCode PCDeflationSetCorrectionFactor(PC pc, PetscScalar fact)
140: {
143: PetscTryMethod(pc, "PCDeflationSetCorrectionFactor_C", (PC, PetscScalar), (pc, fact));
144: return 0;
145: }
147: static PetscErrorCode PCDeflationSetSpaceToCompute_Deflation(PC pc, PCDeflationSpaceType type, PetscInt size)
148: {
149: PC_Deflation *def = (PC_Deflation *)pc->data;
151: if (type) def->spacetype = type;
152: if (size > 0) def->spacesize = size;
153: return 0;
154: }
156: /*@
157: PCDeflationSetSpaceToCompute - Set deflation space type and size to compute.
159: Logically Collective
161: Input Parameters:
162: + pc - the preconditioner context
163: . type - deflation space type to compute (or `PETSC_IGNORE`)
164: - size - size of the space to compute (or `PETSC_DEFAULT`)
166: Options Database Keys:
167: + -pc_deflation_compute_space <haar> - compute `PCDeflationSpaceType` deflation space
168: - -pc_deflation_compute_space_size <1> - size of the deflation space
170: Notes:
171: For wavelet-based deflation, size represents number of levels.
173: The deflation space is computed in `PCSetUp()`.
175: Level: intermediate
177: .seealso: `PCDeflationSetLevels()`, `PCDEFLATION`
178: @*/
179: PetscErrorCode PCDeflationSetSpaceToCompute(PC pc, PCDeflationSpaceType type, PetscInt size)
180: {
184: PetscTryMethod(pc, "PCDeflationSetSpaceToCompute_C", (PC, PCDeflationSpaceType, PetscInt), (pc, type, size));
185: return 0;
186: }
188: static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc, Mat W, PetscBool transpose)
189: {
190: PC_Deflation *def = (PC_Deflation *)pc->data;
192: /* possibly allows W' = Wt (which is valid but not tested) */
193: PetscObjectReference((PetscObject)W);
194: if (transpose) {
195: MatDestroy(&def->Wt);
196: def->Wt = W;
197: } else {
198: MatDestroy(&def->W);
199: def->W = W;
200: }
201: return 0;
202: }
204: /*@
205: PCDeflationSetSpace - Set the deflation space matrix (or its (Hermitian) transpose).
207: Logically Collective
209: Input Parameters:
210: + pc - the preconditioner context
211: . W - deflation matrix
212: - transpose - indicates that W is an explicit transpose of the deflation matrix
214: Notes:
215: Setting W as a multipliplicative `MATCOMPOSITE` enables use of the multilevel
216: deflation. If W = W0*W1*W2*...*Wn, W0 is taken as the first deflation space and
217: the coarse problem (W0'*A*W0)^{-1} is again preconditioned by deflation with
218: W1 as the deflation matrix. This repeats until the maximum level set by
219: PCDeflationSetLevels() is reached or there are no more matrices available.
220: If there are matrices left after reaching the maximum level,
221: they are merged into a deflation matrix ...*W{n-1}*Wn.
223: Level: intermediate
225: .seealso: `PCDeflationSetLevels()`, `PCDEFLATION`, `PCDeflationSetProjectionNullSpaceMat()`
226: @*/
227: PetscErrorCode PCDeflationSetSpace(PC pc, Mat W, PetscBool transpose)
228: {
232: PetscTryMethod(pc, "PCDeflationSetSpace_C", (PC, Mat, PetscBool), (pc, W, transpose));
233: return 0;
234: }
236: static PetscErrorCode PCDeflationSetProjectionNullSpaceMat_Deflation(PC pc, Mat mat)
237: {
238: PC_Deflation *def = (PC_Deflation *)pc->data;
240: PetscObjectReference((PetscObject)mat);
241: MatDestroy(&def->WtA);
242: def->WtA = mat;
243: return 0;
244: }
246: /*@
247: PCDeflationSetProjectionNullSpaceMat - Set the projection null space matrix (W'*A).
249: Collective
251: Input Parameters:
252: + pc - preconditioner context
253: - mat - projection null space matrix
255: Level: developer
257: .seealso: `PCDEFLATION`, `PCDeflationSetSpace()`
258: @*/
259: PetscErrorCode PCDeflationSetProjectionNullSpaceMat(PC pc, Mat mat)
260: {
263: PetscTryMethod(pc, "PCDeflationSetProjectionNullSpaceMat_C", (PC, Mat), (pc, mat));
264: return 0;
265: }
267: static PetscErrorCode PCDeflationSetCoarseMat_Deflation(PC pc, Mat mat)
268: {
269: PC_Deflation *def = (PC_Deflation *)pc->data;
271: PetscObjectReference((PetscObject)mat);
272: MatDestroy(&def->WtAW);
273: def->WtAW = mat;
274: return 0;
275: }
277: /*@
278: PCDeflationSetCoarseMat - Set the coarse problem `Mat`.
280: Collective
282: Input Parameters:
283: + pc - preconditioner context
284: - mat - coarse problem mat
286: Level: developer
288: .seealso: `PCDEFLATION`, `PCDeflationGetCoarseKSP()`
289: @*/
290: PetscErrorCode PCDeflationSetCoarseMat(PC pc, Mat mat)
291: {
294: PetscTryMethod(pc, "PCDeflationSetCoarseMat_C", (PC, Mat), (pc, mat));
295: return 0;
296: }
298: static PetscErrorCode PCDeflationGetCoarseKSP_Deflation(PC pc, KSP *ksp)
299: {
300: PC_Deflation *def = (PC_Deflation *)pc->data;
302: *ksp = def->WtAWinv;
303: return 0;
304: }
306: /*@
307: PCDeflationGetCoarseKSP - Returns the coarse problem `KSP`.
309: Not Collective
311: Input Parameter:
312: . pc - preconditioner context
314: Output Parameter:
315: . ksp - coarse problem `KSP` context
317: Level: advanced
319: .seealso: `PCDEFLATION`, `PCDeflationSetCoarseMat()`
320: @*/
321: PetscErrorCode PCDeflationGetCoarseKSP(PC pc, KSP *ksp)
322: {
325: PetscTryMethod(pc, "PCDeflationGetCoarseKSP_C", (PC, KSP *), (pc, ksp));
326: return 0;
327: }
329: static PetscErrorCode PCDeflationGetPC_Deflation(PC pc, PC *apc)
330: {
331: PC_Deflation *def = (PC_Deflation *)pc->data;
333: *apc = def->pc;
334: return 0;
335: }
337: /*@
338: PCDeflationGetPC - Returns the additional preconditioner M^{-1}.
340: Not Collective
342: Input Parameter:
343: . pc - the preconditioner context
345: Output Parameter:
346: . apc - additional preconditioner
348: Level: advanced
350: .seealso: `PCDEFLATION`, `PCDeflationGetCoarseKSP()`
351: @*/
352: PetscErrorCode PCDeflationGetPC(PC pc, PC *apc)
353: {
356: PetscTryMethod(pc, "PCDeflationGetPC_C", (PC, PC *), (pc, apc));
357: return 0;
358: }
360: /*
361: x <- x + W*(W'*A*W)^{-1}*W'*r = x + Q*r
362: */
363: static PetscErrorCode PCPreSolve_Deflation(PC pc, KSP ksp, Vec b, Vec x)
364: {
365: PC_Deflation *def = (PC_Deflation *)pc->data;
366: Mat A;
367: Vec r, w1, w2;
368: PetscBool nonzero;
370: w1 = def->workcoarse[0];
371: w2 = def->workcoarse[1];
372: r = def->work;
373: PCGetOperators(pc, NULL, &A);
375: KSPGetInitialGuessNonzero(ksp, &nonzero);
376: KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
377: if (nonzero) {
378: MatMult(A, x, r); /* r <- b - Ax */
379: VecAYPX(r, -1.0, b);
380: } else {
381: VecCopy(b, r); /* r <- b (x is 0) */
382: }
384: if (def->Wt) {
385: MatMult(def->Wt, r, w1); /* w1 <- W'*r */
386: } else {
387: MatMultHermitianTranspose(def->W, r, w1); /* w1 <- W'*r */
388: }
389: KSPSolve(def->WtAWinv, w1, w2); /* w2 <- (W'*A*W)^{-1}*w1 */
390: MatMult(def->W, w2, r); /* r <- W*w2 */
391: VecAYPX(x, 1.0, r);
392: return 0;
393: }
395: /*
396: if (def->correct) {
397: z <- M^{-1}r - W*(W'*A*W)^{-1}*(W'*A*M^{-1}r - l*W'*r) = (P*M^{-1} + l*Q)*r
398: } else {
399: z <- M^{-1}*r - W*(W'*A*W)^{-1}*W'*A*M{-1}*r = P*M^{-1}*r
400: }
401: */
402: static PetscErrorCode PCApply_Deflation(PC pc, Vec r, Vec z)
403: {
404: PC_Deflation *def = (PC_Deflation *)pc->data;
405: Mat A;
406: Vec u, w1, w2;
408: w1 = def->workcoarse[0];
409: w2 = def->workcoarse[1];
410: u = def->work;
411: PCGetOperators(pc, NULL, &A);
413: PCApply(def->pc, r, z); /* z <- M^{-1}*r */
414: if (!def->init) {
415: MatMult(def->WtA, z, w1); /* w1 <- W'*A*z */
416: if (def->correct) {
417: if (def->Wt) {
418: MatMult(def->Wt, r, w2); /* w2 <- W'*r */
419: } else {
420: MatMultHermitianTranspose(def->W, r, w2); /* w2 <- W'*r */
421: }
422: VecAXPY(w1, -1.0 * def->correctfact, w2); /* w1 <- w1 - l*w2 */
423: }
424: KSPSolve(def->WtAWinv, w1, w2); /* w2 <- (W'*A*W)^{-1}*w1 */
425: MatMult(def->W, w2, u); /* u <- W*w2 */
426: VecAXPY(z, -1.0, u); /* z <- z - u */
427: }
428: return 0;
429: }
431: static PetscErrorCode PCSetUp_Deflation(PC pc)
432: {
433: PC_Deflation *def = (PC_Deflation *)pc->data;
434: KSP innerksp;
435: PC pcinner;
436: Mat Amat, nextDef = NULL, *mats;
437: PetscInt i, m, red, size;
438: PetscMPIInt commsize;
439: PetscBool match, flgspd, isset, transp = PETSC_FALSE;
440: MatCompositeType ctype;
441: MPI_Comm comm;
442: char prefix[128] = "";
444: if (pc->setupcalled) return 0;
445: PetscObjectGetComm((PetscObject)pc, &comm);
446: PCGetOperators(pc, NULL, &Amat);
447: if (!def->lvl && !def->prefix) PCGetOptionsPrefix(pc, &def->prefix);
448: if (def->lvl) PetscSNPrintf(prefix, sizeof(prefix), "%d_", (int)def->lvl);
450: /* compute a deflation space */
451: if (def->W || def->Wt) {
452: def->spacetype = PC_DEFLATION_SPACE_USER;
453: } else {
454: PCDeflationComputeSpace(pc);
455: }
457: /* nested deflation */
458: if (def->W) {
459: PetscObjectTypeCompare((PetscObject)def->W, MATCOMPOSITE, &match);
460: if (match) {
461: MatCompositeGetType(def->W, &ctype);
462: MatCompositeGetNumberMat(def->W, &size);
463: }
464: } else {
465: MatCreateHermitianTranspose(def->Wt, &def->W);
466: PetscObjectTypeCompare((PetscObject)def->Wt, MATCOMPOSITE, &match);
467: if (match) {
468: MatCompositeGetType(def->Wt, &ctype);
469: MatCompositeGetNumberMat(def->Wt, &size);
470: }
471: transp = PETSC_TRUE;
472: }
473: if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) {
474: if (!transp) {
475: if (def->lvl < def->maxlvl) {
476: PetscMalloc1(size, &mats);
477: for (i = 0; i < size; i++) MatCompositeGetMat(def->W, i, &mats[i]);
478: size -= 1;
479: MatDestroy(&def->W);
480: def->W = mats[size];
481: PetscObjectReference((PetscObject)mats[size]);
482: if (size > 1) {
483: MatCreateComposite(comm, size, mats, &nextDef);
484: MatCompositeSetType(nextDef, MAT_COMPOSITE_MULTIPLICATIVE);
485: } else {
486: nextDef = mats[0];
487: PetscObjectReference((PetscObject)mats[0]);
488: }
489: PetscFree(mats);
490: } else {
491: /* TODO test merge side performance */
492: /* MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT); */
493: MatCompositeMerge(def->W);
494: }
495: } else {
496: if (def->lvl < def->maxlvl) {
497: PetscMalloc1(size, &mats);
498: for (i = 0; i < size; i++) MatCompositeGetMat(def->Wt, i, &mats[i]);
499: size -= 1;
500: MatDestroy(&def->Wt);
501: def->Wt = mats[0];
502: PetscObjectReference((PetscObject)mats[0]);
503: if (size > 1) {
504: MatCreateComposite(comm, size, &mats[1], &nextDef);
505: MatCompositeSetType(nextDef, MAT_COMPOSITE_MULTIPLICATIVE);
506: } else {
507: nextDef = mats[1];
508: PetscObjectReference((PetscObject)mats[1]);
509: }
510: PetscFree(mats);
511: } else {
512: /* MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT); */
513: MatCompositeMerge(def->Wt);
514: }
515: }
516: }
518: if (transp) {
519: MatDestroy(&def->W);
520: MatHermitianTranspose(def->Wt, MAT_INITIAL_MATRIX, &def->W);
521: }
523: /* assemble WtA */
524: if (!def->WtA) {
525: if (def->Wt) {
526: MatMatMult(def->Wt, Amat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &def->WtA);
527: } else {
528: #if defined(PETSC_USE_COMPLEX)
529: MatHermitianTranspose(def->W, MAT_INITIAL_MATRIX, &def->Wt);
530: MatMatMult(def->Wt, Amat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &def->WtA);
531: #else
532: MatTransposeMatMult(def->W, Amat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &def->WtA);
533: #endif
534: }
535: }
536: /* setup coarse problem */
537: if (!def->WtAWinv) {
538: MatGetSize(def->W, NULL, &m);
539: if (!def->WtAW) {
540: MatMatMult(def->WtA, def->W, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &def->WtAW);
541: /* TODO create MatInheritOption(Mat,MatOption) */
542: MatIsSPDKnown(Amat, &isset, &flgspd);
543: if (isset) MatSetOption(def->WtAW, MAT_SPD, flgspd);
544: if (PetscDefined(USE_DEBUG)) {
545: /* Check columns of W are not in kernel of A */
546: PetscReal *norms;
548: PetscMalloc1(m, &norms);
549: MatGetColumnNorms(def->WtAW, NORM_INFINITY, norms);
551: PetscFree(norms);
552: }
553: } else MatIsSPDKnown(def->WtAW, &isset, &flgspd);
555: /* TODO use MATINV ? */
556: KSPCreate(comm, &def->WtAWinv);
557: KSPSetOperators(def->WtAWinv, def->WtAW, def->WtAW);
558: KSPGetPC(def->WtAWinv, &pcinner);
559: /* Setup KSP and PC */
560: if (nextDef) { /* next level for multilevel deflation */
561: innerksp = def->WtAWinv;
562: /* set default KSPtype */
563: if (!def->ksptype) {
564: def->ksptype = KSPFGMRES;
565: if (isset && flgspd) { /* SPD system */
566: def->ksptype = KSPFCG;
567: }
568: }
569: KSPSetType(innerksp, def->ksptype); /* TODO iherit from KSP + tolerances */
570: PCSetType(pcinner, PCDEFLATION); /* TODO create coarse preconditinoner M_c = WtMW ? */
571: PCDeflationSetSpace(pcinner, nextDef, transp);
572: PCDeflationSetLevels_Deflation(pcinner, def->lvl + 1, def->maxlvl);
573: /* inherit options */
574: if (def->prefix) ((PC_Deflation *)(pcinner->data))->prefix = def->prefix;
575: ((PC_Deflation *)(pcinner->data))->init = def->init;
576: ((PC_Deflation *)(pcinner->data))->ksptype = def->ksptype;
577: ((PC_Deflation *)(pcinner->data))->correct = def->correct;
578: ((PC_Deflation *)(pcinner->data))->correctfact = def->correctfact;
579: ((PC_Deflation *)(pcinner->data))->reductionfact = def->reductionfact;
580: MatDestroy(&nextDef);
581: } else { /* the last level */
582: KSPSetType(def->WtAWinv, KSPPREONLY);
583: PCSetType(pcinner, PCTELESCOPE);
584: /* do not overwrite PCTELESCOPE */
585: if (def->prefix) KSPSetOptionsPrefix(def->WtAWinv, def->prefix);
586: KSPAppendOptionsPrefix(def->WtAWinv, "deflation_tel_");
587: PCSetFromOptions(pcinner);
588: PetscObjectTypeCompare((PetscObject)pcinner, PCTELESCOPE, &match);
590: /* Reduction factor choice */
591: red = def->reductionfact;
592: if (red < 0) {
593: MPI_Comm_size(comm, &commsize);
594: red = PetscCeilInt(commsize, PetscCeilInt(m, commsize));
595: PetscObjectTypeCompareAny((PetscObject)(def->WtAW), &match, MATSEQDENSE, MATMPIDENSE, MATDENSE, "");
596: if (match) red = commsize;
597: PetscInfo(pc, "Auto choosing reduction factor %" PetscInt_FMT "\n", red);
598: }
599: PCTelescopeSetReductionFactor(pcinner, red);
600: PCSetUp(pcinner);
601: PCTelescopeGetKSP(pcinner, &innerksp);
602: if (innerksp) {
603: KSPGetPC(innerksp, &pcinner);
604: PCSetType(pcinner, PCLU);
605: #if defined(PETSC_HAVE_SUPERLU)
606: MatGetFactorAvailable(def->WtAW, MATSOLVERSUPERLU, MAT_FACTOR_LU, &match);
607: if (match) PCFactorSetMatSolverType(pcinner, MATSOLVERSUPERLU);
608: #endif
609: #if defined(PETSC_HAVE_SUPERLU_DIST)
610: MatGetFactorAvailable(def->WtAW, MATSOLVERSUPERLU_DIST, MAT_FACTOR_LU, &match);
611: if (match) PCFactorSetMatSolverType(pcinner, MATSOLVERSUPERLU_DIST);
612: #endif
613: }
614: }
616: if (innerksp) {
617: if (def->prefix) {
618: KSPSetOptionsPrefix(innerksp, def->prefix);
619: KSPAppendOptionsPrefix(innerksp, "deflation_");
620: } else {
621: KSPSetOptionsPrefix(innerksp, "deflation_");
622: }
623: KSPAppendOptionsPrefix(innerksp, prefix);
624: KSPSetFromOptions(innerksp);
625: KSPSetUp(innerksp);
626: }
627: }
628: KSPSetFromOptions(def->WtAWinv);
629: KSPSetUp(def->WtAWinv);
631: /* create preconditioner */
632: if (!def->pc) {
633: PCCreate(comm, &def->pc);
634: PCSetOperators(def->pc, Amat, Amat);
635: PCSetType(def->pc, PCNONE);
636: if (def->prefix) PCSetOptionsPrefix(def->pc, def->prefix);
637: PCAppendOptionsPrefix(def->pc, "deflation_");
638: PCAppendOptionsPrefix(def->pc, prefix);
639: PCAppendOptionsPrefix(def->pc, "pc_");
640: PCSetFromOptions(def->pc);
641: PCSetUp(def->pc);
642: }
644: /* create work vecs */
645: MatCreateVecs(Amat, NULL, &def->work);
646: KSPCreateVecs(def->WtAWinv, 2, &def->workcoarse, 0, NULL);
647: return 0;
648: }
650: static PetscErrorCode PCReset_Deflation(PC pc)
651: {
652: PC_Deflation *def = (PC_Deflation *)pc->data;
654: VecDestroy(&def->work);
655: VecDestroyVecs(2, &def->workcoarse);
656: MatDestroy(&def->W);
657: MatDestroy(&def->Wt);
658: MatDestroy(&def->WtA);
659: MatDestroy(&def->WtAW);
660: KSPDestroy(&def->WtAWinv);
661: PCDestroy(&def->pc);
662: return 0;
663: }
665: static PetscErrorCode PCDestroy_Deflation(PC pc)
666: {
667: PCReset_Deflation(pc);
668: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetInitOnly_C", NULL);
669: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetLevels_C", NULL);
670: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetReductionFactor_C", NULL);
671: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetCorrectionFactor_C", NULL);
672: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetSpaceToCompute_C", NULL);
673: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetSpace_C", NULL);
674: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetProjectionNullSpaceMat_C", NULL);
675: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetCoarseMat_C", NULL);
676: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationGetCoarseKSP_C", NULL);
677: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationGetPC_C", NULL);
678: PetscFree(pc->data);
679: return 0;
680: }
682: static PetscErrorCode PCView_Deflation(PC pc, PetscViewer viewer)
683: {
684: PC_Deflation *def = (PC_Deflation *)pc->data;
685: PetscInt its;
686: PetscBool iascii;
688: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
689: if (iascii) {
690: if (def->correct) PetscViewerASCIIPrintf(viewer, "using CP correction, factor = %g+%gi\n", (double)PetscRealPart(def->correctfact), (double)PetscImaginaryPart(def->correctfact));
691: if (!def->lvl) PetscViewerASCIIPrintf(viewer, "deflation space type: %s\n", PCDeflationSpaceTypes[def->spacetype]);
693: PetscViewerASCIIPrintf(viewer, "--- Additional PC:\n");
694: PetscViewerASCIIPushTab(viewer);
695: PCView(def->pc, viewer);
696: PetscViewerASCIIPopTab(viewer);
698: PetscViewerASCIIPrintf(viewer, "--- Coarse problem solver:\n");
699: PetscViewerASCIIPushTab(viewer);
700: KSPGetTotalIterations(def->WtAWinv, &its);
701: PetscViewerASCIIPrintf(viewer, "total number of iterations: %" PetscInt_FMT "\n", its);
702: KSPView(def->WtAWinv, viewer);
703: PetscViewerASCIIPopTab(viewer);
704: }
705: return 0;
706: }
708: static PetscErrorCode PCSetFromOptions_Deflation(PC pc, PetscOptionItems *PetscOptionsObject)
709: {
710: PC_Deflation *def = (PC_Deflation *)pc->data;
712: PetscOptionsHeadBegin(PetscOptionsObject, "Deflation options");
713: PetscOptionsBool("-pc_deflation_init_only", "Use only initialization step - Initdef", "PCDeflationSetInitOnly", def->init, &def->init, NULL);
714: PetscOptionsInt("-pc_deflation_levels", "Maximum of deflation levels", "PCDeflationSetLevels", def->maxlvl, &def->maxlvl, NULL);
715: PetscOptionsInt("-pc_deflation_reduction_factor", "Reduction factor for coarse problem solution using PCTELESCOPE", "PCDeflationSetReductionFactor", def->reductionfact, &def->reductionfact, NULL);
716: PetscOptionsBool("-pc_deflation_correction", "Add coarse problem correction Q to P", "PCDeflationSetCorrectionFactor", def->correct, &def->correct, NULL);
717: PetscOptionsScalar("-pc_deflation_correction_factor", "Set multiple of Q to use as coarse problem correction", "PCDeflationSetCorrectionFactor", def->correctfact, &def->correctfact, NULL);
718: PetscOptionsEnum("-pc_deflation_compute_space", "Compute deflation space", "PCDeflationSetSpace", PCDeflationSpaceTypes, (PetscEnum)def->spacetype, (PetscEnum *)&def->spacetype, NULL);
719: PetscOptionsInt("-pc_deflation_compute_space_size", "Set size of the deflation space to compute", "PCDeflationSetSpace", def->spacesize, &def->spacesize, NULL);
720: PetscOptionsBool("-pc_deflation_space_extend", "Extend deflation space instead of truncating (wavelets)", "PCDeflation", def->extendsp, &def->extendsp, NULL);
721: PetscOptionsHeadEnd();
722: return 0;
723: }
725: /*MC
726: PCDEFLATION - Deflation preconditioner shifts (deflates) part of the spectrum to zero or to a predefined value.
728: Options Database Keys:
729: + -pc_deflation_init_only <false> - if true computes only the special guess
730: . -pc_deflation_max_lvl <0> - maximum number of levels for multilevel deflation
731: . -pc_deflation_reduction_factor <\-1> - reduction factor on bottom level coarse problem for PCTELESCOPE (default based on the size of the coarse problem)
732: . -pc_deflation_correction <false> - if true apply coarse problem correction
733: . -pc_deflation_correction_factor <1.0> - sets coarse problem correction factor
734: . -pc_deflation_compute_space <haar> - compute PCDeflationSpaceType deflation space
735: - -pc_deflation_compute_space_size <1> - size of the deflation space (corresponds to number of levels for wavelet-based deflation)
737: Notes:
738: Given a (complex - transpose is always Hermitian) full rank deflation matrix W, the deflation (introduced in [1,2])
739: preconditioner uses projections Q = W*(W'*A*W)^{-1}*W' and P = I - Q*A, where A is pmat.
741: The deflation computes initial guess x0 = x_{-1} - Q*r_{-1}, which is the solution on the deflation space.
742: If `PCDeflationSetInitOnly()` or -pc_deflation_init_only is set to `PETSC_TRUE` (InitDef scheme), the application of the
743: preconditioner consists only of application of the additional preconditioner M^{-1}. Otherwise, the preconditioner
744: application consists of P*M^{-1} + factor*Q. The first part of the preconditioner (PM^{-1}) shifts some eigenvalues
745: to zero while the addition of the coarse problem correction (factor*Q) makes the preconditioner to shift some
746: eigenvalues to the given factor. The InitDef scheme is recommended for deflation using high accuracy estimates
747: of eigenvectors of A when it exhibits similar convergence to the full deflation but is cheaper.
749: The deflation matrix is by default automatically computed. The type of deflation matrix and its size to compute can
750: be controlled by `PCDeflationSetSpaceToCompute()` or -pc_deflation_compute_space and -pc_deflation_compute_space_size.
751: User can set an arbitrary deflation space matrix with `PCDeflationSetSpace()`. If the deflation matrix
752: is a multiplicative `MATCOMPOSITE`, a multilevel deflation [3] is used. The first matrix in the composite is used as the
753: deflation matrix, and the coarse problem (W'*A*W)^{-1} is solved by `KSPFCG` (if A is `MAT_SPD`) or `KSPFGMRES` preconditioned
754: by deflation with deflation matrix being the next matrix in the `MATCOMPOSITE`. This scheme repeats until the maximum
755: level is reached or there are no more matrices. If the maximum level is reached, the remaining matrices are merged
756: (multiplied) to create the last deflation matrix. The maximum level defaults to 0 and can be set by
757: `PCDeflationSetLevels()` or by -pc_deflation_levels.
759: The coarse problem `KSP` can be controlled from the command line with prefix -deflation_ for the first level and -deflation_[lvl-1]
760: from the second level onward. You can also use
761: `PCDeflationGetCoarseKSP()` to control it from code. The bottom level KSP defaults to
762: `KSPPREONLY` with `PCLU` direct solver (`MATSOLVERSUPERLU`/`MATSOLVERSUPERLU_DIST` if available) wrapped into `PCTELESCOPE`.
763: For convenience, the reduction factor can be set by `PCDeflationSetReductionFactor()`
764: or -pc_deflation_recduction_factor. The default is chosen heuristically based on the coarse problem size.
766: The additional preconditioner can be controlled from command line with prefix -deflation_[lvl]_pc (same rules used for
767: coarse problem `KSP` apply for [lvl]_ part of prefix), e.g., -deflation_1_pc_pc_type bjacobi. You can also use
768: `PCDeflationGetPC()` to control the additional preconditioner from code. It defaults to `PCNONE`.
770: The coarse problem correction term (factor*Q) can be turned on by -pc_deflation_correction and the factor value can
771: be set by pc_deflation_correction_factor or by `PCDeflationSetCorrectionFactor()`. The coarse problem can
772: significantly improve convergence when the deflation coarse problem is not solved with high enough accuracy. We
773: recommend setting factor to some eigenvalue, e.g., the largest eigenvalue so that the preconditioner does not create
774: an isolated eigenvalue.
776: The options are automatically inherited from the previous deflation level.
778: The preconditioner supports `KSPMonitorDynamicTolerance()`. This is useful for the multilevel scheme for which we also
779: recommend limiting the number of iterations for the coarse problems.
781: See section 3 of [4] for additional references and decription of the algorithm when used for conjugate gradients.
782: Section 4 describes some possible choices for the deflation space.
784: Contributed by Jakub Kruzik (PERMON), Institute of Geonics of the Czech
785: Academy of Sciences and VSB - TU Ostrava.
787: Developed from PERMON code used in [4] while on a research stay with
788: Prof. Reinhard Nabben at the Institute of Mathematics, TU Berlin.
790: References:
791: + * - A. Nicolaides. "Deflation of conjugate gradients with applications to boundary value problems", SIAM J. Numer. Anal. 24.2, 1987.
792: . * - Z. Dostal. "Conjugate gradient method with preconditioning by projector", Int J. Comput. Math. 23.3-4, 1988.
793: . * - Y. A. Erlangga and R. Nabben. "Multilevel Projection-Based Nested Krylov Iteration for Boundary Value Problems", SIAM J. Sci. Comput. 30.3, 2008.
794: - * - J. Kruzik "Implementation of the Deflated Variants of the Conjugate Gradient Method", Master's thesis, VSB-TUO, 2018 - http://dspace5.vsb.cz/bitstream/handle/10084/130303/KRU0097_USP_N2658_2612T078_2018.pdf
796: Level: intermediate
798: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
799: `PCDeflationSetInitOnly()`, `PCDeflationSetLevels()`, `PCDeflationSetReductionFactor()`,
800: `PCDeflationSetCorrectionFactor()`, `PCDeflationSetSpaceToCompute()`,
801: `PCDeflationSetSpace()`, `PCDeflationSpaceType`, `PCDeflationSetProjectionNullSpaceMat()`,
802: `PCDeflationSetCoarseMat()`, `PCDeflationGetCoarseKSP()`, `PCDeflationGetPC()`
803: M*/
805: PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc)
806: {
807: PC_Deflation *def;
809: PetscNew(&def);
810: pc->data = (void *)def;
812: def->init = PETSC_FALSE;
813: def->correct = PETSC_FALSE;
814: def->correctfact = 1.0;
815: def->reductionfact = -1;
816: def->spacetype = PC_DEFLATION_SPACE_HAAR;
817: def->spacesize = 1;
818: def->extendsp = PETSC_FALSE;
819: def->lvl = 0;
820: def->maxlvl = 0;
821: def->W = NULL;
822: def->Wt = NULL;
824: pc->ops->apply = PCApply_Deflation;
825: pc->ops->presolve = PCPreSolve_Deflation;
826: pc->ops->setup = PCSetUp_Deflation;
827: pc->ops->reset = PCReset_Deflation;
828: pc->ops->destroy = PCDestroy_Deflation;
829: pc->ops->setfromoptions = PCSetFromOptions_Deflation;
830: pc->ops->view = PCView_Deflation;
832: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetInitOnly_C", PCDeflationSetInitOnly_Deflation);
833: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetLevels_C", PCDeflationSetLevels_Deflation);
834: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetReductionFactor_C", PCDeflationSetReductionFactor_Deflation);
835: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetCorrectionFactor_C", PCDeflationSetCorrectionFactor_Deflation);
836: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetSpaceToCompute_C", PCDeflationSetSpaceToCompute_Deflation);
837: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetSpace_C", PCDeflationSetSpace_Deflation);
838: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetProjectionNullSpaceMat_C", PCDeflationSetProjectionNullSpaceMat_Deflation);
839: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetCoarseMat_C", PCDeflationSetCoarseMat_Deflation);
840: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationGetCoarseKSP_C", PCDeflationGetCoarseKSP_Deflation);
841: PetscObjectComposeFunction((PetscObject)pc, "PCDeflationGetPC_C", PCDeflationGetPC_Deflation);
842: return 0;
843: }