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