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:   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:     graph->commsizelimit = 0; /* don't use the COMM_SELF variant of the graph */
106:     PetscCall(PCBDDCGraphInit(graph, l2g, N, PETSC_MAX_INT));
107:     PetscCall(MatGetRowIJ(lA, 0, PETSC_TRUE, PETSC_FALSE, &graph->nvtxs_csr, (const PetscInt **)&graph->xadj, (const PetscInt **)&graph->adjncy, &flg));
108:     PetscCall(PCBDDCGraphSetUp(graph, vsize, NULL, NULL, 0, NULL, NULL));
109:     PetscCall(MatRestoreRowIJ(lA, 0, PETSC_TRUE, PETSC_FALSE, &graph->nvtxs_csr, (const PetscInt **)&graph->xadj, (const PetscInt **)&graph->adjncy, &flg));
110:     PetscCall(PCBDDCGraphComputeConnectedComponents(graph));
111:   }
112:   l2g = graph->l2gmap;
113:   if (reduced) {
114:     PetscContainer        gcand;
115:     PCBDDCGraphCandidates cand;
116:     PetscErrorCode (*rgdsw)(DM, PetscInt *, IS **);

118:     PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMComputeLocalRGDSWSets", &rgdsw));
119:     PetscCheck(rgdsw, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not supported");
120:     PetscCall(PetscNew(&cand));
121:     PetscCall((*rgdsw)(dm, &cand->nfc, &cand->Faces));
122:     /* filter interior (if any) and guarantee IS are ordered by global numbering */
123:     for (i = 0; i < cand->nfc; i++) {
124:       IS is, is2;

126:       PetscCall(ISLocalToGlobalMappingApplyIS(l2g, cand->Faces[i], &is));
127:       PetscCall(ISDestroy(&cand->Faces[i]));
128:       PetscCall(ISSort(is));
129:       PetscCall(ISGlobalToLocalMappingApplyIS(l2g, IS_GTOLM_DROP, is, &is2));
130:       PetscCall(ISDestroy(&is));
131:       PetscCall(ISGlobalToLocalMappingApplyIS(ipcis->BtoNmap, IS_GTOLM_DROP, is2, &is));
132:       PetscCall(ISDestroy(&is2));
133:       PetscCall(ISLocalToGlobalMappingApplyIS(ipcis->BtoNmap, is, &cand->Faces[i]));
134:       PetscCall(ISDestroy(&is));
135:     }
136:     PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &gcand));
137:     PetscCall(PetscContainerSetPointer(gcand, cand));
138:     PetscCall(PetscContainerSetUserDestroy(gcand, PCBDDCDestroyGraphCandidatesIS));
139:     PetscCall(PetscObjectCompose((PetscObject)l2g, "_PCBDDCGraphCandidatesIS", (PetscObject)gcand));
140:     PetscCall(PetscContainerDestroy(&gcand));
141:   }

