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__ "PCBDDCSubSchursCreate" 8 PetscErrorCode PCBDDCSubSchursCreate(PCBDDCSubSchurs *sub_schurs) 9 { 10 PCBDDCSubSchurs schurs_ctx; 11 PetscErrorCode ierr; 12 13 PetscFunctionBegin; 14 ierr = PetscNew(&schurs_ctx);CHKERRQ(ierr); 15 *sub_schurs = schurs_ctx; 16 PetscFunctionReturn(0); 17 } 18 19 #undef __FUNCT__ 20 #define __FUNCT__ "PCBDDCSubSchursDestroy" 21 PetscErrorCode PCBDDCSubSchursDestroy(PCBDDCSubSchurs *sub_schurs) 22 { 23 PetscErrorCode ierr; 24 25 PetscFunctionBegin; 26 ierr = PCBDDCSubSchursReset(*sub_schurs);CHKERRQ(ierr); 27 ierr = PetscFree(*sub_schurs);CHKERRQ(ierr); 28 PetscFunctionReturn(0); 29 } 30 31 #undef __FUNCT__ 32 #define __FUNCT__ "PCBDDCSubSchursReset" 33 PetscErrorCode PCBDDCSubSchursReset(PCBDDCSubSchurs sub_schurs) 34 { 35 PetscInt i; 36 PetscErrorCode ierr; 37 38 PetscFunctionBegin; 39 for (i=0;i<sub_schurs->n_subs;i++) { 40 ierr = ISDestroy(&sub_schurs->is_AEj_I[i]);CHKERRQ(ierr); 41 ierr = ISDestroy(&sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 42 ierr = MatDestroy(&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 43 ierr = VecDestroy(&sub_schurs->work1[i]);CHKERRQ(ierr); 44 ierr = VecDestroy(&sub_schurs->work2[i]);CHKERRQ(ierr); 45 } 46 ierr = PetscFree5(sub_schurs->is_AEj_I,sub_schurs->is_AEj_B,sub_schurs->S_Ej,sub_schurs->work1,sub_schurs->work2);CHKERRQ(ierr); 47 sub_schurs->n_subs = 0; 48 PetscFunctionReturn(0); 49 } 50 51 #undef __FUNCT__ 52 #define __FUNCT__ "PCBDDCAdjGetNextLayer_Private" 53 PETSC_STATIC_INLINE PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt* queue_tip,PetscInt n_prev,PetscBT touched,PetscInt* xadj,PetscInt* adjncy,PetscInt* n_added) 54 { 55 PetscInt i,j,n; 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 n = 0; 60 for (i=-n_prev;i<0;i++) { 61 PetscInt start_dof = queue_tip[i]; 62 for (j=xadj[start_dof];j<xadj[start_dof+1];j++) { 63 PetscInt dof = adjncy[j]; 64 if (!PetscBTLookup(touched,dof)) { 65 ierr = PetscBTSet(touched,dof);CHKERRQ(ierr); 66 queue_tip[n] = dof; 67 n++; 68 } 69 } 70 } 71 *n_added = n; 72 PetscFunctionReturn(0); 73 } 74 75 #undef __FUNCT__ 76 #define __FUNCT__ "PCBDDCSubSchursSetUp" 77 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) 78 { 79 Mat A_II,A_IB,A_BI,A_BB; 80 Mat AE_II,*AE_IE,*AE_EI,*AE_EE; 81 IS is_I,*is_subset_B; 82 ISLocalToGlobalMapping BtoNmap; 83 PetscInt i; 84 PetscBool is_sorted; 85 PetscErrorCode ierr; 86 87 PetscFunctionBegin; 88 ierr = ISSorted(is_A_I,&is_sorted);CHKERRQ(ierr); 89 if (!is_sorted) { 90 SETERRQ(PetscObjectComm((PetscObject)is_A_I),PETSC_ERR_PLIB,"IS for I dofs should be shorted"); 91 } 92 ierr = ISSorted(is_A_B,&is_sorted);CHKERRQ(ierr); 93 if (!is_sorted) { 94 SETERRQ(PetscObjectComm((PetscObject)is_A_B),PETSC_ERR_PLIB,"IS for B dofs should be shorted"); 95 } 96 97 /* get Schur complement matrices */ 98 ierr = MatSchurComplementGetSubMatrices(S,&A_II,NULL,&A_IB,&A_BI,&A_BB);CHKERRQ(ierr); 99 100 /* allocate space for schur complements */ 101 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); 102 ierr = PetscMalloc4(ncc,&is_subset_B,ncc,&AE_IE,ncc,&AE_EI,ncc,&AE_EE);CHKERRQ(ierr); 103 sub_schurs->n_subs = ncc; 104 105 /* maps */ 106 if (sub_schurs->n_subs && nlayers >= 0 && xadj != NULL && adjncy != NULL) { /* Interior problems can be different from the original one */ 107 ISLocalToGlobalMapping ItoNmap; 108 PetscBT touched; 109 const PetscInt* idx_B; 110 PetscInt n_I,n_B,n_local_dofs,n_prev_added,j,layer,*local_numbering; 111 112 /* get sizes */ 113 ierr = ISGetLocalSize(is_A_I,&n_I);CHKERRQ(ierr); 114 ierr = ISGetLocalSize(is_A_B,&n_B);CHKERRQ(ierr); 115 116 ierr = ISLocalToGlobalMappingCreateIS(is_A_I,&ItoNmap);CHKERRQ(ierr); 117 ierr = PetscMalloc1(n_I+n_B,&local_numbering);CHKERRQ(ierr); 118 ierr = PetscBTCreate(n_I+n_B,&touched);CHKERRQ(ierr); 119 ierr = PetscBTMemzero(n_I+n_B,touched);CHKERRQ(ierr); 120 121 /* all boundary dofs must be skipped when adding layers */ 122 ierr = ISGetIndices(is_A_B,&idx_B);CHKERRQ(ierr); 123 for (j=0;j<n_B;j++) { 124 ierr = PetscBTSet(touched,idx_B[j]);CHKERRQ(ierr); 125 } 126 ierr = PetscMemcpy(local_numbering,idx_B,n_B*sizeof(PetscInt));CHKERRQ(ierr); 127 ierr = ISRestoreIndices(is_A_B,&idx_B);CHKERRQ(ierr); 128 129 /* add next layers of dofs */ 130 n_local_dofs = n_B; 131 n_prev_added = n_B; 132 for (layer=0;layer<nlayers;layer++) { 133 PetscInt n_added; 134 if (n_local_dofs == n_I+n_B) break; 135 if (n_local_dofs > n_I+n_B) { 136 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); 137 } 138 ierr = PCBDDCAdjGetNextLayer_Private(local_numbering+n_local_dofs,n_prev_added,touched,xadj,adjncy,&n_added);CHKERRQ(ierr); 139 n_prev_added = n_added; 140 n_local_dofs += n_added; 141 if (!n_added) break; 142 } 143 ierr = PetscBTDestroy(&touched);CHKERRQ(ierr); 144 145 /* IS for I dofs in original numbering and in I numbering */ 146 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); 147 ierr = PetscFree(local_numbering);CHKERRQ(ierr); 148 ierr = ISSort(sub_schurs->is_AEj_I[0]);CHKERRQ(ierr); 149 ierr = ISGlobalToLocalMappingApplyIS(ItoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_I[0],&is_I);CHKERRQ(ierr); 150 ierr = ISLocalToGlobalMappingDestroy(&ItoNmap);CHKERRQ(ierr); 151 152 /* II block */ 153 ierr = MatGetSubMatrix(A_II,is_I,is_I,MAT_INITIAL_MATRIX,&AE_II);CHKERRQ(ierr); 154 } else { 155 PetscInt n_I; 156 157 /* IS for I dofs in original numbering */ 158 ierr = PetscObjectReference((PetscObject)is_A_I);CHKERRQ(ierr); 159 sub_schurs->is_AEj_I[0] = is_A_I; 160 161 /* IS for I dofs in I numbering (strided 1) */ 162 ierr = ISGetSize(is_A_I,&n_I);CHKERRQ(ierr); 163 ierr = ISCreateStride(PetscObjectComm((PetscObject)is_A_I),n_I,0,1,&is_I);CHKERRQ(ierr); 164 165 /* II block is the same */ 166 ierr = PetscObjectReference((PetscObject)A_II);CHKERRQ(ierr); 167 AE_II = A_II; 168 } 169 170 /* TODO: just for compatibility with the previous version, needs to be fixed */ 171 for (i=1;i<sub_schurs->n_subs;i++) { 172 ierr = PetscObjectReference((PetscObject)sub_schurs->is_AEj_I[0]);CHKERRQ(ierr); 173 sub_schurs->is_AEj_I[i] = sub_schurs->is_AEj_I[0]; 174 } 175 176 /* subsets in original and boundary numbering */ 177 ierr = ISLocalToGlobalMappingCreateIS(is_A_B,&BtoNmap);CHKERRQ(ierr); 178 for (i=0;i<sub_schurs->n_subs;i++) { 179 ierr = ISDuplicate(is_cc[i],&sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 180 ierr = ISSort(sub_schurs->is_AEj_B[i]);CHKERRQ(ierr); 181 ierr = ISGlobalToLocalMappingApplyIS(BtoNmap,IS_GTOLM_DROP,sub_schurs->is_AEj_B[i],&is_subset_B[i]);CHKERRQ(ierr); 182 } 183 ierr = ISLocalToGlobalMappingDestroy(&BtoNmap);CHKERRQ(ierr); 184 185 /* EE block */ 186 for (i=0;i<sub_schurs->n_subs;i++) { 187 ierr = MatGetSubMatrix(A_BB,is_subset_B[i],is_subset_B[i],MAT_INITIAL_MATRIX,&AE_EE[i]);CHKERRQ(ierr); 188 } 189 /* IE block */ 190 for (i=0;i<sub_schurs->n_subs;i++) { 191 ierr = MatGetSubMatrix(A_IB,is_I,is_subset_B[i],MAT_INITIAL_MATRIX,&AE_IE[i]);CHKERRQ(ierr); 192 } 193 /* EI block */ 194 for (i=0;i<sub_schurs->n_subs;i++) { 195 ierr = MatGetSubMatrix(A_BI,is_subset_B[i],is_I,MAT_INITIAL_MATRIX,&AE_EI[i]);CHKERRQ(ierr); 196 } 197 198 /* setup Schur complements on subset */ 199 for (i=0;i<sub_schurs->n_subs;i++) { 200 ierr = MatCreateSchurComplement(AE_II,AE_II,AE_IE[i],AE_EI[i],AE_EE[i],&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 201 ierr = MatCreateVecs(sub_schurs->S_Ej[i],&sub_schurs->work1[i],&sub_schurs->work2[i]);CHKERRQ(ierr); 202 if (AE_II == A_II) { /* we can reuse the same ksp */ 203 KSP ksp; 204 ierr = MatSchurComplementGetKSP(S,&ksp);CHKERRQ(ierr); 205 ierr = MatSchurComplementSetKSP(sub_schurs->S_Ej[i],ksp);CHKERRQ(ierr); 206 } else { /* build new ksp object which inherits ksp and pc types from the original one */ 207 KSP origksp,schurksp; 208 PC origpc,schurpc; 209 KSPType ksp_type; 210 PCType pc_type; 211 PetscInt n_internal; 212 213 ierr = MatSchurComplementGetKSP(S,&origksp);CHKERRQ(ierr); 214 ierr = MatSchurComplementGetKSP(sub_schurs->S_Ej[i],&schurksp);CHKERRQ(ierr); 215 ierr = KSPGetType(origksp,&ksp_type);CHKERRQ(ierr); 216 ierr = KSPSetType(schurksp,ksp_type);CHKERRQ(ierr); 217 ierr = KSPGetPC(schurksp,&schurpc);CHKERRQ(ierr); 218 ierr = KSPGetPC(origksp,&origpc);CHKERRQ(ierr); 219 ierr = PCGetType(origpc,&pc_type);CHKERRQ(ierr); 220 ierr = PCSetType(schurpc,pc_type);CHKERRQ(ierr); 221 ierr = ISGetSize(is_I,&n_internal);CHKERRQ(ierr); 222 if (n_internal) { /* UMFPACK gives error with 0 sized problems */ 223 MatSolverPackage solver=NULL; 224 ierr = PCFactorGetMatSolverPackage(origpc,(const MatSolverPackage*)&solver);CHKERRQ(ierr); 225 if (solver) { 226 ierr = PCFactorSetMatSolverPackage(schurpc,solver);CHKERRQ(ierr); 227 } 228 } 229 ierr = KSPSetUp(schurksp);CHKERRQ(ierr); 230 } 231 } 232 /* free */ 233 ierr = ISDestroy(&is_I);CHKERRQ(ierr); 234 ierr = MatDestroy(&AE_II);CHKERRQ(ierr); 235 for (i=0;i<sub_schurs->n_subs;i++) { 236 ierr = MatDestroy(&AE_EE[i]);CHKERRQ(ierr); 237 ierr = MatDestroy(&AE_IE[i]);CHKERRQ(ierr); 238 ierr = MatDestroy(&AE_EI[i]);CHKERRQ(ierr); 239 ierr = ISDestroy(&is_subset_B[i]);CHKERRQ(ierr); 240 } 241 ierr = PetscFree4(is_subset_B,AE_IE,AE_EI,AE_EE);CHKERRQ(ierr); 242 PetscFunctionReturn(0); 243 } 244