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: }