1 #include <petsc/private/pcmgimpl.h> 2 #include <petsc/private/pcbddcimpl.h> 3 #include <petsc/private/pcbddcprivateimpl.h> 4 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; 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 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 **); 117 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; 125 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 } 142 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; 154 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]; 180 181 PetscCall(ISLocalToGlobalMappingApplyIS(l2g, cref, &sGi[0])); 182 PetscCall(ISDestroy(&cref)); 183 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"); 187 188 PetscCall(PetscMalloc1(ngct + 1, &ccii)); 189 PetscCall(PetscMalloc1(ia[ngct], &cridx)); 190 PetscCall(PetscMalloc1(ia[ngct], &cdata)); 191 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"); 196 197 PetscCall(MatRestoreRowIJ(cmat, 0, PETSC_FALSE, PETSC_FALSE, &i, &ia, &ja, &flg)); 198 PetscCall(MatSeqAIJRestoreArrayRead(cmat, &tdata)); 199 PetscCall(MatDestroy(&cmat)); 200 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)); 208 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 } 216 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; 227 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 } 240 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)); 246 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)); 252 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)); 257 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; 268 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)); 281 282 ptype = MATAIJ; 283 if (PetscDefined(HAVE_DEVICE)) { 284 PetscCall(MatBoundToCPU(A, &flg)); 285 if (!flg) { 286 VecType vtype; 287 char *found = NULL; 288 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)); 302 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; 308 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)); 313 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)); 319 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)); 335 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 } 354