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