143:   /* interface functions */
144:   otol                           = ipcbddc->adaptive_threshold[1];
145:   odbg                           = ipcbddc->dbg_flag;
146:   ominmax[0]                     = ipcbddc->adaptive_nmin;
147:   ominmax[1]                     = ipcbddc->adaptive_nmax;
148:   ipcbddc->adaptive_threshold[1] = tol;
149:   ipcbddc->dbg_flag              = dbg;
150:   ipcbddc->adaptive_nmin         = minmax[0];
151:   ipcbddc->adaptive_nmax         = minmax[1];
152:   if (tol != 0.0) { /* adaptive */
153:     Mat lS;

155:     PetscCall(MatCreateSchurComplement(ipcis->A_II, ipcis->pA_II, ipcis->A_IB, ipcis->A_BI, ipcis->A_BB, &lS));
156:     PetscCall(KSPGetOptionsPrefix(sksp[0], &prefix));
157:     PetscCall(PCBDDCSubSchursCreate(&sub_schurs));
158:     PetscCall(PCBDDCSubSchursInit(sub_schurs, prefix, ipcis->is_I_local, ipcis->is_B_local, graph, ipcis->BtoNmap, PETSC_FALSE, PETSC_TRUE));
159:     if (userdefined) PetscCall(PCBDDCComputeFakeChange(pcbddc, PETSC_FALSE, graph, NULL, &cmat, &cref, NULL, &flg));
160:     else {
161:       cmat = NULL;
162:       cref = NULL;
163:     }
164:     PetscCall(PCBDDCSubSchursSetUp(sub_schurs, lA, lS, PETSC_TRUE, NULL, NULL, -1, NULL, PETSC_TRUE, reuse_solver, PETSC_FALSE, 0, NULL, NULL, cmat, cref));
165:     PetscCall(MatDestroy(&lS));
166:     PetscCall(MatDestroy(&cmat));
167:     PetscCall(ISDestroy(&cref));
168:     if (sub_schurs->reuse_solver) {
169:       PetscCall(KSPSetPC(sksp[0], sub_schurs->reuse_solver->interior_solver));
170:       PetscCall(PCDestroy(&sub_schurs->reuse_solver->interior_solver));
171:       sub_schurs->reuse_solver = NULL;
172:     }
173:   }
174:   PetscCall(PCBDDCComputeFakeChange(pcbddc, PETSC_TRUE, graph, sub_schurs, &cmat, &cref, &sGiM[0], NULL));
175:   PetscCall(PCBDDCSubSchursDestroy(&sub_schurs));
176:   ipcbddc->adaptive_threshold[1] = otol;
177:   ipcbddc->dbg_flag              = odbg;
178:   ipcbddc->adaptive_nmin         = ominmax[0];
179:   ipcbddc->adaptive_nmax         = ominmax[1];

181:   PetscCall(ISLocalToGlobalMappingApplyIS(l2g, cref, &sGi[0]));
182:   PetscCall(ISDestroy(&cref));

184:   PetscCall(MatSeqAIJGetArrayRead(cmat, &tdata));
185:   PetscCall(MatGetRowIJ(cmat, 0, PETSC_FALSE, PETSC_FALSE, &ngct, &ia, &ja, &flg));
186:   PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in MatGetRowIJ");

188:   PetscCall(PetscMalloc1(ngct + 1, &ccii));
189:   PetscCall(PetscMalloc1(ia[ngct], &cridx));
190:   PetscCall(PetscMalloc1(ia[ngct], &cdata));

192:   PetscCall(PetscArraycpy(ccii, ia, ngct + 1));
193:   PetscCall(PetscArraycpy(cdata, tdata, ia[ngct]));
194:   PetscCall(ISGlobalToLocalMappingApply(ipcis->BtoNmap, IS_GTOLM_DROP, ia[ngct], ja, &i, cridx));
195:   PetscCheck(i == ia[ngct], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Error in G2L");

197:   PetscCall(MatRestoreRowIJ(cmat, 0, PETSC_FALSE, PETSC_FALSE, &i, &ia, &ja, &flg));
198:   PetscCall(MatSeqAIJRestoreArrayRead(cmat, &tdata));
199:   PetscCall(MatDestroy(&cmat));

201:   /* populate dense matrix */
202:   PetscCall(ISGetLocalSize(sG[0], &ng));
203:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, ng, ngct, NULL, &sGf[0]));
204:   PetscCall(MatDenseGetArrayWrite(sGf[0], &data));
205:   for (i = 0; i < ngct; i++)
206:     for (j = ccii[i]; j < ccii[i + 1]; j++) data[ng * i + cridx[j]] = cdata[j];
207:   PetscCall(MatDenseRestoreArrayWrite(sGf[0], &data));

209:   PetscCall(PetscFree(cdata));
210:   PetscCall(PetscFree(ccii));
211:   PetscCall(PetscFree(cridx));
212:   PetscCall(PCDestroy(&pcbddc));
213:   PetscCall(MatDestroy(&A));
214:   PetscFunctionReturn(PETSC_SUCCESS);
215: }

