1 #include <../src/ksp/pc/impls/bddc/bddc.h> 2 #include <../src/ksp/pc/impls/bddc/bddcprivate.h> 3 4 static PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt*,PetscInt,PetscBT,PetscInt*,PetscInt*,PetscInt*); 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "PCBDDCSubSchursSetUpNew" 8 PetscErrorCode PCBDDCSubSchursSetUpNew(PCBDDCSubSchurs sub_schurs, PetscInt xadj[], PetscInt adjncy[], PetscInt nlayers) 9 { 10 Mat A_II,A_IB,A_BI,A_BB; 11 Mat AE_II,*AE_IE,*AE_EI,*AE_EE; 12 IS is_I,*is_subset_B; 13 ISLocalToGlobalMapping BtoNmap; 14 PetscInt i; 15 PetscBool implicit_schurs; 16 PetscErrorCode ierr; 17 18 PetscFunctionBegin; 19 ierr = PetscObjectTypeCompare((PetscObject)sub_schurs->S,MATSCHURCOMPLEMENT,&implicit_schurs);CHKERRQ(ierr); 20 21 /* allocate space for schur complements */ 22 ierr = PetscMalloc5(sub_schurs->n_subs,&sub_schurs->is_AEj_I, 23 sub_schurs->n_subs,&sub_schurs->is_AEj_B, 24 sub_schurs->n_subs,&sub_schurs->S_Ej, 25 sub_schurs->n_subs,&sub_schurs->work1, 26 sub_schurs->n_subs,&sub_schurs->work2);CHKERRQ(ierr); 27 28 /* get Schur complement matrices */ 29 if (implicit_schurs) { 30 ierr = MatSchurComplementGetSubMatrices(sub_schurs->S,&A_II,NULL,&A_IB,&A_BI,&A_BB);CHKERRQ(ierr); 31 ierr = PetscMalloc4(sub_schurs->n_subs,&is_subset_B, 32 sub_schurs->n_subs,&AE_IE, 33 sub_schurs->n_subs,&AE_EI, 34 sub_schurs->n_subs,&AE_EE);CHKERRQ(ierr); 35 } 36 37 /* determine interior problems */ 38 if (nlayers >= 0 && xadj != NULL && adjncy != NULL) { /* Interior problems can be different from the original one */ 39 PetscBT touched; 40 const PetscInt* idx_B; 41 PetscInt n_I,n_B,n_local_dofs,n_prev_added,j,layer,*local_numbering; 42 43 /* get sizes */ 44 ierr = ISGetLocalSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr); 45 ierr = ISGetLocalSize(sub_schurs->is_B,&n_B);CHKERRQ(ierr); 46 47 ierr = PetscMalloc1(n_I+n_B,&local_numbering);CHKERRQ(ierr); 48 ierr = PetscBTCreate(n_I+n_B,&touched);CHKERRQ(ierr); 49 ierr = PetscBTMemzero(n_I+n_B,touched);CHKERRQ(ierr); 50 51 /* all boundary dofs must be skipped when adding layers */ 52 ierr = ISGetIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr); 53 for (j=0;j<n_B;j++) { 54 ierr = PetscBTSet(touched,idx_B[j]);CHKERRQ(ierr); 55 } 56 ierr = PetscMemcpy(local_numbering,idx_B,n_B*sizeof(PetscInt));CHKERRQ(ierr); 57 ierr = ISRestoreIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr); 58 59 /* add prescribed number of layers of dofs */ 60 n_local_dofs = n_B; 61 n_prev_added = n_B; 62 for (layer=0;layer<nlayers;layer++) { 63 PetscInt n_added; 64 if (n_local_dofs == n_I+n_B) break; 65 if (n_local_dofs > n_I+n_B) { 66 SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error querying layer %d. Out of bound access (%d > %d)",layer,n_local_dofs,n_I+n_B); 67 } 68 ierr = PCBDDCAdjGetNextLayer_Private(local_numbering+n_local_dofs,n_prev_added,touched,xadj,adjncy,&n_added);CHKERRQ(ierr); 69 n_prev_added = n_added; 70 n_local_dofs += n_added; 71 if (!n_added) break; 72 } 73 ierr = PetscBTDestroy(&touched);CHKERRQ(ierr); 74 75 /* IS for I dofs in original numbering */ 76 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sub_schurs->is_I),n_local_dofs-n_B,local_numbering+n_B,PETSC_COPY_VALUES,&sub_schurs->is_AEj_I[0]);CHKERRQ(ierr); 77 ierr = PetscFree(local_numbering);CHKERRQ(ierr); 78 ierr = ISSort(sub_schurs->is_AEj_I[0]);CHKERRQ(ierr); 79 /* IS for I dofs in boundary numbering */ 80 if (implicit_schurs) { 81 ISLocalToGlobalMapping ItoNmap; 82 ierr = ISLocalToGlobalMappingCreateIS(sub_schurs->is_I,&ItoNmap);CHKERRQ(ierr); 83 ierr = ISGlobalToLocalMappingApplyIS(ItoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_I[0],&is_I);CHKERRQ(ierr); 84 ierr = ISLocalToGlobalMappingDestroy(&ItoNmap);CHKERRQ(ierr); 85 86 /* II block */ 87 ierr = MatGetSubMatrix(A_II,is_I,is_I,MAT_INITIAL_MATRIX,&AE_II);CHKERRQ(ierr); 88 } 89 } else { 90 PetscInt n_I; 91 92 /* IS for I dofs in original numbering */ 93 ierr = PetscObjectReference((PetscObject)sub_schurs->is_I);CHKERRQ(ierr); 94 sub_schurs->is_AEj_I[0] = sub_schurs->is_I; 95 96 /* IS for I dofs in I numbering (strided 1) */ 97 if (implicit_schurs) { 98 ierr = ISGetSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr); 99 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub_schurs->is_I),n_I,0,1,&is_I);CHKERRQ(ierr); 100 101 /* II block is the same */ 102 ierr = PetscObjectReference((PetscObject)A_II);CHKERRQ(ierr); 103 AE_II = A_II; 104 } 105 } 106 107 /* TODO: just for compatibility with the previous version, needs to be fixed */ 108 for (i=1;i<sub_schurs->n_subs;i++) { 109 ierr = PetscObjectReference((PetscObject)sub_schurs->is_AEj_I[0]);CHKERRQ(ierr); 110 sub_schurs->is_AEj_I[i] = sub_schurs->is_AEj_I[0]; 111 } 112 113 if (implicit_schurs) { 114 /* subsets in original and boundary numbering */ 115 ierr = ISLocalToGlobalMappingCreateIS(sub_schurs->is_B,&BtoNmap);CHKERRQ(ierr); 116 for (i=0;i<sub_schurs->n_subs;i++) { 117 ierr = ISDuplicate(sub_schurs->is_subs[i],&sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 118 ierr = ISSort(sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 119 ierr = ISGlobalToLocalMappingApplyIS(BtoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_B[i],&is_subset_B[i]);CHKERRQ(ierr); 120 } 121 ierr = ISLocalToGlobalMappingDestroy(&BtoNmap);CHKERRQ(ierr); 122 123 /* EE block */ 124 for (i=0;i<sub_schurs->n_subs;i++) { 125 ierr = MatGetSubMatrix(A_BB,is_subset_B[i],is_subset_B[i],MAT_INITIAL_MATRIX,&AE_EE[i]);CHKERRQ(ierr); 126 } 127 /* IE block */ 128 for (i=0;i<sub_schurs->n_subs;i++) { 129 ierr = MatGetSubMatrix(A_IB,is_I,is_subset_B[i],MAT_INITIAL_MATRIX,&AE_IE[i]);CHKERRQ(ierr); 130 } 131 /* EI block */ 132 for (i=0;i<sub_schurs->n_subs;i++) { 133 ierr = MatGetSubMatrix(A_BI,is_subset_B[i],is_I,MAT_INITIAL_MATRIX,&AE_EI[i]);CHKERRQ(ierr); 134 } 135 136 /* setup Schur complements on subset */ 137 for (i=0;i<sub_schurs->n_subs;i++) { 138 ierr = MatCreateSchurComplement(AE_II,AE_II,AE_IE[i],AE_EI[i],AE_EE[i],&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 139 ierr = MatCreateVecs(sub_schurs->S_Ej[i],&sub_schurs->work1[i],&sub_schurs->work2[i]);CHKERRQ(ierr); 140 if (AE_II == A_II) { /* we can reuse the same ksp */ 141 KSP ksp; 142 ierr = MatSchurComplementGetKSP(sub_schurs->S,&ksp);CHKERRQ(ierr); 143 ierr = MatSchurComplementSetKSP(sub_schurs->S_Ej[i],ksp);CHKERRQ(ierr); 144 } else { /* build new ksp object which inherits ksp and pc types from the original one */ 145 KSP origksp,schurksp; 146 PC origpc,schurpc; 147 KSPType ksp_type; 148 PCType pc_type; 149 PetscInt n_internal; 150 151 ierr = MatSchurComplementGetKSP(sub_schurs->S,&origksp);CHKERRQ(ierr); 152 ierr = MatSchurComplementGetKSP(sub_schurs->S_Ej[i],&schurksp);CHKERRQ(ierr); 153 ierr = KSPGetType(origksp,&ksp_type);CHKERRQ(ierr); 154 ierr = KSPSetType(schurksp,ksp_type);CHKERRQ(ierr); 155 ierr = KSPGetPC(schurksp,&schurpc);CHKERRQ(ierr); 156 ierr = KSPGetPC(origksp,&origpc);CHKERRQ(ierr); 157 ierr = PCGetType(origpc,&pc_type);CHKERRQ(ierr); 158 ierr = PCSetType(schurpc,pc_type);CHKERRQ(ierr); 159 ierr = ISGetSize(is_I,&n_internal);CHKERRQ(ierr); 160 if (n_internal) { /* UMFPACK gives error with 0 sized problems */ 161 MatSolverPackage solver=NULL; 162 ierr = PCFactorGetMatSolverPackage(origpc,(const MatSolverPackage*)&solver);CHKERRQ(ierr); 163 if (solver) { 164 ierr = PCFactorSetMatSolverPackage(schurpc,solver);CHKERRQ(ierr); 165 } 166 } 167 ierr = KSPSetUp(schurksp);CHKERRQ(ierr); 168 } 169 } 170 /* free */ 171 ierr = ISDestroy(&is_I);CHKERRQ(ierr); 172 ierr = MatDestroy(&AE_II);CHKERRQ(ierr); 173 for (i=0;i<sub_schurs->n_subs;i++) { 174 ierr = MatDestroy(&AE_EE[i]);CHKERRQ(ierr); 175 ierr = MatDestroy(&AE_IE[i]);CHKERRQ(ierr); 176 ierr = MatDestroy(&AE_EI[i]);CHKERRQ(ierr); 177 ierr = ISDestroy(&is_subset_B[i]);CHKERRQ(ierr); 178 } 179 ierr = PetscFree4(is_subset_B,AE_IE,AE_EI,AE_EE);CHKERRQ(ierr); 180 } 181 PetscFunctionReturn(0); 182 } 183 184 #undef __FUNCT__ 185 #define __FUNCT__ "PCBDDCSubSchursInit" 186 PetscErrorCode PCBDDCSubSchursInit(PCBDDCSubSchurs sub_schurs, Mat A, Mat S, IS is_I, IS is_B, PCBDDCGraph graph, PetscInt seqthreshold) 187 { 188 IS *faces,*edges,*all_cc; 189 PetscInt *index_sequential,*index_parallel; 190 PetscInt *auxlocal_sequential,*auxlocal_parallel; 191 PetscInt *auxglobal_sequential,*auxglobal_parallel; 192 PetscInt *auxmapping;//,*idxs; 193 PetscInt i,max_subset_size; 194 PetscInt n_sequential_problems,n_local_sequential_problems,n_parallel_problems,n_local_parallel_problems; 195 PetscInt n_faces,n_edges,n_all_cc; 196 PetscBool is_sorted; 197 PetscErrorCode ierr; 198 199 PetscFunctionBegin; 200 ierr = ISSorted(is_I,&is_sorted);CHKERRQ(ierr); 201 if (!is_sorted) { 202 SETERRQ(PetscObjectComm((PetscObject)is_I),PETSC_ERR_PLIB,"IS for I dofs should be shorted"); 203 } 204 ierr = ISSorted(is_B,&is_sorted);CHKERRQ(ierr); 205 if (!is_sorted) { 206 SETERRQ(PetscObjectComm((PetscObject)is_B),PETSC_ERR_PLIB,"IS for B dofs should be shorted"); 207 } 208 209 /* reset any previous data */ 210 ierr = PCBDDCSubSchursReset(sub_schurs);CHKERRQ(ierr); 211 212 /* get index sets for faces and edges */ 213 ierr = PCBDDCGraphGetCandidatesIS(graph,&n_faces,&faces,&n_edges,&edges,NULL);CHKERRQ(ierr); 214 n_all_cc = n_faces+n_edges; 215 ierr = PetscMalloc1(n_all_cc,&all_cc);CHKERRQ(ierr); 216 for (i=0;i<n_faces;i++) { 217 all_cc[i] = faces[i]; 218 } 219 for (i=0;i<n_edges;i++) { 220 all_cc[n_faces+i] = edges[i]; 221 } 222 ierr = PetscFree(faces);CHKERRQ(ierr); 223 ierr = PetscFree(edges);CHKERRQ(ierr); 224 225 /* map interface's subsets */ 226 max_subset_size = 0; 227 for (i=0;i<n_all_cc;i++) { 228 PetscInt subset_size; 229 ierr = ISGetLocalSize(all_cc[i],&subset_size);CHKERRQ(ierr); 230 max_subset_size = PetscMax(max_subset_size,subset_size); 231 } 232 ierr = PetscMalloc1(max_subset_size,&auxmapping);CHKERRQ(ierr); 233 ierr = PetscMalloc2(graph->ncc,&auxlocal_sequential, 234 graph->ncc,&auxlocal_parallel);CHKERRQ(ierr); 235 ierr = PetscMalloc2(graph->ncc,&index_sequential, 236 graph->ncc,&index_parallel);CHKERRQ(ierr); 237 238 /* if threshold is negative, uses all sequential problems */ 239 if (seqthreshold < 0) seqthreshold = max_subset_size; 240 241 /* determine which problem has to be solved in parallel or sequentially */ 242 n_local_sequential_problems = 0; 243 n_local_parallel_problems = 0; 244 for (i=0;i<n_all_cc;i++) { 245 PetscInt subset_size,j,min_loc = 0; 246 const PetscInt *idxs; 247 248 ierr = ISGetLocalSize(all_cc[i],&subset_size);CHKERRQ(ierr); 249 ierr = ISGetIndices(all_cc[i],&idxs);CHKERRQ(ierr); 250 ierr = ISLocalToGlobalMappingApply(graph->l2gmap,subset_size,idxs,auxmapping);CHKERRQ(ierr); 251 for (j=1;j<subset_size;j++) { 252 if (auxmapping[j]<auxmapping[min_loc]) { 253 min_loc = j; 254 } 255 } 256 if (subset_size > seqthreshold) { 257 index_parallel[n_local_parallel_problems] = i; 258 auxlocal_parallel[n_local_parallel_problems] = idxs[min_loc]; 259 n_local_parallel_problems++; 260 } else { 261 index_sequential[n_local_sequential_problems] = i; 262 auxlocal_sequential[n_local_sequential_problems] = idxs[min_loc]; 263 n_local_sequential_problems++; 264 } 265 ierr = ISRestoreIndices(all_cc[i],&idxs);CHKERRQ(ierr); 266 } 267 268 /* Number parallel problems */ 269 auxglobal_parallel = 0; 270 ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)graph->l2gmap),graph->l2gmap,n_local_parallel_problems,auxlocal_parallel,PETSC_NULL,&n_parallel_problems,&auxglobal_parallel);CHKERRQ(ierr); 271 272 /* Number sequential problems */ 273 auxglobal_sequential = 0; 274 ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)graph->l2gmap),graph->l2gmap,n_local_sequential_problems,auxlocal_sequential,PETSC_NULL,&n_sequential_problems,&auxglobal_sequential);CHKERRQ(ierr); 275 276 /* update info in sub_schurs */ 277 if (A) { 278 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 279 sub_schurs->A = A; 280 } 281 ierr = PetscObjectReference((PetscObject)S);CHKERRQ(ierr); 282 sub_schurs->S = S; 283 ierr = PetscObjectReference((PetscObject)is_I);CHKERRQ(ierr); 284 sub_schurs->is_I = is_I; 285 ierr = PetscObjectReference((PetscObject)is_B);CHKERRQ(ierr); 286 sub_schurs->is_B = is_B; 287 sub_schurs->n_subs_seq = n_local_sequential_problems; 288 sub_schurs->n_subs_par = n_local_parallel_problems; 289 sub_schurs->n_subs_seq_g = n_sequential_problems; 290 sub_schurs->n_subs_par_g = n_parallel_problems; 291 sub_schurs->n_subs = sub_schurs->n_subs_seq + sub_schurs->n_subs_par; 292 sub_schurs->is_subs = all_cc; 293 sub_schurs->index_sequential = index_sequential; 294 sub_schurs->index_parallel = index_parallel; 295 sub_schurs->auxglobal_sequential = auxglobal_sequential; 296 sub_schurs->auxglobal_parallel = auxglobal_parallel; 297 298 /* free workspace */ 299 ierr = PetscFree(auxmapping);CHKERRQ(ierr); 300 ierr = PetscFree2(auxlocal_sequential,auxlocal_parallel);CHKERRQ(ierr); 301 302 PetscFunctionReturn(0); 303 } 304 305 #undef __FUNCT__ 306 #define __FUNCT__ "PCBDDCSubSchursCreate" 307 PetscErrorCode PCBDDCSubSchursCreate(PCBDDCSubSchurs *sub_schurs) 308 { 309 PCBDDCSubSchurs schurs_ctx; 310 PetscErrorCode ierr; 311 312 PetscFunctionBegin; 313 ierr = PetscNew(&schurs_ctx);CHKERRQ(ierr); 314 schurs_ctx->n_subs = 0; 315 *sub_schurs = schurs_ctx; 316 PetscFunctionReturn(0); 317 } 318 319 #undef __FUNCT__ 320 #define __FUNCT__ "PCBDDCSubSchursDestroy" 321 PetscErrorCode PCBDDCSubSchursDestroy(PCBDDCSubSchurs *sub_schurs) 322 { 323 PetscErrorCode ierr; 324 325 PetscFunctionBegin; 326 ierr = PCBDDCSubSchursReset(*sub_schurs);CHKERRQ(ierr); 327 ierr = PetscFree(*sub_schurs);CHKERRQ(ierr); 328 PetscFunctionReturn(0); 329 } 330 331 #undef __FUNCT__ 332 #define __FUNCT__ "PCBDDCSubSchursReset" 333 PetscErrorCode PCBDDCSubSchursReset(PCBDDCSubSchurs sub_schurs) 334 { 335 PetscInt i; 336 PetscErrorCode ierr; 337 338 PetscFunctionBegin; 339 ierr = MatDestroy(&sub_schurs->A);CHKERRQ(ierr); 340 ierr = MatDestroy(&sub_schurs->S);CHKERRQ(ierr); 341 ierr = ISDestroy(&sub_schurs->is_I);CHKERRQ(ierr); 342 ierr = ISDestroy(&sub_schurs->is_B);CHKERRQ(ierr); 343 for (i=0;i<sub_schurs->n_subs;i++) { 344 ierr = ISDestroy(&sub_schurs->is_subs[i]);CHKERRQ(ierr); 345 ierr = ISDestroy(&sub_schurs->is_AEj_I[i]);CHKERRQ(ierr); 346 ierr = ISDestroy(&sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 347 ierr = MatDestroy(&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 348 ierr = VecDestroy(&sub_schurs->work1[i]);CHKERRQ(ierr); 349 ierr = VecDestroy(&sub_schurs->work2[i]);CHKERRQ(ierr); 350 } 351 if (sub_schurs->n_subs) { 352 ierr = PetscFree(sub_schurs->is_subs);CHKERRQ(ierr); 353 ierr = PetscFree5(sub_schurs->is_AEj_I,sub_schurs->is_AEj_B,sub_schurs->S_Ej,sub_schurs->work1,sub_schurs->work2);CHKERRQ(ierr); 354 ierr = PetscFree2(sub_schurs->index_sequential,sub_schurs->index_parallel);CHKERRQ(ierr); 355 ierr = PetscFree(sub_schurs->auxglobal_sequential);CHKERRQ(ierr); 356 ierr = PetscFree(sub_schurs->auxglobal_parallel);CHKERRQ(ierr); 357 } 358 sub_schurs->n_subs = 0; 359 PetscFunctionReturn(0); 360 } 361 362 #undef __FUNCT__ 363 #define __FUNCT__ "PCBDDCAdjGetNextLayer_Private" 364 PETSC_STATIC_INLINE PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt* queue_tip,PetscInt n_prev,PetscBT touched,PetscInt* xadj,PetscInt* adjncy,PetscInt* n_added) 365 { 366 PetscInt i,j,n; 367 PetscErrorCode ierr; 368 369 PetscFunctionBegin; 370 n = 0; 371 for (i=-n_prev;i<0;i++) { 372 PetscInt start_dof = queue_tip[i]; 373 for (j=xadj[start_dof];j<xadj[start_dof+1];j++) { 374 PetscInt dof = adjncy[j]; 375 if (!PetscBTLookup(touched,dof)) { 376 ierr = PetscBTSet(touched,dof);CHKERRQ(ierr); 377 queue_tip[n] = dof; 378 n++; 379 } 380 } 381 } 382 *n_added = n; 383 PetscFunctionReturn(0); 384 } 385 386 #undef __FUNCT__ 387 #define __FUNCT__ "PCBDDCSubSchursSetUp" 388 PetscErrorCode PCBDDCSubSchursSetUp(PCBDDCSubSchurs sub_schurs, Mat S, IS is_A_I, IS is_A_B, PetscInt ncc, IS is_cc[], PetscInt xadj[], PetscInt adjncy[], PetscInt nlayers) 389 { 390 Mat A_II,A_IB,A_BI,A_BB; 391 Mat AE_II,*AE_IE,*AE_EI,*AE_EE; 392 IS is_I,*is_subset_B; 393 ISLocalToGlobalMapping BtoNmap; 394 PetscInt i; 395 PetscBool is_sorted; 396 PetscErrorCode ierr; 397 398 PetscFunctionBegin; 399 ierr = ISSorted(is_A_I,&is_sorted);CHKERRQ(ierr); 400 if (!is_sorted) { 401 SETERRQ(PetscObjectComm((PetscObject)is_A_I),PETSC_ERR_PLIB,"IS for I dofs should be shorted"); 402 } 403 ierr = ISSorted(is_A_B,&is_sorted);CHKERRQ(ierr); 404 if (!is_sorted) { 405 SETERRQ(PetscObjectComm((PetscObject)is_A_B),PETSC_ERR_PLIB,"IS for B dofs should be shorted"); 406 } 407 408 /* get Schur complement matrices */ 409 ierr = MatSchurComplementGetSubMatrices(S,&A_II,NULL,&A_IB,&A_BI,&A_BB);CHKERRQ(ierr); 410 411 /* allocate space for schur complements */ 412 ierr = PetscMalloc5(ncc,&sub_schurs->is_AEj_I,ncc,&sub_schurs->is_AEj_B,ncc,&sub_schurs->S_Ej,ncc,&sub_schurs->work1,ncc,&sub_schurs->work2);CHKERRQ(ierr); 413 ierr = PetscMalloc4(ncc,&is_subset_B,ncc,&AE_IE,ncc,&AE_EI,ncc,&AE_EE);CHKERRQ(ierr); 414 sub_schurs->n_subs = ncc; 415 416 /* maps */ 417 if (sub_schurs->n_subs && nlayers >= 0 && xadj != NULL && adjncy != NULL) { /* Interior problems can be different from the original one */ 418 ISLocalToGlobalMapping ItoNmap; 419 PetscBT touched; 420 const PetscInt* idx_B; 421 PetscInt n_I,n_B,n_local_dofs,n_prev_added,j,layer,*local_numbering; 422 423 /* get sizes */ 424 ierr = ISGetLocalSize(is_A_I,&n_I);CHKERRQ(ierr); 425 ierr = ISGetLocalSize(is_A_B,&n_B);CHKERRQ(ierr); 426 427 ierr = ISLocalToGlobalMappingCreateIS(is_A_I,&ItoNmap);CHKERRQ(ierr); 428 ierr = PetscMalloc1(n_I+n_B,&local_numbering);CHKERRQ(ierr); 429 ierr = PetscBTCreate(n_I+n_B,&touched);CHKERRQ(ierr); 430 ierr = PetscBTMemzero(n_I+n_B,touched);CHKERRQ(ierr); 431 432 /* all boundary dofs must be skipped when adding layers */ 433 ierr = ISGetIndices(is_A_B,&idx_B);CHKERRQ(ierr); 434 for (j=0;j<n_B;j++) { 435 ierr = PetscBTSet(touched,idx_B[j]);CHKERRQ(ierr); 436 } 437 ierr = PetscMemcpy(local_numbering,idx_B,n_B*sizeof(PetscInt));CHKERRQ(ierr); 438 ierr = ISRestoreIndices(is_A_B,&idx_B);CHKERRQ(ierr); 439 440 /* add next layers of dofs */ 441 n_local_dofs = n_B; 442 n_prev_added = n_B; 443 for (layer=0;layer<nlayers;layer++) { 444 PetscInt n_added; 445 if (n_local_dofs == n_I+n_B) break; 446 if (n_local_dofs > n_I+n_B) { 447 SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error querying layer %d. Out of bound access (%d > %d)",layer,n_local_dofs,n_I+n_B); 448 } 449 ierr = PCBDDCAdjGetNextLayer_Private(local_numbering+n_local_dofs,n_prev_added,touched,xadj,adjncy,&n_added);CHKERRQ(ierr); 450 n_prev_added = n_added; 451 n_local_dofs += n_added; 452 if (!n_added) break; 453 } 454 ierr = PetscBTDestroy(&touched);CHKERRQ(ierr); 455 456 /* IS for I dofs in original numbering and in I numbering */ 457 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ItoNmap),n_local_dofs-n_B,local_numbering+n_B,PETSC_COPY_VALUES,&sub_schurs->is_AEj_I[0]);CHKERRQ(ierr); 458 ierr = PetscFree(local_numbering);CHKERRQ(ierr); 459 ierr = ISSort(sub_schurs->is_AEj_I[0]);CHKERRQ(ierr); 460 ierr = ISGlobalToLocalMappingApplyIS(ItoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_I[0],&is_I);CHKERRQ(ierr); 461 ierr = ISLocalToGlobalMappingDestroy(&ItoNmap);CHKERRQ(ierr); 462 463 /* II block */ 464 ierr = MatGetSubMatrix(A_II,is_I,is_I,MAT_INITIAL_MATRIX,&AE_II);CHKERRQ(ierr); 465 } else { 466 PetscInt n_I; 467 468 /* IS for I dofs in original numbering */ 469 ierr = PetscObjectReference((PetscObject)is_A_I);CHKERRQ(ierr); 470 sub_schurs->is_AEj_I[0] = is_A_I; 471 472 /* IS for I dofs in I numbering (strided 1) */ 473 ierr = ISGetSize(is_A_I,&n_I);CHKERRQ(ierr); 474 ierr = ISCreateStride(PetscObjectComm((PetscObject)is_A_I),n_I,0,1,&is_I);CHKERRQ(ierr); 475 476 /* II block is the same */ 477 ierr = PetscObjectReference((PetscObject)A_II);CHKERRQ(ierr); 478 AE_II = A_II; 479 } 480 481 /* TODO: just for compatibility with the previous version, needs to be fixed */ 482 for (i=1;i<sub_schurs->n_subs;i++) { 483 ierr = PetscObjectReference((PetscObject)sub_schurs->is_AEj_I[0]);CHKERRQ(ierr); 484 sub_schurs->is_AEj_I[i] = sub_schurs->is_AEj_I[0]; 485 } 486 487 /* subsets in original and boundary numbering */ 488 ierr = ISLocalToGlobalMappingCreateIS(is_A_B,&BtoNmap);CHKERRQ(ierr); 489 for (i=0;i<sub_schurs->n_subs;i++) { 490 ierr = ISDuplicate(is_cc[i],&sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 491 ierr = ISSort(sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 492 ierr = ISGlobalToLocalMappingApplyIS(BtoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_B[i],&is_subset_B[i]);CHKERRQ(ierr); 493 } 494 ierr = ISLocalToGlobalMappingDestroy(&BtoNmap);CHKERRQ(ierr); 495 496 /* EE block */ 497 for (i=0;i<sub_schurs->n_subs;i++) { 498 ierr = MatGetSubMatrix(A_BB,is_subset_B[i],is_subset_B[i],MAT_INITIAL_MATRIX,&AE_EE[i]);CHKERRQ(ierr); 499 } 500 /* IE block */ 501 for (i=0;i<sub_schurs->n_subs;i++) { 502 ierr = MatGetSubMatrix(A_IB,is_I,is_subset_B[i],MAT_INITIAL_MATRIX,&AE_IE[i]);CHKERRQ(ierr); 503 } 504 /* EI block */ 505 for (i=0;i<sub_schurs->n_subs;i++) { 506 ierr = MatGetSubMatrix(A_BI,is_subset_B[i],is_I,MAT_INITIAL_MATRIX,&AE_EI[i]);CHKERRQ(ierr); 507 } 508 509 /* setup Schur complements on subset */ 510 for (i=0;i<sub_schurs->n_subs;i++) { 511 ierr = MatCreateSchurComplement(AE_II,AE_II,AE_IE[i],AE_EI[i],AE_EE[i],&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 512 ierr = MatCreateVecs(sub_schurs->S_Ej[i],&sub_schurs->work1[i],&sub_schurs->work2[i]);CHKERRQ(ierr); 513 if (AE_II == A_II) { /* we can reuse the same ksp */ 514 KSP ksp; 515 ierr = MatSchurComplementGetKSP(S,&ksp);CHKERRQ(ierr); 516 ierr = MatSchurComplementSetKSP(sub_schurs->S_Ej[i],ksp);CHKERRQ(ierr); 517 } else { /* build new ksp object which inherits ksp and pc types from the original one */ 518 KSP origksp,schurksp; 519 PC origpc,schurpc; 520 KSPType ksp_type; 521 PCType pc_type; 522 PetscInt n_internal; 523 524 ierr = MatSchurComplementGetKSP(S,&origksp);CHKERRQ(ierr); 525 ierr = MatSchurComplementGetKSP(sub_schurs->S_Ej[i],&schurksp);CHKERRQ(ierr); 526 ierr = KSPGetType(origksp,&ksp_type);CHKERRQ(ierr); 527 ierr = KSPSetType(schurksp,ksp_type);CHKERRQ(ierr); 528 ierr = KSPGetPC(schurksp,&schurpc);CHKERRQ(ierr); 529 ierr = KSPGetPC(origksp,&origpc);CHKERRQ(ierr); 530 ierr = PCGetType(origpc,&pc_type);CHKERRQ(ierr); 531 ierr = PCSetType(schurpc,pc_type);CHKERRQ(ierr); 532 ierr = ISGetSize(is_I,&n_internal);CHKERRQ(ierr); 533 if (n_internal) { /* UMFPACK gives error with 0 sized problems */ 534 MatSolverPackage solver=NULL; 535 ierr = PCFactorGetMatSolverPackage(origpc,(const MatSolverPackage*)&solver);CHKERRQ(ierr); 536 if (solver) { 537 ierr = PCFactorSetMatSolverPackage(schurpc,solver);CHKERRQ(ierr); 538 } 539 } 540 ierr = KSPSetUp(schurksp);CHKERRQ(ierr); 541 } 542 } 543 /* free */ 544 ierr = ISDestroy(&is_I);CHKERRQ(ierr); 545 ierr = MatDestroy(&AE_II);CHKERRQ(ierr); 546 for (i=0;i<sub_schurs->n_subs;i++) { 547 ierr = MatDestroy(&AE_EE[i]);CHKERRQ(ierr); 548 ierr = MatDestroy(&AE_IE[i]);CHKERRQ(ierr); 549 ierr = MatDestroy(&AE_EI[i]);CHKERRQ(ierr); 550 ierr = ISDestroy(&is_subset_B[i]);CHKERRQ(ierr); 551 } 552 ierr = PetscFree4(is_subset_B,AE_IE,AE_EI,AE_EE);CHKERRQ(ierr); 553 PetscFunctionReturn(0); 554 } 555