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