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