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