Actual source code: deflationspace.c
1: #include <../src/ksp/pc/impls/deflation/deflation.h>
3: PetscScalar db2[] = {0.7071067811865476, 0.7071067811865476};
5: PetscScalar db4[] = {-0.12940952255092145, 0.22414386804185735, 0.836516303737469, 0.48296291314469025};
7: PetscScalar db8[] = {-0.010597401784997278, 0.032883011666982945, 0.030841381835986965, -0.18703481171888114, -0.02798376941698385, 0.6308807679295904, 0.7148465705525415, 0.23037781330885523};
9: PetscScalar db16[] = {-0.00011747678400228192, 0.0006754494059985568, -0.0003917403729959771, -0.00487035299301066, 0.008746094047015655, 0.013981027917015516, -0.04408825393106472, -0.01736930100202211,
10: 0.128747426620186, 0.00047248457399797254, -0.2840155429624281, -0.015829105256023893, 0.5853546836548691, 0.6756307362980128, 0.3128715909144659, 0.05441584224308161};
12: PetscScalar biorth22[] = {0.0, -0.1767766952966369, 0.3535533905932738, 1.0606601717798214, 0.3535533905932738, -0.1767766952966369};
14: PetscScalar meyer[] = {0.0, -1.009999956941423e-12, 8.519459636796214e-09, -1.111944952595278e-08, -1.0798819539621958e-08, 6.066975741351135e-08, -1.0866516536735883e-07, 8.200680650386481e-08, 1.1783004497663934e-07, -5.506340565252278e-07, 1.1307947017916706e-06, -1.489549216497156e-06, 7.367572885903746e-07, 3.20544191334478e-06, -1.6312699734552807e-05, 6.554305930575149e-05, -0.0006011502343516092, -0.002704672124643725, 0.002202534100911002, 0.006045814097323304, -0.006387718318497156, -0.011061496392513451, 0.015270015130934803, 0.017423434103729693, -0.03213079399021176, -0.024348745906078023, 0.0637390243228016, 0.030655091960824263, -0.13284520043622938, -0.035087555656258346, 0.44459300275757724, 0.7445855923188063, 0.44459300275757724, -0.035087555656258346, -0.13284520043622938, 0.030655091960824263, 0.0637390243228016, -0.024348745906078023, -0.03213079399021176, 0.017423434103729693, 0.015270015130934803, -0.011061496392513451, -0.006387718318497156, 0.006045814097323304, 0.002202534100911002, -0.002704672124643725, -0.0006011502343516092, 6.554305930575149e-05, -1.6312699734552807e-05, 3.20544191334478e-06, 7.367572885903746e-07, -1.489549216497156e-06, 1.1307947017916706e-06, -5.506340565252278e-07, 1.1783004497663934e-07, 8.200680650386481e-08, -1.0866516536735883e-07, 6.066975741351135e-08, -1.0798819539621958e-08, -1.111944952595278e-08, 8.519459636796214e-09, -1.009999956941423e-12};
16: static PetscErrorCode PCDeflationCreateSpaceWave(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt ncoeffs, PetscScalar *coeffs, PetscBool trunc, Mat *H)
17: {
18: Mat defl;
19: PetscInt i, j, k, ilo, ihi, *Iidx;
21: PetscMalloc1(ncoeffs, &Iidx);
23: MatCreate(comm, &defl);
24: MatSetSizes(defl, m, n, M, N);
25: MatSetUp(defl);
26: MatSeqAIJSetPreallocation(defl, ncoeffs, NULL);
27: MatMPIAIJSetPreallocation(defl, ncoeffs, NULL, ncoeffs, NULL);
28: MatSetOption(defl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
29: MatSetOption(defl, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
31: /* Alg 735 Taswell: fvecmat */
32: k = ncoeffs - 2;
33: if (trunc) k = k / 2;
35: MatGetOwnershipRange(defl, &ilo, &ihi);
36: for (i = 0; i < ncoeffs; i++) {
37: Iidx[i] = i + ilo * 2 - k;
38: if (Iidx[i] >= N) Iidx[i] = PETSC_MIN_INT;
39: }
40: for (i = ilo; i < ihi; i++) {
41: MatSetValues(defl, 1, &i, ncoeffs, Iidx, coeffs, INSERT_VALUES);
42: for (j = 0; j < ncoeffs; j++) {
43: Iidx[j] += 2;
44: if (Iidx[j] >= N) Iidx[j] = PETSC_MIN_INT;
45: }
46: }
48: MatAssemblyBegin(defl, MAT_FINAL_ASSEMBLY);
49: MatAssemblyEnd(defl, MAT_FINAL_ASSEMBLY);
51: PetscFree(Iidx);
52: *H = defl;
53: return 0;
54: }
56: PetscErrorCode PCDeflationGetSpaceHaar(PC pc, Mat *W, PetscInt size)
57: {
58: Mat A, defl;
59: PetscInt i, j, len, ilo, ihi, *Iidx, m, M;
60: PetscScalar *col, val;
62: /* Haar basis wavelet, level=size */
63: len = pow(2, size);
64: PetscMalloc2(len, &col, len, &Iidx);
65: val = 1. / pow(2, size / 2.);
66: for (i = 0; i < len; i++) col[i] = val;
68: PCGetOperators(pc, NULL, &A);
69: MatGetLocalSize(A, &m, NULL);
70: MatGetSize(A, &M, NULL);
71: MatCreate(PetscObjectComm((PetscObject)A), &defl);
72: MatSetSizes(defl, m, PETSC_DECIDE, M, PetscCeilInt(M, len));
73: MatSetUp(defl);
74: MatSeqAIJSetPreallocation(defl, size, NULL);
75: MatMPIAIJSetPreallocation(defl, size, NULL, size, NULL);
76: MatSetOption(defl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
78: MatGetOwnershipRangeColumn(defl, &ilo, &ihi);
79: for (i = 0; i < len; i++) Iidx[i] = i + ilo * len;
80: if (M % len && ihi == PetscCeilInt(M, len)) ihi -= 1;
81: for (i = ilo; i < ihi; i++) {
82: MatSetValues(defl, len, Iidx, 1, &i, col, INSERT_VALUES);
83: for (j = 0; j < len; j++) Iidx[j] += len;
84: }
85: if (M % len && ihi + 1 == PetscCeilInt(M, len)) {
86: len = M % len;
87: val = 1. / pow(pow(2, len), 0.5);
88: for (i = 0; i < len; i++) col[i] = val;
89: MatSetValues(defl, len, Iidx, 1, &ihi, col, INSERT_VALUES);
90: }
92: MatAssemblyBegin(defl, MAT_FINAL_ASSEMBLY);
93: MatAssemblyEnd(defl, MAT_FINAL_ASSEMBLY);
95: PetscFree2(col, Iidx);
96: *W = defl;
97: return 0;
98: }
100: PetscErrorCode PCDeflationGetSpaceWave(PC pc, Mat *W, PetscInt size, PetscInt ncoeffs, PetscScalar *coeffs, PetscBool trunc)
101: {
102: Mat A, *H, defl;
103: PetscInt i, m, M, Mdefl, Ndefl;
104: MPI_Comm comm;
106: PetscObjectGetComm((PetscObject)pc, &comm);
107: PetscMalloc1(size, &H);
108: PCGetOperators(pc, &A, NULL);
109: MatGetLocalSize(A, &m, NULL);
110: MatGetSize(A, &M, NULL);
111: Mdefl = M;
112: Ndefl = M;
113: for (i = 0; i < size; i++) {
114: if (Mdefl % 2) {
115: if (trunc) Mdefl = (PetscInt)PetscCeilReal(Mdefl / 2.);
116: else Mdefl = (PetscInt)PetscFloorReal((ncoeffs + Mdefl - 1) / 2.);
117: } else Mdefl = Mdefl / 2;
118: PCDeflationCreateSpaceWave(comm, PETSC_DECIDE, m, Mdefl, Ndefl, ncoeffs, coeffs, trunc, &H[i]);
119: MatGetLocalSize(H[i], &m, NULL);
120: Ndefl = Mdefl;
121: }
122: MatCreateComposite(comm, size, H, &defl);
123: MatCompositeSetType(defl, MAT_COMPOSITE_MULTIPLICATIVE);
124: *W = defl;
126: for (i = 0; i < size; i++) MatDestroy(&H[i]);
127: PetscFree(H);
128: return 0;
129: }
131: PetscErrorCode PCDeflationGetSpaceAggregation(PC pc, Mat *W)
132: {
133: Mat A, defl;
134: PetscInt i, ilo, ihi, *Iidx, M;
135: PetscMPIInt m;
136: PetscScalar *col;
137: MPI_Comm comm;
139: PCGetOperators(pc, &A, NULL);
140: MatGetOwnershipRangeColumn(A, &ilo, &ihi);
141: MatGetSize(A, &M, NULL);
142: PetscObjectGetComm((PetscObject)A, &comm);
143: MPI_Comm_size(comm, &m);
144: MatCreate(comm, &defl);
145: MatSetSizes(defl, ihi - ilo, 1, M, m);
146: MatSetUp(defl);
147: MatSeqAIJSetPreallocation(defl, 1, NULL);
148: MatMPIAIJSetPreallocation(defl, 1, NULL, 0, NULL);
149: MatSetOption(defl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
150: MatSetOption(defl, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
152: PetscMalloc2(ihi - ilo, &col, ihi - ilo, &Iidx);
153: for (i = ilo; i < ihi; i++) {
154: Iidx[i - ilo] = i;
155: col[i - ilo] = 1;
156: }
157: MPI_Comm_rank(comm, &m);
158: i = m;
159: MatSetValues(defl, ihi - ilo, Iidx, 1, &i, col, INSERT_VALUES);
161: MatAssemblyBegin(defl, MAT_FINAL_ASSEMBLY);
162: MatAssemblyEnd(defl, MAT_FINAL_ASSEMBLY);
164: PetscFree2(col, Iidx);
165: *W = defl;
166: return 0;
167: }
169: PetscErrorCode PCDeflationComputeSpace(PC pc)
170: {
171: Mat defl;
172: PetscBool transp = PETSC_TRUE;
173: PC_Deflation *def = (PC_Deflation *)pc->data;
177: switch (def->spacetype) {
178: case PC_DEFLATION_SPACE_HAAR:
179: transp = PETSC_FALSE;
180: PCDeflationGetSpaceHaar(pc, &defl, def->spacesize);
181: break;
182: case PC_DEFLATION_SPACE_DB2:
183: PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 2, db2, PetscNot(def->extendsp));
184: break;
185: case PC_DEFLATION_SPACE_DB4:
186: PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 4, db4, PetscNot(def->extendsp));
187: break;
188: case PC_DEFLATION_SPACE_DB8:
189: PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 8, db8, PetscNot(def->extendsp));
190: break;
191: case PC_DEFLATION_SPACE_DB16:
192: PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 16, db16, PetscNot(def->extendsp));
193: break;
194: case PC_DEFLATION_SPACE_BIORTH22:
195: PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 6, biorth22, PetscNot(def->extendsp));
196: break;
197: case PC_DEFLATION_SPACE_MEYER:
198: PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 62, meyer, PetscNot(def->extendsp));
199: break;
200: case PC_DEFLATION_SPACE_AGGREGATION:
201: transp = PETSC_FALSE;
202: PCDeflationGetSpaceAggregation(pc, &defl);
203: break;
204: default:
205: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Wrong PCDeflationSpaceType specified");
206: }
208: PCDeflationSetSpace(pc, defl, transp);
209: MatDestroy(&defl);
210: return 0;
211: }