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 Mat global_schur_subsets,*submat_global_schur_subsets,work_mat; 13 ISLocalToGlobalMapping l2gmap_subsets; 14 IS is_I,*is_subset_B,temp_is; 15 PetscInt *nnz,*all_local_idx_G,*all_local_idx_B,*all_local_idx_N,*all_permutation_G; 16 PetscInt i,subset_size,max_subset_size; 17 PetscInt local_size,global_size; 18 PetscBool implicit_schurs; 19 PetscErrorCode ierr; 20 21 PetscFunctionBegin; 22 //ierr = PetscObjectTypeCompare((PetscObject)sub_schurs->S,MATSCHURCOMPLEMENT,&implicit_schurs);CHKERRQ(ierr); 23 implicit_schurs = PETSC_TRUE; 24 /* allocate space for schur complements */ 25 ierr = PetscMalloc4(sub_schurs->n_subs,&sub_schurs->is_AEj_B, 26 sub_schurs->n_subs,&sub_schurs->S_Ej, 27 sub_schurs->n_subs,&sub_schurs->work1, 28 sub_schurs->n_subs,&sub_schurs->work2);CHKERRQ(ierr); 29 30 /* get Schur complement matrices */ 31 if (implicit_schurs) { 32 ierr = MatSchurComplementGetSubMatrices(sub_schurs->S,&A_II,NULL,&A_IB,&A_BI,&A_BB);CHKERRQ(ierr); 33 ierr = PetscMalloc4(sub_schurs->n_subs,&is_subset_B, 34 sub_schurs->n_subs,&AE_IE, 35 sub_schurs->n_subs,&AE_EI, 36 sub_schurs->n_subs,&AE_EE);CHKERRQ(ierr); 37 } 38 39 /* determine interior problems */ 40 if (nlayers >= 0 && xadj != NULL && adjncy != NULL) { /* Interior problems can be different from the original one */ 41 PetscBT touched; 42 const PetscInt* idx_B; 43 PetscInt n_I,n_B,n_local_dofs,n_prev_added,j,layer,*local_numbering; 44 45 /* get sizes */ 46 ierr = ISGetLocalSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr); 47 ierr = ISGetLocalSize(sub_schurs->is_B,&n_B);CHKERRQ(ierr); 48 49 ierr = PetscMalloc1(n_I+n_B,&local_numbering);CHKERRQ(ierr); 50 ierr = PetscBTCreate(n_I+n_B,&touched);CHKERRQ(ierr); 51 ierr = PetscBTMemzero(n_I+n_B,touched);CHKERRQ(ierr); 52 53 /* all boundary dofs must be skipped when adding layers */ 54 ierr = ISGetIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr); 55 for (j=0;j<n_B;j++) { 56 ierr = PetscBTSet(touched,idx_B[j]);CHKERRQ(ierr); 57 } 58 ierr = PetscMemcpy(local_numbering,idx_B,n_B*sizeof(PetscInt));CHKERRQ(ierr); 59 ierr = ISRestoreIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr); 60 61 /* add prescribed number of layers of dofs */ 62 n_local_dofs = n_B; 63 n_prev_added = n_B; 64 for (layer=0;layer<nlayers;layer++) { 65 PetscInt n_added; 66 if (n_local_dofs == n_I+n_B) break; 67 if (n_local_dofs > n_I+n_B) { 68 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); 69 } 70 ierr = PCBDDCAdjGetNextLayer_Private(local_numbering+n_local_dofs,n_prev_added,touched,xadj,adjncy,&n_added);CHKERRQ(ierr); 71 n_prev_added = n_added; 72 n_local_dofs += n_added; 73 if (!n_added) break; 74 } 75 ierr = PetscBTDestroy(&touched);CHKERRQ(ierr); 76 77 /* IS for I dofs in original numbering */ 78 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sub_schurs->is_I),n_local_dofs-n_B,local_numbering+n_B,PETSC_COPY_VALUES,&sub_schurs->is_I_layer);CHKERRQ(ierr); 79 ierr = PetscFree(local_numbering);CHKERRQ(ierr); 80 ierr = ISSort(sub_schurs->is_I_layer);CHKERRQ(ierr); 81 /* IS for I dofs in boundary numbering */ 82 if (implicit_schurs) { 83 ISLocalToGlobalMapping ItoNmap; 84 ierr = ISLocalToGlobalMappingCreateIS(sub_schurs->is_I,&ItoNmap);CHKERRQ(ierr); 85 ierr = ISGlobalToLocalMappingApplyIS(ItoNmap,IS_GTOLM_DROP,sub_schurs->is_I_layer,&is_I);CHKERRQ(ierr); 86 ierr = ISLocalToGlobalMappingDestroy(&ItoNmap);CHKERRQ(ierr); 87 88 /* II block */ 89 ierr = MatGetSubMatrix(A_II,is_I,is_I,MAT_INITIAL_MATRIX,&AE_II);CHKERRQ(ierr); 90 } 91 } else { 92 PetscInt n_I; 93 94 /* IS for I dofs in original numbering */ 95 ierr = PetscObjectReference((PetscObject)sub_schurs->is_I);CHKERRQ(ierr); 96 sub_schurs->is_I_layer = sub_schurs->is_I; 97 98 /* IS for I dofs in I numbering (strided 1) */ 99 if (implicit_schurs) { 100 ierr = ISGetSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr); 101 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub_schurs->is_I),n_I,0,1,&is_I);CHKERRQ(ierr); 102 103 /* II block is the same */ 104 ierr = PetscObjectReference((PetscObject)A_II);CHKERRQ(ierr); 105 AE_II = A_II; 106 } 107 } 108 109 if (implicit_schurs) { 110 /* subsets in original and boundary numbering */ 111 for (i=0;i<sub_schurs->n_subs;i++) { 112 ierr = ISDuplicate(sub_schurs->is_subs[i],&sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 113 ierr = ISSort(sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 114 ierr = ISGlobalToLocalMappingApplyIS(sub_schurs->BtoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_B[i],&is_subset_B[i]);CHKERRQ(ierr); 115 } 116 117 /* EE block */ 118 for (i=0;i<sub_schurs->n_subs;i++) { 119 ierr = MatGetSubMatrix(A_BB,is_subset_B[i],is_subset_B[i],MAT_INITIAL_MATRIX,&AE_EE[i]);CHKERRQ(ierr); 120 } 121 /* IE block */ 122 for (i=0;i<sub_schurs->n_subs;i++) { 123 ierr = MatGetSubMatrix(A_IB,is_I,is_subset_B[i],MAT_INITIAL_MATRIX,&AE_IE[i]);CHKERRQ(ierr); 124 } 125 /* EI block */ 126 for (i=0;i<sub_schurs->n_subs;i++) { 127 ierr = MatGetSubMatrix(A_BI,is_subset_B[i],is_I,MAT_INITIAL_MATRIX,&AE_EI[i]);CHKERRQ(ierr); 128 } 129 130 /* setup Schur complements on subset */ 131 for (i=0;i<sub_schurs->n_subs;i++) { 132 ierr = MatCreateSchurComplement(AE_II,AE_II,AE_IE[i],AE_EI[i],AE_EE[i],&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 133 if (AE_II == A_II) { /* we can reuse the same ksp */ 134 KSP ksp; 135 ierr = MatSchurComplementGetKSP(sub_schurs->S,&ksp);CHKERRQ(ierr); 136 ierr = MatSchurComplementSetKSP(sub_schurs->S_Ej[i],ksp);CHKERRQ(ierr); 137 } else { /* build new ksp object which inherits ksp and pc types from the original one */ 138 KSP origksp,schurksp; 139 PC origpc,schurpc; 140 KSPType ksp_type; 141 PCType pc_type; 142 PetscInt n_internal; 143 144 ierr = MatSchurComplementGetKSP(sub_schurs->S,&origksp);CHKERRQ(ierr); 145 ierr = MatSchurComplementGetKSP(sub_schurs->S_Ej[i],&schurksp);CHKERRQ(ierr); 146 ierr = KSPGetType(origksp,&ksp_type);CHKERRQ(ierr); 147 ierr = KSPSetType(schurksp,ksp_type);CHKERRQ(ierr); 148 ierr = KSPGetPC(schurksp,&schurpc);CHKERRQ(ierr); 149 ierr = KSPGetPC(origksp,&origpc);CHKERRQ(ierr); 150 ierr = PCGetType(origpc,&pc_type);CHKERRQ(ierr); 151 ierr = PCSetType(schurpc,pc_type);CHKERRQ(ierr); 152 ierr = ISGetSize(is_I,&n_internal);CHKERRQ(ierr); 153 if (n_internal) { /* UMFPACK gives error with 0 sized problems */ 154 MatSolverPackage solver=NULL; 155 ierr = PCFactorGetMatSolverPackage(origpc,(const MatSolverPackage*)&solver);CHKERRQ(ierr); 156 if (solver) { 157 ierr = PCFactorSetMatSolverPackage(schurpc,solver);CHKERRQ(ierr); 158 } 159 } 160 ierr = KSPSetUp(schurksp);CHKERRQ(ierr); 161 } 162 } 163 /* free */ 164 ierr = ISDestroy(&is_I);CHKERRQ(ierr); 165 ierr = MatDestroy(&AE_II);CHKERRQ(ierr); 166 for (i=0;i<sub_schurs->n_subs;i++) { 167 ierr = MatDestroy(&AE_EE[i]);CHKERRQ(ierr); 168 ierr = MatDestroy(&AE_IE[i]);CHKERRQ(ierr); 169 ierr = MatDestroy(&AE_EI[i]);CHKERRQ(ierr); 170 ierr = ISDestroy(&is_subset_B[i]);CHKERRQ(ierr); 171 } 172 ierr = PetscFree4(is_subset_B,AE_IE,AE_EI,AE_EE);CHKERRQ(ierr); 173 } 174 175 /* TODO: just for compatibility with the previous version, needs to be fixed */ 176 for (i=0;i<sub_schurs->n_subs_par;i++) { 177 PetscInt j = sub_schurs->index_parallel[i]; 178 ierr = MatCreateVecs(sub_schurs->S_Ej[j],&sub_schurs->work1[j],&sub_schurs->work2[j]);CHKERRQ(ierr); 179 } 180 for (i=0;i<sub_schurs->n_subs_seq;i++) { 181 sub_schurs->work1[sub_schurs->index_sequential[i]] = 0; 182 sub_schurs->work2[sub_schurs->index_sequential[i]] = 0; 183 } 184 185 if (!sub_schurs->n_subs_seq_g) { 186 sub_schurs->S_Ej_all = 0; 187 sub_schurs->sum_S_Ej_all = 0; 188 sub_schurs->is_Ej_all = 0; 189 PetscFunctionReturn(0); 190 } 191 192 /* Get info on subset sizes and sum of all subsets sizes */ 193 max_subset_size = 0; 194 local_size = 0; 195 for (i=0;i<sub_schurs->n_subs_seq;i++) { 196 PetscInt j = sub_schurs->index_sequential[i]; 197 ierr = ISGetLocalSize(sub_schurs->is_AEj_B[j],&subset_size);CHKERRQ(ierr); 198 max_subset_size = PetscMax(subset_size,max_subset_size); 199 local_size += subset_size; 200 } 201 202 /* Work arrays for local indices */ 203 ierr = PetscMalloc1(local_size,&all_local_idx_B);CHKERRQ(ierr); 204 ierr = PetscMalloc1(local_size,&all_local_idx_N);CHKERRQ(ierr); 205 ierr = PetscMalloc1(local_size,&nnz);CHKERRQ(ierr); 206 207 /* Get local indices in local whole numbering and local boundary numbering */ 208 local_size = 0; 209 for (i=0;i<sub_schurs->n_subs_seq;i++) { 210 PetscInt j; 211 const PetscInt *idxs; 212 213 PetscInt local_problem_index = sub_schurs->index_sequential[i]; 214 ierr = ISGetLocalSize(sub_schurs->is_AEj_B[local_problem_index],&subset_size);CHKERRQ(ierr); 215 ierr = ISGetIndices(sub_schurs->is_AEj_B[local_problem_index],&idxs);CHKERRQ(ierr); 216 /* subset indices in local numbering */ 217 ierr = PetscMemcpy(all_local_idx_N+local_size,idxs,subset_size*sizeof(PetscInt));CHKERRQ(ierr); 218 ierr = ISRestoreIndices(sub_schurs->is_AEj_B[local_problem_index],&idxs);CHKERRQ(ierr); 219 for (j=0;j<subset_size;j++) nnz[local_size+j] = subset_size; 220 local_size += subset_size; 221 } 222 223 /* subset indices in local boundary numbering */ 224 ierr = ISGlobalToLocalMappingApply(sub_schurs->BtoNmap,IS_GTOLM_DROP,local_size,all_local_idx_N,&subset_size,all_local_idx_B);CHKERRQ(ierr); 225 if (subset_size != local_size) { 226 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in sub_schurs serial (BtoNmap)! %d != %d\n",subset_size,local_size); 227 } 228 ierr = ISCreateGeneral(PETSC_COMM_SELF,local_size,all_local_idx_B,PETSC_OWN_POINTER,&sub_schurs->is_Ej_all);CHKERRQ(ierr); 229 230 /* Number dofs on all subsets (parallel) and sort numbering */ 231 ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)sub_schurs->l2gmap),sub_schurs->l2gmap,local_size,all_local_idx_N,PETSC_NULL,&global_size,&all_local_idx_G);CHKERRQ(ierr); 232 ierr = PetscMalloc1(local_size,&all_permutation_G);CHKERRQ(ierr); 233 for (i=0;i<local_size;i++) { 234 all_permutation_G[i]=i; 235 } 236 ierr = PetscSortIntWithPermutation(local_size,all_local_idx_G,all_permutation_G);CHKERRQ(ierr); 237 238 /* Local matrix of all local Schur on subsets */ 239 ierr = MatCreate(PETSC_COMM_SELF,&sub_schurs->S_Ej_all);CHKERRQ(ierr); 240 ierr = MatSetSizes(sub_schurs->S_Ej_all,PETSC_DECIDE,PETSC_DECIDE,local_size,local_size);CHKERRQ(ierr); 241 ierr = MatSetType(sub_schurs->S_Ej_all,MATAIJ);CHKERRQ(ierr); 242 ierr = MatSeqAIJSetPreallocation(sub_schurs->S_Ej_all,0,nnz);CHKERRQ(ierr); 243 ierr = PetscFree(nnz);CHKERRQ(ierr); 244 245 if (implicit_schurs) { 246 PetscScalar *fill_vals; 247 PetscInt *dummy_idx; 248 249 /* Work arrays */ 250 ierr = PetscMalloc2(max_subset_size,&dummy_idx,max_subset_size*max_subset_size,&fill_vals);CHKERRQ(ierr); 251 252 /* Loop on local problems to compute Schur complements explicitly (TODO; optimize)*/ 253 local_size = 0; 254 for (i=0;i<sub_schurs->n_subs_seq;i++) { 255 Vec work1,work2; 256 PetscInt j,local_problem_index = sub_schurs->index_sequential[i]; 257 258 ierr = MatCreateVecs(sub_schurs->S_Ej[local_problem_index],&work1,&work2);CHKERRQ(ierr); 259 ierr = ISGetLocalSize(sub_schurs->is_AEj_B[local_problem_index],&subset_size);CHKERRQ(ierr); 260 /* local Schur */ 261 for (j=0;j<subset_size;j++) { 262 ierr = VecSet(work1,0.0);CHKERRQ(ierr); 263 ierr = VecSetValue(work1,j,1.0,INSERT_VALUES);CHKERRQ(ierr); 264 ierr = VecPlaceArray(work2,&fill_vals[j*subset_size]);CHKERRQ(ierr); 265 ierr = MatMult(sub_schurs->S_Ej[local_problem_index],work1,work2);CHKERRQ(ierr); 266 ierr = VecResetArray(work2);CHKERRQ(ierr); 267 } 268 for (j=0;j<subset_size;j++) { 269 dummy_idx[j]=local_size+j; 270 } 271 ierr = MatSetValues(sub_schurs->S_Ej_all,subset_size,dummy_idx,subset_size,dummy_idx,fill_vals,INSERT_VALUES);CHKERRQ(ierr); 272 local_size += subset_size; 273 ierr = VecDestroy(&work1);CHKERRQ(ierr); 274 ierr = VecDestroy(&work2);CHKERRQ(ierr); 275 } 276 ierr = PetscFree2(dummy_idx,fill_vals);CHKERRQ(ierr); 277 } 278 ierr = MatAssemblyBegin(sub_schurs->S_Ej_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 279 ierr = MatAssemblyEnd(sub_schurs->S_Ej_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 280 281 /* Global matrix of all assembled Schur on subsets */ 282 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)sub_schurs->l2gmap),1,local_size,all_local_idx_G,PETSC_COPY_VALUES,&l2gmap_subsets);CHKERRQ(ierr); 283 ierr = MatCreateIS(PetscObjectComm((PetscObject)sub_schurs->l2gmap),1,PETSC_DECIDE,PETSC_DECIDE,global_size,global_size,l2gmap_subsets,&work_mat);CHKERRQ(ierr); 284 ierr = ISLocalToGlobalMappingDestroy(&l2gmap_subsets);CHKERRQ(ierr); 285 ierr = MatISSetLocalMat(work_mat,sub_schurs->S_Ej_all);CHKERRQ(ierr); 286 ierr = MatISGetMPIXAIJ(work_mat,MAT_INITIAL_MATRIX,&global_schur_subsets);CHKERRQ(ierr); 287 ierr = MatDestroy(&work_mat);CHKERRQ(ierr); 288 289 /* Get local part of (\sum_j S_Ej) */ 290 for (i=0;i<local_size;i++) { 291 all_local_idx_N[i] = all_local_idx_G[all_permutation_G[i]]; 292 } 293 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sub_schurs->l2gmap),local_size,all_local_idx_N,PETSC_OWN_POINTER,&temp_is);CHKERRQ(ierr); 294 ierr = MatGetSubMatrices(global_schur_subsets,1,&temp_is,&temp_is,MAT_INITIAL_MATRIX,&submat_global_schur_subsets);CHKERRQ(ierr); 295 ierr = MatDestroy(&global_schur_subsets);CHKERRQ(ierr); 296 ierr = ISDestroy(&temp_is);CHKERRQ(ierr); 297 for (i=0;i<local_size;i++) { 298 all_local_idx_G[all_permutation_G[i]] = i; 299 } 300 ierr = ISCreateGeneral(PETSC_COMM_SELF,local_size,all_local_idx_G,PETSC_OWN_POINTER,&temp_is);CHKERRQ(ierr); 301 ierr = ISSetPermutation(temp_is);CHKERRQ(ierr); 302 ierr = MatPermute(submat_global_schur_subsets[0],temp_is,temp_is,&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr); 303 ierr = MatDestroyMatrices(1,&submat_global_schur_subsets);CHKERRQ(ierr); 304 ierr = ISDestroy(&temp_is);CHKERRQ(ierr); 305 ierr = PetscFree(all_permutation_G);CHKERRQ(ierr); 306 PetscFunctionReturn(0); 307 } 308 309 #undef __FUNCT__ 310 #define __FUNCT__ "PCBDDCSubSchursInit" 311 PetscErrorCode PCBDDCSubSchursInit(PCBDDCSubSchurs sub_schurs, Mat A, Mat S, IS is_I, IS is_B, PCBDDCGraph graph, ISLocalToGlobalMapping BtoNmap, PetscInt seqthreshold) 312 { 313 IS *faces,*edges,*all_cc; 314 PetscInt *index_sequential,*index_parallel; 315 PetscInt *auxlocal_sequential,*auxlocal_parallel; 316 PetscInt *auxglobal_sequential,*auxglobal_parallel; 317 PetscInt *auxmapping;//,*idxs; 318 PetscInt i,max_subset_size; 319 PetscInt n_sequential_problems,n_local_sequential_problems,n_parallel_problems,n_local_parallel_problems; 320 PetscInt n_faces,n_edges,n_all_cc; 321 PetscBool is_sorted; 322 PetscErrorCode ierr; 323 324 PetscFunctionBegin; 325 ierr = ISSorted(is_I,&is_sorted);CHKERRQ(ierr); 326 if (!is_sorted) { 327 SETERRQ(PetscObjectComm((PetscObject)is_I),PETSC_ERR_PLIB,"IS for I dofs should be shorted"); 328 } 329 ierr = ISSorted(is_B,&is_sorted);CHKERRQ(ierr); 330 if (!is_sorted) { 331 SETERRQ(PetscObjectComm((PetscObject)is_B),PETSC_ERR_PLIB,"IS for B dofs should be shorted"); 332 } 333 334 /* reset any previous data */ 335 ierr = PCBDDCSubSchursReset(sub_schurs);CHKERRQ(ierr); 336 337 /* get index sets for faces and edges */ 338 ierr = PCBDDCGraphGetCandidatesIS(graph,&n_faces,&faces,&n_edges,&edges,NULL);CHKERRQ(ierr); 339 n_all_cc = n_faces+n_edges; 340 ierr = PetscMalloc1(n_all_cc,&all_cc);CHKERRQ(ierr); 341 for (i=0;i<n_faces;i++) { 342 all_cc[i] = faces[i]; 343 } 344 for (i=0;i<n_edges;i++) { 345 all_cc[n_faces+i] = edges[i]; 346 } 347 ierr = PetscFree(faces);CHKERRQ(ierr); 348 ierr = PetscFree(edges);CHKERRQ(ierr); 349 350 /* map interface's subsets */ 351 max_subset_size = 0; 352 for (i=0;i<n_all_cc;i++) { 353 PetscInt subset_size; 354 ierr = ISGetLocalSize(all_cc[i],&subset_size);CHKERRQ(ierr); 355 max_subset_size = PetscMax(max_subset_size,subset_size); 356 } 357 ierr = PetscMalloc1(max_subset_size,&auxmapping);CHKERRQ(ierr); 358 ierr = PetscMalloc2(graph->ncc,&auxlocal_sequential, 359 graph->ncc,&auxlocal_parallel);CHKERRQ(ierr); 360 ierr = PetscMalloc2(graph->ncc,&index_sequential, 361 graph->ncc,&index_parallel);CHKERRQ(ierr); 362 363 /* if threshold is negative, uses all sequential problems */ 364 if (seqthreshold < 0) seqthreshold = max_subset_size; 365 366 /* determine which problem has to be solved in parallel or sequentially */ 367 n_local_sequential_problems = 0; 368 n_local_parallel_problems = 0; 369 for (i=0;i<n_all_cc;i++) { 370 PetscInt subset_size,j,min_loc = 0; 371 const PetscInt *idxs; 372 373 ierr = ISGetLocalSize(all_cc[i],&subset_size);CHKERRQ(ierr); 374 ierr = ISGetIndices(all_cc[i],&idxs);CHKERRQ(ierr); 375 ierr = ISLocalToGlobalMappingApply(graph->l2gmap,subset_size,idxs,auxmapping);CHKERRQ(ierr); 376 for (j=1;j<subset_size;j++) { 377 if (auxmapping[j]<auxmapping[min_loc]) { 378 min_loc = j; 379 } 380 } 381 if (subset_size > seqthreshold) { 382 index_parallel[n_local_parallel_problems] = i; 383 auxlocal_parallel[n_local_parallel_problems] = idxs[min_loc]; 384 n_local_parallel_problems++; 385 } else { 386 index_sequential[n_local_sequential_problems] = i; 387 auxlocal_sequential[n_local_sequential_problems] = idxs[min_loc]; 388 n_local_sequential_problems++; 389 } 390 ierr = ISRestoreIndices(all_cc[i],&idxs);CHKERRQ(ierr); 391 } 392 393 /* Number parallel problems */ 394 auxglobal_parallel = 0; 395 ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)graph->l2gmap),graph->l2gmap,n_local_parallel_problems,auxlocal_parallel,PETSC_NULL,&n_parallel_problems,&auxglobal_parallel);CHKERRQ(ierr); 396 397 /* Number sequential problems */ 398 auxglobal_sequential = 0; 399 ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)graph->l2gmap),graph->l2gmap,n_local_sequential_problems,auxlocal_sequential,PETSC_NULL,&n_sequential_problems,&auxglobal_sequential);CHKERRQ(ierr); 400 401 /* update info in sub_schurs */ 402 if (A) { 403 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 404 sub_schurs->A = A; 405 } 406 ierr = PetscObjectReference((PetscObject)S);CHKERRQ(ierr); 407 sub_schurs->S = S; 408 ierr = PetscObjectReference((PetscObject)is_I);CHKERRQ(ierr); 409 sub_schurs->is_I = is_I; 410 ierr = PetscObjectReference((PetscObject)is_B);CHKERRQ(ierr); 411 sub_schurs->is_B = is_B; 412 ierr = PetscObjectReference((PetscObject)graph->l2gmap);CHKERRQ(ierr); 413 sub_schurs->l2gmap = graph->l2gmap; 414 ierr = PetscObjectReference((PetscObject)BtoNmap);CHKERRQ(ierr); 415 sub_schurs->BtoNmap = BtoNmap; 416 sub_schurs->n_subs_seq = n_local_sequential_problems; 417 sub_schurs->n_subs_par = n_local_parallel_problems; 418 sub_schurs->n_subs_seq_g = n_sequential_problems; 419 sub_schurs->n_subs_par_g = n_parallel_problems; 420 sub_schurs->n_subs = sub_schurs->n_subs_seq + sub_schurs->n_subs_par; 421 sub_schurs->is_subs = all_cc; 422 sub_schurs->index_sequential = index_sequential; 423 sub_schurs->index_parallel = index_parallel; 424 sub_schurs->auxglobal_sequential = auxglobal_sequential; 425 sub_schurs->auxglobal_parallel = auxglobal_parallel; 426 427 /* free workspace */ 428 ierr = PetscFree(auxmapping);CHKERRQ(ierr); 429 ierr = PetscFree2(auxlocal_sequential,auxlocal_parallel);CHKERRQ(ierr); 430 431 PetscFunctionReturn(0); 432 } 433 434 #undef __FUNCT__ 435 #define __FUNCT__ "PCBDDCSubSchursCreate" 436 PetscErrorCode PCBDDCSubSchursCreate(PCBDDCSubSchurs *sub_schurs) 437 { 438 PCBDDCSubSchurs schurs_ctx; 439 PetscErrorCode ierr; 440 441 PetscFunctionBegin; 442 ierr = PetscNew(&schurs_ctx);CHKERRQ(ierr); 443 schurs_ctx->n_subs = 0; 444 *sub_schurs = schurs_ctx; 445 PetscFunctionReturn(0); 446 } 447 448 #undef __FUNCT__ 449 #define __FUNCT__ "PCBDDCSubSchursDestroy" 450 PetscErrorCode PCBDDCSubSchursDestroy(PCBDDCSubSchurs *sub_schurs) 451 { 452 PetscErrorCode ierr; 453 454 PetscFunctionBegin; 455 ierr = PCBDDCSubSchursReset(*sub_schurs);CHKERRQ(ierr); 456 ierr = PetscFree(*sub_schurs);CHKERRQ(ierr); 457 PetscFunctionReturn(0); 458 } 459 460 #undef __FUNCT__ 461 #define __FUNCT__ "PCBDDCSubSchursReset" 462 PetscErrorCode PCBDDCSubSchursReset(PCBDDCSubSchurs sub_schurs) 463 { 464 PetscInt i; 465 PetscErrorCode ierr; 466 467 PetscFunctionBegin; 468 ierr = MatDestroy(&sub_schurs->A);CHKERRQ(ierr); 469 ierr = MatDestroy(&sub_schurs->S);CHKERRQ(ierr); 470 ierr = ISDestroy(&sub_schurs->is_I);CHKERRQ(ierr); 471 ierr = ISDestroy(&sub_schurs->is_B);CHKERRQ(ierr); 472 ierr = ISLocalToGlobalMappingDestroy(&sub_schurs->l2gmap);CHKERRQ(ierr); 473 ierr = ISLocalToGlobalMappingDestroy(&sub_schurs->BtoNmap);CHKERRQ(ierr); 474 ierr = MatDestroy(&sub_schurs->S_Ej_all);CHKERRQ(ierr); 475 ierr = MatDestroy(&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr); 476 ierr = ISDestroy(&sub_schurs->is_Ej_all);CHKERRQ(ierr); 477 for (i=0;i<sub_schurs->n_subs;i++) { 478 ierr = ISDestroy(&sub_schurs->is_subs[i]);CHKERRQ(ierr); 479 ierr = ISDestroy(&sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 480 ierr = MatDestroy(&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 481 ierr = VecDestroy(&sub_schurs->work1[i]);CHKERRQ(ierr); 482 ierr = VecDestroy(&sub_schurs->work2[i]);CHKERRQ(ierr); 483 } 484 ierr = ISDestroy(&sub_schurs->is_I_layer);CHKERRQ(ierr); 485 if (sub_schurs->n_subs) { 486 ierr = PetscFree(sub_schurs->is_subs);CHKERRQ(ierr); 487 ierr = PetscFree4(sub_schurs->is_AEj_B,sub_schurs->S_Ej,sub_schurs->work1,sub_schurs->work2);CHKERRQ(ierr); 488 ierr = PetscFree2(sub_schurs->index_sequential,sub_schurs->index_parallel);CHKERRQ(ierr); 489 ierr = PetscFree(sub_schurs->auxglobal_sequential);CHKERRQ(ierr); 490 ierr = PetscFree(sub_schurs->auxglobal_parallel);CHKERRQ(ierr); 491 } 492 sub_schurs->n_subs = 0; 493 PetscFunctionReturn(0); 494 } 495 496 #undef __FUNCT__ 497 #define __FUNCT__ "PCBDDCAdjGetNextLayer_Private" 498 PETSC_STATIC_INLINE PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt* queue_tip,PetscInt n_prev,PetscBT touched,PetscInt* xadj,PetscInt* adjncy,PetscInt* n_added) 499 { 500 PetscInt i,j,n; 501 PetscErrorCode ierr; 502 503 PetscFunctionBegin; 504 n = 0; 505 for (i=-n_prev;i<0;i++) { 506 PetscInt start_dof = queue_tip[i]; 507 for (j=xadj[start_dof];j<xadj[start_dof+1];j++) { 508 PetscInt dof = adjncy[j]; 509 if (!PetscBTLookup(touched,dof)) { 510 ierr = PetscBTSet(touched,dof);CHKERRQ(ierr); 511 queue_tip[n] = dof; 512 n++; 513 } 514 } 515 } 516 *n_added = n; 517 PetscFunctionReturn(0); 518 } 519 520 #undef __FUNCT__ 521 #define __FUNCT__ "PCBDDCSubSchursSetUp" 522 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) 523 { 524 Mat A_II,A_IB,A_BI,A_BB; 525 Mat AE_II,*AE_IE,*AE_EI,*AE_EE; 526 IS is_I,*is_subset_B; 527 ISLocalToGlobalMapping BtoNmap; 528 PetscInt i; 529 PetscBool is_sorted; 530 PetscErrorCode ierr; 531 532 PetscFunctionBegin; 533 ierr = ISSorted(is_A_I,&is_sorted);CHKERRQ(ierr); 534 if (!is_sorted) { 535 SETERRQ(PetscObjectComm((PetscObject)is_A_I),PETSC_ERR_PLIB,"IS for I dofs should be shorted"); 536 } 537 ierr = ISSorted(is_A_B,&is_sorted);CHKERRQ(ierr); 538 if (!is_sorted) { 539 SETERRQ(PetscObjectComm((PetscObject)is_A_B),PETSC_ERR_PLIB,"IS for B dofs should be shorted"); 540 } 541 542 /* get Schur complement matrices */ 543 ierr = MatSchurComplementGetSubMatrices(S,&A_II,NULL,&A_IB,&A_BI,&A_BB);CHKERRQ(ierr); 544 545 /* allocate space for schur complements */ 546 sub_schurs->n_subs = ncc; 547 ierr = PetscMalloc4(sub_schurs->n_subs,&sub_schurs->is_AEj_B, 548 sub_schurs->n_subs,&sub_schurs->S_Ej, 549 sub_schurs->n_subs,&sub_schurs->work1, 550 sub_schurs->n_subs,&sub_schurs->work2);CHKERRQ(ierr); 551 ierr = PetscMalloc4(ncc,&is_subset_B,ncc,&AE_IE,ncc,&AE_EI,ncc,&AE_EE);CHKERRQ(ierr); 552 553 /* maps */ 554 if (sub_schurs->n_subs && nlayers >= 0 && xadj != NULL && adjncy != NULL) { /* Interior problems can be different from the original one */ 555 ISLocalToGlobalMapping ItoNmap; 556 PetscBT touched; 557 const PetscInt* idx_B; 558 PetscInt n_I,n_B,n_local_dofs,n_prev_added,j,layer,*local_numbering; 559 560 /* get sizes */ 561 ierr = ISGetLocalSize(is_A_I,&n_I);CHKERRQ(ierr); 562 ierr = ISGetLocalSize(is_A_B,&n_B);CHKERRQ(ierr); 563 564 ierr = ISLocalToGlobalMappingCreateIS(is_A_I,&ItoNmap);CHKERRQ(ierr); 565 ierr = PetscMalloc1(n_I+n_B,&local_numbering);CHKERRQ(ierr); 566 ierr = PetscBTCreate(n_I+n_B,&touched);CHKERRQ(ierr); 567 ierr = PetscBTMemzero(n_I+n_B,touched);CHKERRQ(ierr); 568 569 /* all boundary dofs must be skipped when adding layers */ 570 ierr = ISGetIndices(is_A_B,&idx_B);CHKERRQ(ierr); 571 for (j=0;j<n_B;j++) { 572 ierr = PetscBTSet(touched,idx_B[j]);CHKERRQ(ierr); 573 } 574 ierr = PetscMemcpy(local_numbering,idx_B,n_B*sizeof(PetscInt));CHKERRQ(ierr); 575 ierr = ISRestoreIndices(is_A_B,&idx_B);CHKERRQ(ierr); 576 577 /* add next layers of dofs */ 578 n_local_dofs = n_B; 579 n_prev_added = n_B; 580 for (layer=0;layer<nlayers;layer++) { 581 PetscInt n_added; 582 if (n_local_dofs == n_I+n_B) break; 583 if (n_local_dofs > n_I+n_B) { 584 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); 585 } 586 ierr = PCBDDCAdjGetNextLayer_Private(local_numbering+n_local_dofs,n_prev_added,touched,xadj,adjncy,&n_added);CHKERRQ(ierr); 587 n_prev_added = n_added; 588 n_local_dofs += n_added; 589 if (!n_added) break; 590 } 591 ierr = PetscBTDestroy(&touched);CHKERRQ(ierr); 592 593 /* IS for I dofs in original numbering and in I numbering */ 594 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)ItoNmap),n_local_dofs-n_B,local_numbering+n_B,PETSC_COPY_VALUES,&sub_schurs->is_I_layer);CHKERRQ(ierr); 595 ierr = PetscFree(local_numbering);CHKERRQ(ierr); 596 ierr = ISSort(sub_schurs->is_I_layer);CHKERRQ(ierr); 597 ierr = ISGlobalToLocalMappingApplyIS(ItoNmap,IS_GTOLM_DROP,sub_schurs->is_I_layer,&is_I);CHKERRQ(ierr); 598 ierr = ISLocalToGlobalMappingDestroy(&ItoNmap);CHKERRQ(ierr); 599 600 /* II block */ 601 ierr = MatGetSubMatrix(A_II,is_I,is_I,MAT_INITIAL_MATRIX,&AE_II);CHKERRQ(ierr); 602 } else { 603 PetscInt n_I; 604 605 /* IS for I dofs in original numbering */ 606 ierr = PetscObjectReference((PetscObject)is_A_I);CHKERRQ(ierr); 607 sub_schurs->is_I_layer = is_A_I; 608 609 /* IS for I dofs in I numbering (strided 1) */ 610 ierr = ISGetSize(is_A_I,&n_I);CHKERRQ(ierr); 611 ierr = ISCreateStride(PetscObjectComm((PetscObject)is_A_I),n_I,0,1,&is_I);CHKERRQ(ierr); 612 613 /* II block is the same */ 614 ierr = PetscObjectReference((PetscObject)A_II);CHKERRQ(ierr); 615 AE_II = A_II; 616 } 617 618 /* subsets in original and boundary numbering */ 619 ierr = ISLocalToGlobalMappingCreateIS(is_A_B,&BtoNmap);CHKERRQ(ierr); 620 for (i=0;i<sub_schurs->n_subs;i++) { 621 ierr = ISDuplicate(is_cc[i],&sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 622 ierr = ISSort(sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 623 ierr = ISGlobalToLocalMappingApplyIS(BtoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_B[i],&is_subset_B[i]);CHKERRQ(ierr); 624 } 625 ierr = ISLocalToGlobalMappingDestroy(&BtoNmap);CHKERRQ(ierr); 626 627 /* EE block */ 628 for (i=0;i<sub_schurs->n_subs;i++) { 629 ierr = MatGetSubMatrix(A_BB,is_subset_B[i],is_subset_B[i],MAT_INITIAL_MATRIX,&AE_EE[i]);CHKERRQ(ierr); 630 } 631 /* IE block */ 632 for (i=0;i<sub_schurs->n_subs;i++) { 633 ierr = MatGetSubMatrix(A_IB,is_I,is_subset_B[i],MAT_INITIAL_MATRIX,&AE_IE[i]);CHKERRQ(ierr); 634 } 635 /* EI block */ 636 for (i=0;i<sub_schurs->n_subs;i++) { 637 ierr = MatGetSubMatrix(A_BI,is_subset_B[i],is_I,MAT_INITIAL_MATRIX,&AE_EI[i]);CHKERRQ(ierr); 638 } 639 640 /* setup Schur complements on subset */ 641 for (i=0;i<sub_schurs->n_subs;i++) { 642 ierr = MatCreateSchurComplement(AE_II,AE_II,AE_IE[i],AE_EI[i],AE_EE[i],&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 643 ierr = MatCreateVecs(sub_schurs->S_Ej[i],&sub_schurs->work1[i],&sub_schurs->work2[i]);CHKERRQ(ierr); 644 if (AE_II == A_II) { /* we can reuse the same ksp */ 645 KSP ksp; 646 ierr = MatSchurComplementGetKSP(S,&ksp);CHKERRQ(ierr); 647 ierr = MatSchurComplementSetKSP(sub_schurs->S_Ej[i],ksp);CHKERRQ(ierr); 648 } else { /* build new ksp object which inherits ksp and pc types from the original one */ 649 KSP origksp,schurksp; 650 PC origpc,schurpc; 651 KSPType ksp_type; 652 PCType pc_type; 653 PetscInt n_internal; 654 655 ierr = MatSchurComplementGetKSP(S,&origksp);CHKERRQ(ierr); 656 ierr = MatSchurComplementGetKSP(sub_schurs->S_Ej[i],&schurksp);CHKERRQ(ierr); 657 ierr = KSPGetType(origksp,&ksp_type);CHKERRQ(ierr); 658 ierr = KSPSetType(schurksp,ksp_type);CHKERRQ(ierr); 659 ierr = KSPGetPC(schurksp,&schurpc);CHKERRQ(ierr); 660 ierr = KSPGetPC(origksp,&origpc);CHKERRQ(ierr); 661 ierr = PCGetType(origpc,&pc_type);CHKERRQ(ierr); 662 ierr = PCSetType(schurpc,pc_type);CHKERRQ(ierr); 663 ierr = ISGetSize(is_I,&n_internal);CHKERRQ(ierr); 664 if (n_internal) { /* UMFPACK gives error with 0 sized problems */ 665 MatSolverPackage solver=NULL; 666 ierr = PCFactorGetMatSolverPackage(origpc,(const MatSolverPackage*)&solver);CHKERRQ(ierr); 667 if (solver) { 668 ierr = PCFactorSetMatSolverPackage(schurpc,solver);CHKERRQ(ierr); 669 } 670 } 671 ierr = KSPSetUp(schurksp);CHKERRQ(ierr); 672 } 673 } 674 /* free */ 675 ierr = ISDestroy(&is_I);CHKERRQ(ierr); 676 ierr = MatDestroy(&AE_II);CHKERRQ(ierr); 677 for (i=0;i<sub_schurs->n_subs;i++) { 678 ierr = MatDestroy(&AE_EE[i]);CHKERRQ(ierr); 679 ierr = MatDestroy(&AE_IE[i]);CHKERRQ(ierr); 680 ierr = MatDestroy(&AE_EI[i]);CHKERRQ(ierr); 681 ierr = ISDestroy(&is_subset_B[i]);CHKERRQ(ierr); 682 } 683 ierr = PetscFree4(is_subset_B,AE_IE,AE_EI,AE_EE);CHKERRQ(ierr); 684 PetscFunctionReturn(0); 685 } 686