217: PetscErrorCode PCMGGDSWCreateCoarseSpace_Private(PC pc, PetscInt l, DM dm, KSP smooth, PetscInt Nc, Mat guess, Mat *cspace)
218: {
219:   KSP            *sksp;
220:   Mat             A, *sA_IG, *sGf, preallocator;
221:   IS              Gidx, GidxMult, cG;
222:   IS             *sI, *sG, *sGi, *sGiM;
223:   const PetscInt *cidx;
224:   PetscInt        NG, ns, n, i, c, rbs, cbs[2];
225:   PetscBool       flg;
226:   MatType         ptype;

228:   PetscFunctionBegin;
229:   *cspace = NULL;
230:   if (!l) PetscFunctionReturn(PETSC_SUCCESS);
231:   if (pc->useAmat) {
232:     PetscCall(KSPGetOperatorsSet(smooth, &flg, NULL));
233:     PetscCheck(flg, PetscObjectComm((PetscObject)smooth), PETSC_ERR_ORDER, "Amat not set");
234:     PetscCall(KSPGetOperators(smooth, &A, NULL));
235:   } else {
236:     PetscCall(KSPGetOperatorsSet(smooth, NULL, &flg));
237:     PetscCheck(flg, PetscObjectComm((PetscObject)smooth), PETSC_ERR_ORDER, "Pmat not set");
238:     PetscCall(KSPGetOperators(smooth, NULL, &A));
239:   }

241:   /* Setup (also setup smoother here) */
242:   if (!pc->setupcalled) PetscCall(KSPSetFromOptions(smooth));
243:   PetscCall(KSPSetUp(smooth));
244:   PetscCall(KSPSetUpOnBlocks(smooth));
245:   PetscCall(PCMGGDSWSetUp(pc, l, dm, smooth, Nc, A, &ns, &sA_IG, &sksp, &sI, &sG, &sGf, &sGi, &sGiM));

247:   /* Number GDSW basis functions */
248:   PetscCall(ISConcatenate(PetscObjectComm((PetscObject)A), ns, sGi, &Gidx));
249:   PetscCall(ISConcatenate(PetscObjectComm((PetscObject)A), ns, sGiM, &GidxMult));
250:   PetscCall(ISRenumber(Gidx, GidxMult, &NG, &cG));
251:   PetscCall(ISDestroy(&Gidx));

253:   /* Detect column block size */
254:   PetscCall(ISGetMinMax(GidxMult, &cbs[0], &cbs[1]));
255:   PetscCall(PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)A), cbs, cbs));
256:   PetscCall(ISDestroy(&GidxMult));

258:   /* Construct global interpolation matrix */
259:   PetscCall(MatGetLocalSize(A, NULL, &n));
260:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &preallocator));
261:   PetscCall(MatSetSizes(preallocator, n, PETSC_DECIDE, PETSC_DECIDE, NG));
262:   PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
263:   PetscCall(MatSetUp(preallocator));
264:   PetscCall(ISGetIndices(cG, &cidx));
265:   for (i = 0, c = 0; i < ns; i++) {
266:     const PetscInt *ri, *rg;
267:     PetscInt        nri, nrg, ncg;

269:     PetscCall(ISGetLocalSize(sI[i], &nri));
270:     PetscCall(ISGetLocalSize(sG[i], &nrg));
271:     PetscCall(ISGetIndices(sI[i], &ri));
272:     PetscCall(ISGetIndices(sG[i], &rg));
273:     PetscCall(MatGetSize(sGf[i], NULL, &ncg));
274:     PetscCall(MatSetValues(preallocator, nri, ri, ncg, cidx + c, NULL, INSERT_VALUES));
275:     PetscCall(MatSetValues(preallocator, nrg, rg, ncg, cidx + c, NULL, INSERT_VALUES));
276:     PetscCall(ISRestoreIndices(sI[i], &ri));
277:     PetscCall(ISRestoreIndices(sG[i], &rg));
278:   }
279:   PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
280:   PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));

