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 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(PetscLogObjectParent((PetscObject)pc,(PetscObject)pcbddc)); 75 PetscCall(PCSetType(pcbddc,PCBDDC)); 76 PetscCall(PCSetOperators(pcbddc,smoothA,A)); 77 PetscCall(PCISSetUp(pcbddc,PETSC_TRUE,PETSC_FALSE)); 78 } else { 79 PetscCall(PetscObjectReference((PetscObject)smoothpc)); 80 pcbddc = smoothpc; 81 } 82 ipcis = (PC_IS*)pcbddc->data; 83 ipcbddc = (PC_BDDC*)pcbddc->data; 84 PetscCall(PetscObjectReference((PetscObject)ipcis->A_IB)); 85 PetscCall(PetscObjectReference((PetscObject)ipcis->is_I_global)); 86 PetscCall(PetscObjectReference((PetscObject)ipcis->is_B_global)); 87 sA_IG[0] = ipcis->A_IB; 88 sI[0] = ipcis->is_I_global; 89 sG[0] = ipcis->is_B_global; 90 91 PetscCall(KSPCreate(PetscObjectComm((PetscObject)ipcis->A_II),&sksp[0])); 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 { cmat = NULL; cref = NULL; } 161 PetscCall(PCBDDCSubSchursSetUp(sub_schurs,lA,lS,PETSC_TRUE,NULL,NULL,-1,NULL,PETSC_TRUE,reuse_solver,PETSC_FALSE,0,NULL,NULL,cmat,cref)); 162 PetscCall(MatDestroy(&lS)); 163 PetscCall(MatDestroy(&cmat)); 164 PetscCall(ISDestroy(&cref)); 165 if (sub_schurs->reuse_solver) { 166 PetscCall(KSPSetPC(sksp[0],sub_schurs->reuse_solver->interior_solver)); 167 PetscCall(PCDestroy(&sub_schurs->reuse_solver->interior_solver)); 168 sub_schurs->reuse_solver = NULL; 169 } 170 } 171 PetscCall(PCBDDCComputeFakeChange(pcbddc,PETSC_TRUE,graph,sub_schurs,&cmat,&cref,&sGiM[0],NULL)); 172 PetscCall(PCBDDCSubSchursDestroy(&sub_schurs)); 173 ipcbddc->adaptive_threshold[1] = otol; 174 ipcbddc->dbg_flag = odbg; 175 ipcbddc->adaptive_nmin = ominmax[0]; 176 ipcbddc->adaptive_nmax = ominmax[1]; 177 178 PetscCall(ISLocalToGlobalMappingApplyIS(l2g,cref,&sGi[0])); 179 PetscCall(ISDestroy(&cref)); 180 181 PetscCall(MatSeqAIJGetArrayRead(cmat,&tdata)); 182 PetscCall(MatGetRowIJ(cmat,0,PETSC_FALSE,PETSC_FALSE,&ngct,&ia,&ja,&flg)); 183 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatGetRowIJ"); 184 185 PetscCall(PetscMalloc1(ngct+1,&ccii)); 186 PetscCall(PetscMalloc1(ia[ngct],&cridx)); 187 PetscCall(PetscMalloc1(ia[ngct],&cdata)); 188 189 PetscCall(PetscArraycpy(ccii,ia,ngct+1)); 190 PetscCall(PetscArraycpy(cdata,tdata,ia[ngct])); 191 PetscCall(ISGlobalToLocalMappingApply(ipcis->BtoNmap,IS_GTOLM_DROP,ia[ngct],ja,&i,cridx)); 192 PetscCheck(i == ia[ngct],PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in G2L"); 193 194 PetscCall(MatRestoreRowIJ(cmat,0,PETSC_FALSE,PETSC_FALSE,&i,&ia,&ja,&flg)); 195 PetscCall(MatSeqAIJRestoreArrayRead(cmat,&tdata)); 196 PetscCall(MatDestroy(&cmat)); 197 198 /* populate dense matrix */ 199 PetscCall(ISGetLocalSize(sG[0],&ng)); 200 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ng,ngct,NULL,&sGf[0])); 201 PetscCall(MatDenseGetArrayWrite(sGf[0],&data)); 202 for (i = 0; i < ngct; i++) 203 for (j = ccii[i]; j < ccii[i+1]; j++) 204 data[ng*i + cridx[j]] = cdata[j]; 205 PetscCall(MatDenseRestoreArrayWrite(sGf[0],&data)); 206 207 PetscCall(PetscFree(cdata)); 208 PetscCall(PetscFree(ccii)); 209 PetscCall(PetscFree(cridx)); 210 PetscCall(PCDestroy(&pcbddc)); 211 PetscCall(MatDestroy(&A)); 212 PetscFunctionReturn(0); 213 } 214 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; 225 226 PetscFunctionBegin; 227 *cspace = NULL; 228 if (!l) PetscFunctionReturn(0); 229 if (pc->useAmat) { 230 PetscCall(KSPGetOperatorsSet(smooth,&flg,NULL)); 231 PetscCheck(flg,PetscObjectComm((PetscObject)smooth),PETSC_ERR_ORDER,"Amat not set"); 232 PetscCall(KSPGetOperators(smooth,&A,NULL)); 233 } else { 234 PetscCall(KSPGetOperatorsSet(smooth,NULL,&flg)); 235 PetscCheck(flg,PetscObjectComm((PetscObject)smooth),PETSC_ERR_ORDER,"Pmat not set"); 236 PetscCall(KSPGetOperators(smooth,NULL,&A)); 237 } 238 239 /* Setup (also setup smoother here) */ 240 if (!pc->setupcalled) PetscCall(KSPSetFromOptions(smooth)); 241 PetscCall(KSPSetUp(smooth)); 242 PetscCall(KSPSetUpOnBlocks(smooth)); 243 PetscCall(PCMGGDSWSetUp(pc,l,dm,smooth,Nc,A,&ns,&sA_IG,&sksp,&sI,&sG,&sGf,&sGi,&sGiM)); 244 245 /* Number GDSW basis functions */ 246 PetscCall(ISConcatenate(PetscObjectComm((PetscObject)A),ns,sGi,&Gidx)); 247 PetscCall(ISConcatenate(PetscObjectComm((PetscObject)A),ns,sGiM,&GidxMult)); 248 PetscCall(ISRenumber(Gidx,GidxMult,&NG,&cG)); 249 PetscCall(ISDestroy(&Gidx)); 250 251 /* Detect column block size */ 252 PetscCall(ISGetMinMax(GidxMult,&cbs[0],&cbs[1])); 253 PetscCall(PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)A),cbs,cbs)); 254 PetscCall(ISDestroy(&GidxMult)); 255 256 /* Construct global interpolation matrix */ 257 PetscCall(MatGetLocalSize(A,NULL,&n)); 258 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&preallocator)); 259 PetscCall(MatSetSizes(preallocator,n,PETSC_DECIDE,PETSC_DECIDE,NG)); 260 PetscCall(MatSetType(preallocator,MATPREALLOCATOR)); 261 PetscCall(MatSetUp(preallocator)); 262 PetscCall(ISGetIndices(cG,&cidx)); 263 for (i = 0, c = 0; i < ns; i++) { 264 const PetscInt *ri, *rg; 265 PetscInt nri, nrg, ncg; 266 267 PetscCall(ISGetLocalSize(sI[i],&nri)); 268 PetscCall(ISGetLocalSize(sG[i],&nrg)); 269 PetscCall(ISGetIndices(sI[i],&ri)); 270 PetscCall(ISGetIndices(sG[i],&rg)); 271 PetscCall(MatGetSize(sGf[i],NULL,&ncg)); 272 PetscCall(MatSetValues(preallocator,nri,ri,ncg,cidx+c,NULL,INSERT_VALUES)); 273 PetscCall(MatSetValues(preallocator,nrg,rg,ncg,cidx+c,NULL,INSERT_VALUES)); 274 PetscCall(ISRestoreIndices(sI[i],&ri)); 275 PetscCall(ISRestoreIndices(sG[i],&rg)); 276 } 277 PetscCall(MatAssemblyBegin(preallocator,MAT_FINAL_ASSEMBLY)); 278 PetscCall(MatAssemblyEnd(preallocator,MAT_FINAL_ASSEMBLY)); 279 280 ptype = MATAIJ; 281 if (PetscDefined(HAVE_DEVICE)) { 282 PetscCall(MatBoundToCPU(A,&flg)); 283 if (!flg) { 284 VecType vtype; 285 char *found; 286 287 PetscCall(MatGetVecType(A,&vtype)); 288 PetscCall(PetscStrstr(vtype,"cuda",&found)); 289 if (found) ptype = MATAIJCUSPARSE; 290 } 291 } 292 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),cspace)); 293 PetscCall(MatSetSizes(*cspace,n,PETSC_DECIDE,PETSC_DECIDE,NG)); 294 PetscCall(MatSetType(*cspace,ptype)); 295 PetscCall(MatGetBlockSizes(A,NULL,&rbs)); 296 PetscCall(MatSetBlockSizes(*cspace,rbs,cbs[0] == cbs[1] ? cbs[0] : 1)); 297 PetscCall(MatPreallocatorPreallocate(preallocator,PETSC_FALSE,*cspace)); 298 PetscCall(MatDestroy(&preallocator)); 299 PetscCall(MatSetOption(*cspace,MAT_ROW_ORIENTED,PETSC_FALSE)); 300 301 for (i = 0, c = 0; i < ns; i++) { 302 Mat X, Y; 303 const PetscScalar *v; 304 const PetscInt *ri, *rg; 305 PetscInt nri, nrg, ncg; 306 307 PetscCall(MatMatMult(sA_IG[i],sGf[i],MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Y)); 308 PetscCall(MatScale(Y,-1.0)); 309 PetscCall(MatDuplicate(Y,MAT_DO_NOT_COPY_VALUES,&X)); 310 PetscCall(KSPMatSolve(sksp[i],Y,X)); 311 312 PetscCall(ISGetLocalSize(sI[i],&nri)); 313 PetscCall(ISGetLocalSize(sG[i],&nrg)); 314 PetscCall(ISGetIndices(sI[i],&ri)); 315 PetscCall(ISGetIndices(sG[i],&rg)); 316 PetscCall(MatGetSize(sGf[i],NULL,&ncg)); 317 318 PetscCall(MatDenseGetArrayRead(X,&v)); 319 PetscCall(MatSetValues(*cspace,nri,ri,ncg,cidx+c,v,INSERT_VALUES)); 320 PetscCall(MatDenseRestoreArrayRead(X,&v)); 321 PetscCall(MatDenseGetArrayRead(sGf[i],&v)); 322 PetscCall(MatSetValues(*cspace,nrg,rg,ncg,cidx+c,v,INSERT_VALUES)); 323 PetscCall(MatDenseRestoreArrayRead(sGf[i],&v)); 324 PetscCall(ISRestoreIndices(sI[i],&ri)); 325 PetscCall(ISRestoreIndices(sG[i],&rg)); 326 PetscCall(MatDestroy(&Y)); 327 PetscCall(MatDestroy(&X)); 328 } 329 PetscCall(ISRestoreIndices(cG,&cidx)); 330 PetscCall(ISDestroy(&cG)); 331 PetscCall(MatAssemblyBegin(*cspace,MAT_FINAL_ASSEMBLY)); 332 PetscCall(MatAssemblyEnd(*cspace,MAT_FINAL_ASSEMBLY)); 333 334 for (i = 0; i < ns; i++) { 335 PetscCall(KSPDestroy(&sksp[i])); 336 PetscCall(ISDestroy(&sI[i])); 337 PetscCall(ISDestroy(&sG[i])); 338 PetscCall(ISDestroy(&sGi[i])); 339 PetscCall(ISDestroy(&sGiM[i])); 340 PetscCall(MatDestroy(&sGf[i])); 341 PetscCall(MatDestroy(&sA_IG[i])); 342 } 343 PetscCall(PetscFree(sksp)); 344 PetscCall(PetscFree(sI)); 345 PetscCall(PetscFree(sG)); 346 PetscCall(PetscFree(sGi)); 347 PetscCall(PetscFree(sGiM)); 348 PetscCall(PetscFree(sGf)); 349 PetscCall(PetscFree(sA_IG)); 350 PetscFunctionReturn(0); 351 } 352