1 #include <../src/ksp/pc/impls/bddc/bddc.h> 2 #include <../src/ksp/pc/impls/bddc/bddcprivate.h> 3 #include <petscblaslapack.h> 4 5 PETSC_STATIC_INLINE PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt*,PetscInt,PetscBT,PetscInt*,PetscInt*,PetscInt*); 6 static PetscErrorCode PCBDDCComputeExplicitSchur(Mat,PetscBool,MatReuse,Mat*); 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "PCBDDCComputeExplicitSchur" 10 static PetscErrorCode PCBDDCComputeExplicitSchur(Mat M, PetscBool issym, MatReuse reuse, Mat *S) 11 { 12 Mat B, C, D, Bd, Cd, AinvBd; 13 KSP ksp; 14 PC pc; 15 PetscBool isLU, isILU, isCHOL, Bdense, Cdense; 16 PetscReal fill = 2.0; 17 PetscInt n_I; 18 PetscMPIInt size; 19 PetscErrorCode ierr; 20 21 PetscFunctionBegin; 22 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)M),&size);CHKERRQ(ierr); 23 if (size != 1) { 24 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not for parallel matrices"); 25 } 26 if (reuse == MAT_REUSE_MATRIX) { 27 PetscBool Sdense; 28 29 ierr = PetscObjectTypeCompare((PetscObject)*S, MATSEQDENSE, &Sdense);CHKERRQ(ierr); 30 if (!Sdense) { 31 SETERRQ(PetscObjectComm((PetscObject)M),PETSC_ERR_SUP,"S should dense"); 32 } 33 } 34 ierr = MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D);CHKERRQ(ierr); 35 ierr = MatSchurComplementGetKSP(M, &ksp);CHKERRQ(ierr); 36 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr); 37 ierr = PetscObjectTypeCompare((PetscObject) pc, PCLU, &isLU);CHKERRQ(ierr); 38 ierr = PetscObjectTypeCompare((PetscObject) pc, PCILU, &isILU);CHKERRQ(ierr); 39 ierr = PetscObjectTypeCompare((PetscObject) pc, PCCHOLESKY, &isCHOL);CHKERRQ(ierr); 40 ierr = PetscObjectTypeCompare((PetscObject) B, MATSEQDENSE, &Bdense);CHKERRQ(ierr); 41 ierr = PetscObjectTypeCompare((PetscObject) C, MATSEQDENSE, &Cdense);CHKERRQ(ierr); 42 ierr = MatGetSize(B,&n_I,NULL);CHKERRQ(ierr); 43 if (n_I) { 44 if (!Bdense) { 45 ierr = MatConvert(B, MATSEQDENSE, MAT_INITIAL_MATRIX, &Bd);CHKERRQ(ierr); 46 } else { 47 Bd = B; 48 } 49 50 if (isLU || isILU || isCHOL) { 51 Mat fact; 52 ierr = KSPSetUp(ksp);CHKERRQ(ierr); 53 ierr = PCFactorGetMatrix(pc, &fact);CHKERRQ(ierr); 54 ierr = MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);CHKERRQ(ierr); 55 ierr = MatMatSolve(fact, Bd, AinvBd);CHKERRQ(ierr); 56 } else { 57 PetscBool ex = PETSC_TRUE; 58 59 if (ex) { 60 Mat Ainvd; 61 62 ierr = PCComputeExplicitOperator(pc, &Ainvd);CHKERRQ(ierr); 63 ierr = MatMatMult(Ainvd, Bd, MAT_INITIAL_MATRIX, fill, &AinvBd);CHKERRQ(ierr); 64 ierr = MatDestroy(&Ainvd);CHKERRQ(ierr); 65 } else { 66 Vec sol,rhs; 67 PetscScalar *arrayrhs,*arraysol; 68 PetscInt i,nrhs,n; 69 70 ierr = MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);CHKERRQ(ierr); 71 ierr = MatGetSize(Bd,&n,&nrhs);CHKERRQ(ierr); 72 ierr = MatDenseGetArray(Bd,&arrayrhs);CHKERRQ(ierr); 73 ierr = MatDenseGetArray(AinvBd,&arraysol);CHKERRQ(ierr); 74 ierr = KSPGetSolution(ksp,&sol);CHKERRQ(ierr); 75 ierr = KSPGetRhs(ksp,&rhs);CHKERRQ(ierr); 76 for (i=0;i<nrhs;i++) { 77 ierr = VecPlaceArray(rhs,arrayrhs+i*n);CHKERRQ(ierr); 78 ierr = VecPlaceArray(sol,arraysol+i*n);CHKERRQ(ierr); 79 ierr = KSPSolve(ksp,rhs,sol);CHKERRQ(ierr); 80 ierr = VecResetArray(rhs);CHKERRQ(ierr); 81 ierr = VecResetArray(sol);CHKERRQ(ierr); 82 } 83 ierr = MatDenseRestoreArray(Bd,&arrayrhs);CHKERRQ(ierr); 84 ierr = MatDenseRestoreArray(AinvBd,&arrayrhs);CHKERRQ(ierr); 85 } 86 } 87 if (!Bdense & !issym) { 88 ierr = MatDestroy(&Bd);CHKERRQ(ierr); 89 } 90 91 if (!issym) { 92 if (!Cdense) { 93 ierr = MatConvert(C, MATSEQDENSE, MAT_INITIAL_MATRIX, &Cd);CHKERRQ(ierr); 94 } else { 95 Cd = C; 96 } 97 ierr = MatMatMult(Cd, AinvBd, reuse, fill, S);CHKERRQ(ierr); 98 if (!Cdense) { 99 ierr = MatDestroy(&Cd);CHKERRQ(ierr); 100 } 101 } else { 102 ierr = MatTransposeMatMult(Bd, AinvBd, reuse, fill, S);CHKERRQ(ierr); 103 if (!Bdense) { 104 ierr = MatDestroy(&Bd);CHKERRQ(ierr); 105 } 106 } 107 ierr = MatDestroy(&AinvBd);CHKERRQ(ierr); 108 } 109 110 if (D) { 111 Mat Dd; 112 PetscBool Ddense; 113 114 ierr = PetscObjectTypeCompare((PetscObject)D,MATSEQDENSE,&Ddense);CHKERRQ(ierr); 115 if (!Ddense) { 116 ierr = MatConvert(D, MATSEQDENSE, MAT_INITIAL_MATRIX, &Dd);CHKERRQ(ierr); 117 } else { 118 Dd = D; 119 } 120 if (n_I) { 121 ierr = MatAYPX(*S,-1.0,Dd,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 122 } else { 123 if (reuse == MAT_INITIAL_MATRIX) { 124 ierr = MatDuplicate(Dd,MAT_COPY_VALUES,S);CHKERRQ(ierr); 125 } else { 126 ierr = MatCopy(Dd,*S,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 127 } 128 } 129 if (!Ddense) { 130 ierr = MatDestroy(&Dd);CHKERRQ(ierr); 131 } 132 } else { 133 ierr = MatScale(*S,-1.0);CHKERRQ(ierr); 134 } 135 PetscFunctionReturn(0); 136 } 137 138 #undef __FUNCT__ 139 #define __FUNCT__ "PCBDDCSubSchursSetUp" 140 PetscErrorCode PCBDDCSubSchursSetUp(PCBDDCSubSchurs sub_schurs, PetscInt xadj[], PetscInt adjncy[], PetscInt nlayers, PetscBool compute_Stilda, PetscBool deluxe, PetscBool invert_Stildas, PetscBool use_edges, PetscBool use_faces) 141 { 142 Mat A_II,A_IB,A_BI,A_BB; 143 Mat AE_II,*AE_IE,*AE_EI,*AE_EE; 144 Mat S_all,S_all_inv; 145 Mat global_schur_subsets,work_mat; 146 Mat S_Ej_tilda_all,S_Ej_inv_all; 147 ISLocalToGlobalMapping l2gmap_subsets; 148 IS is_I,*is_subset_B,temp_is; 149 PetscInt *nnz,*all_local_idx_G,*all_local_idx_N; 150 PetscInt i,subset_size,max_subset_size_seq; 151 PetscInt extra,local_size,global_size; 152 PetscBLASInt B_N,B_ierr,B_lwork,*pivots; 153 PetscScalar *Bwork; 154 PetscErrorCode ierr; 155 156 PetscFunctionBegin; 157 /* get Schur complement matrices */ 158 if (!sub_schurs->use_mumps) { 159 PetscBool isseqaij; 160 161 if (compute_Stilda) { 162 SETERRQ(PetscObjectComm((PetscObject)sub_schurs->l2gmap),PETSC_ERR_SUP,"Adaptive selection of constraints requires MUMPS"); 163 } 164 ierr = MatSchurComplementGetSubMatrices(sub_schurs->S,&A_II,NULL,&A_IB,&A_BI,&A_BB);CHKERRQ(ierr); 165 ierr = PetscObjectTypeCompare((PetscObject)A_BB,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 166 if (!isseqaij) { 167 ierr = MatConvert(A_BB,MATSEQAIJ,MAT_REUSE_MATRIX,&A_BB);CHKERRQ(ierr); 168 ierr = MatConvert(A_IB,MATSEQAIJ,MAT_REUSE_MATRIX,&A_IB);CHKERRQ(ierr); 169 ierr = MatConvert(A_BI,MATSEQAIJ,MAT_REUSE_MATRIX,&A_BI);CHKERRQ(ierr); 170 } 171 ierr = PetscMalloc4(sub_schurs->n_subs,&is_subset_B, 172 sub_schurs->n_subs,&AE_IE, 173 sub_schurs->n_subs,&AE_EI, 174 sub_schurs->n_subs,&AE_EE);CHKERRQ(ierr); 175 } else { 176 is_subset_B = NULL; 177 AE_IE = NULL; 178 AE_EI = NULL; 179 AE_EE = NULL; 180 } 181 182 /* determine interior problems */ 183 ierr = ISDestroy(&sub_schurs->is_I_layer);CHKERRQ(ierr); 184 ierr = ISGetLocalSize(sub_schurs->is_I,&i);CHKERRQ(ierr); 185 if (nlayers >= 0 && i) { /* Interior problems can be different from the original one */ 186 PetscBT touched; 187 const PetscInt* idx_B; 188 PetscInt n_I,n_B,n_local_dofs,n_prev_added,j,layer,*local_numbering; 189 190 if (xadj == NULL || adjncy == NULL) { 191 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Cannot request layering without adjacency"); 192 } 193 /* get sizes */ 194 ierr = ISGetLocalSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr); 195 ierr = ISGetLocalSize(sub_schurs->is_B,&n_B);CHKERRQ(ierr); 196 197 ierr = PetscMalloc1(n_I+n_B,&local_numbering);CHKERRQ(ierr); 198 ierr = PetscBTCreate(n_I+n_B,&touched);CHKERRQ(ierr); 199 ierr = PetscBTMemzero(n_I+n_B,touched);CHKERRQ(ierr); 200 201 /* all boundary dofs must be skipped when adding layers */ 202 ierr = ISGetIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr); 203 for (j=0;j<n_B;j++) { 204 ierr = PetscBTSet(touched,idx_B[j]);CHKERRQ(ierr); 205 } 206 ierr = PetscMemcpy(local_numbering,idx_B,n_B*sizeof(PetscInt));CHKERRQ(ierr); 207 ierr = ISRestoreIndices(sub_schurs->is_B,&idx_B);CHKERRQ(ierr); 208 209 /* add prescribed number of layers of dofs */ 210 n_local_dofs = n_B; 211 n_prev_added = n_B; 212 for (layer=0;layer<nlayers;layer++) { 213 PetscInt n_added; 214 if (n_local_dofs == n_I+n_B) break; 215 if (n_local_dofs > n_I+n_B) { 216 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); 217 } 218 ierr = PCBDDCAdjGetNextLayer_Private(local_numbering+n_local_dofs,n_prev_added,touched,xadj,adjncy,&n_added);CHKERRQ(ierr); 219 n_prev_added = n_added; 220 n_local_dofs += n_added; 221 if (!n_added) break; 222 } 223 ierr = PetscBTDestroy(&touched);CHKERRQ(ierr); 224 225 /* IS for I layer dofs in original numbering */ 226 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); 227 ierr = PetscFree(local_numbering);CHKERRQ(ierr); 228 ierr = ISSort(sub_schurs->is_I_layer);CHKERRQ(ierr); 229 /* IS for I layer dofs in I numbering */ 230 if (!sub_schurs->use_mumps) { 231 ISLocalToGlobalMapping ItoNmap; 232 ierr = ISLocalToGlobalMappingCreateIS(sub_schurs->is_I,&ItoNmap);CHKERRQ(ierr); 233 ierr = ISGlobalToLocalMappingApplyIS(ItoNmap,IS_GTOLM_DROP,sub_schurs->is_I_layer,&is_I);CHKERRQ(ierr); 234 ierr = ISLocalToGlobalMappingDestroy(&ItoNmap);CHKERRQ(ierr); 235 236 /* II block */ 237 ierr = MatGetSubMatrix(A_II,is_I,is_I,MAT_INITIAL_MATRIX,&AE_II);CHKERRQ(ierr); 238 } 239 } else { 240 PetscInt n_I; 241 242 /* IS for I dofs in original numbering */ 243 ierr = PetscObjectReference((PetscObject)sub_schurs->is_I);CHKERRQ(ierr); 244 sub_schurs->is_I_layer = sub_schurs->is_I; 245 246 /* IS for I dofs in I numbering (strided 1) */ 247 if (!sub_schurs->use_mumps) { 248 ierr = ISGetSize(sub_schurs->is_I,&n_I);CHKERRQ(ierr); 249 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub_schurs->is_I),n_I,0,1,&is_I);CHKERRQ(ierr); 250 251 /* II block is the same */ 252 ierr = PetscObjectReference((PetscObject)A_II);CHKERRQ(ierr); 253 AE_II = A_II; 254 } 255 } 256 /* Get info on subset sizes and sum of all subsets sizes */ 257 max_subset_size_seq = 0; 258 local_size = 0; 259 for (i=0;i<sub_schurs->n_subs_seq;i++) { 260 PetscInt j = sub_schurs->index_sequential[i]; 261 ierr = ISGetLocalSize(sub_schurs->is_subs[j],&subset_size);CHKERRQ(ierr); 262 max_subset_size_seq = PetscMax(subset_size,max_subset_size_seq); 263 local_size += subset_size; 264 } 265 266 /* Work arrays for local indices */ 267 extra = 0; 268 if (sub_schurs->use_mumps) { 269 ierr = ISGetLocalSize(sub_schurs->is_I_layer,&extra);CHKERRQ(ierr); 270 } 271 ierr = PetscMalloc1(local_size+extra,&all_local_idx_N);CHKERRQ(ierr); 272 if (extra) { 273 const PetscInt *idxs; 274 ierr = ISGetIndices(sub_schurs->is_I_layer,&idxs);CHKERRQ(ierr); 275 ierr = PetscMemcpy(all_local_idx_N,idxs,extra*sizeof(PetscInt));CHKERRQ(ierr); 276 ierr = ISRestoreIndices(sub_schurs->is_I_layer,&idxs);CHKERRQ(ierr); 277 } 278 ierr = PetscMalloc1(local_size,&nnz);CHKERRQ(ierr); 279 280 /* Get local indices in local numbering */ 281 local_size = 0; 282 for (i=0;i<sub_schurs->n_subs_seq;i++) { 283 PetscInt j; 284 const PetscInt *idxs; 285 286 PetscInt local_problem_index = sub_schurs->index_sequential[i]; 287 ierr = ISGetLocalSize(sub_schurs->is_subs[local_problem_index],&subset_size);CHKERRQ(ierr); 288 ierr = ISGetIndices(sub_schurs->is_subs[local_problem_index],&idxs);CHKERRQ(ierr); 289 /* subset indices in local numbering */ 290 ierr = PetscMemcpy(all_local_idx_N+local_size+extra,idxs,subset_size*sizeof(PetscInt));CHKERRQ(ierr); 291 ierr = ISRestoreIndices(sub_schurs->is_subs[local_problem_index],&idxs);CHKERRQ(ierr); 292 for (j=0;j<subset_size;j++) nnz[local_size+j] = subset_size; 293 local_size += subset_size; 294 } 295 296 S_all = NULL; 297 S_all_inv = NULL; 298 S_Ej_tilda_all = NULL; 299 S_Ej_inv_all = NULL; 300 sub_schurs->is_hermitian = PETSC_FALSE; 301 sub_schurs->is_posdef = PETSC_FALSE; 302 if (sub_schurs->A) { 303 ierr = MatIsHermitian(sub_schurs->A,0.0,&sub_schurs->is_hermitian);CHKERRQ(ierr); 304 if (sub_schurs->is_hermitian) { 305 PetscScalar val; 306 Vec vec1,vec2; 307 308 ierr = MatCreateVecs(sub_schurs->A,&vec1,&vec2);CHKERRQ(ierr); 309 ierr = VecSetRandom(vec1,NULL); 310 ierr = VecCopy(vec1,vec2);CHKERRQ(ierr); 311 ierr = MatMult(sub_schurs->A,vec2,vec1);CHKERRQ(ierr); 312 ierr = VecDot(vec1,vec2,&val);CHKERRQ(ierr); 313 if (PetscRealPart(val) > 0. && PetscImaginaryPart(val) == 0.) sub_schurs->is_posdef = PETSC_TRUE; 314 ierr = VecDestroy(&vec1);CHKERRQ(ierr); 315 ierr = VecDestroy(&vec2);CHKERRQ(ierr); 316 } 317 } 318 Bwork = NULL; 319 pivots = NULL; 320 if (sub_schurs->n_subs && deluxe && compute_Stilda && !sub_schurs->is_hermitian) { /* workspace needed only for GETRI */ 321 PetscScalar lwork; 322 323 B_lwork = -1; 324 ierr = PetscBLASIntCast(local_size,&B_N);CHKERRQ(ierr); 325 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 326 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,Bwork,&B_N,pivots,&lwork,&B_lwork,&B_ierr)); 327 ierr = PetscFPTrapPop();CHKERRQ(ierr); 328 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to GETRI Lapack routine %d",(int)B_ierr); 329 ierr = PetscBLASIntCast((PetscInt)PetscRealPart(lwork),&B_lwork);CHKERRQ(ierr); 330 ierr = PetscMalloc2(B_lwork,&Bwork,B_N,&pivots);CHKERRQ(ierr); 331 } 332 333 /* prepare parallel matrices for summing up properly schurs on subsets */ 334 ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)sub_schurs->l2gmap),sub_schurs->l2gmap,local_size,all_local_idx_N+extra,PETSC_NULL,&global_size,&all_local_idx_G);CHKERRQ(ierr); 335 ierr = ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)sub_schurs->l2gmap),1,local_size,all_local_idx_G,PETSC_COPY_VALUES,&l2gmap_subsets);CHKERRQ(ierr); 336 ierr = MatCreateIS(PetscObjectComm((PetscObject)sub_schurs->l2gmap),1,PETSC_DECIDE,PETSC_DECIDE,global_size,global_size,l2gmap_subsets,&work_mat);CHKERRQ(ierr); 337 ierr = ISLocalToGlobalMappingDestroy(&l2gmap_subsets);CHKERRQ(ierr); 338 ierr = MatCreate(PetscObjectComm((PetscObject)work_mat),&global_schur_subsets);CHKERRQ(ierr); 339 ierr = MatSetSizes(global_schur_subsets,PETSC_DECIDE,PETSC_DECIDE,global_size,global_size);CHKERRQ(ierr); 340 ierr = MatSetType(global_schur_subsets,MATMPIAIJ);CHKERRQ(ierr); 341 342 if (sub_schurs->n_subs) { 343 if (!sub_schurs->use_mumps) { 344 /* subsets in original and boundary numbering */ 345 for (i=0;i<sub_schurs->n_subs;i++) { 346 ierr = ISGlobalToLocalMappingApplyIS(sub_schurs->BtoNmap,IS_GTOLM_DROP,sub_schurs->is_subs[i],&is_subset_B[i]);CHKERRQ(ierr); 347 } 348 349 /* EE block */ 350 for (i=0;i<sub_schurs->n_subs;i++) { 351 ierr = MatGetSubMatrix(A_BB,is_subset_B[i],is_subset_B[i],MAT_INITIAL_MATRIX,&AE_EE[i]);CHKERRQ(ierr); 352 } 353 354 /* IE block */ 355 for (i=0;i<sub_schurs->n_subs;i++) { 356 ierr = MatGetSubMatrix(A_IB,is_I,is_subset_B[i],MAT_INITIAL_MATRIX,&AE_IE[i]);CHKERRQ(ierr); 357 } 358 359 /* EI block */ 360 for (i=0;i<sub_schurs->n_subs;i++) { 361 if (sub_schurs->is_hermitian) { 362 ierr = MatCreateTranspose(AE_IE[i],&AE_EI[i]);CHKERRQ(ierr); 363 } else { 364 ierr = MatGetSubMatrix(A_BI,is_subset_B[i],is_I,MAT_INITIAL_MATRIX,&AE_EI[i]);CHKERRQ(ierr); 365 } 366 } 367 368 /* setup Schur complements on subset */ 369 for (i=0;i<sub_schurs->n_subs;i++) { 370 ierr = MatDestroy(&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 371 ierr = MatCreateSchurComplement(AE_II,AE_II,AE_IE[i],AE_EI[i],AE_EE[i],&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 372 ierr = MatDestroy(&AE_EE[i]);CHKERRQ(ierr); 373 ierr = MatDestroy(&AE_IE[i]);CHKERRQ(ierr); 374 ierr = MatDestroy(&AE_EI[i]);CHKERRQ(ierr); 375 if (AE_II == A_II) { /* we can reuse the same ksp */ 376 KSP ksp; 377 ierr = MatSchurComplementGetKSP(sub_schurs->S,&ksp);CHKERRQ(ierr); 378 ierr = MatSchurComplementSetKSP(sub_schurs->S_Ej[i],ksp);CHKERRQ(ierr); 379 } else { /* build new ksp object which inherits ksp and pc types from the original one */ 380 KSP origksp,schurksp; 381 PC origpc,schurpc; 382 KSPType ksp_type; 383 PCType pc_type; 384 PetscInt n_internal; 385 386 ierr = MatSchurComplementGetKSP(sub_schurs->S,&origksp);CHKERRQ(ierr); 387 ierr = MatSchurComplementGetKSP(sub_schurs->S_Ej[i],&schurksp);CHKERRQ(ierr); 388 ierr = KSPGetType(origksp,&ksp_type);CHKERRQ(ierr); 389 ierr = KSPSetType(schurksp,ksp_type);CHKERRQ(ierr); 390 ierr = KSPGetPC(schurksp,&schurpc);CHKERRQ(ierr); 391 ierr = KSPGetPC(origksp,&origpc);CHKERRQ(ierr); 392 ierr = PCGetType(origpc,&pc_type);CHKERRQ(ierr); 393 ierr = PCSetType(schurpc,pc_type);CHKERRQ(ierr); 394 ierr = ISGetSize(is_I,&n_internal);CHKERRQ(ierr); 395 if (n_internal) { /* UMFPACK gives error with 0 sized problems */ 396 MatSolverPackage solver=NULL; 397 ierr = PCFactorGetMatSolverPackage(origpc,(const MatSolverPackage*)&solver);CHKERRQ(ierr); 398 if (solver) { 399 ierr = PCFactorSetMatSolverPackage(schurpc,solver);CHKERRQ(ierr); 400 } 401 } 402 ierr = KSPSetUp(schurksp);CHKERRQ(ierr); 403 } 404 } 405 /* free */ 406 ierr = ISDestroy(&is_I);CHKERRQ(ierr); 407 ierr = MatDestroy(&AE_II);CHKERRQ(ierr); 408 for (i=0;i<sub_schurs->n_subs;i++) { 409 ierr = ISDestroy(&is_subset_B[i]);CHKERRQ(ierr); 410 } 411 ierr = PetscFree4(is_subset_B,AE_IE,AE_EI,AE_EE);CHKERRQ(ierr); 412 #if defined(PETSC_HAVE_MUMPS) 413 } else { 414 Mat A,F; 415 IS is_A_all; 416 PetscInt *idxs_schur,n_I; 417 418 /* get working mat */ 419 ierr = ISGetLocalSize(sub_schurs->is_I_layer,&n_I);CHKERRQ(ierr); 420 ierr = ISCreateGeneral(PETSC_COMM_SELF,local_size+n_I,all_local_idx_N,PETSC_COPY_VALUES,&is_A_all);CHKERRQ(ierr); 421 ierr = MatGetSubMatrixUnsorted(sub_schurs->A,is_A_all,is_A_all,&A);CHKERRQ(ierr); 422 ierr = ISDestroy(&is_A_all);CHKERRQ(ierr); 423 424 if (n_I) { 425 if (sub_schurs->is_hermitian && sub_schurs->is_posdef) { 426 ierr = MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&F);CHKERRQ(ierr); 427 } else { 428 ierr = MatGetFactor(A,MATSOLVERMUMPS,MAT_FACTOR_LU,&F);CHKERRQ(ierr); 429 } 430 431 /* subsets ordered last */ 432 ierr = PetscMalloc1(local_size,&idxs_schur);CHKERRQ(ierr); 433 for (i=0;i<local_size;i++) { 434 idxs_schur[i] = n_I+i+1; 435 } 436 ierr = MatMumpsSetSchurIndices(F,local_size,idxs_schur);CHKERRQ(ierr); 437 ierr = PetscFree(idxs_schur);CHKERRQ(ierr); 438 439 /* factorization step */ 440 if (sub_schurs->is_hermitian && sub_schurs->is_posdef) { 441 ierr = MatCholeskyFactorSymbolic(F,A,NULL,NULL);CHKERRQ(ierr); 442 ierr = MatCholeskyFactorNumeric(F,A,NULL);CHKERRQ(ierr); 443 } else { 444 ierr = MatLUFactorSymbolic(F,A,NULL,NULL,NULL);CHKERRQ(ierr); 445 ierr = MatLUFactorNumeric(F,A,NULL);CHKERRQ(ierr); 446 } 447 448 /* get explicit Schur Complement computed during numeric factorization */ 449 ierr = MatMumpsGetSchurComplement(F,&S_all);CHKERRQ(ierr); 450 ierr = MatDestroy(&F);CHKERRQ(ierr); 451 } else { 452 ierr = MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&S_all);CHKERRQ(ierr); 453 } 454 ierr = MatDestroy(&A);CHKERRQ(ierr); 455 #endif 456 } 457 } else { 458 ierr = PetscFree(nnz);CHKERRQ(ierr); 459 } 460 461 if (!sub_schurs->n_subs_seq_g) { /* need to get rid of the old parallel deluxe */ 462 ierr = PetscFree(all_local_idx_G);CHKERRQ(ierr); 463 ierr = PetscFree(all_local_idx_N);CHKERRQ(ierr); 464 ierr = PetscFree2(Bwork,pivots);CHKERRQ(ierr); 465 ierr = MatDestroy(&global_schur_subsets);CHKERRQ(ierr); 466 ierr = MatDestroy(&work_mat);CHKERRQ(ierr); 467 PetscFunctionReturn(0); 468 } 469 470 /* subset indices in local boundary numbering */ 471 if (!sub_schurs->is_Ej_all) { 472 PetscInt *all_local_idx_B; 473 474 ierr = PetscMalloc1(local_size,&all_local_idx_B);CHKERRQ(ierr); 475 ierr = ISGlobalToLocalMappingApply(sub_schurs->BtoNmap,IS_GTOLM_DROP,local_size,all_local_idx_N+extra,&subset_size,all_local_idx_B);CHKERRQ(ierr); 476 if (subset_size != local_size) { 477 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in sub_schurs serial (BtoNmap)! %d != %d\n",subset_size,local_size); 478 } 479 ierr = ISCreateGeneral(PETSC_COMM_SELF,local_size,all_local_idx_B,PETSC_OWN_POINTER,&sub_schurs->is_Ej_all);CHKERRQ(ierr); 480 } 481 ierr = PetscFree(all_local_idx_N);CHKERRQ(ierr); 482 483 /* Local matrix of all local Schur on subsets (transposed) */ 484 if (!sub_schurs->S_Ej_all) { 485 ierr = MatCreate(PETSC_COMM_SELF,&sub_schurs->S_Ej_all);CHKERRQ(ierr); 486 ierr = MatSetSizes(sub_schurs->S_Ej_all,PETSC_DECIDE,PETSC_DECIDE,local_size,local_size);CHKERRQ(ierr); 487 ierr = MatSetType(sub_schurs->S_Ej_all,MATAIJ);CHKERRQ(ierr); 488 ierr = MatSeqAIJSetPreallocation(sub_schurs->S_Ej_all,0,nnz);CHKERRQ(ierr); 489 } else { 490 ierr = MatZeroEntries(sub_schurs->S_Ej_all);CHKERRQ(ierr); 491 } 492 493 ierr = PetscBTMemzero(sub_schurs->n_subs,sub_schurs->computed_Stilda_subs);CHKERRQ(ierr); 494 if (!sub_schurs->use_mumps) { 495 if (sub_schurs->n_subs_seq) { 496 PetscScalar *work; 497 PetscInt *dummy_idx; 498 499 /* Work arrays */ 500 ierr = PetscMalloc2(max_subset_size_seq,&dummy_idx,max_subset_size_seq*max_subset_size_seq,&work);CHKERRQ(ierr); 501 502 /* Loop on local problems end compute Schur complements explicitly */ 503 local_size = 0; 504 for (i=0;i<sub_schurs->n_subs_seq;i++) { 505 Mat S_Ej_expl; 506 PetscInt j,lpi = sub_schurs->index_sequential[i]; 507 PetscBool Sdense; 508 509 ierr = ISGetLocalSize(sub_schurs->is_subs[lpi],&subset_size);CHKERRQ(ierr); 510 ierr = MatCreateSeqDense(PETSC_COMM_SELF,subset_size,subset_size,work,&S_Ej_expl);CHKERRQ(ierr); 511 ierr = PCBDDCComputeExplicitSchur(sub_schurs->S_Ej[lpi],sub_schurs->is_hermitian,MAT_REUSE_MATRIX,&S_Ej_expl);CHKERRQ(ierr); 512 ierr = PetscObjectTypeCompare((PetscObject)S_Ej_expl,MATSEQDENSE,&Sdense);CHKERRQ(ierr); 513 if (Sdense) { 514 for (j=0;j<subset_size;j++) { 515 dummy_idx[j]=local_size+j; 516 } 517 ierr = MatSetValues(sub_schurs->S_Ej_all,subset_size,dummy_idx,subset_size,dummy_idx,work,INSERT_VALUES);CHKERRQ(ierr); 518 } else { 519 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet implemented for sparse matrices"); 520 } 521 local_size += subset_size; 522 ierr = MatDestroy(&sub_schurs->S_Ej[lpi]);CHKERRQ(ierr); 523 sub_schurs->S_Ej[lpi] = S_Ej_expl; 524 } 525 ierr = PetscFree2(dummy_idx,work);CHKERRQ(ierr); 526 } 527 } else { 528 PetscInt *dummy_idx,n_all; 529 PetscScalar *work; 530 531 if (compute_Stilda) { 532 ierr = MatCreate(PETSC_COMM_SELF,&S_Ej_tilda_all);CHKERRQ(ierr); 533 ierr = MatSetSizes(S_Ej_tilda_all,PETSC_DECIDE,PETSC_DECIDE,local_size,local_size);CHKERRQ(ierr); 534 ierr = MatSetType(S_Ej_tilda_all,MATAIJ);CHKERRQ(ierr); 535 ierr = MatSeqAIJSetPreallocation(S_Ej_tilda_all,0,nnz);CHKERRQ(ierr); 536 if (deluxe) { 537 ierr = MatCreate(PETSC_COMM_SELF,&S_Ej_inv_all);CHKERRQ(ierr); 538 ierr = MatSetSizes(S_Ej_inv_all,PETSC_DECIDE,PETSC_DECIDE,local_size,local_size);CHKERRQ(ierr); 539 ierr = MatSetType(S_Ej_inv_all,MATAIJ);CHKERRQ(ierr); 540 ierr = MatSeqAIJSetPreallocation(S_Ej_inv_all,0,nnz);CHKERRQ(ierr); 541 } 542 } 543 544 if (S_all) { /* multilevel */ 545 ierr = MatGetSize(S_all,&n_all,NULL);CHKERRQ(ierr); 546 } else { 547 n_all = 0; 548 } 549 local_size = 0; 550 551 /* Work arrays */ 552 if (compute_Stilda) { 553 ierr = PetscMalloc2(max_subset_size_seq,&dummy_idx,n_all*n_all,&work);CHKERRQ(ierr); 554 } else { 555 ierr = PetscMalloc2(max_subset_size_seq,&dummy_idx,0,&work);CHKERRQ(ierr); 556 } 557 558 if (compute_Stilda && S_all) { 559 PetscScalar *vals; 560 ierr = PetscBLASIntCast(n_all,&B_N);CHKERRQ(ierr); 561 ierr = MatDuplicate(S_all,MAT_COPY_VALUES,&S_all_inv);CHKERRQ(ierr); 562 ierr = MatDenseGetArray(S_all_inv,&vals);CHKERRQ(ierr); 563 if (!sub_schurs->is_hermitian) { 564 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,vals,&B_N,pivots,&B_ierr)); 565 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr); 566 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,vals,&B_N,pivots,Bwork,&B_lwork,&B_ierr)); 567 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr); 568 } else { 569 PetscInt j,k; 570 571 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&B_N,vals,&B_N,&B_ierr)); 572 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRF Lapack routine %d",(int)B_ierr); 573 PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&B_N,vals,&B_N,&B_ierr)); 574 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRI Lapack routine %d",(int)B_ierr); 575 for (j=0;j<B_N;j++) { 576 for (k=j+1;k<B_N;k++) { 577 vals[k*B_N+j] = vals[j*B_N+k]; 578 } 579 } 580 } 581 ierr = MatDenseRestoreArray(S_all_inv,&vals);CHKERRQ(ierr); 582 } 583 584 for (i=0;i<sub_schurs->n_subs;i++) { 585 IS is_E; 586 PetscScalar *vals; 587 PetscInt j; 588 589 ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr); 590 for (j=0;j<subset_size;j++) { 591 dummy_idx[j]=local_size+j; 592 } 593 ierr = ISCreateStride(PetscObjectComm((PetscObject)sub_schurs->is_I),subset_size,local_size,1,&is_E);CHKERRQ(ierr); 594 if (sub_schurs->S_Ej[i]) { 595 ierr = MatGetSubMatrix(S_all,is_E,is_E,MAT_REUSE_MATRIX,&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 596 } else { 597 ierr = MatGetSubMatrix(S_all,is_E,is_E,MAT_INITIAL_MATRIX,&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 598 } 599 ierr = MatDenseGetArray(sub_schurs->S_Ej[i],&vals);CHKERRQ(ierr); 600 ierr = MatSetValues(sub_schurs->S_Ej_all,subset_size,dummy_idx,subset_size,dummy_idx,vals,INSERT_VALUES);CHKERRQ(ierr); 601 ierr = MatDenseRestoreArray(sub_schurs->S_Ej[i],&vals);CHKERRQ(ierr); 602 603 if (compute_Stilda && ((PetscBTLookup(sub_schurs->is_edge,i) && use_edges) || (!PetscBTLookup(sub_schurs->is_edge,i) && use_faces))) { 604 Mat Stilda; 605 606 ierr = PetscBTSet(sub_schurs->computed_Stilda_subs,i);CHKERRQ(ierr); 607 ierr = MatGetSubMatrix(S_all_inv,is_E,is_E,MAT_INITIAL_MATRIX,&Stilda);CHKERRQ(ierr); 608 ierr = MatDenseGetArray(Stilda,&vals);CHKERRQ(ierr); 609 if (deluxe) { /* when using deluxe scaling, we need (S_1^-1+S_2^-1)^-1 */ 610 PetscScalar *vals2; 611 612 ierr = PetscBLASIntCast(subset_size,&B_N);CHKERRQ(ierr); 613 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 614 ierr = MatDenseGetArray(sub_schurs->S_Ej[i],&vals2);CHKERRQ(ierr); 615 if (!sub_schurs->is_hermitian) { 616 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,vals2,&B_N,pivots,&B_ierr)); 617 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr); 618 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,vals2,&B_N,pivots,Bwork,&B_lwork,&B_ierr)); 619 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr); 620 } else { 621 PetscInt j,k; 622 623 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&B_N,vals2,&B_N,&B_ierr)); 624 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRF Lapack routine %d",(int)B_ierr); 625 PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&B_N,vals2,&B_N,&B_ierr)); 626 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRI Lapack routine %d",(int)B_ierr); 627 for (j=0;j<B_N;j++) { 628 for (k=j+1;k<B_N;k++) { 629 vals2[k*B_N+j] = vals2[j*B_N+k]; 630 } 631 } 632 } 633 ierr = PetscFPTrapPop();CHKERRQ(ierr); 634 ierr = MatSetValues(S_Ej_inv_all,subset_size,dummy_idx,subset_size,dummy_idx,vals2,INSERT_VALUES);CHKERRQ(ierr); 635 ierr = MatDenseRestoreArray(sub_schurs->S_Ej[i],&vals2);CHKERRQ(ierr); 636 } 637 ierr = MatSetValues(S_Ej_tilda_all,subset_size,dummy_idx,subset_size,dummy_idx,vals,INSERT_VALUES);CHKERRQ(ierr); 638 ierr = MatDenseRestoreArray(Stilda,&vals);CHKERRQ(ierr); 639 ierr = MatDestroy(&Stilda);CHKERRQ(ierr); 640 } 641 ierr = ISDestroy(&is_E);CHKERRQ(ierr); 642 local_size += subset_size; 643 } 644 ierr = PetscFree2(dummy_idx,work);CHKERRQ(ierr); 645 } 646 ierr = PetscFree(nnz);CHKERRQ(ierr); 647 ierr = MatDestroy(&S_all);CHKERRQ(ierr); 648 ierr = MatDestroy(&S_all_inv);CHKERRQ(ierr); 649 ierr = MatAssemblyBegin(sub_schurs->S_Ej_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 650 ierr = MatAssemblyEnd(sub_schurs->S_Ej_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 651 if (S_Ej_tilda_all) { 652 ierr = MatAssemblyBegin(S_Ej_tilda_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 653 ierr = MatAssemblyEnd(S_Ej_tilda_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 654 } 655 if (S_Ej_inv_all) { 656 ierr = MatAssemblyBegin(S_Ej_inv_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 657 ierr = MatAssemblyEnd(S_Ej_inv_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 658 } 659 660 /* Global matrix of all assembled Schur on subsets */ 661 ierr = MatISSetLocalMat(work_mat,sub_schurs->S_Ej_all);CHKERRQ(ierr); 662 ierr = MatISSetMPIXAIJPreallocation_Private(work_mat,global_schur_subsets,PETSC_TRUE);CHKERRQ(ierr); 663 ierr = MatISGetMPIXAIJ(work_mat,MAT_REUSE_MATRIX,&global_schur_subsets);CHKERRQ(ierr); 664 /* Get local part of (\sum_j S_Ej) */ 665 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)sub_schurs->l2gmap),local_size,all_local_idx_G,PETSC_OWN_POINTER,&temp_is);CHKERRQ(ierr); 666 ierr = MatDestroy(&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr); 667 ierr = MatGetSubMatrixUnsorted(global_schur_subsets,temp_is,temp_is,&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr); 668 669 if (S_Ej_tilda_all) { 670 ierr = MatISSetLocalMat(work_mat,S_Ej_tilda_all);CHKERRQ(ierr); 671 ierr = MatISGetMPIXAIJ(work_mat,MAT_REUSE_MATRIX,&global_schur_subsets);CHKERRQ(ierr); 672 ierr = MatDestroy(&sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr); 673 ierr = MatGetSubMatrixUnsorted(global_schur_subsets,temp_is,temp_is,&sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr); 674 } 675 if (S_Ej_inv_all) { 676 PetscInt cum; 677 PetscScalar *array,*array2; 678 ierr = MatISSetLocalMat(work_mat,S_Ej_inv_all);CHKERRQ(ierr); 679 ierr = MatISGetMPIXAIJ(work_mat,MAT_REUSE_MATRIX,&global_schur_subsets);CHKERRQ(ierr); 680 ierr = MatDestroy(&sub_schurs->sum_S_Ej_inv_all);CHKERRQ(ierr); 681 ierr = MatGetSubMatrixUnsorted(global_schur_subsets,temp_is,temp_is,&sub_schurs->sum_S_Ej_inv_all);CHKERRQ(ierr); 682 /* invert blocks of sum_S_Ej_inv_all */ 683 ierr = MatSeqAIJGetArray(sub_schurs->sum_S_Ej_inv_all,&array);CHKERRQ(ierr); 684 ierr = MatSeqAIJGetArray(sub_schurs->sum_S_Ej_tilda_all,&array2);CHKERRQ(ierr); 685 cum = 0; 686 for (i=0;i<sub_schurs->n_subs;i++) { 687 ierr = ISGetLocalSize(sub_schurs->is_subs[i],&subset_size);CHKERRQ(ierr); 688 if (PetscBTLookup(sub_schurs->computed_Stilda_subs,i)) { 689 ierr = PetscBLASIntCast(subset_size,&B_N);CHKERRQ(ierr); 690 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 691 if (!sub_schurs->is_hermitian) { 692 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,array+cum,&B_N,pivots,&B_ierr)); 693 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr); 694 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,array+cum,&B_N,pivots,Bwork,&B_lwork,&B_ierr)); 695 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr); 696 } else { 697 PetscInt j,k; 698 699 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&B_N,array+cum,&B_N,&B_ierr)); 700 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRF Lapack routine %d",(int)B_ierr); 701 PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&B_N,array+cum,&B_N,&B_ierr)); 702 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRI Lapack routine %d",(int)B_ierr); 703 for (j=0;j<B_N;j++) { 704 for (k=j+1;k<B_N;k++) { 705 array[k*B_N+j+cum] = array[j*B_N+k+cum]; 706 } 707 } 708 } 709 if (invert_Stildas) { /* Stildas can be singular */ 710 if (!sub_schurs->is_hermitian) { 711 PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&B_N,&B_N,array2+cum,&B_N,pivots,&B_ierr)); 712 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRF Lapack routine %d",(int)B_ierr); 713 PetscStackCallBLAS("LAPACKgetri",LAPACKgetri_(&B_N,array2+cum,&B_N,pivots,Bwork,&B_lwork,&B_ierr)); 714 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GETRI Lapack routine %d",(int)B_ierr); 715 } else { 716 PetscInt j,k; 717 718 PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("L",&B_N,array2+cum,&B_N,&B_ierr)); 719 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRF Lapack routine %d",(int)B_ierr); 720 PetscStackCallBLAS("LAPACKpotri",LAPACKpotri_("L",&B_N,array2+cum,&B_N,&B_ierr)); 721 if (B_ierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in POTRI Lapack routine %d",(int)B_ierr); 722 for (j=0;j<B_N;j++) { 723 for (k=j+1;k<B_N;k++) { 724 array2[k*B_N+j+cum] = array2[j*B_N+k+cum]; 725 } 726 } 727 } 728 } 729 ierr = PetscFPTrapPop();CHKERRQ(ierr); 730 } 731 cum += subset_size*subset_size; 732 } 733 ierr = MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_inv_all,&array);CHKERRQ(ierr); 734 ierr = MatSeqAIJRestoreArray(sub_schurs->sum_S_Ej_tilda_all,&array2);CHKERRQ(ierr); 735 } 736 737 ierr = PetscFree2(Bwork,pivots);CHKERRQ(ierr); 738 ierr = MatDestroy(&global_schur_subsets);CHKERRQ(ierr); 739 ierr = MatDestroy(&S_Ej_tilda_all);CHKERRQ(ierr); 740 ierr = MatDestroy(&S_Ej_inv_all);CHKERRQ(ierr); 741 ierr = MatDestroy(&work_mat);CHKERRQ(ierr); 742 ierr = ISDestroy(&temp_is);CHKERRQ(ierr); 743 PetscFunctionReturn(0); 744 } 745 746 #undef __FUNCT__ 747 #define __FUNCT__ "PCBDDCSubSchursInit" 748 PetscErrorCode PCBDDCSubSchursInit(PCBDDCSubSchurs sub_schurs, Mat A, Mat S, IS is_I, IS is_B, PCBDDCGraph graph, ISLocalToGlobalMapping BtoNmap, PetscInt seqthreshold) 749 { 750 IS *faces,*edges,*all_cc,vertices; 751 PetscInt *index_sequential,*index_parallel; 752 PetscInt *auxlocal_sequential,*auxlocal_parallel; 753 PetscInt *auxglobal_sequential,*auxglobal_parallel; 754 PetscInt *auxmapping; 755 PetscInt i,max_subset_size; 756 PetscInt n_sequential_problems,n_local_sequential_problems,n_parallel_problems,n_local_parallel_problems; 757 PetscInt n_faces,n_edges,n_all_cc; 758 PetscBool is_sorted; 759 PetscErrorCode ierr; 760 761 PetscFunctionBegin; 762 ierr = ISSorted(is_I,&is_sorted);CHKERRQ(ierr); 763 if (!is_sorted) { 764 SETERRQ(PetscObjectComm((PetscObject)is_I),PETSC_ERR_PLIB,"IS for I dofs should be shorted"); 765 } 766 ierr = ISSorted(is_B,&is_sorted);CHKERRQ(ierr); 767 if (!is_sorted) { 768 SETERRQ(PetscObjectComm((PetscObject)is_B),PETSC_ERR_PLIB,"IS for B dofs should be shorted"); 769 } 770 771 /* reset any previous data */ 772 ierr = PCBDDCSubSchursReset(sub_schurs);CHKERRQ(ierr); 773 774 /* get index sets for faces and edges */ 775 ierr = PCBDDCGraphGetCandidatesIS(graph,&n_faces,&faces,&n_edges,&edges,&vertices);CHKERRQ(ierr); 776 n_all_cc = n_faces+n_edges; 777 ierr = PetscBTCreate(n_all_cc,&sub_schurs->is_edge);CHKERRQ(ierr); 778 ierr = PetscBTCreate(n_all_cc,&sub_schurs->computed_Stilda_subs);CHKERRQ(ierr); 779 ierr = PetscMalloc1(n_all_cc,&all_cc);CHKERRQ(ierr); 780 for (i=0;i<n_faces;i++) { 781 all_cc[i] = faces[i]; 782 } 783 for (i=0;i<n_edges;i++) { 784 all_cc[n_faces+i] = edges[i]; 785 ierr = PetscBTSet(sub_schurs->is_edge,n_faces+i);CHKERRQ(ierr); 786 } 787 ierr = PetscFree(faces);CHKERRQ(ierr); 788 ierr = PetscFree(edges);CHKERRQ(ierr); 789 790 /* map interface's subsets */ 791 max_subset_size = 0; 792 for (i=0;i<n_all_cc;i++) { 793 PetscInt subset_size; 794 ierr = ISGetLocalSize(all_cc[i],&subset_size);CHKERRQ(ierr); 795 max_subset_size = PetscMax(max_subset_size,subset_size); 796 } 797 ierr = PetscMalloc1(max_subset_size,&auxmapping);CHKERRQ(ierr); 798 ierr = PetscMalloc2(graph->ncc,&auxlocal_sequential, 799 graph->ncc,&auxlocal_parallel);CHKERRQ(ierr); 800 ierr = PetscMalloc2(graph->ncc,&index_sequential, 801 graph->ncc,&index_parallel);CHKERRQ(ierr); 802 803 /* if threshold is negative uses all sequential problems (possibly using MUMPS) */ 804 sub_schurs->use_mumps = PETSC_FALSE; 805 if (seqthreshold < 0) { 806 seqthreshold = max_subset_size; 807 #if defined(PETSC_HAVE_MUMPS) 808 sub_schurs->use_mumps = (PetscBool)(!!A); 809 #endif 810 } 811 812 813 /* determine which problem has to be solved in parallel or sequentially */ 814 n_local_sequential_problems = 0; 815 n_local_parallel_problems = 0; 816 for (i=0;i<n_all_cc;i++) { 817 PetscInt subset_size,j,min_loc = 0; 818 const PetscInt *idxs; 819 820 ierr = ISGetLocalSize(all_cc[i],&subset_size);CHKERRQ(ierr); 821 ierr = ISGetIndices(all_cc[i],&idxs);CHKERRQ(ierr); 822 ierr = ISLocalToGlobalMappingApply(graph->l2gmap,subset_size,idxs,auxmapping);CHKERRQ(ierr); 823 for (j=1;j<subset_size;j++) { 824 if (auxmapping[j]<auxmapping[min_loc]) { 825 min_loc = j; 826 } 827 } 828 if (subset_size > seqthreshold) { 829 index_parallel[n_local_parallel_problems] = i; 830 auxlocal_parallel[n_local_parallel_problems] = idxs[min_loc]; 831 n_local_parallel_problems++; 832 } else { 833 index_sequential[n_local_sequential_problems] = i; 834 auxlocal_sequential[n_local_sequential_problems] = idxs[min_loc]; 835 n_local_sequential_problems++; 836 } 837 ierr = ISRestoreIndices(all_cc[i],&idxs);CHKERRQ(ierr); 838 } 839 840 /* Number parallel problems */ 841 auxglobal_parallel = 0; 842 ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)graph->l2gmap),graph->l2gmap,n_local_parallel_problems,auxlocal_parallel,PETSC_NULL,&n_parallel_problems,&auxglobal_parallel);CHKERRQ(ierr); 843 844 /* Number sequential problems */ 845 auxglobal_sequential = 0; 846 ierr = PCBDDCSubsetNumbering(PetscObjectComm((PetscObject)graph->l2gmap),graph->l2gmap,n_local_sequential_problems,auxlocal_sequential,PETSC_NULL,&n_sequential_problems,&auxglobal_sequential);CHKERRQ(ierr); 847 848 /* update info in sub_schurs */ 849 if (A) { 850 PetscBool isseqaij; 851 852 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 853 if (isseqaij) { 854 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 855 sub_schurs->A = A; 856 } else { /* SeqBAIJ matrices does not support symmetry checking, SeqSBAIJ does not support MatPermute */ 857 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&sub_schurs->A);CHKERRQ(ierr); 858 } 859 } 860 ierr = PetscObjectReference((PetscObject)S);CHKERRQ(ierr); 861 sub_schurs->S = S; 862 ierr = PetscObjectReference((PetscObject)is_I);CHKERRQ(ierr); 863 sub_schurs->is_I = is_I; 864 ierr = PetscObjectReference((PetscObject)is_B);CHKERRQ(ierr); 865 sub_schurs->is_B = is_B; 866 ierr = PetscObjectReference((PetscObject)graph->l2gmap);CHKERRQ(ierr); 867 sub_schurs->l2gmap = graph->l2gmap; 868 ierr = PetscObjectReference((PetscObject)BtoNmap);CHKERRQ(ierr); 869 sub_schurs->BtoNmap = BtoNmap; 870 sub_schurs->n_subs_seq = n_local_sequential_problems; 871 sub_schurs->n_subs_par = n_local_parallel_problems; 872 sub_schurs->n_subs_seq_g = n_sequential_problems; 873 sub_schurs->n_subs_par_g = n_parallel_problems; 874 sub_schurs->n_subs = sub_schurs->n_subs_seq + sub_schurs->n_subs_par; 875 sub_schurs->is_subs = all_cc; 876 if (!sub_schurs->use_mumps) { /* for adaptive selection */ 877 for (i=0;i<sub_schurs->n_subs;i++) { 878 ierr = ISSort(sub_schurs->is_subs[i]);CHKERRQ(ierr); 879 } 880 } 881 sub_schurs->is_Ej_com = vertices; 882 sub_schurs->index_sequential = index_sequential; 883 sub_schurs->index_parallel = index_parallel; 884 sub_schurs->auxglobal_sequential = auxglobal_sequential; 885 sub_schurs->auxglobal_parallel = auxglobal_parallel; 886 887 888 /* allocate space for schur complements */ 889 ierr = PetscCalloc1(sub_schurs->n_subs,&sub_schurs->S_Ej);CHKERRQ(ierr); 890 sub_schurs->S_Ej_all = NULL; 891 sub_schurs->sum_S_Ej_all = NULL; 892 sub_schurs->sum_S_Ej_inv_all = NULL; 893 sub_schurs->sum_S_Ej_tilda_all = NULL; 894 sub_schurs->is_Ej_all = NULL; 895 896 /* free workspace */ 897 ierr = PetscFree(auxmapping);CHKERRQ(ierr); 898 ierr = PetscFree2(auxlocal_sequential,auxlocal_parallel);CHKERRQ(ierr); 899 900 PetscFunctionReturn(0); 901 } 902 903 #undef __FUNCT__ 904 #define __FUNCT__ "PCBDDCSubSchursCreate" 905 PetscErrorCode PCBDDCSubSchursCreate(PCBDDCSubSchurs *sub_schurs) 906 { 907 PCBDDCSubSchurs schurs_ctx; 908 PetscErrorCode ierr; 909 910 PetscFunctionBegin; 911 ierr = PetscNew(&schurs_ctx);CHKERRQ(ierr); 912 schurs_ctx->n_subs = 0; 913 *sub_schurs = schurs_ctx; 914 PetscFunctionReturn(0); 915 } 916 917 #undef __FUNCT__ 918 #define __FUNCT__ "PCBDDCSubSchursDestroy" 919 PetscErrorCode PCBDDCSubSchursDestroy(PCBDDCSubSchurs *sub_schurs) 920 { 921 PetscErrorCode ierr; 922 923 PetscFunctionBegin; 924 ierr = PCBDDCSubSchursReset(*sub_schurs);CHKERRQ(ierr); 925 ierr = PetscFree(*sub_schurs);CHKERRQ(ierr); 926 PetscFunctionReturn(0); 927 } 928 929 #undef __FUNCT__ 930 #define __FUNCT__ "PCBDDCSubSchursReset" 931 PetscErrorCode PCBDDCSubSchursReset(PCBDDCSubSchurs sub_schurs) 932 { 933 PetscInt i; 934 PetscErrorCode ierr; 935 936 PetscFunctionBegin; 937 ierr = MatDestroy(&sub_schurs->A);CHKERRQ(ierr); 938 ierr = MatDestroy(&sub_schurs->S);CHKERRQ(ierr); 939 ierr = ISDestroy(&sub_schurs->is_I);CHKERRQ(ierr); 940 ierr = ISDestroy(&sub_schurs->is_B);CHKERRQ(ierr); 941 ierr = ISLocalToGlobalMappingDestroy(&sub_schurs->l2gmap);CHKERRQ(ierr); 942 ierr = ISLocalToGlobalMappingDestroy(&sub_schurs->BtoNmap);CHKERRQ(ierr); 943 ierr = MatDestroy(&sub_schurs->S_Ej_all);CHKERRQ(ierr); 944 ierr = MatDestroy(&sub_schurs->sum_S_Ej_all);CHKERRQ(ierr); 945 ierr = MatDestroy(&sub_schurs->sum_S_Ej_inv_all);CHKERRQ(ierr); 946 ierr = MatDestroy(&sub_schurs->sum_S_Ej_tilda_all);CHKERRQ(ierr); 947 ierr = ISDestroy(&sub_schurs->is_Ej_all);CHKERRQ(ierr); 948 ierr = ISDestroy(&sub_schurs->is_Ej_com);CHKERRQ(ierr); 949 ierr = PetscBTDestroy(&sub_schurs->is_edge);CHKERRQ(ierr); 950 ierr = PetscBTDestroy(&sub_schurs->computed_Stilda_subs);CHKERRQ(ierr); 951 for (i=0;i<sub_schurs->n_subs;i++) { 952 ierr = ISDestroy(&sub_schurs->is_subs[i]);CHKERRQ(ierr); 953 ierr = MatDestroy(&sub_schurs->S_Ej[i]);CHKERRQ(ierr); 954 } 955 ierr = ISDestroy(&sub_schurs->is_I_layer);CHKERRQ(ierr); 956 if (sub_schurs->n_subs) { 957 ierr = PetscFree(sub_schurs->is_subs);CHKERRQ(ierr); 958 ierr = PetscFree(sub_schurs->S_Ej);CHKERRQ(ierr); 959 } 960 ierr = PetscFree2(sub_schurs->index_sequential,sub_schurs->index_parallel);CHKERRQ(ierr); 961 ierr = PetscFree(sub_schurs->auxglobal_sequential);CHKERRQ(ierr); 962 ierr = PetscFree(sub_schurs->auxglobal_parallel);CHKERRQ(ierr); 963 sub_schurs->n_subs = 0; 964 PetscFunctionReturn(0); 965 } 966 967 #undef __FUNCT__ 968 #define __FUNCT__ "PCBDDCAdjGetNextLayer_Private" 969 PETSC_STATIC_INLINE PetscErrorCode PCBDDCAdjGetNextLayer_Private(PetscInt* queue_tip,PetscInt n_prev,PetscBT touched,PetscInt* xadj,PetscInt* adjncy,PetscInt* n_added) 970 { 971 PetscInt i,j,n; 972 PetscErrorCode ierr; 973 974 PetscFunctionBegin; 975 n = 0; 976 for (i=-n_prev;i<0;i++) { 977 PetscInt start_dof = queue_tip[i]; 978 for (j=xadj[start_dof];j<xadj[start_dof+1];j++) { 979 PetscInt dof = adjncy[j]; 980 if (!PetscBTLookup(touched,dof)) { 981 ierr = PetscBTSet(touched,dof);CHKERRQ(ierr); 982 queue_tip[n] = dof; 983 n++; 984 } 985 } 986 } 987 *n_added = n; 988 PetscFunctionReturn(0); 989 } 990