282:   ptype = MATAIJ;
283:   if (PetscDefined(HAVE_DEVICE)) {
284:     PetscCall(MatBoundToCPU(A, &flg));
285:     if (!flg) {
286:       VecType vtype;
287:       char   *found = NULL;

289:       PetscCall(MatGetVecType(A, &vtype));
290:       PetscCall(PetscStrstr(vtype, "cuda", &found));
291:       if (found) ptype = MATAIJCUSPARSE;
292:     }
293:   }
294:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), cspace));
295:   PetscCall(MatSetSizes(*cspace, n, PETSC_DECIDE, PETSC_DECIDE, NG));
296:   PetscCall(MatSetType(*cspace, ptype));
297:   PetscCall(MatGetBlockSizes(A, NULL, &rbs));
298:   PetscCall(MatSetBlockSizes(*cspace, rbs, cbs[0] == cbs[1] ? cbs[0] : 1));
299:   PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_FALSE, *cspace));
300:   PetscCall(MatDestroy(&preallocator));
301:   PetscCall(MatSetOption(*cspace, MAT_ROW_ORIENTED, PETSC_FALSE));

303:   for (i = 0, c = 0; i < ns; i++) {
304:     Mat                X, Y;
305:     const PetscScalar *v;
306:     const PetscInt    *ri, *rg;
307:     PetscInt           nri, nrg, ncg;

309:     PetscCall(MatMatMult(sA_IG[i], sGf[i], MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Y));
310:     PetscCall(MatScale(Y, -1.0));
311:     PetscCall(MatDuplicate(Y, MAT_DO_NOT_COPY_VALUES, &X));
312:     PetscCall(KSPMatSolve(sksp[i], Y, X));

314:     PetscCall(ISGetLocalSize(sI[i], &nri));
315:     PetscCall(ISGetLocalSize(sG[i], &nrg));
316:     PetscCall(ISGetIndices(sI[i], &ri));
317:     PetscCall(ISGetIndices(sG[i], &rg));
318:     PetscCall(MatGetSize(sGf[i], NULL, &ncg));

320:     PetscCall(MatDenseGetArrayRead(X, &v));
321:     PetscCall(MatSetValues(*cspace, nri, ri, ncg, cidx + c, v, INSERT_VALUES));
322:     PetscCall(MatDenseRestoreArrayRead(X, &v));
323:     PetscCall(MatDenseGetArrayRead(sGf[i], &v));
324:     PetscCall(MatSetValues(*cspace, nrg, rg, ncg, cidx + c, v, INSERT_VALUES));
325:     PetscCall(MatDenseRestoreArrayRead(sGf[i], &v));
326:     PetscCall(ISRestoreIndices(sI[i], &ri));
327:     PetscCall(ISRestoreIndices(sG[i], &rg));
328:     PetscCall(MatDestroy(&Y));
329:     PetscCall(MatDestroy(&X));
330:   }
331:   PetscCall(ISRestoreIndices(cG, &cidx));
332:   PetscCall(ISDestroy(&cG));
333:   PetscCall(MatAssemblyBegin(*cspace, MAT_FINAL_ASSEMBLY));
334:   PetscCall(MatAssemblyEnd(*cspace, MAT_FINAL_ASSEMBLY));

336:   for (i = 0; i < ns; i++) {
337:     PetscCall(KSPDestroy(&sksp[i]));
338:     PetscCall(ISDestroy(&sI[i]));
339:     PetscCall(ISDestroy(&sG[i]));
340:     PetscCall(ISDestroy(&sGi[i]));
341:     PetscCall(ISDestroy(&sGiM[i]));
342:     PetscCall(MatDestroy(&sGf[i]));
343:     PetscCall(MatDestroy(&sA_IG[i]));
344:   }
345:   PetscCall(PetscFree(sksp));
346:   PetscCall(PetscFree(sI));
347:   PetscCall(PetscFree(sG));
348:   PetscCall(PetscFree(sGi));
349:   PetscCall(PetscFree(sGiM));
350:   PetscCall(PetscFree(sGf));
351:   PetscCall(PetscFree(sA_IG));
352:   PetscFunctionReturn(PETSC_SUCCESS);
353: }