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