1 #include <petsc/private/pcmgimpl.h>
2 #include <petsc/private/pcbddcimpl.h>
3 #include <petsc/private/pcbddcprivateimpl.h>
4
PCMGGDSWSetUp(PC pc,PetscInt l,DM dm,KSP smooth,PetscInt Nc,Mat A,PetscInt * ns,Mat ** sA_IG_n,KSP ** sksp_n,IS ** sI_n,IS ** sG_n,Mat ** sGf_n,IS ** sGi_n,IS ** sGiM_n)5 static PetscErrorCode PCMGGDSWSetUp(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, Mat A, PetscInt *ns, Mat **sA_IG_n, KSP **sksp_n, IS **sI_n, IS **sG_n, Mat **sGf_n, IS **sGi_n, IS **sGiM_n)
6 {
7 KSP *sksp;
8 PC pcbddc = NULL, smoothpc;
9 PC_BDDC *ipcbddc;
10 PC_IS *ipcis;
11 Mat *sA_IG, *sGf, cmat, lA;
12 ISLocalToGlobalMapping l2g;
13 IS *sI, *sG, *sGi, *sGiM, cref;
14 PCBDDCSubSchurs sub_schurs = NULL;
15 PCBDDCGraph graph;
16 const char *prefix;
17 const PetscScalar *tdata;
18 PetscScalar *data, *cdata;
19 PetscReal tol = 0.0, otol;
20 const PetscInt *ia, *ja;
21 PetscInt *ccii, *cridx;
22 PetscInt i, j, ngct, ng, dbg = 0, odbg, minmax[2] = {0, PETSC_INT_MAX}, ominmax[2], vsize;
23 PetscBool flg, userdefined = PETSC_TRUE, reuse_solver = PETSC_TRUE, reduced = PETSC_FALSE;
24
25 PetscFunctionBegin;
26 PetscCall(MatGetBlockSize(A, &vsize));
27 PetscCall(KSPGetOptionsPrefix(smooth, &prefix));
28 PetscOptionsBegin(PetscObjectComm((PetscObject)smooth), prefix, "GDSW options", "PC");
29 PetscCall(PetscOptionsReal("-gdsw_tolerance", "Tolerance for eigenvalue problem", NULL, tol, &tol, NULL));
30 PetscCall(PetscOptionsBool("-gdsw_userdefined", "Use user-defined functions in addition to those adaptively generated", NULL, userdefined, &userdefined, NULL));
31 PetscCall(PetscOptionsIntArray("-gdsw_minmax", "Minimum and maximum number of basis functions per connected component for adaptive GDSW", NULL, minmax, (i = 2, &i), NULL));
32 PetscCall(PetscOptionsInt("-gdsw_vertex_size", "Connected components smaller or equal to vertex size will be considered as vertices", NULL, vsize, &vsize, NULL));
33 PetscCall(PetscOptionsBool("-gdsw_reuse", "Reuse interior solver from Schur complement computations", NULL, reuse_solver, &reuse_solver, NULL));
34 PetscCall(PetscOptionsBool("-gdsw_reduced", "Reduced GDSW", NULL, reduced, &reduced, NULL));
35 PetscCall(PetscOptionsInt("-gdsw_debug", "Debug output", NULL, dbg, &dbg, NULL));
36 PetscOptionsEnd();
37
38 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &flg));
39 if (!flg) {
40 MatNullSpace nnsp;
41
42 PetscCall(MatGetNearNullSpace(A, &nnsp));
43 PetscCall(PetscObjectReference((PetscObject)nnsp));
44 PetscCall(MatConvert(A, MATIS, MAT_INITIAL_MATRIX, &A));
45 PetscCall(MatSetNearNullSpace(A, nnsp));
46 PetscCall(MatNullSpaceDestroy(&nnsp));
47 } else PetscCall(PetscObjectReference((PetscObject)A));
48
49 /* TODO Multi sub */
50 *ns = 1;
51 PetscCall(PetscMalloc1(*ns, &sA_IG));
52 PetscCall(PetscMalloc1(*ns, &sksp));
53 PetscCall(PetscMalloc1(*ns, &sI));
54 PetscCall(PetscMalloc1(*ns, &sG));
55 PetscCall(PetscMalloc1(*ns, &sGf));
56 PetscCall(PetscMalloc1(*ns, &sGi));
57 PetscCall(PetscMalloc1(*ns, &sGiM));
58 *sA_IG_n = sA_IG;
59 *sksp_n = sksp;
60 *sI_n = sI;
61 *sG_n = sG;
62 *sGf_n = sGf;
63 *sGi_n = sGi;
64 *sGiM_n = sGiM;
65
66 /* submatrices and solvers */
67 PetscCall(KSPGetPC(smooth, &smoothpc));
68 PetscCall(PetscObjectTypeCompareAny((PetscObject)smoothpc, &flg, PCBDDC, ""));
69 if (!flg) {
70 Mat smoothA;
71
72 PetscCall(PCGetOperators(smoothpc, &smoothA, NULL));
73 PetscCall(PCCreate(PetscObjectComm((PetscObject)A), &pcbddc));
74 PetscCall(PCSetType(pcbddc, PCBDDC));
75 PetscCall(PCSetOperators(pcbddc, smoothA, A));
76 PetscCall(PCISSetUp(pcbddc, PETSC_TRUE, PETSC_FALSE));
77 } else {
78 PetscCall(PetscObjectReference((PetscObject)smoothpc));
79 pcbddc = smoothpc;
80 }
81 ipcis = (PC_IS *)pcbddc->data;
82 ipcbddc = (PC_BDDC *)pcbddc->data;
83 PetscCall(PetscObjectReference((PetscObject)ipcis->A_IB));
84 PetscCall(PetscObjectReference((PetscObject)ipcis->is_I_global));
85 PetscCall(PetscObjectReference((PetscObject)ipcis->is_B_global));
86 sA_IG[0] = ipcis->A_IB;
87 sI[0] = ipcis->is_I_global;
88 sG[0] = ipcis->is_B_global;
89
90 PetscCall(KSPCreate(PetscObjectComm((PetscObject)ipcis->A_II), &sksp[0]));
91 PetscCall(KSPSetNestLevel(sksp[0], pc->kspnestlevel));
92 PetscCall(KSPSetOperators(sksp[0], ipcis->A_II, ipcis->pA_II));
93 PetscCall(KSPSetOptionsPrefix(sksp[0], prefix));
94 PetscCall(KSPAppendOptionsPrefix(sksp[0], "gdsw_"));
95 PetscCall(KSPSetFromOptions(sksp[0]));
96
97 /* analyze interface */
98 PetscCall(MatISGetLocalMat(A, &lA));
99 graph = ipcbddc->mat_graph;
100 if (!flg) {
101 PetscInt N;
102
103 PetscCall(MatISGetLocalToGlobalMapping(A, &l2g, NULL));
104 PetscCall(MatGetSize(A, &N, NULL));
105 PetscCall(PCBDDCGraphInit(graph, l2g, N, PETSC_INT_MAX));
106 PetscCall(MatGetRowIJ(lA, 0, PETSC_TRUE, PETSC_FALSE, &graph->nvtxs_csr, (const PetscInt **)&graph->xadj, (const PetscInt **)&graph->adjncy, &flg));
107 PetscCall(PCBDDCGraphSetUp(graph, vsize, NULL, NULL, 0, NULL, NULL));
108 PetscCall(MatRestoreRowIJ(lA, 0, PETSC_TRUE, PETSC_FALSE, &graph->nvtxs_csr, (const PetscInt **)&graph->xadj, (const PetscInt **)&graph->adjncy, &flg));
109 PetscCall(PCBDDCGraphComputeConnectedComponents(graph));
110 }
111 l2g = graph->l2gmap;
112 if (reduced) {
113 PetscContainer gcand;
114 PCBDDCGraphCandidates cand;
115 PetscErrorCode (*rgdsw)(DM, PetscInt *, IS **);
116
117 PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMComputeLocalRGDSWSets", &rgdsw));
118 PetscCheck(rgdsw, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not supported");
119 PetscCall(PetscNew(&cand));
120 PetscCall((*rgdsw)(dm, &cand->nfc, &cand->Faces));
121 /* filter interior (if any) and guarantee IS are ordered by global numbering */
122 for (i = 0; i < cand->nfc; i++) {
123 IS is, is2;
124
125 PetscCall(ISLocalToGlobalMappingApplyIS(l2g, cand->Faces[i], &is));
126 PetscCall(ISDestroy(&cand->Faces[i]));
127 PetscCall(ISSort(is));
128 PetscCall(ISGlobalToLocalMappingApplyIS(l2g, IS_GTOLM_DROP, is, &is2));
129 PetscCall(ISDestroy(&is));
130 PetscCall(ISGlobalToLocalMappingApplyIS(ipcis->BtoNmap, IS_GTOLM_DROP, is2, &is));
131 PetscCall(ISDestroy(&is2));
132 PetscCall(ISLocalToGlobalMappingApplyIS(ipcis->BtoNmap, is, &cand->Faces[i]));
133 PetscCall(ISDestroy(&is));
134 }
135 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &gcand));
136 PetscCall(PetscContainerSetPointer(gcand, cand));
137 PetscCall(PetscContainerSetCtxDestroy(gcand, PCBDDCDestroyGraphCandidatesIS));
138 PetscCall(PetscObjectCompose((PetscObject)l2g, "_PCBDDCGraphCandidatesIS", (PetscObject)gcand));
139 PetscCall(PetscContainerDestroy(&gcand));
140 }
141
142 /* interface functions */
143 otol = ipcbddc->adaptive_threshold[1];
144 odbg = ipcbddc->dbg_flag;
145 ominmax[0] = ipcbddc->adaptive_nmin;
146 ominmax[1] = ipcbddc->adaptive_nmax;
147 ipcbddc->adaptive_threshold[1] = tol;
148 ipcbddc->dbg_flag = dbg;
149 ipcbddc->adaptive_nmin = minmax[0];
150 ipcbddc->adaptive_nmax = minmax[1];
151 if (tol != 0.0) { /* adaptive */
152 Mat lS;
153
154 PetscCall(MatCreateSchurComplement(ipcis->A_II, ipcis->pA_II, ipcis->A_IB, ipcis->A_BI, ipcis->A_BB, &lS));
155 PetscCall(KSPGetOptionsPrefix(sksp[0], &prefix));
156 PetscCall(PCBDDCSubSchursCreate(&sub_schurs));
157 PetscCall(PCBDDCSubSchursInit(sub_schurs, prefix, ipcis->is_I_local, ipcis->is_B_local, graph, ipcis->BtoNmap, PETSC_FALSE, PETSC_TRUE));
158 if (userdefined) PetscCall(PCBDDCComputeFakeChange(pcbddc, PETSC_FALSE, graph, NULL, &cmat, &cref, NULL, &flg));
159 else {
160 cmat = NULL;
161 cref = NULL;
162 }
163 PetscCall(PCBDDCSubSchursSetUp(sub_schurs, lA, lS, PETSC_TRUE, NULL, NULL, -1, NULL, PETSC_TRUE, reuse_solver, PETSC_FALSE, 0, NULL, NULL, cmat, cref));
164 PetscCall(MatDestroy(&lS));
165 PetscCall(MatDestroy(&cmat));
166 PetscCall(ISDestroy(&cref));
167 if (sub_schurs->reuse_solver) {
168 PetscCall(KSPSetPC(sksp[0], sub_schurs->reuse_solver->interior_solver));
169 PetscCall(PCDestroy(&sub_schurs->reuse_solver->interior_solver));
170 sub_schurs->reuse_solver = NULL;
171 }
172 }
173 PetscCall(PCBDDCComputeFakeChange(pcbddc, PETSC_TRUE, graph, sub_schurs, &cmat, &cref, &sGiM[0], NULL));
174 PetscCall(PCBDDCSubSchursDestroy(&sub_schurs));
175 ipcbddc->adaptive_threshold[1] = otol;
176 ipcbddc->dbg_flag = odbg;
177 ipcbddc->adaptive_nmin = ominmax[0];
178 ipcbddc->adaptive_nmax = ominmax[1];
179
180 PetscCall(ISLocalToGlobalMappingApplyIS(l2g, cref, &sGi[0]));
181 PetscCall(ISDestroy(&cref));
182
183 PetscCall(MatSeqAIJGetArrayRead(cmat, &tdata));
184 PetscCall(MatGetRowIJ(cmat, 0, PETSC_FALSE, PETSC_FALSE, &ngct, &ia, &ja, &flg));
185 PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in MatGetRowIJ");
186
187 PetscCall(PetscMalloc1(ngct + 1, &ccii));
188 PetscCall(PetscMalloc1(ia[ngct], &cridx));
189 PetscCall(PetscMalloc1(ia[ngct], &cdata));
190
191 PetscCall(PetscArraycpy(ccii, ia, ngct + 1));
192 PetscCall(PetscArraycpy(cdata, tdata, ia[ngct]));
193 PetscCall(ISGlobalToLocalMappingApply(ipcis->BtoNmap, IS_GTOLM_DROP, ia[ngct], ja, &i, cridx));
194 PetscCheck(i == ia[ngct], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in G2L");
195
196 PetscCall(MatRestoreRowIJ(cmat, 0, PETSC_FALSE, PETSC_FALSE, &i, &ia, &ja, &flg));
197 PetscCall(MatSeqAIJRestoreArrayRead(cmat, &tdata));
198 PetscCall(MatDestroy(&cmat));
199
200 /* populate dense matrix */
201 PetscCall(ISGetLocalSize(sG[0], &ng));
202 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, ng, ngct, NULL, &sGf[0]));
203 PetscCall(MatDenseGetArrayWrite(sGf[0], &data));
204 for (i = 0; i < ngct; i++)
205 for (j = ccii[i]; j < ccii[i + 1]; j++) data[ng * i + cridx[j]] = cdata[j];
206 PetscCall(MatDenseRestoreArrayWrite(sGf[0], &data));
207
208 PetscCall(PetscFree(cdata));
209 PetscCall(PetscFree(ccii));
210 PetscCall(PetscFree(cridx));
211 PetscCall(PCDestroy(&pcbddc));
212 PetscCall(MatDestroy(&A));
213 PetscFunctionReturn(PETSC_SUCCESS);
214 }
215
PCMGGDSWCreateCoarseSpace_Private(PC pc,PetscInt l,DM dm,KSP smooth,PetscInt Nc,Mat guess,Mat * cspace)216 PetscErrorCode PCMGGDSWCreateCoarseSpace_Private(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, Mat guess, Mat *cspace)
217 {
218 KSP *sksp;
219 Mat A, *sA_IG, *sGf, preallocator;
220 IS Gidx, GidxMult, cG;
221 IS *sI, *sG, *sGi, *sGiM;
222 const PetscInt *cidx;
223 PetscInt NG, ns, n, i, c, rbs, cbs[2];
224 PetscBool flg;
225 MatType ptype;
226
227 PetscFunctionBegin;
228 *cspace = NULL;
229 if (!l) PetscFunctionReturn(PETSC_SUCCESS);
230 if (pc->useAmat) {
231 PetscCall(KSPGetOperatorsSet(smooth, &flg, NULL));
232 PetscCheck(flg, PetscObjectComm((PetscObject)smooth), PETSC_ERR_ORDER, "Amat not set");
233 PetscCall(KSPGetOperators(smooth, &A, NULL));
234 } else {
235 PetscCall(KSPGetOperatorsSet(smooth, NULL, &flg));
236 PetscCheck(flg, PetscObjectComm((PetscObject)smooth), PETSC_ERR_ORDER, "Pmat not set");
237 PetscCall(KSPGetOperators(smooth, NULL, &A));
238 }
239
240 /* Setup (also setup smoother here) */
241 if (!pc->setupcalled) PetscCall(KSPSetFromOptions(smooth));
242 PetscCall(KSPSetUp(smooth));
243 PetscCall(KSPSetUpOnBlocks(smooth));
244 PetscCall(PCMGGDSWSetUp(pc, l, dm, smooth, Nc, A, &ns, &sA_IG, &sksp, &sI, &sG, &sGf, &sGi, &sGiM));
245
246 /* Number GDSW basis functions */
247 PetscCall(ISConcatenate(PetscObjectComm((PetscObject)A), ns, sGi, &Gidx));
248 PetscCall(ISConcatenate(PetscObjectComm((PetscObject)A), ns, sGiM, &GidxMult));
249 PetscCall(ISRenumber(Gidx, GidxMult, &NG, &cG));
250 PetscCall(ISDestroy(&Gidx));
251
252 /* Detect column block size */
253 PetscCall(ISGetMinMax(GidxMult, &cbs[0], &cbs[1]));
254 PetscCall(PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)A), cbs, cbs));
255 PetscCall(ISDestroy(&GidxMult));
256
257 /* Construct global interpolation matrix */
258 PetscCall(MatGetLocalSize(A, NULL, &n));
259 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &preallocator));
260 PetscCall(MatSetSizes(preallocator, n, PETSC_DECIDE, PETSC_DECIDE, NG));
261 PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
262 PetscCall(MatSetUp(preallocator));
263 PetscCall(ISGetIndices(cG, &cidx));
264 for (i = 0, c = 0; i < ns; i++) {
265 const PetscInt *ri, *rg;
266 PetscInt nri, nrg, ncg;
267
268 PetscCall(ISGetLocalSize(sI[i], &nri));
269 PetscCall(ISGetLocalSize(sG[i], &nrg));
270 PetscCall(ISGetIndices(sI[i], &ri));
271 PetscCall(ISGetIndices(sG[i], &rg));
272 PetscCall(MatGetSize(sGf[i], NULL, &ncg));
273 PetscCall(MatSetValues(preallocator, nri, ri, ncg, cidx + c, NULL, INSERT_VALUES));
274 PetscCall(MatSetValues(preallocator, nrg, rg, ncg, cidx + c, NULL, INSERT_VALUES));
275 PetscCall(ISRestoreIndices(sI[i], &ri));
276 PetscCall(ISRestoreIndices(sG[i], &rg));
277 }
278 PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
279 PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
280
281 ptype = MATAIJ;
282 if (PetscDefined(HAVE_DEVICE)) {
283 PetscCall(MatBoundToCPU(A, &flg));
284 if (!flg) {
285 VecType vtype;
286 char *found = NULL;
287
288 PetscCall(MatGetVecType(A, &vtype));
289 PetscCall(PetscStrstr(vtype, "cuda", &found));
290 if (found) ptype = MATAIJCUSPARSE;
291 }
292 }
293 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), cspace));
294 PetscCall(MatSetSizes(*cspace, n, PETSC_DECIDE, PETSC_DECIDE, NG));
295 PetscCall(MatSetType(*cspace, ptype));
296 PetscCall(MatGetBlockSizes(A, NULL, &rbs));
297 PetscCall(MatSetBlockSizes(*cspace, rbs, cbs[0] == cbs[1] ? cbs[0] : 1));
298 PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_FALSE, *cspace));
299 PetscCall(MatDestroy(&preallocator));
300 PetscCall(MatSetOption(*cspace, MAT_ROW_ORIENTED, PETSC_FALSE));
301
302 for (i = 0, c = 0; i < ns; i++) {
303 Mat X, Y;
304 const PetscScalar *v;
305 const PetscInt *ri, *rg;
306 PetscInt nri, nrg, ncg;
307
308 PetscCall(MatMatMult(sA_IG[i], sGf[i], MAT_INITIAL_MATRIX, PETSC_CURRENT, &Y));
309 PetscCall(MatScale(Y, -1.0));
310 PetscCall(MatDuplicate(Y, MAT_DO_NOT_COPY_VALUES, &X));
311 PetscCall(KSPMatSolve(sksp[i], Y, X));
312
313 PetscCall(ISGetLocalSize(sI[i], &nri));
314 PetscCall(ISGetLocalSize(sG[i], &nrg));
315 PetscCall(ISGetIndices(sI[i], &ri));
316 PetscCall(ISGetIndices(sG[i], &rg));
317 PetscCall(MatGetSize(sGf[i], NULL, &ncg));
318
319 PetscCall(MatDenseGetArrayRead(X, &v));
320 PetscCall(MatSetValues(*cspace, nri, ri, ncg, cidx + c, v, INSERT_VALUES));
321 PetscCall(MatDenseRestoreArrayRead(X, &v));
322 PetscCall(MatDenseGetArrayRead(sGf[i], &v));
323 PetscCall(MatSetValues(*cspace, nrg, rg, ncg, cidx + c, v, INSERT_VALUES));
324 PetscCall(MatDenseRestoreArrayRead(sGf[i], &v));
325 PetscCall(ISRestoreIndices(sI[i], &ri));
326 PetscCall(ISRestoreIndices(sG[i], &rg));
327 PetscCall(MatDestroy(&Y));
328 PetscCall(MatDestroy(&X));
329 }
330 PetscCall(ISRestoreIndices(cG, &cidx));
331 PetscCall(ISDestroy(&cG));
332 PetscCall(MatAssemblyBegin(*cspace, MAT_FINAL_ASSEMBLY));
333 PetscCall(MatAssemblyEnd(*cspace, MAT_FINAL_ASSEMBLY));
334
335 for (i = 0; i < ns; i++) {
336 PetscCall(KSPDestroy(&sksp[i]));
337 PetscCall(ISDestroy(&sI[i]));
338 PetscCall(ISDestroy(&sG[i]));
339 PetscCall(ISDestroy(&sGi[i]));
340 PetscCall(ISDestroy(&sGiM[i]));
341 PetscCall(MatDestroy(&sGf[i]));
342 PetscCall(MatDestroy(&sA_IG[i]));
343 }
344 PetscCall(PetscFree(sksp));
345 PetscCall(PetscFree(sI));
346 PetscCall(PetscFree(sG));
347 PetscCall(PetscFree(sGi));
348 PetscCall(PetscFree(sGiM));
349 PetscCall(PetscFree(sGf));
350 PetscCall(PetscFree(sA_IG));
351 PetscFunctionReturn(PETSC_SUCCESS);
352 }
353