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 ierr = MatDestroy(&sub_schurs->S_Ej_all);CHKERRQ(ierr); 344 ierr = MatDestroy(&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr); 345 for (i=0;i<sub_schurs->n_subs;i++) { 346 ierr = ISDestroy(&sub_schurs->is_subs[i]);CHKERRQ(ierr); 347 ierr = ISDestroy(&sub_schurs->is_AEj_I[i]);CHKERRQ(ierr); 348 ierr = ISDestroy(&sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 349 ierr = MatDestroy(&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 350 ierr = VecDestroy(&sub_schurs->work1[i]);CHKERRQ(ierr); 351 ierr = VecDestroy(&sub_schurs->work2[i]);CHKERRQ(ierr); 352 } 353 if (sub_schurs->n_subs) { 354 ierr = PetscFree(sub_schurs->is_subs);CHKERRQ(ierr); 355 ierr = PetscFree5(sub_schurs->is_AEj_I,sub_schurs->is_AEj_B,sub_schurs->S_Ej,sub_schurs->work1,sub_schurs->work2);CHKERRQ(ierr); 356 ierr = PetscFree2(sub_schurs->index_sequential,sub_schurs->index_parallel);CHKERRQ(ierr); 357 ierr = PetscFree(sub_schurs->auxglobal_sequential);CHKERRQ(ierr); 358 ierr = PetscFree(sub_schurs->auxglobal_parallel);CHKERRQ(ierr); 359 } 360 sub_schurs->n_subs = 0; 361 PetscFunctionReturn(0); 362 } 363 364 #undef __FUNCT__ 365 #define __FUNCT__ "PCBDDCAdjGetNextLayer_Private" 366 PETSC_STATIC_INLINE PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt* queue_tip,PetscInt n_prev,PetscBT touched,PetscInt* xadj,PetscInt* adjncy,PetscInt* n_added) 367 { 368 PetscInt i,j,n; 369 PetscErrorCode ierr; 370 371 PetscFunctionBegin; 372 n = 0; 373 for (i=-n_prev;i<0;i++) { 374 PetscInt start_dof = queue_tip[i]; 375 for (j=xadj[start_dof];j<xadj[start_dof+1];j++) { 376 PetscInt dof = adjncy[j]; 377 if (!PetscBTLookup(touched,dof)) { 378 ierr = PetscBTSet(touched,dof);CHKERRQ(ierr); 379 queue_tip[n] = dof; 380 n++; 381 } 382 } 383 } 384 *n_added = n; 385 PetscFunctionReturn(0); 386 } 387 388 #undef __FUNCT__ 389 #define __FUNCT__ "PCBDDCSubSchursSetUp" 390 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) 391 { 392 Mat A_II,A_IB,A_BI,A_BB; 393 Mat AE_II,*AE_IE,*AE_EI,*AE_EE; 394 IS is_I,*is_subset_B; 395 ISLocalToGlobalMapping BtoNmap; 396 PetscInt i; 397 PetscBool is_sorted; 398 PetscErrorCode ierr; 399 400 PetscFunctionBegin; 401 ierr = ISSorted(is_A_I,&is_sorted);CHKERRQ(ierr); 402 if (!is_sorted) { 403 SETERRQ(PetscObjectComm((PetscObject)is_A_I),PETSC_ERR_PLIB,"IS for I dofs should be shorted"); 404 } 405 ierr = ISSorted(is_A_B,&is_sorted);CHKERRQ(ierr); 406 if (!is_sorted) { 407 SETERRQ(PetscObjectComm((PetscObject)is_A_B),PETSC_ERR_PLIB,"IS for B dofs should be shorted"); 408 } 409 410 /* get Schur complement matrices */ 411 ierr = MatSchurComplementGetSubMatrices(S,&A_II,NULL,&A_IB,&A_BI,&A_BB);CHKERRQ(ierr); 412 413 /* allocate space for schur complements */ 414 ierr = PetscMalloc5(sub_schurs->n_subs,&sub_schurs->is_AEj_I, 415 sub_schurs->n_subs,&sub_schurs->is_AEj_B, 416 sub_schurs->n_subs,&sub_schurs->S_Ej, 417 sub_schurs->n_subs,&sub_schurs->work1, 418 sub_schurs->n_subs,&sub_schurs->work2);CHKERRQ(ierr); 419 ierr = PetscMalloc4(ncc,&is_subset_B,ncc,&AE_IE,ncc,&AE_EI,ncc,&AE_EE);CHKERRQ(ierr); 420 sub_schurs->n_subs = ncc; 421 422 /* maps */ 423 if (sub_schurs->n_subs && nlayers >= 0 && xadj != NULL && adjncy != NULL) { /* Interior problems can be different from the original one */ 424 ISLocalToGlobalMapping ItoNmap; 425 PetscBT touched; 426 const PetscInt* idx_B; 427 PetscInt n_I,n_B,n_local_dofs,n_prev_added,j,layer,*local_numbering; 428 429 /* get sizes */ 430 ierr = ISGetLocalSize(is_A_I,&n_I);CHKERRQ(ierr); 431 ierr = ISGetLocalSize(is_A_B,&n_B);CHKERRQ(ierr); 432 433 ierr = ISLocalToGlobalMappingCreateIS(is_A_I,&ItoNmap);CHKERRQ(ierr); 434 ierr = PetscMalloc1(n_I+n_B,&local_numbering);CHKERRQ(ierr); 435 ierr = PetscBTCreate(n_I+n_B,&touched);CHKERRQ(ierr); 436 ierr = PetscBTMemzero(n_I+n_B,touched);CHKERRQ(ierr); 437 438 /* all boundary dofs must be skipped when adding layers */ 439 ierr = ISGetIndices(is_A_B,&idx_B);CHKERRQ(ierr); 440 for (j=0;j<n_B;j++) { 441 ierr = PetscBTSet(touched,idx_B[j]);CHKERRQ(ierr); 442 } 443 ierr = PetscMemcpy(local_numbering,idx_B,n_B*sizeof(PetscInt));CHKERRQ(ierr); 444 ierr = ISRestoreIndices(is_A_B,&idx_B);CHKERRQ(ierr); 445 446 /* add next layers of dofs */ 447 n_local_dofs = n_B; 448 n_prev_added = n_B; 449 for (layer=0;layer<nlayers;layer++) { 450 PetscInt n_added; 451 if (n_local_dofs == n_I+n_B) break; 452 if (n_local_dofs > n_I+n_B) { 453 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); 454 } 455 ierr = PCBDDCAdjGetNextLayer_Private(local_numbering+n_local_dofs,n_prev_added,touched,xadj,adjncy,&n_added);CHKERRQ(ierr); 456 n_prev_added = n_added; 457 n_local_dofs += n_added; 458 if (!n_added) break; 459 } 460 ierr = PetscBTDestroy(&touched);CHKERRQ(ierr); 461 462 /* IS for I dofs in original numbering and in I numbering */ 463 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); 464 ierr = PetscFree(local_numbering);CHKERRQ(ierr); 465 ierr = ISSort(sub_schurs->is_AEj_I[0]);CHKERRQ(ierr); 466 ierr = ISGlobalToLocalMappingApplyIS(ItoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_I[0],&is_I);CHKERRQ(ierr); 467 ierr = ISLocalToGlobalMappingDestroy(&ItoNmap);CHKERRQ(ierr); 468 469 /* II block */ 470 ierr = MatGetSubMatrix(A_II,is_I,is_I,MAT_INITIAL_MATRIX,&AE_II);CHKERRQ(ierr); 471 } else { 472 PetscInt n_I; 473 474 /* IS for I dofs in original numbering */ 475 ierr = PetscObjectReference((PetscObject)is_A_I);CHKERRQ(ierr); 476 sub_schurs->is_AEj_I[0] = is_A_I; 477 478 /* IS for I dofs in I numbering (strided 1) */ 479 ierr = ISGetSize(is_A_I,&n_I);CHKERRQ(ierr); 480 ierr = ISCreateStride(PetscObjectComm((PetscObject)is_A_I),n_I,0,1,&is_I);CHKERRQ(ierr); 481 482 /* II block is the same */ 483 ierr = PetscObjectReference((PetscObject)A_II);CHKERRQ(ierr); 484 AE_II = A_II; 485 } 486 487 /* TODO: just for compatibility with the previous version, needs to be fixed */ 488 for (i=1;i<sub_schurs->n_subs;i++) { 489 ierr = PetscObjectReference((PetscObject)sub_schurs->is_AEj_I[0]);CHKERRQ(ierr); 490 sub_schurs->is_AEj_I[i] = sub_schurs->is_AEj_I[0]; 491 } 492 493 /* subsets in original and boundary numbering */ 494 ierr = ISLocalToGlobalMappingCreateIS(is_A_B,&BtoNmap);CHKERRQ(ierr); 495 for (i=0;i<sub_schurs->n_subs;i++) { 496 ierr = ISDuplicate(is_cc[i],&sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 497 ierr = ISSort(sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 498 ierr = ISGlobalToLocalMappingApplyIS(BtoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_B[i],&is_subset_B[i]);CHKERRQ(ierr); 499 } 500 ierr = ISLocalToGlobalMappingDestroy(&BtoNmap);CHKERRQ(ierr); 501 502 /* EE block */ 503 for (i=0;i<sub_schurs->n_subs;i++) { 504 ierr = MatGetSubMatrix(A_BB,is_subset_B[i],is_subset_B[i],MAT_INITIAL_MATRIX,&AE_EE[i]);CHKERRQ(ierr); 505 } 506 /* IE block */ 507 for (i=0;i<sub_schurs->n_subs;i++) { 508 ierr = MatGetSubMatrix(A_IB,is_I,is_subset_B[i],MAT_INITIAL_MATRIX,&AE_IE[i]);CHKERRQ(ierr); 509 } 510 /* EI block */ 511 for (i=0;i<sub_schurs->n_subs;i++) { 512 ierr = MatGetSubMatrix(A_BI,is_subset_B[i],is_I,MAT_INITIAL_MATRIX,&AE_EI[i]);CHKERRQ(ierr); 513 } 514 515 /* setup Schur complements on subset */ 516 for (i=0;i<sub_schurs->n_subs;i++) { 517 ierr = MatCreateSchurComplement(AE_II,AE_II,AE_IE[i],AE_EI[i],AE_EE[i],&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 518 ierr = MatCreateVecs(sub_schurs->S_Ej[i],&sub_schurs->work1[i],&sub_schurs->work2[i]);CHKERRQ(ierr); 519 if (AE_II == A_II) { /* we can reuse the same ksp */ 520 KSP ksp; 521 ierr = MatSchurComplementGetKSP(S,&ksp);CHKERRQ(ierr); 522 ierr = MatSchurComplementSetKSP(sub_schurs->S_Ej[i],ksp);CHKERRQ(ierr); 523 } else { /* build new ksp object which inherits ksp and pc types from the original one */ 524 KSP origksp,schurksp; 525 PC origpc,schurpc; 526 KSPType ksp_type; 527 PCType pc_type; 528 PetscInt n_internal; 529 530 ierr = MatSchurComplementGetKSP(S,&origksp);CHKERRQ(ierr); 531 ierr = MatSchurComplementGetKSP(sub_schurs->S_Ej[i],&schurksp);CHKERRQ(ierr); 532 ierr = KSPGetType(origksp,&ksp_type);CHKERRQ(ierr); 533 ierr = KSPSetType(schurksp,ksp_type);CHKERRQ(ierr); 534 ierr = KSPGetPC(schurksp,&schurpc);CHKERRQ(ierr); 535 ierr = KSPGetPC(origksp,&origpc);CHKERRQ(ierr); 536 ierr = PCGetType(origpc,&pc_type);CHKERRQ(ierr); 537 ierr = PCSetType(schurpc,pc_type);CHKERRQ(ierr); 538 ierr = ISGetSize(is_I,&n_internal);CHKERRQ(ierr); 539 if (n_internal) { /* UMFPACK gives error with 0 sized problems */ 540 MatSolverPackage solver=NULL; 541 ierr = PCFactorGetMatSolverPackage(origpc,(const MatSolverPackage*)&solver);CHKERRQ(ierr); 542 if (solver) { 543 ierr = PCFactorSetMatSolverPackage(schurpc,solver);CHKERRQ(ierr); 544 } 545 } 546 ierr = KSPSetUp(schurksp);CHKERRQ(ierr); 547 } 548 } 549 /* free */ 550 ierr = ISDestroy(&is_I);CHKERRQ(ierr); 551 ierr = MatDestroy(&AE_II);CHKERRQ(ierr); 552 for (i=0;i<sub_schurs->n_subs;i++) { 553 ierr = MatDestroy(&AE_EE[i]);CHKERRQ(ierr); 554 ierr = MatDestroy(&AE_IE[i]);CHKERRQ(ierr); 555 ierr = MatDestroy(&AE_EI[i]);CHKERRQ(ierr); 556 ierr = ISDestroy(&is_subset_B[i]);CHKERRQ(ierr); 557 } 558 ierr = PetscFree4(is_subset_B,AE_IE,AE_EI,AE_EE);CHKERRQ(ierr); 559 PetscFunctionReturn(0); 560 } 561