Actual source code: gdsw.c
1: #include <petsc/private/pcmgimpl.h>
2: #include <petsc/private/pcbddcimpl.h>
3: #include <petsc/private/pcbddcprivateimpl.h>
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;
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();
38: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &flg));
39: if (!flg) {
40: MatNullSpace nnsp;
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));
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;
66: /* submatrices and solvers */
67: PetscCall(KSPGetPC(smooth, &smoothpc));
68: PetscCall(PetscObjectTypeCompareAny((PetscObject)smoothpc, &flg, PCBDDC, ""));
69: if (!flg) {
70: Mat smoothA;
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;
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]));
97: /* analyze interface */
98: PetscCall(MatISGetLocalMat(A, &lA));
99: graph = ipcbddc->mat_graph;
100: if (!flg) {
101: PetscInt N;
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 **);
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;
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(PetscContainerSetUserDestroy(gcand, PCBDDCDestroyGraphCandidatesIS));
138: PetscCall(PetscObjectCompose((PetscObject)l2g, "_PCBDDCGraphCandidatesIS", (PetscObject)gcand));
139: PetscCall(PetscContainerDestroy(&gcand));
140: }
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;
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];
180: PetscCall(ISLocalToGlobalMappingApplyIS(l2g, cref, &sGi[0]));
181: PetscCall(ISDestroy(&cref));
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");
187: PetscCall(PetscMalloc1(ngct + 1, &ccii));
188: PetscCall(PetscMalloc1(ia[ngct], &cridx));
189: PetscCall(PetscMalloc1(ia[ngct], &cdata));
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");
196: PetscCall(MatRestoreRowIJ(cmat, 0, PETSC_FALSE, PETSC_FALSE, &i, &ia, &ja, &flg));
197: PetscCall(MatSeqAIJRestoreArrayRead(cmat, &tdata));
198: PetscCall(MatDestroy(&cmat));
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));
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: }
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;
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: }
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));
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));
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));
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;
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));
281: ptype = MATAIJ;
282: if (PetscDefined(HAVE_DEVICE)) {
283: PetscCall(MatBoundToCPU(A, &flg));
284: if (!flg) {
285: VecType vtype;
286: char *found = NULL;
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));
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;
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));
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));
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));
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: }