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