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: PetscFunctionBegin;
10: def->init = flg;
11: PetscFunctionReturn(PETSC_SUCCESS);
12: }
14: /*@
15: PCDeflationSetInitOnly - Do only initialization step.
16: Sets initial guess to the solution on the deflation space but does not apply
17: the deflation preconditioner. The additional preconditioner is still applied.
19: Logically Collective
21: Input Parameters:
22: + pc - the preconditioner context
23: - flg - default `PETSC_FALSE`
25: Options Database Key:
26: . -pc_deflation_init_only <false> - if true computes only the special guess
28: Level: intermediate
30: .seealso: [](ch_ksp), `PCDEFLATION`
31: @*/
32: PetscErrorCode PCDeflationSetInitOnly(PC pc, PetscBool flg)
33: {
34: PetscFunctionBegin;
37: PetscTryMethod(pc, "PCDeflationSetInitOnly_C", (PC, PetscBool), (pc, flg));
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: static PetscErrorCode PCDeflationSetLevels_Deflation(PC pc, PetscInt current, PetscInt max)
42: {
43: PC_Deflation *def = (PC_Deflation *)pc->data;
45: PetscFunctionBegin;
46: if (current) def->lvl = current;
47: def->maxlvl = max;
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: /*@
52: PCDeflationSetLevels - Set the maximum level of deflation nesting.
54: Logically Collective
56: Input Parameters:
57: + pc - the preconditioner context
58: - max - maximum deflation level
60: Options Database Key:
61: . -pc_deflation_max_lvl <0> - maximum number of levels for multilevel deflation
63: Level: intermediate
65: .seealso: [](ch_ksp), `PCDeflationSetSpaceToCompute()`, `PCDeflationSetSpace()`, `PCDEFLATION`
66: @*/
67: PetscErrorCode PCDeflationSetLevels(PC pc, PetscInt max)
68: {
69: PetscFunctionBegin;
72: PetscTryMethod(pc, "PCDeflationSetLevels_C", (PC, PetscInt, PetscInt), (pc, 0, max));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: static PetscErrorCode PCDeflationSetReductionFactor_Deflation(PC pc, PetscInt red)
77: {
78: PC_Deflation *def = (PC_Deflation *)pc->data;
80: PetscFunctionBegin;
81: def->reductionfact = red;
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: /*@
86: PCDeflationSetReductionFactor - Set reduction factor for the `PCDEFLATION`
88: Logically Collective
90: Input Parameters:
91: + pc - the preconditioner context
92: - red - reduction factor (or `PETSC_DETERMINE`)
94: Options Database Key:
95: . -pc_deflation_reduction_factor <\-1> - reduction factor on bottom level coarse problem for `PCDEFLATION`
97: Note:
98: Default is computed based on the size of the coarse problem.
100: Level: intermediate
102: .seealso: [](ch_ksp), `PCTELESCOPE`, `PCDEFLATION`, `PCDeflationSetLevels()`
103: @*/
104: PetscErrorCode PCDeflationSetReductionFactor(PC pc, PetscInt red)
105: {
106: PetscFunctionBegin;
109: PetscTryMethod(pc, "PCDeflationSetReductionFactor_C", (PC, PetscInt), (pc, red));
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: static PetscErrorCode PCDeflationSetCorrectionFactor_Deflation(PC pc, PetscScalar fact)
114: {
115: PC_Deflation *def = (PC_Deflation *)pc->data;
117: PetscFunctionBegin;
118: /* TODO PETSC_DETERMINE -> compute max eigenvalue with power method */
119: def->correct = PETSC_TRUE;
120: def->correctfact = fact;
121: if (def->correctfact == 0.0) def->correct = PETSC_FALSE;
122: PetscFunctionReturn(PETSC_SUCCESS);
123: }
125: /*@
126: PCDeflationSetCorrectionFactor - Set coarse problem correction factor.
127: The preconditioner becomes P*M^{-1} + fact*Q.
129: Logically Collective
131: Input Parameters:
132: + pc - the preconditioner context
133: - fact - correction factor
135: Options Database Keys:
136: + -pc_deflation_correction <false> - if true apply coarse problem correction
137: - -pc_deflation_correction_factor <1.0> - sets coarse problem correction factor
139: Note:
140: Any non-zero fact enables the coarse problem correction.
142: Level: intermediate
144: .seealso: [](ch_ksp), `PCDEFLATION`, `PCDeflationSetLevels()`, `PCDeflationSetReductionFactor()`
145: @*/
146: PetscErrorCode PCDeflationSetCorrectionFactor(PC pc, PetscScalar fact)
147: {
148: PetscFunctionBegin;
151: PetscTryMethod(pc, "PCDeflationSetCorrectionFactor_C", (PC, PetscScalar), (pc, fact));
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: static PetscErrorCode PCDeflationSetSpaceToCompute_Deflation(PC pc, PCDeflationSpaceType type, PetscInt size)
156: {
157: PC_Deflation *def = (PC_Deflation *)pc->data;
159: PetscFunctionBegin;
160: if (type) def->spacetype = type;
161: if (size > 0) def->spacesize = size;
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: /*@
166: PCDeflationSetSpaceToCompute - Set deflation space type and size to compute.
168: Logically Collective
170: Input Parameters:
171: + pc - the preconditioner context
172: . type - deflation space type to compute (or `PETSC_IGNORE`)
173: - size - size of the space to compute (or `PETSC_DEFAULT`)
175: Options Database Keys:
176: + -pc_deflation_compute_space <haar> - compute `PCDeflationSpaceType` deflation space
177: - -pc_deflation_compute_space_size <1> - size of the deflation space
179: Notes:
180: For wavelet-based deflation, size represents number of levels.
182: The deflation space is computed in `PCSetUp()`.
184: Level: intermediate
186: .seealso: [](ch_ksp), `PCDeflationSetLevels()`, `PCDEFLATION`
187: @*/
188: PetscErrorCode PCDeflationSetSpaceToCompute(PC pc, PCDeflationSpaceType type, PetscInt size)
189: {
190: PetscFunctionBegin;
194: PetscTryMethod(pc, "PCDeflationSetSpaceToCompute_C", (PC, PCDeflationSpaceType, PetscInt), (pc, type, size));
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }
198: static PetscErrorCode PCDeflationSetSpace_Deflation(PC pc, Mat W, PetscBool transpose)
199: {
200: PC_Deflation *def = (PC_Deflation *)pc->data;
202: PetscFunctionBegin;
203: /* possibly allows W' = Wt (which is valid but not tested) */
204: PetscCall(PetscObjectReference((PetscObject)W));
205: if (transpose) {
206: PetscCall(MatDestroy(&def->Wt));
207: def->Wt = W;
208: } else {
209: PetscCall(MatDestroy(&def->W));
210: def->W = W;
211: }
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: /*@
216: PCDeflationSetSpace - Set the deflation space matrix (or its (Hermitian) transpose).
218: Logically Collective
220: Input Parameters:
221: + pc - the preconditioner context
222: . W - deflation matrix
223: - transpose - indicates that W is an explicit transpose of the deflation matrix
225: Level: intermediate
227: Notes:
228: Setting W as a multipliplicative `MATCOMPOSITE` enables use of the multilevel
229: deflation. If W = W0*W1*W2*...*Wn, W0 is taken as the first deflation space and
230: the coarse problem (W0'*A*W0)^{-1} is again preconditioned by deflation with
231: W1 as the deflation matrix. This repeats until the maximum level set by
232: PCDeflationSetLevels() is reached or there are no more matrices available.
233: If there are matrices left after reaching the maximum level,
234: they are merged into a deflation matrix ...*W{n-1}*Wn.
236: .seealso: [](ch_ksp), `PCDeflationSetLevels()`, `PCDEFLATION`, `PCDeflationSetProjectionNullSpaceMat()`
237: @*/
238: PetscErrorCode PCDeflationSetSpace(PC pc, Mat W, PetscBool transpose)
239: {
240: PetscFunctionBegin;
244: PetscTryMethod(pc, "PCDeflationSetSpace_C", (PC, Mat, PetscBool), (pc, W, transpose));
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: static PetscErrorCode PCDeflationSetProjectionNullSpaceMat_Deflation(PC pc, Mat mat)
249: {
250: PC_Deflation *def = (PC_Deflation *)pc->data;
252: PetscFunctionBegin;
253: PetscCall(PetscObjectReference((PetscObject)mat));
254: PetscCall(MatDestroy(&def->WtA));
255: def->WtA = mat;
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: /*@
260: PCDeflationSetProjectionNullSpaceMat - Set the projection null space matrix (W'*A).
262: Collective
264: Input Parameters:
265: + pc - preconditioner context
266: - mat - projection null space matrix
268: Level: developer
270: .seealso: [](ch_ksp), `PCDEFLATION`, `PCDeflationSetSpace()`
271: @*/
272: PetscErrorCode PCDeflationSetProjectionNullSpaceMat(PC pc, Mat mat)
273: {
274: PetscFunctionBegin;
277: PetscTryMethod(pc, "PCDeflationSetProjectionNullSpaceMat_C", (PC, Mat), (pc, mat));
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: static PetscErrorCode PCDeflationSetCoarseMat_Deflation(PC pc, Mat mat)
282: {
283: PC_Deflation *def = (PC_Deflation *)pc->data;
285: PetscFunctionBegin;
286: PetscCall(PetscObjectReference((PetscObject)mat));
287: PetscCall(MatDestroy(&def->WtAW));
288: def->WtAW = mat;
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: /*@
293: PCDeflationSetCoarseMat - Set the coarse problem `Mat`.
295: Collective
297: Input Parameters:
298: + pc - preconditioner context
299: - mat - coarse problem mat
301: Level: developer
303: .seealso: [](ch_ksp), `PCDEFLATION`, `PCDeflationGetCoarseKSP()`
304: @*/
305: PetscErrorCode PCDeflationSetCoarseMat(PC pc, Mat mat)
306: {
307: PetscFunctionBegin;
310: PetscTryMethod(pc, "PCDeflationSetCoarseMat_C", (PC, Mat), (pc, mat));
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: static PetscErrorCode PCDeflationGetCoarseKSP_Deflation(PC pc, KSP *ksp)
315: {
316: PC_Deflation *def = (PC_Deflation *)pc->data;
318: PetscFunctionBegin;
319: *ksp = def->WtAWinv;
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: /*@
324: PCDeflationGetCoarseKSP - Returns the coarse problem `KSP`.
326: Not Collective
328: Input Parameter:
329: . pc - preconditioner context
331: Output Parameter:
332: . ksp - coarse problem `KSP` context
334: Level: advanced
336: .seealso: [](ch_ksp), `PCDEFLATION`, `PCDeflationSetCoarseMat()`
337: @*/
338: PetscErrorCode PCDeflationGetCoarseKSP(PC pc, KSP *ksp)
339: {
340: PetscFunctionBegin;
342: PetscAssertPointer(ksp, 2);
343: PetscTryMethod(pc, "PCDeflationGetCoarseKSP_C", (PC, KSP *), (pc, ksp));
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: static PetscErrorCode PCDeflationGetPC_Deflation(PC pc, PC *apc)
348: {
349: PC_Deflation *def = (PC_Deflation *)pc->data;
351: PetscFunctionBegin;
352: *apc = def->pc;
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: /*@
357: PCDeflationGetPC - Returns the additional preconditioner M^{-1}.
359: Not Collective
361: Input Parameter:
362: . pc - the preconditioner context
364: Output Parameter:
365: . apc - additional preconditioner
367: Level: advanced
369: .seealso: [](ch_ksp), `PCDEFLATION`, `PCDeflationGetCoarseKSP()`
370: @*/
371: PetscErrorCode PCDeflationGetPC(PC pc, PC *apc)
372: {
373: PetscFunctionBegin;
375: PetscAssertPointer(pc, 1);
376: PetscTryMethod(pc, "PCDeflationGetPC_C", (PC, PC *), (pc, apc));
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: /*
381: x <- x + W*(W'*A*W)^{-1}*W'*r = x + Q*r
382: */
383: static PetscErrorCode PCPreSolve_Deflation(PC pc, KSP ksp, Vec b, Vec x)
384: {
385: PC_Deflation *def = (PC_Deflation *)pc->data;
386: Mat A;
387: Vec r, w1, w2;
388: PetscBool nonzero;
390: PetscFunctionBegin;
391: w1 = def->workcoarse[0];
392: w2 = def->workcoarse[1];
393: r = def->work;
394: PetscCall(PCGetOperators(pc, NULL, &A));
396: PetscCall(KSPGetInitialGuessNonzero(ksp, &nonzero));
397: PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
398: if (nonzero) {
399: PetscCall(MatMult(A, x, r)); /* r <- b - Ax */
400: PetscCall(VecAYPX(r, -1.0, b));
401: } else {
402: PetscCall(VecCopy(b, r)); /* r <- b (x is 0) */
403: }
405: if (def->Wt) {
406: PetscCall(MatMult(def->Wt, r, w1)); /* w1 <- W'*r */
407: } else {
408: PetscCall(MatMultHermitianTranspose(def->W, r, w1)); /* w1 <- W'*r */
409: }
410: PetscCall(KSPSolve(def->WtAWinv, w1, w2)); /* w2 <- (W'*A*W)^{-1}*w1 */
411: PetscCall(MatMult(def->W, w2, r)); /* r <- W*w2 */
412: PetscCall(VecAYPX(x, 1.0, r));
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: /*
417: if (def->correct) {
418: z <- M^{-1}r - W*(W'*A*W)^{-1}*(W'*A*M^{-1}r - l*W'*r) = (P*M^{-1} + l*Q)*r
419: } else {
420: z <- M^{-1}*r - W*(W'*A*W)^{-1}*W'*A*M{-1}*r = P*M^{-1}*r
421: }
422: */
423: static PetscErrorCode PCApply_Deflation(PC pc, Vec r, Vec z)
424: {
425: PC_Deflation *def = (PC_Deflation *)pc->data;
426: Mat A;
427: Vec u, w1, w2;
429: PetscFunctionBegin;
430: w1 = def->workcoarse[0];
431: w2 = def->workcoarse[1];
432: u = def->work;
433: PetscCall(PCGetOperators(pc, NULL, &A));
435: PetscCall(PCApply(def->pc, r, z)); /* z <- M^{-1}*r */
436: if (!def->init) {
437: PetscCall(MatMult(def->WtA, z, w1)); /* w1 <- W'*A*z */
438: if (def->correct) {
439: if (def->Wt) {
440: PetscCall(MatMult(def->Wt, r, w2)); /* w2 <- W'*r */
441: } else {
442: PetscCall(MatMultHermitianTranspose(def->W, r, w2)); /* w2 <- W'*r */
443: }
444: PetscCall(VecAXPY(w1, -1.0 * def->correctfact, w2)); /* w1 <- w1 - l*w2 */
445: }
446: PetscCall(KSPSolve(def->WtAWinv, w1, w2)); /* w2 <- (W'*A*W)^{-1}*w1 */
447: PetscCall(MatMult(def->W, w2, u)); /* u <- W*w2 */
448: PetscCall(VecAXPY(z, -1.0, u)); /* z <- z - u */
449: }
450: PetscFunctionReturn(PETSC_SUCCESS);
451: }
453: static PetscErrorCode PCSetUp_Deflation(PC pc)
454: {
455: PC_Deflation *def = (PC_Deflation *)pc->data;
456: DM dm;
457: KSP innerksp;
458: PC pcinner;
459: Mat Amat, nextDef = NULL, *mats;
460: PetscInt i, m, red, size;
461: PetscMPIInt commsize;
462: PetscBool match, flgspd, isset, transp = PETSC_FALSE;
463: MatCompositeType ctype;
464: MPI_Comm comm;
465: char prefix[128] = "";
467: PetscFunctionBegin;
468: if (pc->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
469: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
470: PetscCall(PCGetOperators(pc, NULL, &Amat));
471: if (!def->lvl && !def->prefix) PetscCall(PCGetOptionsPrefix(pc, &def->prefix));
472: if (def->lvl) PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%d_", (int)def->lvl));
474: /* compute a deflation space */
475: if (def->W || def->Wt) {
476: def->spacetype = PC_DEFLATION_SPACE_USER;
477: } else {
478: PetscCall(PCDeflationComputeSpace(pc));
479: }
481: /* nested deflation */
482: if (def->W) {
483: PetscCall(PetscObjectTypeCompare((PetscObject)def->W, MATCOMPOSITE, &match));
484: if (match) {
485: PetscCall(MatCompositeGetType(def->W, &ctype));
486: PetscCall(MatCompositeGetNumberMat(def->W, &size));
487: }
488: } else {
489: PetscCall(MatCreateHermitianTranspose(def->Wt, &def->W));
490: PetscCall(PetscObjectTypeCompare((PetscObject)def->Wt, MATCOMPOSITE, &match));
491: if (match) {
492: PetscCall(MatCompositeGetType(def->Wt, &ctype));
493: PetscCall(MatCompositeGetNumberMat(def->Wt, &size));
494: }
495: transp = PETSC_TRUE;
496: }
497: if (match && ctype == MAT_COMPOSITE_MULTIPLICATIVE) {
498: if (!transp) {
499: if (def->lvl < def->maxlvl) {
500: PetscCall(PetscMalloc1(size, &mats));
501: for (i = 0; i < size; i++) PetscCall(MatCompositeGetMat(def->W, i, &mats[i]));
502: size -= 1;
503: PetscCall(MatDestroy(&def->W));
504: def->W = mats[size];
505: PetscCall(PetscObjectReference((PetscObject)mats[size]));
506: if (size > 1) {
507: PetscCall(MatCreateComposite(comm, size, mats, &nextDef));
508: PetscCall(MatCompositeSetType(nextDef, MAT_COMPOSITE_MULTIPLICATIVE));
509: } else {
510: nextDef = mats[0];
511: PetscCall(PetscObjectReference((PetscObject)mats[0]));
512: }
513: PetscCall(PetscFree(mats));
514: } else {
515: /* TODO test merge side performance */
516: /* PetscCall(MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT)); */
517: PetscCall(MatCompositeMerge(def->W));
518: }
519: } else {
520: if (def->lvl < def->maxlvl) {
521: PetscCall(PetscMalloc1(size, &mats));
522: for (i = 0; i < size; i++) PetscCall(MatCompositeGetMat(def->Wt, i, &mats[i]));
523: size -= 1;
524: PetscCall(MatDestroy(&def->Wt));
525: def->Wt = mats[0];
526: PetscCall(PetscObjectReference((PetscObject)mats[0]));
527: if (size > 1) {
528: PetscCall(MatCreateComposite(comm, size, &mats[1], &nextDef));
529: PetscCall(MatCompositeSetType(nextDef, MAT_COMPOSITE_MULTIPLICATIVE));
530: } else {
531: nextDef = mats[1];
532: PetscCall(PetscObjectReference((PetscObject)mats[1]));
533: }
534: PetscCall(PetscFree(mats));
535: } else {
536: /* PetscCall(MatCompositeSetMergeType(def->W,MAT_COMPOSITE_MERGE_LEFT)); */
537: PetscCall(MatCompositeMerge(def->Wt));
538: }
539: }
540: }
542: if (transp) {
543: PetscCall(MatDestroy(&def->W));
544: PetscCall(MatHermitianTranspose(def->Wt, MAT_INITIAL_MATRIX, &def->W));
545: }
547: /* assemble WtA */
548: if (!def->WtA) {
549: if (def->Wt) {
550: PetscCall(MatMatMult(def->Wt, Amat, MAT_INITIAL_MATRIX, PETSC_CURRENT, &def->WtA));
551: } else {
552: #if defined(PETSC_USE_COMPLEX)
553: PetscCall(MatHermitianTranspose(def->W, MAT_INITIAL_MATRIX, &def->Wt));
554: PetscCall(MatMatMult(def->Wt, Amat, MAT_INITIAL_MATRIX, PETSC_CURRENT, &def->WtA));
555: #else
556: PetscCall(MatTransposeMatMult(def->W, Amat, MAT_INITIAL_MATRIX, PETSC_CURRENT, &def->WtA));
557: #endif
558: }
559: }
560: /* setup coarse problem */
561: if (!def->WtAWinv) {
562: PetscCall(MatGetSize(def->W, NULL, &m));
563: if (!def->WtAW) {
564: PetscCall(MatMatMult(def->WtA, def->W, MAT_INITIAL_MATRIX, PETSC_CURRENT, &def->WtAW));
565: /* TODO create MatInheritOption(Mat,MatOption) */
566: PetscCall(MatIsSPDKnown(Amat, &isset, &flgspd));
567: if (isset) PetscCall(MatSetOption(def->WtAW, MAT_SPD, flgspd));
568: if (PetscDefined(USE_DEBUG)) {
569: /* Check columns of W are not in kernel of A */
570: PetscReal *norms;
572: PetscCall(PetscMalloc1(m, &norms));
573: PetscCall(MatGetColumnNorms(def->WtAW, NORM_INFINITY, norms));
574: for (i = 0; i < m; i++) PetscCheck(norms[i] > 100 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)def->WtAW), PETSC_ERR_SUP, "Column %" PetscInt_FMT " of W is in kernel of A.", i);
575: PetscCall(PetscFree(norms));
576: }
577: } else PetscCall(MatIsSPDKnown(def->WtAW, &isset, &flgspd));
579: /* TODO use MATINV ? */
580: PetscCall(KSPCreate(comm, &def->WtAWinv));
581: PetscCall(KSPSetNestLevel(def->WtAWinv, pc->kspnestlevel));
582: PetscCall(KSPSetOperators(def->WtAWinv, def->WtAW, def->WtAW));
583: PetscCall(KSPGetPC(def->WtAWinv, &pcinner));
584: /* Setup KSP and PC */
585: if (nextDef) { /* next level for multilevel deflation */
586: innerksp = def->WtAWinv;
587: /* set default KSPtype */
588: if (!def->ksptype) {
589: def->ksptype = KSPFGMRES;
590: if (isset && flgspd) { /* SPD system */
591: def->ksptype = KSPFCG;
592: }
593: }
594: PetscCall(KSPSetType(innerksp, def->ksptype)); /* TODO iherit from KSP + tolerances */
595: PetscCall(PCSetType(pcinner, PCDEFLATION)); /* TODO create coarse preconditinoner M_c = WtMW ? */
596: PetscCall(PCDeflationSetSpace(pcinner, nextDef, transp));
597: PetscCall(PCDeflationSetLevels_Deflation(pcinner, def->lvl + 1, def->maxlvl));
598: /* inherit options */
599: if (def->prefix) ((PC_Deflation *)pcinner->data)->prefix = def->prefix;
600: ((PC_Deflation *)pcinner->data)->init = def->init;
601: ((PC_Deflation *)pcinner->data)->ksptype = def->ksptype;
602: ((PC_Deflation *)pcinner->data)->correct = def->correct;
603: ((PC_Deflation *)pcinner->data)->correctfact = def->correctfact;
604: ((PC_Deflation *)pcinner->data)->reductionfact = def->reductionfact;
605: PetscCall(MatDestroy(&nextDef));
606: } else { /* the last level */
607: PetscCall(KSPSetType(def->WtAWinv, KSPPREONLY));
608: PetscCall(PCSetType(pcinner, PCTELESCOPE));
609: /* do not overwrite PCTELESCOPE */
610: if (def->prefix) PetscCall(KSPSetOptionsPrefix(def->WtAWinv, def->prefix));
611: PetscCall(KSPAppendOptionsPrefix(def->WtAWinv, "deflation_tel_"));
612: PetscCall(PCSetFromOptions(pcinner));
613: PetscCall(PetscObjectTypeCompare((PetscObject)pcinner, PCTELESCOPE, &match));
614: PetscCheck(match, comm, PETSC_ERR_SUP, "User can not overwrite PCTELESCOPE on bottom level, use reduction factor = 1 instead.");
615: /* Reduction factor choice */
616: red = def->reductionfact;
617: if (red < 0) {
618: PetscCallMPI(MPI_Comm_size(comm, &commsize));
619: red = PetscCeilInt(commsize, PetscCeilInt(m, commsize));
620: PetscCall(PetscObjectTypeCompareAny((PetscObject)def->WtAW, &match, MATSEQDENSE, MATMPIDENSE, MATDENSE, ""));
621: if (match) red = commsize;
622: PetscCall(PetscInfo(pc, "Auto choosing reduction factor %" PetscInt_FMT "\n", red));
623: }
624: PetscCall(PCTelescopeSetReductionFactor(pcinner, red));
625: PetscCall(PCSetUp(pcinner));
626: PetscCall(PCTelescopeGetKSP(pcinner, &innerksp));
627: if (innerksp) {
628: PetscCall(KSPGetPC(innerksp, &pcinner));
629: PetscCall(PCSetType(pcinner, PCLU));
630: #if defined(PETSC_HAVE_SUPERLU)
631: PetscCall(MatGetFactorAvailable(def->WtAW, MATSOLVERSUPERLU, MAT_FACTOR_LU, &match));
632: if (match) PetscCall(PCFactorSetMatSolverType(pcinner, MATSOLVERSUPERLU));
633: #endif
634: #if defined(PETSC_HAVE_SUPERLU_DIST)
635: PetscCall(MatGetFactorAvailable(def->WtAW, MATSOLVERSUPERLU_DIST, MAT_FACTOR_LU, &match));
636: if (match) PetscCall(PCFactorSetMatSolverType(pcinner, MATSOLVERSUPERLU_DIST));
637: #endif
638: }
639: }
641: if (innerksp) {
642: if (def->prefix) {
643: PetscCall(KSPSetOptionsPrefix(innerksp, def->prefix));
644: PetscCall(KSPAppendOptionsPrefix(innerksp, "deflation_"));
645: } else {
646: PetscCall(KSPSetOptionsPrefix(innerksp, "deflation_"));
647: }
648: PetscCall(KSPAppendOptionsPrefix(innerksp, prefix));
649: PetscCall(KSPSetFromOptions(innerksp));
650: PetscCall(KSPSetUp(innerksp));
651: }
652: }
653: PetscCall(KSPSetFromOptions(def->WtAWinv));
654: PetscCall(KSPSetUp(def->WtAWinv));
656: /* create preconditioner */
657: if (!def->pc) {
658: PetscCall(PCCreate(comm, &def->pc));
659: PetscCall(PCSetOperators(def->pc, Amat, Amat));
660: PetscCall(PCSetType(def->pc, PCNONE));
661: if (def->prefix) PetscCall(PCSetOptionsPrefix(def->pc, def->prefix));
662: PetscCall(PCAppendOptionsPrefix(def->pc, "deflation_"));
663: PetscCall(PCAppendOptionsPrefix(def->pc, prefix));
664: PetscCall(PCAppendOptionsPrefix(def->pc, "pc_"));
665: PetscCall(PCGetDM(pc, &dm));
666: PetscCall(PCSetDM(def->pc, dm));
667: PetscCall(PCSetFromOptions(def->pc));
668: PetscCall(PCSetUp(def->pc));
669: }
671: /* create work vecs */
672: PetscCall(MatCreateVecs(Amat, NULL, &def->work));
673: PetscCall(KSPCreateVecs(def->WtAWinv, 2, &def->workcoarse, 0, NULL));
674: PetscFunctionReturn(PETSC_SUCCESS);
675: }
677: static PetscErrorCode PCReset_Deflation(PC pc)
678: {
679: PC_Deflation *def = (PC_Deflation *)pc->data;
681: PetscFunctionBegin;
682: PetscCall(VecDestroy(&def->work));
683: PetscCall(VecDestroyVecs(2, &def->workcoarse));
684: PetscCall(MatDestroy(&def->W));
685: PetscCall(MatDestroy(&def->Wt));
686: PetscCall(MatDestroy(&def->WtA));
687: PetscCall(MatDestroy(&def->WtAW));
688: PetscCall(KSPDestroy(&def->WtAWinv));
689: PetscCall(PCDestroy(&def->pc));
690: PetscFunctionReturn(PETSC_SUCCESS);
691: }
693: static PetscErrorCode PCDestroy_Deflation(PC pc)
694: {
695: PetscFunctionBegin;
696: PetscCall(PCReset_Deflation(pc));
697: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetInitOnly_C", NULL));
698: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetLevels_C", NULL));
699: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetReductionFactor_C", NULL));
700: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetCorrectionFactor_C", NULL));
701: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetSpaceToCompute_C", NULL));
702: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetSpace_C", NULL));
703: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetProjectionNullSpaceMat_C", NULL));
704: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetCoarseMat_C", NULL));
705: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationGetCoarseKSP_C", NULL));
706: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationGetPC_C", NULL));
707: PetscCall(PetscFree(pc->data));
708: PetscFunctionReturn(PETSC_SUCCESS);
709: }
711: static PetscErrorCode PCView_Deflation(PC pc, PetscViewer viewer)
712: {
713: PC_Deflation *def = (PC_Deflation *)pc->data;
714: PetscInt its;
715: PetscBool iascii;
717: PetscFunctionBegin;
718: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
719: if (iascii) {
720: if (def->correct) PetscCall(PetscViewerASCIIPrintf(viewer, "using CP correction, factor = %g+%gi\n", (double)PetscRealPart(def->correctfact), (double)PetscImaginaryPart(def->correctfact)));
721: if (!def->lvl) PetscCall(PetscViewerASCIIPrintf(viewer, "deflation space type: %s\n", PCDeflationSpaceTypes[def->spacetype]));
723: PetscCall(PetscViewerASCIIPrintf(viewer, "--- Additional PC:\n"));
724: PetscCall(PetscViewerASCIIPushTab(viewer));
725: PetscCall(PCView(def->pc, viewer));
726: PetscCall(PetscViewerASCIIPopTab(viewer));
728: PetscCall(PetscViewerASCIIPrintf(viewer, "--- Coarse problem solver:\n"));
729: PetscCall(PetscViewerASCIIPushTab(viewer));
730: PetscCall(KSPGetTotalIterations(def->WtAWinv, &its));
731: PetscCall(PetscViewerASCIIPrintf(viewer, "total number of iterations: %" PetscInt_FMT "\n", its));
732: PetscCall(KSPView(def->WtAWinv, viewer));
733: PetscCall(PetscViewerASCIIPopTab(viewer));
734: }
735: PetscFunctionReturn(PETSC_SUCCESS);
736: }
738: static PetscErrorCode PCSetFromOptions_Deflation(PC pc, PetscOptionItems *PetscOptionsObject)
739: {
740: PC_Deflation *def = (PC_Deflation *)pc->data;
742: PetscFunctionBegin;
743: PetscOptionsHeadBegin(PetscOptionsObject, "Deflation options");
744: PetscCall(PetscOptionsBool("-pc_deflation_init_only", "Use only initialization step - Initdef", "PCDeflationSetInitOnly", def->init, &def->init, NULL));
745: PetscCall(PetscOptionsInt("-pc_deflation_levels", "Maximum of deflation levels", "PCDeflationSetLevels", def->maxlvl, &def->maxlvl, NULL));
746: PetscCall(PetscOptionsInt("-pc_deflation_reduction_factor", "Reduction factor for coarse problem solution using PCTELESCOPE", "PCDeflationSetReductionFactor", def->reductionfact, &def->reductionfact, NULL));
747: PetscCall(PetscOptionsBool("-pc_deflation_correction", "Add coarse problem correction Q to P", "PCDeflationSetCorrectionFactor", def->correct, &def->correct, NULL));
748: PetscCall(PetscOptionsScalar("-pc_deflation_correction_factor", "Set multiple of Q to use as coarse problem correction", "PCDeflationSetCorrectionFactor", def->correctfact, &def->correctfact, NULL));
749: PetscCall(PetscOptionsEnum("-pc_deflation_compute_space", "Compute deflation space", "PCDeflationSetSpace", PCDeflationSpaceTypes, (PetscEnum)def->spacetype, (PetscEnum *)&def->spacetype, NULL));
750: PetscCall(PetscOptionsInt("-pc_deflation_compute_space_size", "Set size of the deflation space to compute", "PCDeflationSetSpace", def->spacesize, &def->spacesize, NULL));
751: PetscCall(PetscOptionsBool("-pc_deflation_space_extend", "Extend deflation space instead of truncating (wavelets)", "PCDeflation", def->extendsp, &def->extendsp, NULL));
752: PetscOptionsHeadEnd();
753: PetscFunctionReturn(PETSC_SUCCESS);
754: }
756: /*MC
757: PCDEFLATION - Deflation preconditioner shifts (deflates) part of the spectrum to zero or to a predefined value.
759: Options Database Keys:
760: + -pc_deflation_init_only <false> - if true computes only the special guess
761: . -pc_deflation_max_lvl <0> - maximum number of levels for multilevel deflation
762: . -pc_deflation_reduction_factor <\-1> - reduction factor on bottom level coarse problem for PCTELESCOPE (default based on the size of the coarse problem)
763: . -pc_deflation_correction <false> - if true apply coarse problem correction
764: . -pc_deflation_correction_factor <1.0> - sets coarse problem correction factor
765: . -pc_deflation_compute_space <haar> - compute PCDeflationSpaceType deflation space
766: - -pc_deflation_compute_space_size <1> - size of the deflation space (corresponds to number of levels for wavelet-based deflation)
768: Notes:
769: Given a (complex - transpose is always Hermitian) full rank deflation matrix W, the deflation (introduced in [1,2])
770: preconditioner uses projections Q = W*(W'*A*W)^{-1}*W' and P = I - Q*A, where A is pmat.
772: The deflation computes initial guess x0 = x_{-1} - Q*r_{-1}, which is the solution on the deflation space.
773: If `PCDeflationSetInitOnly()` or -pc_deflation_init_only is set to `PETSC_TRUE` (InitDef scheme), the application of the
774: preconditioner consists only of application of the additional preconditioner M^{-1}. Otherwise, the preconditioner
775: application consists of P*M^{-1} + factor*Q. The first part of the preconditioner (PM^{-1}) shifts some eigenvalues
776: to zero while the addition of the coarse problem correction (factor*Q) makes the preconditioner to shift some
777: eigenvalues to the given factor. The InitDef scheme is recommended for deflation using high accuracy estimates
778: of eigenvectors of A when it exhibits similar convergence to the full deflation but is cheaper.
780: The deflation matrix is by default automatically computed. The type of deflation matrix and its size to compute can
781: be controlled by `PCDeflationSetSpaceToCompute()` or -pc_deflation_compute_space and -pc_deflation_compute_space_size.
782: User can set an arbitrary deflation space matrix with `PCDeflationSetSpace()`. If the deflation matrix
783: is a multiplicative `MATCOMPOSITE`, a multilevel deflation [3] is used. The first matrix in the composite is used as the
784: deflation matrix, and the coarse problem (W'*A*W)^{-1} is solved by `KSPFCG` (if A is `MAT_SPD`) or `KSPFGMRES` preconditioned
785: by deflation with deflation matrix being the next matrix in the `MATCOMPOSITE`. This scheme repeats until the maximum
786: level is reached or there are no more matrices. If the maximum level is reached, the remaining matrices are merged
787: (multiplied) to create the last deflation matrix. The maximum level defaults to 0 and can be set by
788: `PCDeflationSetLevels()` or by -pc_deflation_levels.
790: The coarse problem `KSP` can be controlled from the command line with prefix -deflation_ for the first level and -deflation_[lvl-1]
791: from the second level onward. You can also use
792: `PCDeflationGetCoarseKSP()` to control it from code. The bottom level KSP defaults to
793: `KSPPREONLY` with `PCLU` direct solver (`MATSOLVERSUPERLU`/`MATSOLVERSUPERLU_DIST` if available) wrapped into `PCTELESCOPE`.
794: For convenience, the reduction factor can be set by `PCDeflationSetReductionFactor()`
795: or -pc_deflation_recduction_factor. The default is chosen heuristically based on the coarse problem size.
797: The additional preconditioner can be controlled from command line with prefix -deflation_[lvl]_pc (same rules used for
798: coarse problem `KSP` apply for [lvl]_ part of prefix), e.g., -deflation_1_pc_pc_type bjacobi. You can also use
799: `PCDeflationGetPC()` to control the additional preconditioner from code. It defaults to `PCNONE`.
801: The coarse problem correction term (factor*Q) can be turned on by -pc_deflation_correction and the factor value can
802: be set by pc_deflation_correction_factor or by `PCDeflationSetCorrectionFactor()`. The coarse problem can
803: significantly improve convergence when the deflation coarse problem is not solved with high enough accuracy. We
804: recommend setting factor to some eigenvalue, e.g., the largest eigenvalue so that the preconditioner does not create
805: an isolated eigenvalue.
807: The options are automatically inherited from the previous deflation level.
809: The preconditioner supports `KSPMonitorDynamicTolerance()`. This is useful for the multilevel scheme for which we also
810: recommend limiting the number of iterations for the coarse problems.
812: See section 3 of [4] for additional references and description of the algorithm when used for conjugate gradients.
813: Section 4 describes some possible choices for the deflation space.
815: Contributed by Jakub Kruzik (PERMON), Institute of Geonics of the Czech
816: Academy of Sciences and VSB - TU Ostrava.
818: Developed from PERMON code used in [4] while on a research stay with
819: Prof. Reinhard Nabben at the Institute of Mathematics, TU Berlin.
821: References:
822: + * - A. Nicolaides. "Deflation of conjugate gradients with applications to boundary value problems", SIAM J. Numer. Anal. 24.2, 1987.
823: . * - Z. Dostal. "Conjugate gradient method with preconditioning by projector", Int J. Comput. Math. 23.3-4, 1988.
824: . * - Y. A. Erlangga and R. Nabben. "Multilevel Projection-Based Nested Krylov Iteration for Boundary Value Problems", SIAM J. Sci. Comput. 30.3, 2008.
825: - * - 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
827: Level: intermediate
829: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
830: `PCDeflationSetInitOnly()`, `PCDeflationSetLevels()`, `PCDeflationSetReductionFactor()`,
831: `PCDeflationSetCorrectionFactor()`, `PCDeflationSetSpaceToCompute()`,
832: `PCDeflationSetSpace()`, `PCDeflationSpaceType`, `PCDeflationSetProjectionNullSpaceMat()`,
833: `PCDeflationSetCoarseMat()`, `PCDeflationGetCoarseKSP()`, `PCDeflationGetPC()`
834: M*/
836: PETSC_EXTERN PetscErrorCode PCCreate_Deflation(PC pc)
837: {
838: PC_Deflation *def;
840: PetscFunctionBegin;
841: PetscCall(PetscNew(&def));
842: pc->data = (void *)def;
844: def->init = PETSC_FALSE;
845: def->correct = PETSC_FALSE;
846: def->correctfact = 1.0;
847: def->reductionfact = -1;
848: def->spacetype = PC_DEFLATION_SPACE_HAAR;
849: def->spacesize = 1;
850: def->extendsp = PETSC_FALSE;
851: def->lvl = 0;
852: def->maxlvl = 0;
853: def->W = NULL;
854: def->Wt = NULL;
856: pc->ops->apply = PCApply_Deflation;
857: pc->ops->presolve = PCPreSolve_Deflation;
858: pc->ops->setup = PCSetUp_Deflation;
859: pc->ops->reset = PCReset_Deflation;
860: pc->ops->destroy = PCDestroy_Deflation;
861: pc->ops->setfromoptions = PCSetFromOptions_Deflation;
862: pc->ops->view = PCView_Deflation;
864: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetInitOnly_C", PCDeflationSetInitOnly_Deflation));
865: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetLevels_C", PCDeflationSetLevels_Deflation));
866: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetReductionFactor_C", PCDeflationSetReductionFactor_Deflation));
867: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetCorrectionFactor_C", PCDeflationSetCorrectionFactor_Deflation));
868: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetSpaceToCompute_C", PCDeflationSetSpaceToCompute_Deflation));
869: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetSpace_C", PCDeflationSetSpace_Deflation));
870: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetProjectionNullSpaceMat_C", PCDeflationSetProjectionNullSpaceMat_Deflation));
871: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationSetCoarseMat_C", PCDeflationSetCoarseMat_Deflation));
872: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationGetCoarseKSP_C", PCDeflationGetCoarseKSP_Deflation));
873: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCDeflationGetPC_C", PCDeflationGetPC_Deflation));
874: PetscFunctionReturn(PETSC_SUCCESS);
875: }