1 #include <../src/ksp/pc/impls/deflation/deflation.h> /*I "petscksp.h" I*/
2
3 PetscScalar db2[] = {0.7071067811865476, 0.7071067811865476};
4
5 PetscScalar db4[] = {-0.12940952255092145, 0.22414386804185735, 0.836516303737469, 0.48296291314469025};
6
7 PetscScalar db8[] = {-0.010597401784997278, 0.032883011666982945, 0.030841381835986965, -0.18703481171888114, -0.02798376941698385, 0.6308807679295904, 0.7148465705525415, 0.23037781330885523};
8
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};
11
12 PetscScalar biorth22[] = {0.0, -0.1767766952966369, 0.3535533905932738, 1.0606601717798214, 0.3535533905932738, -0.1767766952966369};
13
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};
15
PCDeflationCreateSpaceWave(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt ncoeffs,PetscScalar * coeffs,PetscBool trunc,Mat * H)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;
20
21 PetscFunctionBegin;
22 PetscCall(PetscMalloc1(ncoeffs, &Iidx));
23
24 PetscCall(MatCreate(comm, &defl));
25 PetscCall(MatSetSizes(defl, m, n, M, N));
26 PetscCall(MatSetUp(defl));
27 PetscCall(MatSeqAIJSetPreallocation(defl, ncoeffs, NULL));
28 PetscCall(MatMPIAIJSetPreallocation(defl, ncoeffs, NULL, ncoeffs, NULL));
29 PetscCall(MatSetOption(defl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
30 PetscCall(MatSetOption(defl, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
31
32 /* Alg 735 Taswell: fvecmat */
33 k = ncoeffs - 2;
34 if (trunc) k = k / 2;
35
36 PetscCall(MatGetOwnershipRange(defl, &ilo, &ihi));
37 for (i = 0; i < ncoeffs; i++) {
38 Iidx[i] = i + ilo * 2 - k;
39 if (Iidx[i] >= N) Iidx[i] = PETSC_INT_MIN;
40 }
41 for (i = ilo; i < ihi; i++) {
42 PetscCall(MatSetValues(defl, 1, &i, ncoeffs, Iidx, coeffs, INSERT_VALUES));
43 for (j = 0; j < ncoeffs; j++) {
44 Iidx[j] += 2;
45 if (Iidx[j] >= N) Iidx[j] = PETSC_INT_MIN;
46 }
47 }
48
49 PetscCall(MatAssemblyBegin(defl, MAT_FINAL_ASSEMBLY));
50 PetscCall(MatAssemblyEnd(defl, MAT_FINAL_ASSEMBLY));
51
52 PetscCall(PetscFree(Iidx));
53 *H = defl;
54 PetscFunctionReturn(PETSC_SUCCESS);
55 }
56
PCDeflationGetSpaceHaar(PC pc,Mat * W,PetscInt size)57 static PetscErrorCode PCDeflationGetSpaceHaar(PC pc, Mat *W, PetscInt size)
58 {
59 Mat A, defl;
60 PetscInt i, j, len, ilo, ihi, *Iidx, m, M;
61 PetscScalar *col, val;
62
63 PetscFunctionBegin;
64 /* Haar basis wavelet, level=size */
65 len = pow(2, size);
66 PetscCall(PetscMalloc2(len, &col, len, &Iidx));
67 val = 1. / pow(2, size / 2.);
68 for (i = 0; i < len; i++) col[i] = val;
69
70 PetscCall(PCGetOperators(pc, NULL, &A));
71 PetscCall(MatGetLocalSize(A, &m, NULL));
72 PetscCall(MatGetSize(A, &M, NULL));
73 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &defl));
74 PetscCall(MatSetSizes(defl, m, PETSC_DECIDE, M, PetscCeilInt(M, len)));
75 PetscCall(MatSetUp(defl));
76 PetscCall(MatSeqAIJSetPreallocation(defl, size, NULL));
77 PetscCall(MatMPIAIJSetPreallocation(defl, size, NULL, size, NULL));
78 PetscCall(MatSetOption(defl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
79
80 PetscCall(MatGetOwnershipRangeColumn(defl, &ilo, &ihi));
81 for (i = 0; i < len; i++) Iidx[i] = i + ilo * len;
82 if (M % len && ihi == PetscCeilInt(M, len)) ihi -= 1;
83 for (i = ilo; i < ihi; i++) {
84 PetscCall(MatSetValues(defl, len, Iidx, 1, &i, col, INSERT_VALUES));
85 for (j = 0; j < len; j++) Iidx[j] += len;
86 }
87 if (M % len && ihi + 1 == PetscCeilInt(M, len)) {
88 len = M % len;
89 val = 1. / pow(pow(2, len), 0.5);
90 for (i = 0; i < len; i++) col[i] = val;
91 PetscCall(MatSetValues(defl, len, Iidx, 1, &ihi, col, INSERT_VALUES));
92 }
93
94 PetscCall(MatAssemblyBegin(defl, MAT_FINAL_ASSEMBLY));
95 PetscCall(MatAssemblyEnd(defl, MAT_FINAL_ASSEMBLY));
96
97 PetscCall(PetscFree2(col, Iidx));
98 *W = defl;
99 PetscFunctionReturn(PETSC_SUCCESS);
100 }
101
PCDeflationGetSpaceWave(PC pc,Mat * W,PetscInt size,PetscInt ncoeffs,PetscScalar * coeffs,PetscBool trunc)102 static PetscErrorCode PCDeflationGetSpaceWave(PC pc, Mat *W, PetscInt size, PetscInt ncoeffs, PetscScalar *coeffs, PetscBool trunc)
103 {
104 Mat A, *H, defl;
105 PetscInt i, m, M, Mdefl, Ndefl;
106 MPI_Comm comm;
107
108 PetscFunctionBegin;
109 PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
110 PetscCall(PetscMalloc1(size, &H));
111 PetscCall(PCGetOperators(pc, &A, NULL));
112 PetscCall(MatGetLocalSize(A, &m, NULL));
113 PetscCall(MatGetSize(A, &M, NULL));
114 Mdefl = M;
115 Ndefl = M;
116 for (i = 0; i < size; i++) {
117 if (Mdefl % 2) {
118 if (trunc) Mdefl = (PetscInt)PetscCeilReal(Mdefl / 2.);
119 else Mdefl = (PetscInt)PetscFloorReal((ncoeffs + Mdefl - 1) / 2.);
120 } else Mdefl = Mdefl / 2;
121 PetscCall(PCDeflationCreateSpaceWave(comm, PETSC_DECIDE, m, Mdefl, Ndefl, ncoeffs, coeffs, trunc, &H[i]));
122 PetscCall(MatGetLocalSize(H[i], &m, NULL));
123 Ndefl = Mdefl;
124 }
125 PetscCall(MatCreateComposite(comm, size, H, &defl));
126 PetscCall(MatCompositeSetType(defl, MAT_COMPOSITE_MULTIPLICATIVE));
127 *W = defl;
128
129 for (i = 0; i < size; i++) PetscCall(MatDestroy(&H[i]));
130 PetscCall(PetscFree(H));
131 PetscFunctionReturn(PETSC_SUCCESS);
132 }
133
PCDeflationGetSpaceAggregation(PC pc,Mat * W)134 static PetscErrorCode PCDeflationGetSpaceAggregation(PC pc, Mat *W)
135 {
136 Mat A, defl;
137 PetscInt i, ilo, ihi, *Iidx, M;
138 PetscMPIInt m;
139 PetscScalar *col;
140 MPI_Comm comm;
141
142 PetscFunctionBegin;
143 PetscCall(PCGetOperators(pc, &A, NULL));
144 PetscCall(MatGetOwnershipRangeColumn(A, &ilo, &ihi));
145 PetscCall(MatGetSize(A, &M, NULL));
146 PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
147 PetscCallMPI(MPI_Comm_size(comm, &m));
148 PetscCall(MatCreate(comm, &defl));
149 PetscCall(MatSetSizes(defl, ihi - ilo, 1, M, m));
150 PetscCall(MatSetUp(defl));
151 PetscCall(MatSeqAIJSetPreallocation(defl, 1, NULL));
152 PetscCall(MatMPIAIJSetPreallocation(defl, 1, NULL, 0, NULL));
153 PetscCall(MatSetOption(defl, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
154 PetscCall(MatSetOption(defl, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
155
156 PetscCall(PetscMalloc2(ihi - ilo, &col, ihi - ilo, &Iidx));
157 for (i = ilo; i < ihi; i++) {
158 Iidx[i - ilo] = i;
159 col[i - ilo] = 1;
160 }
161 PetscCallMPI(MPI_Comm_rank(comm, &m));
162 i = m;
163 PetscCall(MatSetValues(defl, ihi - ilo, Iidx, 1, &i, col, INSERT_VALUES));
164
165 PetscCall(MatAssemblyBegin(defl, MAT_FINAL_ASSEMBLY));
166 PetscCall(MatAssemblyEnd(defl, MAT_FINAL_ASSEMBLY));
167
168 PetscCall(PetscFree2(col, Iidx));
169 *W = defl;
170 PetscFunctionReturn(PETSC_SUCCESS);
171 }
172
PCDeflationComputeSpace(PC pc)173 PetscErrorCode PCDeflationComputeSpace(PC pc)
174 {
175 Mat defl;
176 PetscBool transp = PETSC_TRUE;
177 PC_Deflation *def = (PC_Deflation *)pc->data;
178
179 PetscFunctionBegin;
180 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
181 PetscCheck(def->spacesize >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Wrong PCDeflation space size specified: %" PetscInt_FMT, def->spacesize);
182 switch (def->spacetype) {
183 case PC_DEFLATION_SPACE_HAAR:
184 transp = PETSC_FALSE;
185 PetscCall(PCDeflationGetSpaceHaar(pc, &defl, def->spacesize));
186 break;
187 case PC_DEFLATION_SPACE_DB2:
188 PetscCall(PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 2, db2, PetscNot(def->extendsp)));
189 break;
190 case PC_DEFLATION_SPACE_DB4:
191 PetscCall(PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 4, db4, PetscNot(def->extendsp)));
192 break;
193 case PC_DEFLATION_SPACE_DB8:
194 PetscCall(PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 8, db8, PetscNot(def->extendsp)));
195 break;
196 case PC_DEFLATION_SPACE_DB16:
197 PetscCall(PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 16, db16, PetscNot(def->extendsp)));
198 break;
199 case PC_DEFLATION_SPACE_BIORTH22:
200 PetscCall(PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 6, biorth22, PetscNot(def->extendsp)));
201 break;
202 case PC_DEFLATION_SPACE_MEYER:
203 PetscCall(PCDeflationGetSpaceWave(pc, &defl, def->spacesize, 62, meyer, PetscNot(def->extendsp)));
204 break;
205 case PC_DEFLATION_SPACE_AGGREGATION:
206 transp = PETSC_FALSE;
207 PetscCall(PCDeflationGetSpaceAggregation(pc, &defl));
208 break;
209 default:
210 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Wrong PCDeflationSpaceType specified");
211 }
212
213 PetscCall(PCDeflationSetSpace(pc, defl, transp));
214 PetscCall(MatDestroy(&defl));
215 PetscFunctionReturn(PETSC_SUCCESS);
216 }
217