1 /* TODOLIST 2 DofSplitting and DM attached to pc? 3 Change SetNeumannBoundaries to SetNeumannBoundariesLocal and provide new SetNeumannBoundaries (same Dirichlet) 4 change how to deal with the coarse problem (PCBDDCSetCoarseEnvironment): 5 - simplify coarse problem structure -> PCBDDC or PCREDUDANT, nothing else -> same comm for all levels? 6 - remove coarse enums and allow use of PCBDDCGetCoarseKSP 7 - remove metis dependency -> use MatPartitioning for multilevel -> Assemble serial adjacency in PCBDDCAnalyzeInterface? 8 code refactoring: 9 - pick up better names for static functions 10 change options structure: 11 - insert BDDC into MG framework? 12 provide other ops? Ask to developers 13 remove all unused printf 14 man pages 15 */ 16 17 /* ---------------------------------------------------------------------------------------------------------------------------------------------- 18 Implementation of BDDC preconditioner based on: 19 C. Dohrmann "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 20 ---------------------------------------------------------------------------------------------------------------------------------------------- */ 21 22 #include "bddc.h" /*I "petscpc.h" I*/ /* includes for fortran wrappers */ 23 #include "bddcprivate.h" 24 #include <petscblaslapack.h> 25 26 /* prototypes for static functions contained in bddc.c */ 27 static PetscErrorCode PCBDDCSetUseExactDirichlet(PC,PetscBool); 28 static PetscErrorCode PCBDDCSetLevel(PC,PetscInt); 29 static PetscErrorCode PCBDDCCoarseSetUp(PC); 30 static PetscErrorCode PCBDDCSetUpCoarseEnvironment(PC,PetscScalar*); 31 32 /* -------------------------------------------------------------------------- */ 33 #undef __FUNCT__ 34 #define __FUNCT__ "PCSetFromOptions_BDDC" 35 PetscErrorCode PCSetFromOptions_BDDC(PC pc) 36 { 37 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 ierr = PetscOptionsHead("BDDC options");CHKERRQ(ierr); 42 /* Verbose debugging of main data structures */ 43 ierr = PetscOptionsInt("-pc_bddc_check_level" ,"Verbose (debugging) output for PCBDDC" ,"none",pcbddc->dbg_flag ,&pcbddc->dbg_flag ,NULL);CHKERRQ(ierr); 44 /* Some customization for default primal space */ 45 ierr = PetscOptionsBool("-pc_bddc_vertices_only" ,"Use only vertices in coarse space (i.e. discard constraints)","none",pcbddc->vertices_flag ,&pcbddc->vertices_flag ,NULL);CHKERRQ(ierr); 46 ierr = PetscOptionsBool("-pc_bddc_constraints_only","Use only constraints in coarse space (i.e. discard vertices)","none",pcbddc->constraints_flag,&pcbddc->constraints_flag,NULL);CHKERRQ(ierr); 47 ierr = PetscOptionsBool("-pc_bddc_faces_only" ,"Use only faces among constraints of coarse space (i.e. discard edges)" ,"none",pcbddc->faces_flag ,&pcbddc->faces_flag ,NULL);CHKERRQ(ierr); 48 ierr = PetscOptionsBool("-pc_bddc_edges_only" ,"Use only edges among constraints of coarse space (i.e. discard faces)" ,"none",pcbddc->edges_flag ,&pcbddc->edges_flag ,NULL);CHKERRQ(ierr); 49 /* Coarse solver context */ 50 static const char * const avail_coarse_problems[] = {"sequential","replicated","parallel","multilevel","CoarseProblemType","PC_BDDC_",0}; /*order of choiches depends on ENUM defined in bddc.h */ 51 ierr = PetscOptionsEnum("-pc_bddc_coarse_problem_type","Set coarse problem type","none",avail_coarse_problems,(PetscEnum)pcbddc->coarse_problem_type,(PetscEnum*)&pcbddc->coarse_problem_type,NULL);CHKERRQ(ierr); 52 /* Two different application of BDDC to the whole set of dofs, internal and interface */ 53 ierr = PetscOptionsBool("-pc_bddc_switch_preconditioning_type","Switch between M_2 (default) and M_3 preconditioners (as defined by Dohrmann)","none",pcbddc->inexact_prec_type,&pcbddc->inexact_prec_type,NULL);CHKERRQ(ierr); 54 ierr = PetscOptionsBool("-pc_bddc_use_change_of_basis","Use change of basis approach for primal space","none",pcbddc->use_change_of_basis,&pcbddc->use_change_of_basis,NULL);CHKERRQ(ierr); 55 ierr = PetscOptionsBool("-pc_bddc_use_change_on_faces","Use change of basis approach for face constraints","none",pcbddc->use_change_on_faces,&pcbddc->use_change_on_faces,NULL);CHKERRQ(ierr); 56 if (!pcbddc->use_change_of_basis) { 57 pcbddc->use_change_on_faces = PETSC_FALSE; 58 } 59 ierr = PetscOptionsInt("-pc_bddc_coarsening_ratio","Set coarsening ratio used in multilevel coarsening","none",pcbddc->coarsening_ratio,&pcbddc->coarsening_ratio,NULL);CHKERRQ(ierr); 60 ierr = PetscOptionsInt("-pc_bddc_max_levels","Set maximum number of levels for multilevel","none",pcbddc->max_levels,&pcbddc->max_levels,NULL);CHKERRQ(ierr); 61 ierr = PetscOptionsBool("-pc_bddc_use_deluxe_scaling","Use deluxe scaling for BDDC","none",pcbddc->use_deluxe_scaling,&pcbddc->use_deluxe_scaling,NULL);CHKERRQ(ierr); 62 ierr = PetscOptionsTail();CHKERRQ(ierr); 63 PetscFunctionReturn(0); 64 } 65 /* -------------------------------------------------------------------------- */ 66 #undef __FUNCT__ 67 #define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS_BDDC" 68 static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices) 69 { 70 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 71 PetscErrorCode ierr; 72 73 PetscFunctionBegin; 74 ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr); 75 ierr = PetscObjectReference((PetscObject)PrimalVertices);CHKERRQ(ierr); 76 pcbddc->user_primal_vertices = PrimalVertices; 77 PetscFunctionReturn(0); 78 } 79 #undef __FUNCT__ 80 #define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS" 81 /*@ 82 PCBDDCSetPrimalVerticesLocalIS - Set user defined primal vertices in PCBDDC. 83 84 Not collective 85 86 Input Parameters: 87 + pc - the preconditioning context 88 - PrimalVertices - index sets of primal vertices in local numbering 89 90 Level: intermediate 91 92 Notes: 93 94 .seealso: PCBDDC 95 @*/ 96 PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices) 97 { 98 PetscErrorCode ierr; 99 100 PetscFunctionBegin; 101 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 102 PetscValidHeaderSpecific(PrimalVertices,IS_CLASSID,2); 103 ierr = PetscTryMethod(pc,"PCBDDCSetPrimalVerticesLocalIS_C",(PC,IS),(pc,PrimalVertices));CHKERRQ(ierr); 104 PetscFunctionReturn(0); 105 } 106 /* -------------------------------------------------------------------------- */ 107 #undef __FUNCT__ 108 #define __FUNCT__ "PCBDDCSetCoarseProblemType_BDDC" 109 static PetscErrorCode PCBDDCSetCoarseProblemType_BDDC(PC pc, CoarseProblemType CPT) 110 { 111 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 112 113 PetscFunctionBegin; 114 pcbddc->coarse_problem_type = CPT; 115 PetscFunctionReturn(0); 116 } 117 118 #undef __FUNCT__ 119 #define __FUNCT__ "PCBDDCSetCoarseProblemType" 120 /*@ 121 PCBDDCSetCoarseProblemType - Set coarse problem type in PCBDDC. 122 123 Not collective 124 125 Input Parameters: 126 + pc - the preconditioning context 127 - CoarseProblemType - pick a better name and explain what this is 128 129 Level: intermediate 130 131 Notes: 132 Not collective but all procs must call with same arguments. 133 134 .seealso: PCBDDC 135 @*/ 136 PetscErrorCode PCBDDCSetCoarseProblemType(PC pc, CoarseProblemType CPT) 137 { 138 PetscErrorCode ierr; 139 140 PetscFunctionBegin; 141 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 142 ierr = PetscTryMethod(pc,"PCBDDCSetCoarseProblemType_C",(PC,CoarseProblemType),(pc,CPT));CHKERRQ(ierr); 143 PetscFunctionReturn(0); 144 } 145 /* -------------------------------------------------------------------------- */ 146 #undef __FUNCT__ 147 #define __FUNCT__ "PCBDDCSetCoarseningRatio_BDDC" 148 static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc,PetscInt k) 149 { 150 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 151 152 PetscFunctionBegin; 153 pcbddc->coarsening_ratio=k; 154 PetscFunctionReturn(0); 155 } 156 157 #undef __FUNCT__ 158 #define __FUNCT__ "PCBDDCSetCoarseningRatio" 159 /*@ 160 PCBDDCSetCoarseningRatio - Set coarsening ratio used in multilevel coarsening 161 162 Logically collective on PC 163 164 Input Parameters: 165 + pc - the preconditioning context 166 - k - coarsening ratio 167 168 Approximatively k subdomains at the finer level will be aggregated into a single subdomain at the coarser level. 169 170 Level: intermediate 171 172 Notes: 173 174 .seealso: PCBDDC 175 @*/ 176 PetscErrorCode PCBDDCSetCoarseningRatio(PC pc,PetscInt k) 177 { 178 PetscErrorCode ierr; 179 180 PetscFunctionBegin; 181 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 182 ierr = PetscTryMethod(pc,"PCBDDCSetCoarseningRatio_C",(PC,PetscInt),(pc,k));CHKERRQ(ierr); 183 PetscFunctionReturn(0); 184 } 185 /* -------------------------------------------------------------------------- */ 186 187 #undef __FUNCT__ 188 #define __FUNCT__ "PCBDDCSetMaxLevels_BDDC" 189 static PetscErrorCode PCBDDCSetMaxLevels_BDDC(PC pc,PetscInt max_levels) 190 { 191 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 192 193 PetscFunctionBegin; 194 pcbddc->max_levels=max_levels; 195 PetscFunctionReturn(0); 196 } 197 198 #undef __FUNCT__ 199 #define __FUNCT__ "PCBDDCSetMaxLevels" 200 /*@ 201 PCBDDCSetMaxLevels - Sets the maximum number of levels within the multilevel approach. 202 203 Logically collective on PC 204 205 Input Parameters: 206 + pc - the preconditioning context 207 - max_levels - the maximum number of levels 208 209 Default value is 1, i.e. coarse problem will be solved inexactly with one application 210 of PCBDDC preconditioner if the multilevel approach is requested. 211 212 Level: intermediate 213 214 Notes: 215 216 .seealso: PCBDDC 217 @*/ 218 PetscErrorCode PCBDDCSetMaxLevels(PC pc,PetscInt max_levels) 219 { 220 PetscErrorCode ierr; 221 222 PetscFunctionBegin; 223 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 224 ierr = PetscTryMethod(pc,"PCBDDCSetMaxLevels_C",(PC,PetscInt),(pc,max_levels));CHKERRQ(ierr); 225 PetscFunctionReturn(0); 226 } 227 /* -------------------------------------------------------------------------- */ 228 229 #undef __FUNCT__ 230 #define __FUNCT__ "PCBDDCSetNullSpace_BDDC" 231 static PetscErrorCode PCBDDCSetNullSpace_BDDC(PC pc,MatNullSpace NullSpace) 232 { 233 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 234 PetscErrorCode ierr; 235 236 PetscFunctionBegin; 237 ierr = PetscObjectReference((PetscObject)NullSpace);CHKERRQ(ierr); 238 ierr = MatNullSpaceDestroy(&pcbddc->NullSpace);CHKERRQ(ierr); 239 pcbddc->NullSpace=NullSpace; 240 PetscFunctionReturn(0); 241 } 242 243 #undef __FUNCT__ 244 #define __FUNCT__ "PCBDDCSetNullSpace" 245 /*@ 246 PCBDDCSetNullSpace - Set NullSpace of global operator of BDDC preconditioned mat. 247 248 Logically collective on PC and MatNullSpace 249 250 Input Parameters: 251 + pc - the preconditioning context 252 - NullSpace - Null space of the linear operator to be preconditioned. 253 254 Level: intermediate 255 256 Notes: 257 258 .seealso: PCBDDC 259 @*/ 260 PetscErrorCode PCBDDCSetNullSpace(PC pc,MatNullSpace NullSpace) 261 { 262 PetscErrorCode ierr; 263 264 PetscFunctionBegin; 265 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 266 PetscValidHeaderSpecific(NullSpace,MAT_NULLSPACE_CLASSID,2); 267 ierr = PetscTryMethod(pc,"PCBDDCSetNullSpace_C",(PC,MatNullSpace),(pc,NullSpace));CHKERRQ(ierr); 268 PetscFunctionReturn(0); 269 } 270 /* -------------------------------------------------------------------------- */ 271 272 #undef __FUNCT__ 273 #define __FUNCT__ "PCBDDCSetDirichletBoundaries_BDDC" 274 static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc,IS DirichletBoundaries) 275 { 276 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 277 PetscErrorCode ierr; 278 279 PetscFunctionBegin; 280 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 281 ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr); 282 pcbddc->DirichletBoundaries=DirichletBoundaries; 283 PetscFunctionReturn(0); 284 } 285 286 #undef __FUNCT__ 287 #define __FUNCT__ "PCBDDCSetDirichletBoundaries" 288 /*@ 289 PCBDDCSetDirichletBoundaries - Set index set defining subdomain part (in local ordering) 290 of Dirichlet boundaries for the global problem. 291 292 Not collective 293 294 Input Parameters: 295 + pc - the preconditioning context 296 - DirichletBoundaries - sequential index set defining the subdomain part of Dirichlet boundaries (can be NULL) 297 298 Level: intermediate 299 300 Notes: 301 302 .seealso: PCBDDC 303 @*/ 304 PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc,IS DirichletBoundaries) 305 { 306 PetscErrorCode ierr; 307 308 PetscFunctionBegin; 309 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 310 PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2); 311 ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundaries_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr); 312 PetscFunctionReturn(0); 313 } 314 /* -------------------------------------------------------------------------- */ 315 316 #undef __FUNCT__ 317 #define __FUNCT__ "PCBDDCSetNeumannBoundaries_BDDC" 318 static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc,IS NeumannBoundaries) 319 { 320 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 321 PetscErrorCode ierr; 322 323 PetscFunctionBegin; 324 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 325 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 326 pcbddc->NeumannBoundaries=NeumannBoundaries; 327 PetscFunctionReturn(0); 328 } 329 330 #undef __FUNCT__ 331 #define __FUNCT__ "PCBDDCSetNeumannBoundaries" 332 /*@ 333 PCBDDCSetNeumannBoundaries - Set index set defining subdomain part (in local ordering) 334 of Neumann boundaries for the global problem. 335 336 Not collective 337 338 Input Parameters: 339 + pc - the preconditioning context 340 - NeumannBoundaries - sequential index set defining the subdomain part of Neumann boundaries (can be NULL) 341 342 Level: intermediate 343 344 Notes: 345 346 .seealso: PCBDDC 347 @*/ 348 PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc,IS NeumannBoundaries) 349 { 350 PetscErrorCode ierr; 351 352 PetscFunctionBegin; 353 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 354 PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2); 355 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundaries_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 356 PetscFunctionReturn(0); 357 } 358 /* -------------------------------------------------------------------------- */ 359 360 #undef __FUNCT__ 361 #define __FUNCT__ "PCBDDCGetDirichletBoundaries_BDDC" 362 static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc,IS *DirichletBoundaries) 363 { 364 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 365 366 PetscFunctionBegin; 367 *DirichletBoundaries = pcbddc->DirichletBoundaries; 368 PetscFunctionReturn(0); 369 } 370 371 #undef __FUNCT__ 372 #define __FUNCT__ "PCBDDCGetDirichletBoundaries" 373 /*@ 374 PCBDDCGetDirichletBoundaries - Get index set defining subdomain part (in local ordering) 375 of Dirichlet boundaries for the global problem. 376 377 Not collective 378 379 Input Parameters: 380 + pc - the preconditioning context 381 382 Output Parameters: 383 + DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries 384 385 Level: intermediate 386 387 Notes: 388 389 .seealso: PCBDDC 390 @*/ 391 PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc,IS *DirichletBoundaries) 392 { 393 PetscErrorCode ierr; 394 395 PetscFunctionBegin; 396 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 397 ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundaries_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr); 398 PetscFunctionReturn(0); 399 } 400 /* -------------------------------------------------------------------------- */ 401 402 #undef __FUNCT__ 403 #define __FUNCT__ "PCBDDCGetNeumannBoundaries_BDDC" 404 static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc,IS *NeumannBoundaries) 405 { 406 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 407 408 PetscFunctionBegin; 409 *NeumannBoundaries = pcbddc->NeumannBoundaries; 410 PetscFunctionReturn(0); 411 } 412 413 #undef __FUNCT__ 414 #define __FUNCT__ "PCBDDCGetNeumannBoundaries" 415 /*@ 416 PCBDDCGetNeumannBoundaries - Get index set defining subdomain part (in local ordering) 417 of Neumann boundaries for the global problem. 418 419 Not collective 420 421 Input Parameters: 422 + pc - the preconditioning context 423 424 Output Parameters: 425 + NeumannBoundaries - index set defining the subdomain part of Neumann boundaries 426 427 Level: intermediate 428 429 Notes: 430 431 .seealso: PCBDDC 432 @*/ 433 PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc,IS *NeumannBoundaries) 434 { 435 PetscErrorCode ierr; 436 437 PetscFunctionBegin; 438 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 439 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundaries_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 440 PetscFunctionReturn(0); 441 } 442 /* -------------------------------------------------------------------------- */ 443 444 #undef __FUNCT__ 445 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph_BDDC" 446 static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 447 { 448 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 449 PCBDDCGraph mat_graph = pcbddc->mat_graph; 450 PetscErrorCode ierr; 451 452 PetscFunctionBegin; 453 /* free old CSR */ 454 ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr); 455 /* TODO: PCBDDCGraphSetAdjacency */ 456 /* get CSR into graph structure */ 457 if (copymode == PETSC_COPY_VALUES) { 458 ierr = PetscMalloc((nvtxs+1)*sizeof(PetscInt),&mat_graph->xadj);CHKERRQ(ierr); 459 ierr = PetscMalloc(xadj[nvtxs]*sizeof(PetscInt),&mat_graph->adjncy);CHKERRQ(ierr); 460 ierr = PetscMemcpy(mat_graph->xadj,xadj,(nvtxs+1)*sizeof(PetscInt));CHKERRQ(ierr); 461 ierr = PetscMemcpy(mat_graph->adjncy,adjncy,xadj[nvtxs]*sizeof(PetscInt));CHKERRQ(ierr); 462 } else if (copymode == PETSC_OWN_POINTER) { 463 mat_graph->xadj = (PetscInt*)xadj; 464 mat_graph->adjncy = (PetscInt*)adjncy; 465 } 466 mat_graph->nvtxs_csr = nvtxs; 467 PetscFunctionReturn(0); 468 } 469 470 #undef __FUNCT__ 471 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph" 472 /*@ 473 PCBDDCSetLocalAdjacencyGraph - Set CSR graph of local matrix for use of PCBDDC. 474 475 Not collective 476 477 Input Parameters: 478 + pc - the preconditioning context 479 - nvtxs - number of local vertices of the graph 480 - xadj, adjncy - the CSR graph 481 - copymode - either PETSC_COPY_VALUES or PETSC_OWN_POINTER. In the former case the user must free the array passed in; 482 in the latter case, memory must be obtained with PetscMalloc. 483 484 Level: intermediate 485 486 Notes: 487 488 .seealso: PCBDDC 489 @*/ 490 PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc,PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 491 { 492 void (*f)(void) = 0; 493 PetscErrorCode ierr; 494 495 PetscFunctionBegin; 496 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 497 PetscValidIntPointer(xadj,3); 498 PetscValidIntPointer(xadj,4); 499 if (copymode != PETSC_COPY_VALUES && copymode != PETSC_OWN_POINTER) { 500 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported copy mode %d in %s\n",copymode,__FUNCT__); 501 } 502 ierr = PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));CHKERRQ(ierr); 503 /* free arrays if PCBDDC is not the PC type */ 504 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);CHKERRQ(ierr); 505 if (!f && copymode == PETSC_OWN_POINTER) { 506 ierr = PetscFree(xadj);CHKERRQ(ierr); 507 ierr = PetscFree(adjncy);CHKERRQ(ierr); 508 } 509 PetscFunctionReturn(0); 510 } 511 /* -------------------------------------------------------------------------- */ 512 513 #undef __FUNCT__ 514 #define __FUNCT__ "PCBDDCSetDofsSplitting_BDDC" 515 static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 516 { 517 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 518 PetscInt i; 519 PetscErrorCode ierr; 520 521 PetscFunctionBegin; 522 /* Destroy ISes if they were already set */ 523 for (i=0;i<pcbddc->n_ISForDofs;i++) { 524 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 525 } 526 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 527 /* allocate space then set */ 528 ierr = PetscMalloc(n_is*sizeof(IS),&pcbddc->ISForDofs);CHKERRQ(ierr); 529 for (i=0;i<n_is;i++) { 530 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 531 pcbddc->ISForDofs[i]=ISForDofs[i]; 532 } 533 pcbddc->n_ISForDofs=n_is; 534 PetscFunctionReturn(0); 535 } 536 537 #undef __FUNCT__ 538 #define __FUNCT__ "PCBDDCSetDofsSplitting" 539 /*@ 540 PCBDDCSetDofsSplitting - Set index sets defining fields of local mat. 541 542 Not collective 543 544 Input Parameters: 545 + pc - the preconditioning context 546 - n - number of index sets defining the fields 547 - IS[] - array of IS describing the fields 548 549 Level: intermediate 550 551 Notes: 552 553 .seealso: PCBDDC 554 @*/ 555 PetscErrorCode PCBDDCSetDofsSplitting(PC pc,PetscInt n_is, IS ISForDofs[]) 556 { 557 PetscErrorCode ierr; 558 559 PetscFunctionBegin; 560 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 561 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 562 PetscFunctionReturn(0); 563 } 564 /* -------------------------------------------------------------------------- */ 565 #undef __FUNCT__ 566 #define __FUNCT__ "PCPreSolve_BDDC" 567 /* -------------------------------------------------------------------------- */ 568 /* 569 PCPreSolve_BDDC - Changes the right hand side and (if necessary) the initial 570 guess if a transformation of basis approach has been selected. 571 572 Input Parameter: 573 + pc - the preconditioner contex 574 575 Application Interface Routine: PCPreSolve() 576 577 Notes: 578 The interface routine PCPreSolve() is not usually called directly by 579 the user, but instead is called by KSPSolve(). 580 */ 581 static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 582 { 583 PetscErrorCode ierr; 584 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 585 PC_IS *pcis = (PC_IS*)(pc->data); 586 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 587 Mat temp_mat; 588 IS dirIS; 589 PetscInt dirsize,i,*is_indices; 590 PetscScalar *array_x,*array_diagonal; 591 Vec used_vec; 592 PetscBool guess_nonzero; 593 594 PetscFunctionBegin; 595 if (x) { 596 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr); 597 used_vec = x; 598 } else { 599 ierr = PetscObjectReference((PetscObject)pcbddc->temp_solution);CHKERRQ(ierr); 600 used_vec = pcbddc->temp_solution; 601 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 602 } 603 /* hack into ksp data structure PCPreSolve comes earlier in src/ksp/ksp/interface/itfunc.c */ 604 if (ksp) { 605 ierr = KSPGetInitialGuessNonzero(ksp,&guess_nonzero);CHKERRQ(ierr); 606 if ( !guess_nonzero ) { 607 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 608 } 609 } 610 611 if (rhs) { 612 /* store the original rhs */ 613 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 614 615 /* Take into account zeroed rows -> change rhs and store solution removed */ 616 ierr = MatGetDiagonal(pc->pmat,pcis->vec1_global);CHKERRQ(ierr); 617 ierr = VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);CHKERRQ(ierr); 618 ierr = VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 619 ierr = VecScatterEnd(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 620 ierr = VecScatterBegin(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 621 ierr = VecScatterEnd(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 622 ierr = PCBDDCGetDirichletBoundaries(pc,&dirIS);CHKERRQ(ierr); 623 if (dirIS) { 624 ierr = ISGetSize(dirIS,&dirsize);CHKERRQ(ierr); 625 ierr = VecGetArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 626 ierr = VecGetArray(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 627 ierr = ISGetIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 628 for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]]; 629 ierr = ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 630 ierr = VecRestoreArray(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 631 ierr = VecRestoreArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 632 } 633 ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 634 ierr = VecScatterEnd(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 635 636 /* remove the computed solution from the rhs */ 637 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 638 ierr = MatMultAdd(pc->pmat,used_vec,rhs,rhs);CHKERRQ(ierr); 639 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 640 } 641 642 /* store partially computed solution and set initial guess */ 643 if (x) { 644 ierr = VecCopy(used_vec,pcbddc->temp_solution);CHKERRQ(ierr); 645 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 646 if (pcbddc->use_exact_dirichlet && !pcbddc->coarse_psi_B) { 647 ierr = VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 648 ierr = VecScatterEnd (pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 649 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 650 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 651 ierr = VecScatterEnd (pcis->global_to_D,pcis->vec2_D,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 652 if (ksp) { 653 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 654 } 655 } 656 } 657 658 if (pcbddc->use_change_of_basis) { 659 /* swap pointers for local matrices */ 660 temp_mat = matis->A; 661 matis->A = pcbddc->local_mat; 662 pcbddc->local_mat = temp_mat; 663 } 664 if (pcbddc->use_change_of_basis && rhs) { 665 /* Get local rhs and apply transformation of basis */ 666 ierr = VecScatterBegin(pcis->global_to_B,rhs,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 667 ierr = VecScatterEnd (pcis->global_to_B,rhs,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 668 /* from original basis to modified basis */ 669 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,pcis->vec1_B,pcis->vec2_B);CHKERRQ(ierr); 670 /* put back modified values into the global vec using INSERT_VALUES copy mode */ 671 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec2_B,rhs,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 672 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec2_B,rhs,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 673 } 674 if (ksp && pcbddc->NullSpace) { 675 ierr = MatNullSpaceRemove(pcbddc->NullSpace,used_vec);CHKERRQ(ierr); 676 ierr = MatNullSpaceRemove(pcbddc->NullSpace,rhs);CHKERRQ(ierr); 677 } 678 ierr = VecDestroy(&used_vec);CHKERRQ(ierr); 679 PetscFunctionReturn(0); 680 } 681 /* -------------------------------------------------------------------------- */ 682 #undef __FUNCT__ 683 #define __FUNCT__ "PCPostSolve_BDDC" 684 /* -------------------------------------------------------------------------- */ 685 /* 686 PCPostSolve_BDDC - Changes the computed solution if a transformation of basis 687 approach has been selected. Also, restores rhs to its original state. 688 689 Input Parameter: 690 + pc - the preconditioner contex 691 692 Application Interface Routine: PCPostSolve() 693 694 Notes: 695 The interface routine PCPostSolve() is not usually called directly by 696 the user, but instead is called by KSPSolve(). 697 */ 698 static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 699 { 700 PetscErrorCode ierr; 701 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 702 PC_IS *pcis = (PC_IS*)(pc->data); 703 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 704 Mat temp_mat; 705 706 PetscFunctionBegin; 707 if (pcbddc->use_change_of_basis) { 708 /* swap pointers for local matrices */ 709 temp_mat = matis->A; 710 matis->A = pcbddc->local_mat; 711 pcbddc->local_mat = temp_mat; 712 } 713 if (pcbddc->use_change_of_basis && x) { 714 /* Get Local boundary and apply transformation of basis to solution vector */ 715 ierr = VecScatterBegin(pcis->global_to_B,x,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 716 ierr = VecScatterEnd (pcis->global_to_B,x,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 717 /* from modified basis to original basis */ 718 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcis->vec1_B,pcis->vec2_B);CHKERRQ(ierr); 719 /* put back modified values into the global vec using INSERT_VALUES copy mode */ 720 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec2_B,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 721 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec2_B,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 722 } 723 /* add solution removed in presolve */ 724 if (x) { 725 ierr = VecAXPY(x,1.0,pcbddc->temp_solution);CHKERRQ(ierr); 726 } 727 /* restore rhs to its original state */ 728 if (rhs) { 729 ierr = VecCopy(pcbddc->original_rhs,rhs);CHKERRQ(ierr); 730 } 731 PetscFunctionReturn(0); 732 } 733 /* -------------------------------------------------------------------------- */ 734 #undef __FUNCT__ 735 #define __FUNCT__ "PCSetUp_BDDC" 736 /* -------------------------------------------------------------------------- */ 737 /* 738 PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner 739 by setting data structures and options. 740 741 Input Parameter: 742 + pc - the preconditioner context 743 744 Application Interface Routine: PCSetUp() 745 746 Notes: 747 The interface routine PCSetUp() is not usually called directly by 748 the user, but instead is called by PCApply() if necessary. 749 */ 750 PetscErrorCode PCSetUp_BDDC(PC pc) 751 { 752 PetscErrorCode ierr; 753 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 754 MatStructure flag; 755 PetscBool computeis,computetopography,computesolvers; 756 757 PetscFunctionBegin; 758 /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other nonoverlapping preconditioners */ 759 /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup 760 So, we set to pcnone the Neumann problem of pcis in order to avoid unneeded computation 761 Also, we decide to directly build the (same) Dirichlet problem */ 762 ierr = PetscOptionsSetValue("-is_localN_pc_type","none");CHKERRQ(ierr); 763 ierr = PetscOptionsSetValue("-is_localD_pc_type","none");CHKERRQ(ierr); 764 /* Get stdout for dbg */ 765 if (pcbddc->dbg_flag && !pcbddc->dbg_viewer) { 766 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&pcbddc->dbg_viewer);CHKERRQ(ierr); 767 ierr = PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr); 768 } 769 /* first attempt to split work */ 770 if (pc->setupcalled) { 771 computeis = PETSC_FALSE; 772 ierr = PCGetOperators(pc,NULL,NULL,&flag);CHKERRQ(ierr); 773 if (flag == SAME_PRECONDITIONER) { 774 computetopography = PETSC_FALSE; 775 computesolvers = PETSC_FALSE; 776 } else if (flag == SAME_NONZERO_PATTERN) { 777 computetopography = PETSC_FALSE; 778 computesolvers = PETSC_TRUE; 779 } else { /* DIFFERENT_NONZERO_PATTERN */ 780 computetopography = PETSC_TRUE; 781 computesolvers = PETSC_TRUE; 782 } 783 } else { 784 computeis = PETSC_TRUE; 785 computetopography = PETSC_TRUE; 786 computesolvers = PETSC_TRUE; 787 } 788 /* Set up all the "iterative substructuring" common block */ 789 if (computeis) { 790 ierr = PCISSetUp(pc);CHKERRQ(ierr); 791 } 792 /* Analyze interface and set up local constraint and change of basis matrices */ 793 if (computetopography) { 794 /* reset data */ 795 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 796 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 797 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 798 } 799 if (computesolvers) { 800 /* reset data */ 801 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 802 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 803 /* Create coarse and local stuffs used for evaluating action of preconditioner */ 804 ierr = PCBDDCCoarseSetUp(pc);CHKERRQ(ierr); 805 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 806 } 807 PetscFunctionReturn(0); 808 } 809 810 /* -------------------------------------------------------------------------- */ 811 /* 812 PCApply_BDDC - Applies the BDDC preconditioner to a vector. 813 814 Input Parameters: 815 . pc - the preconditioner context 816 . r - input vector (global) 817 818 Output Parameter: 819 . z - output vector (global) 820 821 Application Interface Routine: PCApply() 822 */ 823 #undef __FUNCT__ 824 #define __FUNCT__ "PCApply_BDDC" 825 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 826 { 827 PC_IS *pcis = (PC_IS*)(pc->data); 828 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 829 PetscErrorCode ierr; 830 const PetscScalar one = 1.0; 831 const PetscScalar m_one = -1.0; 832 const PetscScalar zero = 0.0; 833 834 /* This code is similar to that provided in nn.c for PCNN 835 NN interface preconditioner changed to BDDC 836 Added support for M_3 preconditioner in the reference article (code is active if pcbddc->inexact_prec_type = PETSC_TRUE) */ 837 838 PetscFunctionBegin; 839 if (!pcbddc->use_exact_dirichlet || pcbddc->coarse_psi_B) { 840 /* First Dirichlet solve */ 841 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 842 ierr = VecScatterEnd (pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 843 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 844 /* 845 Assembling right hand side for BDDC operator 846 - pcis->vec1_D for the Dirichlet part (if needed, i.e. prec_flag=PETSC_TRUE) 847 - pcis->vec1_B the interface part of the global vector z 848 */ 849 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 850 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 851 if (pcbddc->inexact_prec_type) { ierr = MatMultAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 852 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 853 ierr = VecCopy(r,z);CHKERRQ(ierr); 854 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 855 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 856 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 857 } else { 858 ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr); 859 ierr = VecSet(pcis->vec2_D,zero);CHKERRQ(ierr); 860 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 861 } 862 863 /* Apply interface preconditioner 864 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 865 ierr = PCBDDCApplyInterfacePreconditioner(pc);CHKERRQ(ierr); 866 867 /* Apply transpose of partition of unity operator */ 868 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 869 870 /* Second Dirichlet solve and assembling of output */ 871 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 872 ierr = VecScatterEnd (pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 873 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 874 if (pcbddc->inexact_prec_type) { ierr = MatMultAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 875 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcbddc->vec4_D);CHKERRQ(ierr); 876 ierr = VecScale(pcbddc->vec4_D,m_one);CHKERRQ(ierr); 877 if (pcbddc->inexact_prec_type) { ierr = VecAXPY (pcbddc->vec4_D,one,pcis->vec1_D);CHKERRQ(ierr); } 878 ierr = VecAXPY (pcis->vec2_D,one,pcbddc->vec4_D);CHKERRQ(ierr); 879 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 880 ierr = VecScatterEnd (pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 881 PetscFunctionReturn(0); 882 } 883 /* -------------------------------------------------------------------------- */ 884 885 #undef __FUNCT__ 886 #define __FUNCT__ "PCDestroy_BDDC" 887 PetscErrorCode PCDestroy_BDDC(PC pc) 888 { 889 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 890 PetscErrorCode ierr; 891 892 PetscFunctionBegin; 893 /* free data created by PCIS */ 894 ierr = PCISDestroy(pc);CHKERRQ(ierr); 895 /* free BDDC custom data */ 896 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 897 /* destroy objects related to topography */ 898 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 899 /* free allocated graph structure */ 900 ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr); 901 /* free data for scaling operator */ 902 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 903 /* free solvers stuff */ 904 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 905 /* remove functions */ 906 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 907 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 908 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetMaxLevels_C",NULL);CHKERRQ(ierr); 909 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",NULL);CHKERRQ(ierr); 910 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 911 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 912 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 913 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 914 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseProblemType_C",NULL);CHKERRQ(ierr); 915 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 916 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 917 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 918 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 919 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 920 /* Free the private data structure */ 921 ierr = PetscFree(pc->data);CHKERRQ(ierr); 922 PetscFunctionReturn(0); 923 } 924 /* -------------------------------------------------------------------------- */ 925 926 #undef __FUNCT__ 927 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC" 928 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 929 { 930 FETIDPMat_ctx mat_ctx; 931 PC_IS* pcis; 932 PC_BDDC* pcbddc; 933 PetscErrorCode ierr; 934 935 PetscFunctionBegin; 936 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 937 pcis = (PC_IS*)mat_ctx->pc->data; 938 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 939 940 /* change of basis for physical rhs if needed 941 It also changes the rhs in case of dirichlet boundaries */ 942 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,standard_rhs,NULL);CHKERRQ(ierr); 943 /* store vectors for computation of fetidp final solution */ 944 ierr = VecScatterBegin(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 945 ierr = VecScatterEnd(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 946 /* scale rhs since it should be unassembled */ 947 /* TODO use counter scaling? (also below) */ 948 ierr = VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 949 ierr = VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 950 /* Apply partition of unity */ 951 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 952 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 953 if (!pcbddc->inexact_prec_type) { 954 /* compute partially subassembled Schur complement right-hand side */ 955 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 956 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 957 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 958 ierr = VecSet(standard_rhs,0.0);CHKERRQ(ierr); 959 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 960 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 961 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 962 ierr = VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 963 ierr = VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 964 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 965 } 966 /* BDDC rhs */ 967 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 968 if (pcbddc->inexact_prec_type) { 969 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 970 } 971 /* apply BDDC */ 972 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc);CHKERRQ(ierr); 973 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 974 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 975 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 976 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 977 ierr = VecScatterEnd (mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 978 /* restore original rhs */ 979 ierr = VecCopy(pcbddc->original_rhs,standard_rhs);CHKERRQ(ierr); 980 PetscFunctionReturn(0); 981 } 982 983 #undef __FUNCT__ 984 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS" 985 /*@ 986 PCBDDCMatFETIDPGetRHS - Get rhs for FETIDP linear system. 987 988 Collective 989 990 Input Parameters: 991 + fetidp_mat - the FETIDP mat obtained by a call to PCBDDCCreateFETIDPOperators 992 + standard_rhs - the rhs of your linear system 993 994 Output Parameters: 995 + fetidp_flux_rhs - the rhs of the FETIDP linear system 996 997 Level: developer 998 999 Notes: 1000 1001 .seealso: PCBDDC 1002 @*/ 1003 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 1004 { 1005 FETIDPMat_ctx mat_ctx; 1006 PetscErrorCode ierr; 1007 1008 PetscFunctionBegin; 1009 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1010 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 1011 PetscFunctionReturn(0); 1012 } 1013 /* -------------------------------------------------------------------------- */ 1014 1015 #undef __FUNCT__ 1016 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC" 1017 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 1018 { 1019 FETIDPMat_ctx mat_ctx; 1020 PC_IS* pcis; 1021 PC_BDDC* pcbddc; 1022 PetscErrorCode ierr; 1023 1024 PetscFunctionBegin; 1025 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1026 pcis = (PC_IS*)mat_ctx->pc->data; 1027 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 1028 1029 /* apply B_delta^T */ 1030 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1031 ierr = VecScatterEnd (mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1032 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 1033 /* compute rhs for BDDC application */ 1034 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1035 if (pcbddc->inexact_prec_type) { 1036 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1037 } 1038 /* apply BDDC */ 1039 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc);CHKERRQ(ierr); 1040 /* put values into standard global vector */ 1041 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1042 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1043 if (!pcbddc->inexact_prec_type) { 1044 /* compute values into the interior if solved for the partially subassembled Schur complement */ 1045 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 1046 ierr = VecAXPY(mat_ctx->temp_solution_D,-1.0,pcis->vec1_D);CHKERRQ(ierr); 1047 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1048 } 1049 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1050 ierr = VecScatterEnd (pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1051 /* final change of basis if needed 1052 Is also sums the dirichlet part removed during RHS assembling */ 1053 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 1054 PetscFunctionReturn(0); 1055 } 1056 1057 #undef __FUNCT__ 1058 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution" 1059 /*@ 1060 PCBDDCMatFETIDPGetSolution - Get Solution for FETIDP linear system. 1061 1062 Collective 1063 1064 Input Parameters: 1065 + fetidp_mat - the FETIDP mat obtained by a call to PCBDDCCreateFETIDPOperators 1066 + fetidp_flux_sol - the solution of the FETIDP linear system 1067 1068 Output Parameters: 1069 + standard_sol - the solution on the global domain 1070 1071 Level: developer 1072 1073 Notes: 1074 1075 .seealso: PCBDDC 1076 @*/ 1077 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 1078 { 1079 FETIDPMat_ctx mat_ctx; 1080 PetscErrorCode ierr; 1081 1082 PetscFunctionBegin; 1083 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1084 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 1085 PetscFunctionReturn(0); 1086 } 1087 /* -------------------------------------------------------------------------- */ 1088 1089 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec); 1090 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat); 1091 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec); 1092 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC); 1093 1094 #undef __FUNCT__ 1095 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC" 1096 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 1097 { 1098 1099 FETIDPMat_ctx fetidpmat_ctx; 1100 Mat newmat; 1101 FETIDPPC_ctx fetidppc_ctx; 1102 PC newpc; 1103 MPI_Comm comm; 1104 PetscErrorCode ierr; 1105 1106 PetscFunctionBegin; 1107 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 1108 /* FETIDP linear matrix */ 1109 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 1110 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 1111 ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 1112 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 1113 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 1114 ierr = MatSetUp(newmat);CHKERRQ(ierr); 1115 /* FETIDP preconditioner */ 1116 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 1117 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 1118 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 1119 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 1120 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 1121 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 1122 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 1123 ierr = PCSetOperators(newpc,newmat,newmat,SAME_PRECONDITIONER);CHKERRQ(ierr); 1124 ierr = PCSetUp(newpc);CHKERRQ(ierr); 1125 /* return pointers for objects created */ 1126 *fetidp_mat=newmat; 1127 *fetidp_pc=newpc; 1128 PetscFunctionReturn(0); 1129 } 1130 1131 #undef __FUNCT__ 1132 #define __FUNCT__ "PCBDDCCreateFETIDPOperators" 1133 /*@ 1134 PCBDDCCreateFETIDPOperators - Create operators for FETIDP. 1135 1136 Collective 1137 1138 Input Parameters: 1139 + pc - the BDDC preconditioning context (setup must be already called) 1140 1141 Level: developer 1142 1143 Notes: 1144 1145 .seealso: PCBDDC 1146 @*/ 1147 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 1148 { 1149 PetscErrorCode ierr; 1150 1151 PetscFunctionBegin; 1152 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1153 if (pc->setupcalled) { 1154 ierr = PetscTryMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,Mat*,PC*),(pc,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 1155 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n"); 1156 PetscFunctionReturn(0); 1157 } 1158 /* -------------------------------------------------------------------------- */ 1159 /*MC 1160 PCBDDC - Balancing Domain Decomposition by Constraints. 1161 1162 Options Database Keys: 1163 . -pcbddc ??? - 1164 1165 Level: intermediate 1166 1167 Notes: The matrix used with this preconditioner must be of type MATIS 1168 1169 Unlike more 'conventional' interface preconditioners, this iterates over ALL the 1170 degrees of freedom, NOT just those on the interface (this allows the use of approximate solvers 1171 on the subdomains). 1172 1173 Options for the coarse grid preconditioner can be set with - 1174 Options for the Dirichlet subproblem can be set with - 1175 Options for the Neumann subproblem can be set with - 1176 1177 Contributed by Stefano Zampini 1178 1179 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 1180 M*/ 1181 1182 #undef __FUNCT__ 1183 #define __FUNCT__ "PCCreate_BDDC" 1184 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 1185 { 1186 PetscErrorCode ierr; 1187 PC_BDDC *pcbddc; 1188 1189 PetscFunctionBegin; 1190 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 1191 ierr = PetscNewLog(pc,PC_BDDC,&pcbddc);CHKERRQ(ierr); 1192 pc->data = (void*)pcbddc; 1193 1194 /* create PCIS data structure */ 1195 ierr = PCISCreate(pc);CHKERRQ(ierr); 1196 1197 /* BDDC specific */ 1198 pcbddc->user_primal_vertices = 0; 1199 pcbddc->NullSpace = 0; 1200 pcbddc->temp_solution = 0; 1201 pcbddc->original_rhs = 0; 1202 pcbddc->local_mat = 0; 1203 pcbddc->ChangeOfBasisMatrix = 0; 1204 pcbddc->use_change_of_basis = PETSC_TRUE; 1205 pcbddc->use_change_on_faces = PETSC_FALSE; 1206 pcbddc->coarse_vec = 0; 1207 pcbddc->coarse_rhs = 0; 1208 pcbddc->coarse_ksp = 0; 1209 pcbddc->coarse_phi_B = 0; 1210 pcbddc->coarse_phi_D = 0; 1211 pcbddc->coarse_psi_B = 0; 1212 pcbddc->coarse_psi_D = 0; 1213 pcbddc->vec1_P = 0; 1214 pcbddc->vec1_R = 0; 1215 pcbddc->vec2_R = 0; 1216 pcbddc->local_auxmat1 = 0; 1217 pcbddc->local_auxmat2 = 0; 1218 pcbddc->R_to_B = 0; 1219 pcbddc->R_to_D = 0; 1220 pcbddc->ksp_D = 0; 1221 pcbddc->ksp_R = 0; 1222 pcbddc->local_primal_indices = 0; 1223 pcbddc->inexact_prec_type = PETSC_FALSE; 1224 pcbddc->NeumannBoundaries = 0; 1225 pcbddc->ISForDofs = 0; 1226 pcbddc->ConstraintMatrix = 0; 1227 pcbddc->use_nnsp_true = PETSC_FALSE; 1228 pcbddc->local_primal_sizes = 0; 1229 pcbddc->local_primal_displacements = 0; 1230 pcbddc->coarse_loc_to_glob = 0; 1231 pcbddc->dbg_flag = 0; 1232 pcbddc->coarsening_ratio = 8; 1233 pcbddc->use_exact_dirichlet = PETSC_TRUE; 1234 pcbddc->current_level = 0; 1235 pcbddc->max_levels = 1; 1236 pcbddc->replicated_local_primal_indices = 0; 1237 pcbddc->replicated_local_primal_values = 0; 1238 1239 /* create local graph structure */ 1240 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 1241 1242 /* scaling */ 1243 pcbddc->use_deluxe_scaling = PETSC_FALSE; 1244 pcbddc->work_scaling = 0; 1245 1246 /* function pointers */ 1247 pc->ops->apply = PCApply_BDDC; 1248 pc->ops->applytranspose = 0; 1249 pc->ops->setup = PCSetUp_BDDC; 1250 pc->ops->destroy = PCDestroy_BDDC; 1251 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 1252 pc->ops->view = 0; 1253 pc->ops->applyrichardson = 0; 1254 pc->ops->applysymmetricleft = 0; 1255 pc->ops->applysymmetricright = 0; 1256 pc->ops->presolve = PCPreSolve_BDDC; 1257 pc->ops->postsolve = PCPostSolve_BDDC; 1258 1259 /* composing function */ 1260 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 1261 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 1262 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetMaxLevels_C",PCBDDCSetMaxLevels_BDDC);CHKERRQ(ierr); 1263 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",PCBDDCSetNullSpace_BDDC);CHKERRQ(ierr); 1264 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 1265 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 1266 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 1267 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 1268 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseProblemType_C",PCBDDCSetCoarseProblemType_BDDC);CHKERRQ(ierr); 1269 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 1270 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 1271 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 1272 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 1273 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 1274 PetscFunctionReturn(0); 1275 } 1276 1277 /* -------------------------------------------------------------------------- */ 1278 /* All static functions from now on */ 1279 /* -------------------------------------------------------------------------- */ 1280 1281 #undef __FUNCT__ 1282 #define __FUNCT__ "PCBDDCSetUseExactDirichlet" 1283 static PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool use) 1284 { 1285 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1286 1287 PetscFunctionBegin; 1288 pcbddc->use_exact_dirichlet=use; 1289 PetscFunctionReturn(0); 1290 } 1291 1292 #undef __FUNCT__ 1293 #define __FUNCT__ "PCBDDCSetLevel" 1294 static PetscErrorCode PCBDDCSetLevel(PC pc,PetscInt level) 1295 { 1296 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1297 1298 PetscFunctionBegin; 1299 pcbddc->current_level=level; 1300 PetscFunctionReturn(0); 1301 } 1302 1303 /* -------------------------------------------------------------------------- */ 1304 #undef __FUNCT__ 1305 #define __FUNCT__ "PCBDDCCoarseSetUp" 1306 static PetscErrorCode PCBDDCCoarseSetUp(PC pc) 1307 { 1308 PetscErrorCode ierr; 1309 1310 PC_IS* pcis = (PC_IS*)(pc->data); 1311 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 1312 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 1313 IS is_R_local; 1314 VecType impVecType; 1315 MatType impMatType; 1316 PetscInt n_R=0; 1317 PetscInt n_D=0; 1318 PetscInt n_B=0; 1319 PetscScalar zero=0.0; 1320 PetscScalar one=1.0; 1321 PetscScalar m_one=-1.0; 1322 PetscScalar* array; 1323 PetscScalar *coarse_submat_vals; 1324 PetscInt *idx_R_local; 1325 PetscReal *coarsefunctions_errors,*constraints_errors; 1326 /* auxiliary indices */ 1327 PetscInt i,j,k; 1328 /* for verbose output of bddc */ 1329 PetscViewer viewer=pcbddc->dbg_viewer; 1330 PetscInt dbg_flag=pcbddc->dbg_flag; 1331 /* for counting coarse dofs */ 1332 PetscInt n_vertices,n_constraints; 1333 PetscInt size_of_constraint; 1334 PetscInt *row_cmat_indices; 1335 PetscScalar *row_cmat_values; 1336 PetscInt *vertices; 1337 1338 PetscFunctionBegin; 1339 /* Set Non-overlapping dimensions */ 1340 n_B = pcis->n_B; n_D = pcis->n - n_B; 1341 1342 /* transform local matrices if needed */ 1343 if (pcbddc->use_change_of_basis) { 1344 Mat change_mat_all; 1345 PetscInt *nnz,*is_indices,*temp_indices; 1346 1347 ierr = PetscMalloc(pcis->n*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 1348 ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1349 for (i=0;i<n_D;i++) nnz[is_indices[i]] = 1; 1350 ierr = ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1351 ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1352 k=1; 1353 for (i=0;i<n_B;i++) { 1354 ierr = MatGetRow(pcbddc->ChangeOfBasisMatrix,i,&j,NULL,NULL);CHKERRQ(ierr); 1355 nnz[is_indices[i]]=j; 1356 if (k < j) k = j; 1357 ierr = MatRestoreRow(pcbddc->ChangeOfBasisMatrix,i,&j,NULL,NULL);CHKERRQ(ierr); 1358 } 1359 ierr = ISRestoreIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1360 /* assemble change of basis matrix on the whole set of local dofs */ 1361 ierr = PetscMalloc(k*sizeof(PetscInt),&temp_indices);CHKERRQ(ierr); 1362 ierr = MatCreate(PETSC_COMM_SELF,&change_mat_all);CHKERRQ(ierr); 1363 ierr = MatSetSizes(change_mat_all,pcis->n,pcis->n,pcis->n,pcis->n);CHKERRQ(ierr); 1364 ierr = MatSetType(change_mat_all,MATSEQAIJ);CHKERRQ(ierr); 1365 ierr = MatSeqAIJSetPreallocation(change_mat_all,0,nnz);CHKERRQ(ierr); 1366 ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1367 for (i=0;i<n_D;i++) { 1368 ierr = MatSetValue(change_mat_all,is_indices[i],is_indices[i],1.0,INSERT_VALUES);CHKERRQ(ierr); 1369 } 1370 ierr = ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1371 ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1372 for (i=0;i<n_B;i++) { 1373 ierr = MatGetRow(pcbddc->ChangeOfBasisMatrix,i,&j,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);CHKERRQ(ierr); 1374 for (k=0; k<j; k++) temp_indices[k]=is_indices[row_cmat_indices[k]]; 1375 ierr = MatSetValues(change_mat_all,1,&is_indices[i],j,temp_indices,row_cmat_values,INSERT_VALUES);CHKERRQ(ierr); 1376 ierr = MatRestoreRow(pcbddc->ChangeOfBasisMatrix,i,&j,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);CHKERRQ(ierr); 1377 } 1378 ierr = MatAssemblyBegin(change_mat_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1379 ierr = MatAssemblyEnd(change_mat_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1380 /* TODO: HOW TO WORK WITH BAIJ? PtAP not provided */ 1381 ierr = MatGetBlockSize(matis->A,&i);CHKERRQ(ierr); 1382 if (i==1) { 1383 ierr = MatPtAP(matis->A,change_mat_all,MAT_INITIAL_MATRIX,1.0,&pcbddc->local_mat);CHKERRQ(ierr); 1384 } else { 1385 Mat work_mat; 1386 ierr = MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&work_mat);CHKERRQ(ierr); 1387 ierr = MatPtAP(work_mat,change_mat_all,MAT_INITIAL_MATRIX,1.0,&pcbddc->local_mat);CHKERRQ(ierr); 1388 ierr = MatDestroy(&work_mat);CHKERRQ(ierr); 1389 } 1390 ierr = MatDestroy(&change_mat_all);CHKERRQ(ierr); 1391 ierr = PetscFree(nnz);CHKERRQ(ierr); 1392 ierr = PetscFree(temp_indices);CHKERRQ(ierr); 1393 } else { 1394 /* without change of basis, the local matrix is unchanged */ 1395 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1396 pcbddc->local_mat = matis->A; 1397 } 1398 /* need to rebuild PCIS matrices during SNES or TS -> TODO move this to PCIS code */ 1399 ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr); 1400 ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr); 1401 ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr); 1402 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);CHKERRQ(ierr); 1403 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);CHKERRQ(ierr); 1404 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);CHKERRQ(ierr); 1405 /* Change global null space passed in by the user if change of basis has been requested */ 1406 if (pcbddc->NullSpace && pcbddc->use_change_of_basis) { 1407 ierr = PCBDDCNullSpaceAdaptGlobal(pc);CHKERRQ(ierr); 1408 } 1409 1410 /* Set types for local objects needed by BDDC precondtioner */ 1411 impMatType = MATSEQDENSE; 1412 impVecType = VECSEQ; 1413 /* get vertex indices from constraint matrix */ 1414 ierr = PCBDDCGetPrimalVerticesLocalIdx(pc,&n_vertices,&vertices);CHKERRQ(ierr); 1415 /* Set number of constraints */ 1416 n_constraints = pcbddc->local_primal_size-n_vertices; 1417 /* Dohrmann's notation: dofs splitted in R (Remaining: all dofs but the vertices) and V (Vertices) */ 1418 ierr = VecSet(pcis->vec1_N,one);CHKERRQ(ierr); 1419 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1420 for (i=0;i<n_vertices;i++) array[vertices[i]] = zero; 1421 ierr = PetscMalloc((pcis->n-n_vertices)*sizeof(PetscInt),&idx_R_local);CHKERRQ(ierr); 1422 for (i=0, n_R=0; i<pcis->n; i++) { 1423 if (array[i] == one) { 1424 idx_R_local[n_R] = i; 1425 n_R++; 1426 } 1427 } 1428 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1429 ierr = PetscFree(vertices);CHKERRQ(ierr); 1430 if (dbg_flag) { 1431 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 1432 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1433 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d local dimensions\n",PetscGlobalRank);CHKERRQ(ierr); 1434 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local_size = %d, dirichlet_size = %d, boundary_size = %d\n",pcis->n,n_D,n_B);CHKERRQ(ierr); 1435 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"r_size = %d, v_size = %d, constraints = %d, local_primal_size = %d\n",n_R,n_vertices,n_constraints,pcbddc->local_primal_size);CHKERRQ(ierr); 1436 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"pcbddc->n_vertices = %d, pcbddc->n_constraints = %d\n",pcbddc->n_vertices,pcbddc->n_constraints);CHKERRQ(ierr); 1437 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1438 } 1439 1440 /* Allocate needed vectors */ 1441 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 1442 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);CHKERRQ(ierr); 1443 ierr = VecDuplicate(pcis->vec1_D,&pcbddc->vec4_D);CHKERRQ(ierr); 1444 ierr = VecCreate(PETSC_COMM_SELF,&pcbddc->vec1_R);CHKERRQ(ierr); 1445 ierr = VecSetSizes(pcbddc->vec1_R,n_R,n_R);CHKERRQ(ierr); 1446 ierr = VecSetType(pcbddc->vec1_R,impVecType);CHKERRQ(ierr); 1447 ierr = VecDuplicate(pcbddc->vec1_R,&pcbddc->vec2_R);CHKERRQ(ierr); 1448 ierr = VecCreate(PETSC_COMM_SELF,&pcbddc->vec1_P);CHKERRQ(ierr); 1449 ierr = VecSetSizes(pcbddc->vec1_P,pcbddc->local_primal_size,pcbddc->local_primal_size);CHKERRQ(ierr); 1450 ierr = VecSetType(pcbddc->vec1_P,impVecType);CHKERRQ(ierr); 1451 1452 /* Creating some index sets needed */ 1453 /* For submatrices */ 1454 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_R,idx_R_local,PETSC_OWN_POINTER,&is_R_local);CHKERRQ(ierr); 1455 1456 /* For VecScatters pcbddc->R_to_B and (optionally) pcbddc->R_to_D */ 1457 { 1458 IS is_aux1,is_aux2; 1459 PetscInt *aux_array1; 1460 PetscInt *aux_array2; 1461 PetscInt *idx_I_local; 1462 1463 ierr = PetscMalloc((pcis->n_B-n_vertices)*sizeof(PetscInt),&aux_array1);CHKERRQ(ierr); 1464 ierr = PetscMalloc((pcis->n_B-n_vertices)*sizeof(PetscInt),&aux_array2);CHKERRQ(ierr); 1465 1466 ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&idx_I_local);CHKERRQ(ierr); 1467 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1468 for (i=0; i<n_D; i++) array[idx_I_local[i]] = 0; 1469 ierr = ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&idx_I_local);CHKERRQ(ierr); 1470 for (i=0, j=0; i<n_R; i++) { 1471 if (array[idx_R_local[i]] == one) { 1472 aux_array1[j] = i; 1473 j++; 1474 } 1475 } 1476 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1477 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,aux_array1,PETSC_COPY_VALUES,&is_aux1);CHKERRQ(ierr); 1478 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1479 ierr = VecScatterEnd (pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1480 ierr = VecGetArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1481 for (i=0, j=0; i<n_B; i++) { 1482 if (array[i] == one) { 1483 aux_array2[j] = i; j++; 1484 } 1485 } 1486 ierr = VecRestoreArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1487 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,aux_array2,PETSC_COPY_VALUES,&is_aux2);CHKERRQ(ierr); 1488 ierr = VecScatterCreate(pcbddc->vec1_R,is_aux1,pcis->vec1_B,is_aux2,&pcbddc->R_to_B);CHKERRQ(ierr); 1489 ierr = PetscFree(aux_array1);CHKERRQ(ierr); 1490 ierr = PetscFree(aux_array2);CHKERRQ(ierr); 1491 ierr = ISDestroy(&is_aux1);CHKERRQ(ierr); 1492 ierr = ISDestroy(&is_aux2);CHKERRQ(ierr); 1493 1494 if (pcbddc->inexact_prec_type || dbg_flag ) { 1495 ierr = PetscMalloc(n_D*sizeof(PetscInt),&aux_array1);CHKERRQ(ierr); 1496 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1497 for (i=0, j=0; i<n_R; i++) { 1498 if (array[idx_R_local[i]] == zero) { 1499 aux_array1[j] = i; 1500 j++; 1501 } 1502 } 1503 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1504 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,aux_array1,PETSC_COPY_VALUES,&is_aux1);CHKERRQ(ierr); 1505 ierr = VecScatterCreate(pcbddc->vec1_R,is_aux1,pcis->vec1_D,(IS)0,&pcbddc->R_to_D);CHKERRQ(ierr); 1506 ierr = PetscFree(aux_array1);CHKERRQ(ierr); 1507 ierr = ISDestroy(&is_aux1);CHKERRQ(ierr); 1508 } 1509 } 1510 1511 /* Creating PC contexts for local Dirichlet and Neumann problems */ 1512 { 1513 Mat A_RR; 1514 PC pc_temp; 1515 MatStructure matstruct; 1516 /* Matrix for Dirichlet problem is A_II */ 1517 /* HACK (TODO) A_II can be changed between nonlinear iterations */ 1518 if (pc->setupcalled) { /* we dont need to rebuild dirichlet problem the first time we build BDDC */ 1519 ierr = PCGetOperators(pc,NULL,NULL,&matstruct);CHKERRQ(ierr); 1520 if (matstruct == SAME_NONZERO_PATTERN) { 1521 ierr = MatGetSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,MAT_REUSE_MATRIX,&pcis->A_II);CHKERRQ(ierr); 1522 } else { 1523 ierr = MatDestroy(&pcis->A_II);CHKERRQ(ierr); 1524 ierr = MatGetSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_II);CHKERRQ(ierr); 1525 } 1526 } 1527 ierr = KSPCreate(PETSC_COMM_SELF,&pcbddc->ksp_D);CHKERRQ(ierr); 1528 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_D,(PetscObject)pc,1);CHKERRQ(ierr); 1529 ierr = KSPSetOperators(pcbddc->ksp_D,pcis->A_II,pcis->A_II,SAME_PRECONDITIONER);CHKERRQ(ierr); 1530 ierr = KSPSetType(pcbddc->ksp_D,KSPPREONLY);CHKERRQ(ierr); 1531 ierr = KSPSetOptionsPrefix(pcbddc->ksp_D,"dirichlet_");CHKERRQ(ierr); 1532 /* default */ 1533 ierr = KSPGetPC(pcbddc->ksp_D,&pc_temp);CHKERRQ(ierr); 1534 ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr); 1535 /* Allow user's customization */ 1536 ierr = KSPSetFromOptions(pcbddc->ksp_D);CHKERRQ(ierr); 1537 /* umfpack interface has a bug when matrix dimension is zero */ 1538 if (!n_D) { 1539 ierr = PCSetType(pc_temp,PCNONE);CHKERRQ(ierr); 1540 } 1541 /* Set Up KSP for Dirichlet problem of BDDC */ 1542 ierr = KSPSetUp(pcbddc->ksp_D);CHKERRQ(ierr); 1543 /* set ksp_D into pcis data */ 1544 ierr = KSPDestroy(&pcis->ksp_D);CHKERRQ(ierr); 1545 ierr = PetscObjectReference((PetscObject)pcbddc->ksp_D);CHKERRQ(ierr); 1546 pcis->ksp_D = pcbddc->ksp_D; 1547 /* Matrix for Neumann problem is A_RR -> we need to create it */ 1548 ierr = MatGetSubMatrix(pcbddc->local_mat,is_R_local,is_R_local,MAT_INITIAL_MATRIX,&A_RR);CHKERRQ(ierr); 1549 ierr = KSPCreate(PETSC_COMM_SELF,&pcbddc->ksp_R);CHKERRQ(ierr); 1550 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_R,(PetscObject)pc,1);CHKERRQ(ierr); 1551 ierr = KSPSetOperators(pcbddc->ksp_R,A_RR,A_RR,SAME_PRECONDITIONER);CHKERRQ(ierr); 1552 ierr = KSPSetType(pcbddc->ksp_R,KSPPREONLY);CHKERRQ(ierr); 1553 ierr = KSPSetOptionsPrefix(pcbddc->ksp_R,"neumann_");CHKERRQ(ierr); 1554 /* default */ 1555 ierr = KSPGetPC(pcbddc->ksp_R,&pc_temp);CHKERRQ(ierr); 1556 ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr); 1557 /* Allow user's customization */ 1558 ierr = KSPSetFromOptions(pcbddc->ksp_R);CHKERRQ(ierr); 1559 /* umfpack interface has a bug when matrix dimension is zero */ 1560 if (!n_R) { 1561 ierr = PCSetType(pc_temp,PCNONE);CHKERRQ(ierr); 1562 } 1563 /* Set Up KSP for Neumann problem of BDDC */ 1564 ierr = KSPSetUp(pcbddc->ksp_R);CHKERRQ(ierr); 1565 /* check Dirichlet and Neumann solvers and adapt them if a nullspace correction is needed */ 1566 { 1567 Vec temp_vec; 1568 PetscReal value; 1569 PetscInt use_exact,use_exact_reduced; 1570 1571 ierr = VecDuplicate(pcis->vec1_D,&temp_vec);CHKERRQ(ierr); 1572 ierr = VecSetRandom(pcis->vec1_D,NULL);CHKERRQ(ierr); 1573 ierr = MatMult(pcis->A_II,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1574 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec2_D,temp_vec);CHKERRQ(ierr); 1575 ierr = VecAXPY(temp_vec,m_one,pcis->vec1_D);CHKERRQ(ierr); 1576 ierr = VecNorm(temp_vec,NORM_INFINITY,&value);CHKERRQ(ierr); 1577 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 1578 use_exact = 1; 1579 if (PetscAbsReal(value) > 1.e-4) use_exact = 0; 1580 ierr = MPI_Allreduce(&use_exact,&use_exact_reduced,1,MPIU_INT,MPI_LAND,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 1581 ierr = PCBDDCSetUseExactDirichlet(pc,(PetscBool)use_exact_reduced);CHKERRQ(ierr); 1582 if (dbg_flag) { 1583 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1584 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 1585 ierr = PetscViewerASCIIPrintf(viewer,"Checking solution of Dirichlet and Neumann problems\n");CHKERRQ(ierr); 1586 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d infinity error for Dirichlet solve = % 1.14e \n",PetscGlobalRank,value);CHKERRQ(ierr); 1587 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1588 } 1589 if (n_D && pcbddc->NullSpace && !use_exact_reduced && !pcbddc->inexact_prec_type) { 1590 ierr = PCBDDCNullSpaceAssembleCorrection(pc,pcis->is_I_local); 1591 } 1592 ierr = VecDuplicate(pcbddc->vec1_R,&temp_vec);CHKERRQ(ierr); 1593 ierr = VecSetRandom(pcbddc->vec1_R,NULL);CHKERRQ(ierr); 1594 ierr = MatMult(A_RR,pcbddc->vec1_R,pcbddc->vec2_R);CHKERRQ(ierr); 1595 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec2_R,temp_vec);CHKERRQ(ierr); 1596 ierr = VecAXPY(temp_vec,m_one,pcbddc->vec1_R);CHKERRQ(ierr); 1597 ierr = VecNorm(temp_vec,NORM_INFINITY,&value);CHKERRQ(ierr); 1598 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 1599 use_exact = 1; 1600 if (PetscAbsReal(value) > 1.e-4) use_exact = 0; 1601 ierr = MPI_Allreduce(&use_exact,&use_exact_reduced,1,MPIU_INT,MPI_LAND,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 1602 if (dbg_flag) { 1603 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d infinity error for Neumann solve = % 1.14e \n",PetscGlobalRank,value);CHKERRQ(ierr); 1604 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1605 } 1606 if (n_R && pcbddc->NullSpace && !use_exact_reduced) { /* is it the right logic? */ 1607 ierr = PCBDDCNullSpaceAssembleCorrection(pc,is_R_local); 1608 } 1609 } 1610 /* free Neumann problem's matrix */ 1611 ierr = MatDestroy(&A_RR);CHKERRQ(ierr); 1612 } 1613 1614 /* Assemble all remaining stuff needed to apply BDDC */ 1615 { 1616 Mat A_RV,A_VR,A_VV; 1617 Mat M1; 1618 Mat C_CR; 1619 Mat AUXMAT; 1620 Vec vec1_C; 1621 Vec vec2_C; 1622 Vec vec1_V; 1623 Vec vec2_V; 1624 IS is_C_local,is_V_local,is_aux1; 1625 ISLocalToGlobalMapping BtoNmap; 1626 PetscInt *nnz; 1627 PetscInt *idx_V_B; 1628 PetscInt *auxindices; 1629 PetscInt index; 1630 PetscScalar* array2; 1631 MatFactorInfo matinfo; 1632 PetscBool setsym=PETSC_FALSE,issym=PETSC_FALSE; 1633 1634 /* Allocating some extra storage just to be safe */ 1635 ierr = PetscMalloc (pcis->n*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 1636 ierr = PetscMalloc (pcis->n*sizeof(PetscInt),&auxindices);CHKERRQ(ierr); 1637 for (i=0;i<pcis->n;i++) auxindices[i]=i; 1638 1639 ierr = PCBDDCGetPrimalVerticesLocalIdx(pc,&n_vertices,&vertices);CHKERRQ(ierr); 1640 /* vertices in boundary numbering */ 1641 ierr = PetscMalloc(n_vertices*sizeof(PetscInt),&idx_V_B);CHKERRQ(ierr); 1642 ierr = ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&BtoNmap);CHKERRQ(ierr); 1643 ierr = ISGlobalToLocalMappingApply(BtoNmap,IS_GTOLM_DROP,n_vertices,vertices,&i,idx_V_B);CHKERRQ(ierr); 1644 ierr = ISLocalToGlobalMappingDestroy(&BtoNmap);CHKERRQ(ierr); 1645 if (i != n_vertices) { 1646 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Error in boundary numbering for BDDC vertices! %d != %d\n",n_vertices,i); 1647 } 1648 1649 /* some work vectors on vertices and/or constraints */ 1650 if (n_vertices) { 1651 ierr = VecCreate(PETSC_COMM_SELF,&vec1_V);CHKERRQ(ierr); 1652 ierr = VecSetSizes(vec1_V,n_vertices,n_vertices);CHKERRQ(ierr); 1653 ierr = VecSetType(vec1_V,impVecType);CHKERRQ(ierr); 1654 ierr = VecDuplicate(vec1_V,&vec2_V);CHKERRQ(ierr); 1655 } 1656 if (n_constraints) { 1657 ierr = VecCreate(PETSC_COMM_SELF,&vec1_C);CHKERRQ(ierr); 1658 ierr = VecSetSizes(vec1_C,n_constraints,n_constraints);CHKERRQ(ierr); 1659 ierr = VecSetType(vec1_C,impVecType);CHKERRQ(ierr); 1660 ierr = VecDuplicate(vec1_C,&vec2_C);CHKERRQ(ierr); 1661 ierr = VecDuplicate(vec1_C,&pcbddc->vec1_C);CHKERRQ(ierr); 1662 } 1663 /* Precompute stuffs needed for preprocessing and application of BDDC*/ 1664 if (n_constraints) { 1665 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->local_auxmat2);CHKERRQ(ierr); 1666 ierr = MatSetSizes(pcbddc->local_auxmat2,n_R,n_constraints,n_R,n_constraints);CHKERRQ(ierr); 1667 ierr = MatSetType(pcbddc->local_auxmat2,impMatType);CHKERRQ(ierr); 1668 ierr = MatSeqDenseSetPreallocation(pcbddc->local_auxmat2,NULL);CHKERRQ(ierr); 1669 1670 /* Create Constraint matrix on R nodes: C_{CR} */ 1671 ierr = ISCreateStride(PETSC_COMM_SELF,n_constraints,n_vertices,1,&is_C_local);CHKERRQ(ierr); 1672 ierr = MatGetSubMatrix(pcbddc->ConstraintMatrix,is_C_local,is_R_local,MAT_INITIAL_MATRIX,&C_CR);CHKERRQ(ierr); 1673 ierr = ISDestroy(&is_C_local);CHKERRQ(ierr); 1674 1675 /* Assemble local_auxmat2 = - A_{RR}^{-1} C^T_{CR} needed by BDDC application */ 1676 for (i=0;i<n_constraints;i++) { 1677 ierr = VecSet(pcbddc->vec1_R,zero);CHKERRQ(ierr); 1678 /* Get row of constraint matrix in R numbering */ 1679 ierr = VecGetArray(pcbddc->vec1_R,&array);CHKERRQ(ierr); 1680 ierr = MatGetRow(C_CR,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);CHKERRQ(ierr); 1681 for (j=0;j<size_of_constraint;j++) array[row_cmat_indices[j]] = -row_cmat_values[j]; 1682 ierr = MatRestoreRow(C_CR,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);CHKERRQ(ierr); 1683 ierr = VecRestoreArray(pcbddc->vec1_R,&array);CHKERRQ(ierr); 1684 1685 /* Solve for row of constraint matrix in R numbering */ 1686 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec2_R);CHKERRQ(ierr); 1687 1688 /* Set values */ 1689 ierr = VecGetArray(pcbddc->vec2_R,&array);CHKERRQ(ierr); 1690 ierr = MatSetValues(pcbddc->local_auxmat2,n_R,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1691 ierr = VecRestoreArray(pcbddc->vec2_R,&array);CHKERRQ(ierr); 1692 } 1693 ierr = MatAssemblyBegin(pcbddc->local_auxmat2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1694 ierr = MatAssemblyEnd(pcbddc->local_auxmat2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1695 1696 /* Assemble AUXMAT = ( LUFactor )( -C_{CR} A_{RR}^{-1} C^T_{CR} )^{-1} */ 1697 ierr = MatMatMult(C_CR,pcbddc->local_auxmat2,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AUXMAT);CHKERRQ(ierr); 1698 ierr = MatFactorInfoInitialize(&matinfo);CHKERRQ(ierr); 1699 ierr = ISCreateStride(PETSC_COMM_SELF,n_constraints,0,1,&is_aux1);CHKERRQ(ierr); 1700 ierr = MatLUFactor(AUXMAT,is_aux1,is_aux1,&matinfo);CHKERRQ(ierr); 1701 ierr = ISDestroy(&is_aux1);CHKERRQ(ierr); 1702 1703 /* Assemble explicitly M1 = ( C_{CR} A_{RR}^{-1} C^T_{CR} )^{-1} needed in preproc */ 1704 ierr = MatCreate(PETSC_COMM_SELF,&M1);CHKERRQ(ierr); 1705 ierr = MatSetSizes(M1,n_constraints,n_constraints,n_constraints,n_constraints);CHKERRQ(ierr); 1706 ierr = MatSetType(M1,impMatType);CHKERRQ(ierr); 1707 ierr = MatSeqDenseSetPreallocation(M1,NULL);CHKERRQ(ierr); 1708 for (i=0;i<n_constraints;i++) { 1709 ierr = VecSet(vec1_C,zero);CHKERRQ(ierr); 1710 ierr = VecSetValue(vec1_C,i,one,INSERT_VALUES);CHKERRQ(ierr); 1711 ierr = VecAssemblyBegin(vec1_C);CHKERRQ(ierr); 1712 ierr = VecAssemblyEnd(vec1_C);CHKERRQ(ierr); 1713 ierr = MatSolve(AUXMAT,vec1_C,vec2_C);CHKERRQ(ierr); 1714 ierr = VecScale(vec2_C,m_one);CHKERRQ(ierr); 1715 ierr = VecGetArray(vec2_C,&array);CHKERRQ(ierr); 1716 ierr = MatSetValues(M1,n_constraints,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1717 ierr = VecRestoreArray(vec2_C,&array);CHKERRQ(ierr); 1718 } 1719 ierr = MatAssemblyBegin(M1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1720 ierr = MatAssemblyEnd(M1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1721 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 1722 /* Assemble local_auxmat1 = M1*C_{CR} needed by BDDC application in KSP and in preproc */ 1723 ierr = MatMatMult(M1,C_CR,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pcbddc->local_auxmat1);CHKERRQ(ierr); 1724 1725 } 1726 1727 /* Get submatrices from subdomain matrix */ 1728 if (n_vertices) { 1729 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_vertices,vertices,PETSC_COPY_VALUES,&is_V_local);CHKERRQ(ierr); 1730 ierr = MatGetSubMatrix(pcbddc->local_mat,is_R_local,is_V_local,MAT_INITIAL_MATRIX,&A_RV);CHKERRQ(ierr); 1731 ierr = MatGetSubMatrix(pcbddc->local_mat,is_V_local,is_R_local,MAT_INITIAL_MATRIX,&A_VR);CHKERRQ(ierr); 1732 ierr = MatGetSubMatrix(pcbddc->local_mat,is_V_local,is_V_local,MAT_INITIAL_MATRIX,&A_VV);CHKERRQ(ierr); 1733 ierr = ISDestroy(&is_V_local);CHKERRQ(ierr); 1734 } 1735 1736 /* Matrix of coarse basis functions (local) */ 1737 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_phi_B);CHKERRQ(ierr); 1738 ierr = MatSetSizes(pcbddc->coarse_phi_B,n_B,pcbddc->local_primal_size,n_B,pcbddc->local_primal_size);CHKERRQ(ierr); 1739 ierr = MatSetType(pcbddc->coarse_phi_B,impMatType);CHKERRQ(ierr); 1740 ierr = MatSeqDenseSetPreallocation(pcbddc->coarse_phi_B,NULL);CHKERRQ(ierr); 1741 if (pcbddc->inexact_prec_type || dbg_flag ) { 1742 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_phi_D);CHKERRQ(ierr); 1743 ierr = MatSetSizes(pcbddc->coarse_phi_D,n_D,pcbddc->local_primal_size,n_D,pcbddc->local_primal_size);CHKERRQ(ierr); 1744 ierr = MatSetType(pcbddc->coarse_phi_D,impMatType);CHKERRQ(ierr); 1745 ierr = MatSeqDenseSetPreallocation(pcbddc->coarse_phi_D,NULL);CHKERRQ(ierr); 1746 } 1747 1748 if (dbg_flag) { 1749 ierr = PetscMalloc(2*pcbddc->local_primal_size*sizeof(*coarsefunctions_errors),&coarsefunctions_errors);CHKERRQ(ierr); 1750 ierr = PetscMalloc(2*pcbddc->local_primal_size*sizeof(*constraints_errors),&constraints_errors);CHKERRQ(ierr); 1751 } 1752 /* Subdomain contribution (Non-overlapping) to coarse matrix */ 1753 ierr = PetscMalloc ((pcbddc->local_primal_size)*(pcbddc->local_primal_size)*sizeof(PetscScalar),&coarse_submat_vals);CHKERRQ(ierr); 1754 1755 /* We are now ready to evaluate coarse basis functions and subdomain contribution to coarse problem */ 1756 for (i=0;i<n_vertices;i++){ 1757 ierr = VecSet(vec1_V,zero);CHKERRQ(ierr); 1758 ierr = VecSetValue(vec1_V,i,one,INSERT_VALUES);CHKERRQ(ierr); 1759 ierr = VecAssemblyBegin(vec1_V);CHKERRQ(ierr); 1760 ierr = VecAssemblyEnd(vec1_V);CHKERRQ(ierr); 1761 /* solution of saddle point problem */ 1762 ierr = MatMult(A_RV,vec1_V,pcbddc->vec1_R);CHKERRQ(ierr); 1763 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec1_R);CHKERRQ(ierr); 1764 ierr = VecScale(pcbddc->vec1_R,m_one);CHKERRQ(ierr); 1765 if (n_constraints) { 1766 ierr = MatMult(pcbddc->local_auxmat1,pcbddc->vec1_R,vec1_C);CHKERRQ(ierr); 1767 ierr = MatMultAdd(pcbddc->local_auxmat2,vec1_C,pcbddc->vec1_R,pcbddc->vec1_R);CHKERRQ(ierr); 1768 ierr = VecScale(vec1_C,m_one);CHKERRQ(ierr); 1769 } 1770 ierr = MatMult(A_VR,pcbddc->vec1_R,vec2_V);CHKERRQ(ierr); 1771 ierr = MatMultAdd(A_VV,vec1_V,vec2_V,vec2_V);CHKERRQ(ierr); 1772 1773 /* Set values in coarse basis function and subdomain part of coarse_mat */ 1774 /* coarse basis functions */ 1775 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1776 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1777 ierr = VecScatterEnd (pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1778 ierr = VecGetArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1779 ierr = MatSetValues(pcbddc->coarse_phi_B,n_B,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1780 ierr = VecRestoreArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1781 ierr = MatSetValue(pcbddc->coarse_phi_B,idx_V_B[i],i,one,INSERT_VALUES);CHKERRQ(ierr); 1782 if ( pcbddc->inexact_prec_type || dbg_flag ) { 1783 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1784 ierr = VecScatterEnd(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1785 ierr = VecGetArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1786 ierr = MatSetValues(pcbddc->coarse_phi_D,n_D,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1787 ierr = VecRestoreArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1788 } 1789 /* subdomain contribution to coarse matrix */ 1790 ierr = VecGetArray(vec2_V,&array);CHKERRQ(ierr); 1791 for (j=0; j<n_vertices; j++) coarse_submat_vals[i*pcbddc->local_primal_size+j] = array[j]; /* WARNING -> column major ordering */ 1792 ierr = VecRestoreArray(vec2_V,&array);CHKERRQ(ierr); 1793 if (n_constraints) { 1794 ierr = VecGetArray(vec1_C,&array);CHKERRQ(ierr); 1795 for (j=0; j<n_constraints; j++) coarse_submat_vals[i*pcbddc->local_primal_size+j+n_vertices] = array[j]; /* WARNING -> column major ordering */ 1796 ierr = VecRestoreArray(vec1_C,&array);CHKERRQ(ierr); 1797 } 1798 1799 if ( dbg_flag ) { 1800 /* assemble subdomain vector on nodes */ 1801 ierr = VecSet(pcis->vec1_N,zero);CHKERRQ(ierr); 1802 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1803 ierr = VecGetArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1804 for (j=0;j<n_R;j++) array[idx_R_local[j]] = array2[j]; 1805 array[ vertices[i] ] = one; 1806 ierr = VecRestoreArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1807 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1808 /* assemble subdomain vector of lagrange multipliers (i.e. primal nodes) */ 1809 ierr = VecSet(pcbddc->vec1_P,zero);CHKERRQ(ierr); 1810 ierr = VecGetArray(pcbddc->vec1_P,&array2);CHKERRQ(ierr); 1811 ierr = VecGetArray(vec2_V,&array);CHKERRQ(ierr); 1812 for (j=0;j<n_vertices;j++) array2[j]=array[j]; 1813 ierr = VecRestoreArray(vec2_V,&array);CHKERRQ(ierr); 1814 if (n_constraints) { 1815 ierr = VecGetArray(vec1_C,&array);CHKERRQ(ierr); 1816 for (j=0;j<n_constraints;j++) array2[j+n_vertices]=array[j]; 1817 ierr = VecRestoreArray(vec1_C,&array);CHKERRQ(ierr); 1818 } 1819 ierr = VecRestoreArray(pcbddc->vec1_P,&array2);CHKERRQ(ierr); 1820 ierr = VecScale(pcbddc->vec1_P,m_one);CHKERRQ(ierr); 1821 /* check saddle point solution */ 1822 ierr = MatMult(pcbddc->local_mat,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1823 ierr = MatMultTransposeAdd(pcbddc->ConstraintMatrix,pcbddc->vec1_P,pcis->vec2_N,pcis->vec2_N);CHKERRQ(ierr); 1824 ierr = VecNorm(pcis->vec2_N,NORM_INFINITY,&coarsefunctions_errors[i]);CHKERRQ(ierr); 1825 ierr = MatMult(pcbddc->ConstraintMatrix,pcis->vec1_N,pcbddc->vec1_P);CHKERRQ(ierr); 1826 ierr = VecGetArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1827 array[i]=array[i]+m_one; /* shift by the identity matrix */ 1828 ierr = VecRestoreArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1829 ierr = VecNorm(pcbddc->vec1_P,NORM_INFINITY,&constraints_errors[i]);CHKERRQ(ierr); 1830 } 1831 } 1832 1833 for (i=0;i<n_constraints;i++){ 1834 ierr = VecSet(vec2_C,zero);CHKERRQ(ierr); 1835 ierr = VecSetValue(vec2_C,i,m_one,INSERT_VALUES);CHKERRQ(ierr); 1836 ierr = VecAssemblyBegin(vec2_C);CHKERRQ(ierr); 1837 ierr = VecAssemblyEnd(vec2_C);CHKERRQ(ierr); 1838 /* solution of saddle point problem */ 1839 ierr = MatMult(M1,vec2_C,vec1_C);CHKERRQ(ierr); 1840 ierr = MatMult(pcbddc->local_auxmat2,vec1_C,pcbddc->vec1_R);CHKERRQ(ierr); 1841 ierr = VecScale(vec1_C,m_one);CHKERRQ(ierr); 1842 if (n_vertices) { ierr = MatMult(A_VR,pcbddc->vec1_R,vec2_V);CHKERRQ(ierr); } 1843 /* Set values in coarse basis function and subdomain part of coarse_mat */ 1844 /* coarse basis functions */ 1845 index=i+n_vertices; 1846 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1847 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1848 ierr = VecScatterEnd(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1849 ierr = VecGetArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1850 ierr = MatSetValues(pcbddc->coarse_phi_B,n_B,auxindices,1,&index,array,INSERT_VALUES);CHKERRQ(ierr); 1851 ierr = VecRestoreArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1852 if ( pcbddc->inexact_prec_type || dbg_flag ) { 1853 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1854 ierr = VecScatterEnd(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1855 ierr = VecGetArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1856 ierr = MatSetValues(pcbddc->coarse_phi_D,n_D,auxindices,1,&index,array,INSERT_VALUES);CHKERRQ(ierr); 1857 ierr = VecRestoreArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1858 } 1859 /* subdomain contribution to coarse matrix */ 1860 if (n_vertices) { 1861 ierr = VecGetArray(vec2_V,&array);CHKERRQ(ierr); 1862 for (j=0; j<n_vertices; j++) coarse_submat_vals[index*pcbddc->local_primal_size+j]=array[j]; /* WARNING -> column major ordering */ 1863 ierr = VecRestoreArray(vec2_V,&array);CHKERRQ(ierr); 1864 } 1865 ierr = VecGetArray(vec1_C,&array);CHKERRQ(ierr); 1866 for (j=0; j<n_constraints; j++) coarse_submat_vals[index*pcbddc->local_primal_size+j+n_vertices]=array[j]; /* WARNING -> column major ordering */ 1867 ierr = VecRestoreArray(vec1_C,&array);CHKERRQ(ierr); 1868 1869 if ( dbg_flag ) { 1870 /* assemble subdomain vector on nodes */ 1871 ierr = VecSet(pcis->vec1_N,zero);CHKERRQ(ierr); 1872 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1873 ierr = VecGetArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1874 for (j=0;j<n_R;j++) array[idx_R_local[j]] = array2[j]; 1875 ierr = VecRestoreArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1876 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1877 /* assemble subdomain vector of lagrange multipliers */ 1878 ierr = VecSet(pcbddc->vec1_P,zero);CHKERRQ(ierr); 1879 ierr = VecGetArray(pcbddc->vec1_P,&array2);CHKERRQ(ierr); 1880 if ( n_vertices) { 1881 ierr = VecGetArray(vec2_V,&array);CHKERRQ(ierr); 1882 for (j=0;j<n_vertices;j++) array2[j]=-array[j]; 1883 ierr = VecRestoreArray(vec2_V,&array);CHKERRQ(ierr); 1884 } 1885 ierr = VecGetArray(vec1_C,&array);CHKERRQ(ierr); 1886 for (j=0;j<n_constraints;j++) {array2[j+n_vertices]=-array[j];} 1887 ierr = VecRestoreArray(vec1_C,&array);CHKERRQ(ierr); 1888 ierr = VecRestoreArray(pcbddc->vec1_P,&array2);CHKERRQ(ierr); 1889 /* check saddle point solution */ 1890 ierr = MatMult(pcbddc->local_mat,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1891 ierr = MatMultTransposeAdd(pcbddc->ConstraintMatrix,pcbddc->vec1_P,pcis->vec2_N,pcis->vec2_N);CHKERRQ(ierr); 1892 ierr = VecNorm(pcis->vec2_N,NORM_INFINITY,&coarsefunctions_errors[index]);CHKERRQ(ierr); 1893 ierr = MatMult(pcbddc->ConstraintMatrix,pcis->vec1_N,pcbddc->vec1_P);CHKERRQ(ierr); 1894 ierr = VecGetArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1895 array[index]=array[index]+m_one; /* shift by the identity matrix */ 1896 ierr = VecRestoreArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1897 ierr = VecNorm(pcbddc->vec1_P,NORM_INFINITY,&constraints_errors[index]);CHKERRQ(ierr); 1898 } 1899 } 1900 ierr = MatAssemblyBegin(pcbddc->coarse_phi_B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1901 ierr = MatAssemblyEnd(pcbddc->coarse_phi_B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1902 if ( pcbddc->inexact_prec_type || dbg_flag ) { 1903 ierr = MatAssemblyBegin(pcbddc->coarse_phi_D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1904 ierr = MatAssemblyEnd(pcbddc->coarse_phi_D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1905 } 1906 /* compute other basis functions for non-symmetric problems */ 1907 ierr = MatIsSymmetricKnown(pc->pmat,&setsym,&issym);CHKERRQ(ierr); 1908 if ( !setsym || (setsym && !issym) ) { 1909 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_psi_B);CHKERRQ(ierr); 1910 ierr = MatSetSizes(pcbddc->coarse_psi_B,n_B,pcbddc->local_primal_size,n_B,pcbddc->local_primal_size);CHKERRQ(ierr); 1911 ierr = MatSetType(pcbddc->coarse_psi_B,impMatType);CHKERRQ(ierr); 1912 ierr = MatSeqDenseSetPreallocation(pcbddc->coarse_psi_B,NULL);CHKERRQ(ierr); 1913 if (pcbddc->inexact_prec_type || dbg_flag ) { 1914 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_psi_D);CHKERRQ(ierr); 1915 ierr = MatSetSizes(pcbddc->coarse_psi_D,n_D,pcbddc->local_primal_size,n_D,pcbddc->local_primal_size);CHKERRQ(ierr); 1916 ierr = MatSetType(pcbddc->coarse_psi_D,impMatType);CHKERRQ(ierr); 1917 ierr = MatSeqDenseSetPreallocation(pcbddc->coarse_psi_D,NULL);CHKERRQ(ierr); 1918 } 1919 for (i=0;i<pcbddc->local_primal_size;i++) { 1920 if (n_constraints) { 1921 ierr = VecSet(vec1_C,zero);CHKERRQ(ierr); 1922 ierr = VecGetArray(vec1_C,&array);CHKERRQ(ierr); 1923 for (j=0;j<n_constraints;j++) { 1924 array[j]=coarse_submat_vals[(j+n_vertices)*pcbddc->local_primal_size+i]; 1925 } 1926 ierr = VecRestoreArray(vec1_C,&array);CHKERRQ(ierr); 1927 } 1928 if (i<n_vertices) { 1929 ierr = VecSet(vec1_V,zero);CHKERRQ(ierr); 1930 ierr = VecSetValue(vec1_V,i,m_one,INSERT_VALUES);CHKERRQ(ierr); 1931 ierr = VecAssemblyBegin(vec1_V);CHKERRQ(ierr); 1932 ierr = VecAssemblyEnd(vec1_V);CHKERRQ(ierr); 1933 ierr = MatMultTranspose(A_VR,vec1_V,pcbddc->vec1_R);CHKERRQ(ierr); 1934 if (n_constraints) { 1935 ierr = MatMultTransposeAdd(C_CR,vec1_C,pcbddc->vec1_R,pcbddc->vec1_R);CHKERRQ(ierr); 1936 } 1937 } else { 1938 ierr = MatMultTranspose(C_CR,vec1_C,pcbddc->vec1_R);CHKERRQ(ierr); 1939 } 1940 ierr = KSPSolveTranspose(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec1_R);CHKERRQ(ierr); 1941 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1942 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1943 ierr = VecScatterEnd(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1944 ierr = VecGetArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1945 ierr = MatSetValues(pcbddc->coarse_psi_B,n_B,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1946 ierr = VecRestoreArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1947 if (i<n_vertices) { 1948 ierr = MatSetValue(pcbddc->coarse_psi_B,idx_V_B[i],i,one,INSERT_VALUES);CHKERRQ(ierr); 1949 } 1950 if ( pcbddc->inexact_prec_type || dbg_flag ) { 1951 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1952 ierr = VecScatterEnd(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1953 ierr = VecGetArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1954 ierr = MatSetValues(pcbddc->coarse_psi_D,n_D,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1955 ierr = VecRestoreArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1956 } 1957 1958 if ( dbg_flag ) { 1959 /* assemble subdomain vector on nodes */ 1960 ierr = VecSet(pcis->vec1_N,zero);CHKERRQ(ierr); 1961 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1962 ierr = VecGetArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1963 for (j=0;j<n_R;j++) array[idx_R_local[j]] = array2[j]; 1964 if (i<n_vertices) array[vertices[i]] = one; 1965 ierr = VecRestoreArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1966 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1967 /* assemble subdomain vector of lagrange multipliers */ 1968 ierr = VecGetArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1969 for (j=0;j<pcbddc->local_primal_size;j++) { 1970 array[j]=-coarse_submat_vals[j*pcbddc->local_primal_size+i]; 1971 } 1972 ierr = VecRestoreArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1973 /* check saddle point solution */ 1974 ierr = MatMultTranspose(pcbddc->local_mat,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1975 ierr = MatMultTransposeAdd(pcbddc->ConstraintMatrix,pcbddc->vec1_P,pcis->vec2_N,pcis->vec2_N);CHKERRQ(ierr); 1976 ierr = VecNorm(pcis->vec2_N,NORM_INFINITY,&coarsefunctions_errors[i+pcbddc->local_primal_size]);CHKERRQ(ierr); 1977 ierr = MatMult(pcbddc->ConstraintMatrix,pcis->vec1_N,pcbddc->vec1_P);CHKERRQ(ierr); 1978 ierr = VecGetArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1979 array[i]=array[i]+m_one; /* shift by the identity matrix */ 1980 ierr = VecRestoreArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1981 ierr = VecNorm(pcbddc->vec1_P,NORM_INFINITY,&constraints_errors[i+pcbddc->local_primal_size]);CHKERRQ(ierr); 1982 } 1983 } 1984 ierr = MatAssemblyBegin(pcbddc->coarse_psi_B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1985 ierr = MatAssemblyEnd(pcbddc->coarse_psi_B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1986 if ( pcbddc->inexact_prec_type || dbg_flag ) { 1987 ierr = MatAssemblyBegin(pcbddc->coarse_psi_D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1988 ierr = MatAssemblyEnd(pcbddc->coarse_psi_D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1989 } 1990 } 1991 ierr = PetscFree(idx_V_B);CHKERRQ(ierr); 1992 /* Checking coarse_sub_mat and coarse basis functios */ 1993 /* Symmetric case : It should be \Phi^{(j)^T} A^{(j)} \Phi^{(j)}=coarse_sub_mat */ 1994 /* Non-symmetric case : It should be \Psi^{(j)^T} A^{(j)} \Phi^{(j)}=coarse_sub_mat */ 1995 if (dbg_flag) { 1996 Mat coarse_sub_mat; 1997 Mat TM1,TM2,TM3,TM4; 1998 Mat coarse_phi_D,coarse_phi_B; 1999 Mat coarse_psi_D,coarse_psi_B; 2000 Mat A_II,A_BB,A_IB,A_BI; 2001 MatType checkmattype=MATSEQAIJ; 2002 PetscReal real_value; 2003 2004 ierr = MatConvert(pcis->A_II,checkmattype,MAT_INITIAL_MATRIX,&A_II);CHKERRQ(ierr); 2005 ierr = MatConvert(pcis->A_IB,checkmattype,MAT_INITIAL_MATRIX,&A_IB);CHKERRQ(ierr); 2006 ierr = MatConvert(pcis->A_BI,checkmattype,MAT_INITIAL_MATRIX,&A_BI);CHKERRQ(ierr); 2007 ierr = MatConvert(pcis->A_BB,checkmattype,MAT_INITIAL_MATRIX,&A_BB);CHKERRQ(ierr); 2008 ierr = MatConvert(pcbddc->coarse_phi_D,checkmattype,MAT_INITIAL_MATRIX,&coarse_phi_D);CHKERRQ(ierr); 2009 ierr = MatConvert(pcbddc->coarse_phi_B,checkmattype,MAT_INITIAL_MATRIX,&coarse_phi_B);CHKERRQ(ierr); 2010 if (pcbddc->coarse_psi_B) { 2011 ierr = MatConvert(pcbddc->coarse_psi_D,checkmattype,MAT_INITIAL_MATRIX,&coarse_psi_D);CHKERRQ(ierr); 2012 ierr = MatConvert(pcbddc->coarse_psi_B,checkmattype,MAT_INITIAL_MATRIX,&coarse_psi_B);CHKERRQ(ierr); 2013 } 2014 ierr = MatCreateSeqDense(PETSC_COMM_SELF,pcbddc->local_primal_size,pcbddc->local_primal_size,coarse_submat_vals,&coarse_sub_mat);CHKERRQ(ierr); 2015 2016 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 2017 ierr = PetscViewerASCIIPrintf(viewer,"Check coarse sub mat and local basis functions\n");CHKERRQ(ierr); 2018 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2019 if (pcbddc->coarse_psi_B) { 2020 ierr = MatMatMult(A_II,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 2021 ierr = MatTransposeMatMult(coarse_psi_D,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM1);CHKERRQ(ierr); 2022 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 2023 ierr = MatMatMult(A_BB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 2024 ierr = MatTransposeMatMult(coarse_psi_B,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM2);CHKERRQ(ierr); 2025 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 2026 ierr = MatMatMult(A_IB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 2027 ierr = MatTransposeMatMult(coarse_psi_D,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM3);CHKERRQ(ierr); 2028 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 2029 ierr = MatMatMult(A_BI,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 2030 ierr = MatTransposeMatMult(coarse_psi_B,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM4);CHKERRQ(ierr); 2031 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 2032 } else { 2033 ierr = MatPtAP(A_II,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&TM1);CHKERRQ(ierr); 2034 ierr = MatPtAP(A_BB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&TM2);CHKERRQ(ierr); 2035 ierr = MatMatMult(A_IB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 2036 ierr = MatTransposeMatMult(coarse_phi_D,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM3);CHKERRQ(ierr); 2037 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 2038 ierr = MatMatMult(A_BI,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 2039 ierr = MatTransposeMatMult(coarse_phi_B,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM4);CHKERRQ(ierr); 2040 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 2041 } 2042 ierr = MatAXPY(TM1,one,TM2,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 2043 ierr = MatAXPY(TM1,one,TM3,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 2044 ierr = MatAXPY(TM1,one,TM4,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 2045 ierr = MatConvert(TM1,MATSEQDENSE,MAT_REUSE_MATRIX,&TM1);CHKERRQ(ierr); 2046 ierr = MatAXPY(TM1,m_one,coarse_sub_mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2047 ierr = MatNorm(TM1,NORM_INFINITY,&real_value);CHKERRQ(ierr); 2048 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"----------------------------------\n");CHKERRQ(ierr); 2049 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d \n",PetscGlobalRank);CHKERRQ(ierr); 2050 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"matrix error = % 1.14e\n",real_value);CHKERRQ(ierr); 2051 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"coarse functions (phi) errors\n");CHKERRQ(ierr); 2052 for (i=0;i<pcbddc->local_primal_size;i++) { 2053 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local %02d-th function error = % 1.14e\n",i,coarsefunctions_errors[i]);CHKERRQ(ierr); 2054 } 2055 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"constraints (phi) errors\n");CHKERRQ(ierr); 2056 for (i=0;i<pcbddc->local_primal_size;i++) { 2057 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local %02d-th function error = % 1.14e\n",i,constraints_errors[i]);CHKERRQ(ierr); 2058 } 2059 if (pcbddc->coarse_psi_B) { 2060 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"coarse functions (psi) errors\n");CHKERRQ(ierr); 2061 for (i=pcbddc->local_primal_size;i<2*pcbddc->local_primal_size;i++) { 2062 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local %02d-th function error = % 1.14e\n",i-pcbddc->local_primal_size,coarsefunctions_errors[i]);CHKERRQ(ierr); 2063 } 2064 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"constraints (psi) errors\n");CHKERRQ(ierr); 2065 for (i=pcbddc->local_primal_size;i<2*pcbddc->local_primal_size;i++) { 2066 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local %02d-th function error = % 1.14e\n",i-pcbddc->local_primal_size,constraints_errors[i]);CHKERRQ(ierr); 2067 } 2068 } 2069 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2070 ierr = MatDestroy(&A_II);CHKERRQ(ierr); 2071 ierr = MatDestroy(&A_BB);CHKERRQ(ierr); 2072 ierr = MatDestroy(&A_IB);CHKERRQ(ierr); 2073 ierr = MatDestroy(&A_BI);CHKERRQ(ierr); 2074 ierr = MatDestroy(&TM1);CHKERRQ(ierr); 2075 ierr = MatDestroy(&TM2);CHKERRQ(ierr); 2076 ierr = MatDestroy(&TM3);CHKERRQ(ierr); 2077 ierr = MatDestroy(&TM4);CHKERRQ(ierr); 2078 ierr = MatDestroy(&coarse_phi_D);CHKERRQ(ierr); 2079 ierr = MatDestroy(&coarse_phi_B);CHKERRQ(ierr); 2080 if (pcbddc->coarse_psi_B) { 2081 ierr = MatDestroy(&coarse_psi_D);CHKERRQ(ierr); 2082 ierr = MatDestroy(&coarse_psi_B);CHKERRQ(ierr); 2083 } 2084 ierr = MatDestroy(&coarse_sub_mat);CHKERRQ(ierr); 2085 ierr = PetscFree(coarsefunctions_errors);CHKERRQ(ierr); 2086 ierr = PetscFree(constraints_errors);CHKERRQ(ierr); 2087 } 2088 /* free memory */ 2089 if (n_vertices) { 2090 ierr = PetscFree(vertices);CHKERRQ(ierr); 2091 ierr = VecDestroy(&vec1_V);CHKERRQ(ierr); 2092 ierr = VecDestroy(&vec2_V);CHKERRQ(ierr); 2093 ierr = MatDestroy(&A_RV);CHKERRQ(ierr); 2094 ierr = MatDestroy(&A_VR);CHKERRQ(ierr); 2095 ierr = MatDestroy(&A_VV);CHKERRQ(ierr); 2096 } 2097 if (n_constraints) { 2098 ierr = VecDestroy(&vec1_C);CHKERRQ(ierr); 2099 ierr = VecDestroy(&vec2_C);CHKERRQ(ierr); 2100 ierr = MatDestroy(&M1);CHKERRQ(ierr); 2101 ierr = MatDestroy(&C_CR);CHKERRQ(ierr); 2102 } 2103 ierr = PetscFree(auxindices);CHKERRQ(ierr); 2104 ierr = PetscFree(nnz);CHKERRQ(ierr); 2105 /* create coarse matrix and data structures for message passing associated actual choice of coarse problem type */ 2106 ierr = PCBDDCSetUpCoarseEnvironment(pc,coarse_submat_vals);CHKERRQ(ierr); 2107 ierr = PetscFree(coarse_submat_vals);CHKERRQ(ierr); 2108 } 2109 /* free memory */ 2110 ierr = ISDestroy(&is_R_local);CHKERRQ(ierr); 2111 2112 PetscFunctionReturn(0); 2113 } 2114 2115 /* -------------------------------------------------------------------------- */ 2116 2117 /* BDDC requires metis 5.0.1 for multilevel */ 2118 #if defined(PETSC_HAVE_METIS) 2119 #include "metis.h" 2120 #define MetisInt idx_t 2121 #define MetisScalar real_t 2122 #endif 2123 2124 #undef __FUNCT__ 2125 #define __FUNCT__ "PCBDDCSetUpCoarseEnvironment" 2126 static PetscErrorCode PCBDDCSetUpCoarseEnvironment(PC pc,PetscScalar* coarse_submat_vals) 2127 { 2128 2129 2130 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 2131 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 2132 PC_IS *pcis = (PC_IS*)pc->data; 2133 MPI_Comm prec_comm; 2134 MPI_Comm coarse_comm; 2135 2136 MatNullSpace CoarseNullSpace; 2137 2138 /* common to all choiches */ 2139 PetscScalar *temp_coarse_mat_vals; 2140 PetscScalar *ins_coarse_mat_vals; 2141 PetscInt *ins_local_primal_indices; 2142 PetscMPIInt *localsizes2,*localdispl2; 2143 PetscMPIInt size_prec_comm; 2144 PetscMPIInt rank_prec_comm; 2145 PetscMPIInt active_rank=MPI_PROC_NULL; 2146 PetscMPIInt master_proc=0; 2147 PetscInt ins_local_primal_size; 2148 /* specific to MULTILEVEL_BDDC */ 2149 PetscMPIInt *ranks_recv=0; 2150 PetscMPIInt count_recv=0; 2151 PetscMPIInt rank_coarse_proc_send_to=-1; 2152 PetscMPIInt coarse_color = MPI_UNDEFINED; 2153 ISLocalToGlobalMapping coarse_ISLG; 2154 /* some other variables */ 2155 PetscErrorCode ierr; 2156 MatType coarse_mat_type; 2157 PCType coarse_pc_type; 2158 KSPType coarse_ksp_type; 2159 PC pc_temp; 2160 PetscInt i,j,k; 2161 PetscInt max_it_coarse_ksp=1; /* don't increase this value */ 2162 /* verbose output viewer */ 2163 PetscViewer viewer=pcbddc->dbg_viewer; 2164 PetscInt dbg_flag=pcbddc->dbg_flag; 2165 2166 PetscInt offset,offset2; 2167 PetscMPIInt im_active,active_procs; 2168 PetscInt *dnz,*onz; 2169 2170 PetscBool setsym,issym=PETSC_FALSE; 2171 2172 PetscFunctionBegin; 2173 ierr = PetscObjectGetComm((PetscObject)pc,&prec_comm);CHKERRQ(ierr); 2174 ins_local_primal_indices = 0; 2175 ins_coarse_mat_vals = 0; 2176 localsizes2 = 0; 2177 localdispl2 = 0; 2178 temp_coarse_mat_vals = 0; 2179 coarse_ISLG = 0; 2180 2181 ierr = MPI_Comm_size(prec_comm,&size_prec_comm);CHKERRQ(ierr); 2182 ierr = MPI_Comm_rank(prec_comm,&rank_prec_comm);CHKERRQ(ierr); 2183 ierr = MatIsSymmetricKnown(pc->pmat,&setsym,&issym);CHKERRQ(ierr); 2184 2185 /* Assign global numbering to coarse dofs */ 2186 { 2187 PetscInt *auxlocal_primal,*aux_idx; 2188 PetscMPIInt mpi_local_primal_size; 2189 PetscScalar coarsesum,*array; 2190 2191 mpi_local_primal_size = (PetscMPIInt)pcbddc->local_primal_size; 2192 2193 /* Construct needed data structures for message passing */ 2194 j = 0; 2195 if (rank_prec_comm == 0 || pcbddc->coarse_problem_type == REPLICATED_BDDC || pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 2196 j = size_prec_comm; 2197 } 2198 ierr = PetscMalloc(j*sizeof(PetscMPIInt),&pcbddc->local_primal_sizes);CHKERRQ(ierr); 2199 ierr = PetscMalloc(j*sizeof(PetscMPIInt),&pcbddc->local_primal_displacements);CHKERRQ(ierr); 2200 /* Gather local_primal_size information for all processes */ 2201 if (pcbddc->coarse_problem_type == REPLICATED_BDDC || pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 2202 ierr = MPI_Allgather(&mpi_local_primal_size,1,MPIU_INT,&pcbddc->local_primal_sizes[0],1,MPIU_INT,prec_comm);CHKERRQ(ierr); 2203 } else { 2204 ierr = MPI_Gather(&mpi_local_primal_size,1,MPIU_INT,&pcbddc->local_primal_sizes[0],1,MPIU_INT,0,prec_comm);CHKERRQ(ierr); 2205 } 2206 pcbddc->replicated_primal_size = 0; 2207 for (i=0; i<j; i++) { 2208 pcbddc->local_primal_displacements[i] = pcbddc->replicated_primal_size ; 2209 pcbddc->replicated_primal_size += pcbddc->local_primal_sizes[i]; 2210 } 2211 2212 /* First let's count coarse dofs. 2213 This code fragment assumes that the number of local constraints per connected component 2214 is not greater than the number of nodes defined for the connected component 2215 (otherwise we will surely have linear dependence between constraints and thus a singular coarse problem) */ 2216 ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(PetscInt),&auxlocal_primal);CHKERRQ(ierr); 2217 ierr = PCBDDCGetPrimalVerticesLocalIdx(pc,&i,&aux_idx);CHKERRQ(ierr); 2218 ierr = PetscMemcpy(auxlocal_primal,aux_idx,i*sizeof(PetscInt));CHKERRQ(ierr); 2219 ierr = PetscFree(aux_idx);CHKERRQ(ierr); 2220 ierr = PCBDDCGetPrimalConstraintsLocalIdx(pc,&j,&aux_idx);CHKERRQ(ierr); 2221 ierr = PetscMemcpy(&auxlocal_primal[i],aux_idx,j*sizeof(PetscInt));CHKERRQ(ierr); 2222 ierr = PetscFree(aux_idx);CHKERRQ(ierr); 2223 /* Compute number of coarse dofs */ 2224 ierr = PCBDDCSubsetNumbering(prec_comm,matis->mapping,pcbddc->local_primal_size,auxlocal_primal,NULL,&pcbddc->coarse_size,&pcbddc->local_primal_indices);CHKERRQ(ierr); 2225 2226 if (dbg_flag) { 2227 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2228 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 2229 ierr = PetscViewerASCIIPrintf(viewer,"Check coarse indices\n");CHKERRQ(ierr); 2230 ierr = VecSet(pcis->vec1_N,0.0);CHKERRQ(ierr); 2231 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 2232 for (i=0;i<pcbddc->local_primal_size;i++) array[auxlocal_primal[i]]=1.0; 2233 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 2234 ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr); 2235 ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2236 ierr = VecScatterEnd (matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2237 ierr = VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2238 ierr = VecScatterEnd (matis->ctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2239 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 2240 for (i=0;i<pcis->n;i++) { 2241 if (array[i] == 1.0) { 2242 ierr = ISLocalToGlobalMappingApply(matis->mapping,1,&i,&j);CHKERRQ(ierr); 2243 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d: WRONG COARSE INDEX %d (local %d)\n",PetscGlobalRank,j,i);CHKERRQ(ierr); 2244 } 2245 } 2246 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2247 for (i=0;i<pcis->n;i++) { 2248 if (PetscRealPart(array[i]) > 0.0) array[i] = 1.0/PetscRealPart(array[i]); 2249 } 2250 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 2251 ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr); 2252 ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2253 ierr = VecScatterEnd (matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2254 ierr = VecSum(pcis->vec1_global,&coarsesum);CHKERRQ(ierr); 2255 ierr = PetscViewerASCIIPrintf(viewer,"Size of coarse problem SHOULD be %lf\n",coarsesum);CHKERRQ(ierr); 2256 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2257 } 2258 ierr = PetscFree(auxlocal_primal);CHKERRQ(ierr); 2259 } 2260 2261 if (dbg_flag) { 2262 ierr = PetscViewerASCIIPrintf(viewer,"Size of coarse problem is %d\n",pcbddc->coarse_size);CHKERRQ(ierr); 2263 if (dbg_flag > 1) { 2264 ierr = PetscViewerASCIIPrintf(viewer,"Distribution of local primal indices\n");CHKERRQ(ierr); 2265 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2266 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d\n",PetscGlobalRank);CHKERRQ(ierr); 2267 for (i=0;i<pcbddc->local_primal_size;i++) { 2268 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local_primal_indices[%d]=%d \n",i,pcbddc->local_primal_indices[i]); 2269 } 2270 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2271 } 2272 } 2273 2274 im_active = 0; 2275 if (pcis->n) im_active = 1; 2276 ierr = MPI_Allreduce(&im_active,&active_procs,1,MPIU_INT,MPI_SUM,prec_comm);CHKERRQ(ierr); 2277 2278 /* adapt coarse problem type */ 2279 #if defined(PETSC_HAVE_METIS) 2280 if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 2281 if (pcbddc->current_level < pcbddc->max_levels) { 2282 if ( (active_procs/pcbddc->coarsening_ratio) < 2 ) { 2283 if (dbg_flag) { 2284 ierr = PetscViewerASCIIPrintf(viewer,"Not enough active processes on level %d (active %d,ratio %d). Parallel direct solve for coarse problem\n",pcbddc->current_level,active_procs,pcbddc->coarsening_ratio);CHKERRQ(ierr); 2285 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2286 } 2287 pcbddc->coarse_problem_type = PARALLEL_BDDC; 2288 } 2289 } else { 2290 if (dbg_flag) { 2291 ierr = PetscViewerASCIIPrintf(viewer,"Max number of levels reached. Using parallel direct solve for coarse problem\n",pcbddc->max_levels,active_procs,pcbddc->coarsening_ratio);CHKERRQ(ierr); 2292 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2293 } 2294 pcbddc->coarse_problem_type = PARALLEL_BDDC; 2295 } 2296 } 2297 #else 2298 pcbddc->coarse_problem_type = PARALLEL_BDDC; 2299 #endif 2300 2301 switch(pcbddc->coarse_problem_type){ 2302 2303 case(MULTILEVEL_BDDC): /* we define a coarse mesh where subdomains are elements */ 2304 { 2305 #if defined(PETSC_HAVE_METIS) 2306 /* we need additional variables */ 2307 MetisInt n_subdomains,n_parts,objval,ncon,faces_nvtxs; 2308 MetisInt *metis_coarse_subdivision; 2309 MetisInt options[METIS_NOPTIONS]; 2310 PetscMPIInt size_coarse_comm,rank_coarse_comm; 2311 PetscMPIInt procs_jumps_coarse_comm; 2312 PetscMPIInt *coarse_subdivision; 2313 PetscMPIInt *total_count_recv; 2314 PetscMPIInt *total_ranks_recv; 2315 PetscMPIInt *displacements_recv; 2316 PetscMPIInt *my_faces_connectivity; 2317 PetscMPIInt *petsc_faces_adjncy; 2318 MetisInt *faces_adjncy; 2319 MetisInt *faces_xadj; 2320 PetscMPIInt *number_of_faces; 2321 PetscMPIInt *faces_displacements; 2322 PetscInt *array_int; 2323 PetscMPIInt my_faces=0; 2324 PetscMPIInt total_faces=0; 2325 PetscInt ranks_stretching_ratio; 2326 2327 /* define some quantities */ 2328 pcbddc->coarse_communications_type = SCATTERS_BDDC; 2329 coarse_mat_type = MATIS; 2330 coarse_pc_type = PCBDDC; 2331 coarse_ksp_type = KSPRICHARDSON; 2332 2333 /* details of coarse decomposition */ 2334 n_subdomains = active_procs; 2335 n_parts = n_subdomains/pcbddc->coarsening_ratio; 2336 ranks_stretching_ratio = size_prec_comm/active_procs; 2337 procs_jumps_coarse_comm = pcbddc->coarsening_ratio*ranks_stretching_ratio; 2338 2339 #if 0 2340 PetscMPIInt *old_ranks; 2341 PetscInt *new_ranks,*jj,*ii; 2342 MatPartitioning mat_part; 2343 IS coarse_new_decomposition,is_numbering; 2344 PetscViewer viewer_test; 2345 MPI_Comm test_coarse_comm; 2346 PetscMPIInt test_coarse_color; 2347 Mat mat_adj; 2348 /* Create new communicator for coarse problem splitting the old one */ 2349 /* procs with coarse_color = MPI_UNDEFINED will have coarse_comm = MPI_COMM_NULL (from mpi standards) 2350 key = rank_prec_comm -> keep same ordering of ranks from the old to the new communicator */ 2351 test_coarse_color = ( im_active ? 0 : MPI_UNDEFINED ); 2352 test_coarse_comm = MPI_COMM_NULL; 2353 ierr = MPI_Comm_split(prec_comm,test_coarse_color,rank_prec_comm,&test_coarse_comm);CHKERRQ(ierr); 2354 if (im_active) { 2355 ierr = PetscMalloc(n_subdomains*sizeof(PetscMPIInt),&old_ranks); 2356 ierr = PetscMalloc(size_prec_comm*sizeof(PetscInt),&new_ranks); 2357 ierr = MPI_Comm_rank(test_coarse_comm,&rank_coarse_comm);CHKERRQ(ierr); 2358 ierr = MPI_Comm_size(test_coarse_comm,&j);CHKERRQ(ierr); 2359 ierr = MPI_Allgather(&rank_prec_comm,1,MPIU_INT,old_ranks,1,MPIU_INT,test_coarse_comm);CHKERRQ(ierr); 2360 for (i=0; i<size_prec_comm; i++) new_ranks[i] = -1; 2361 for (i=0; i<n_subdomains; i++) new_ranks[old_ranks[i]] = i; 2362 ierr = PetscViewerASCIIOpen(test_coarse_comm,"test_mat_part.out",&viewer_test);CHKERRQ(ierr); 2363 k = pcis->n_neigh-1; 2364 ierr = PetscMalloc(2*sizeof(PetscInt),&ii); 2365 ii[0]=0; 2366 ii[1]=k; 2367 ierr = PetscMalloc(k*sizeof(PetscInt),&jj); 2368 for (i=0; i<k; i++) jj[i]=new_ranks[pcis->neigh[i+1]]; 2369 ierr = PetscSortInt(k,jj);CHKERRQ(ierr); 2370 ierr = MatCreateMPIAdj(test_coarse_comm,1,n_subdomains,ii,jj,NULL,&mat_adj);CHKERRQ(ierr); 2371 ierr = MatView(mat_adj,viewer_test);CHKERRQ(ierr); 2372 ierr = MatPartitioningCreate(test_coarse_comm,&mat_part);CHKERRQ(ierr); 2373 ierr = MatPartitioningSetAdjacency(mat_part,mat_adj);CHKERRQ(ierr); 2374 ierr = MatPartitioningSetFromOptions(mat_part);CHKERRQ(ierr); 2375 printf("Setting Nparts %d\n",n_parts); 2376 ierr = MatPartitioningSetNParts(mat_part,n_parts);CHKERRQ(ierr); 2377 ierr = MatPartitioningView(mat_part,viewer_test);CHKERRQ(ierr); 2378 ierr = MatPartitioningApply(mat_part,&coarse_new_decomposition);CHKERRQ(ierr); 2379 ierr = ISView(coarse_new_decomposition,viewer_test);CHKERRQ(ierr); 2380 ierr = ISPartitioningToNumbering(coarse_new_decomposition,&is_numbering);CHKERRQ(ierr); 2381 ierr = ISView(is_numbering,viewer_test);CHKERRQ(ierr); 2382 ierr = PetscViewerDestroy(&viewer_test);CHKERRQ(ierr); 2383 ierr = ISDestroy(&coarse_new_decomposition);CHKERRQ(ierr); 2384 ierr = ISDestroy(&is_numbering);CHKERRQ(ierr); 2385 ierr = MatPartitioningDestroy(&mat_part);CHKERRQ(ierr); 2386 ierr = PetscFree(old_ranks);CHKERRQ(ierr); 2387 ierr = PetscFree(new_ranks);CHKERRQ(ierr); 2388 ierr = MPI_Comm_free(&test_coarse_comm);CHKERRQ(ierr); 2389 } 2390 #endif 2391 2392 /* build CSR graph of subdomains' connectivity */ 2393 ierr = PetscMalloc (pcis->n*sizeof(PetscInt),&array_int);CHKERRQ(ierr); 2394 ierr = PetscMemzero(array_int,pcis->n*sizeof(PetscInt));CHKERRQ(ierr); 2395 for (i=1;i<pcis->n_neigh;i++){/* i=1 so I don't count myself -> faces nodes counts to 1 */ 2396 for (j=0;j<pcis->n_shared[i];j++){ 2397 array_int[ pcis->shared[i][j] ]+=1; 2398 } 2399 } 2400 for (i=1;i<pcis->n_neigh;i++){ 2401 for (j=0;j<pcis->n_shared[i];j++){ 2402 if (array_int[ pcis->shared[i][j] ] > 0 ){ 2403 my_faces++; 2404 break; 2405 } 2406 } 2407 } 2408 2409 ierr = MPI_Reduce(&my_faces,&total_faces,1,MPIU_INT,MPI_SUM,master_proc,prec_comm);CHKERRQ(ierr); 2410 ierr = PetscMalloc (my_faces*sizeof(PetscInt),&my_faces_connectivity);CHKERRQ(ierr); 2411 my_faces=0; 2412 for (i=1;i<pcis->n_neigh;i++){ 2413 for (j=0;j<pcis->n_shared[i];j++){ 2414 if (array_int[ pcis->shared[i][j] ] > 0 ){ 2415 my_faces_connectivity[my_faces]=pcis->neigh[i]; 2416 my_faces++; 2417 break; 2418 } 2419 } 2420 } 2421 if (rank_prec_comm == master_proc) { 2422 ierr = PetscMalloc (total_faces*sizeof(PetscMPIInt),&petsc_faces_adjncy);CHKERRQ(ierr); 2423 ierr = PetscMalloc (size_prec_comm*sizeof(PetscMPIInt),&number_of_faces);CHKERRQ(ierr); 2424 ierr = PetscMalloc (total_faces*sizeof(MetisInt),&faces_adjncy);CHKERRQ(ierr); 2425 ierr = PetscMalloc ((n_subdomains+1)*sizeof(MetisInt),&faces_xadj);CHKERRQ(ierr); 2426 ierr = PetscMalloc ((size_prec_comm+1)*sizeof(PetscMPIInt),&faces_displacements);CHKERRQ(ierr); 2427 } 2428 ierr = MPI_Gather(&my_faces,1,MPIU_INT,&number_of_faces[0],1,MPIU_INT,master_proc,prec_comm);CHKERRQ(ierr); 2429 if (rank_prec_comm == master_proc) { 2430 faces_xadj[0]=0; 2431 faces_displacements[0]=0; 2432 j=0; 2433 for (i=1;i<size_prec_comm+1;i++) { 2434 faces_displacements[i]=faces_displacements[i-1]+number_of_faces[i-1]; 2435 if (number_of_faces[i-1]) { 2436 j++; 2437 faces_xadj[j]=faces_xadj[j-1]+number_of_faces[i-1]; 2438 } 2439 } 2440 } 2441 ierr = MPI_Gatherv(&my_faces_connectivity[0],my_faces,MPIU_INT,&petsc_faces_adjncy[0],number_of_faces,faces_displacements,MPIU_INT,master_proc,prec_comm);CHKERRQ(ierr); 2442 ierr = PetscFree(my_faces_connectivity);CHKERRQ(ierr); 2443 ierr = PetscFree(array_int);CHKERRQ(ierr); 2444 if (rank_prec_comm == master_proc) { 2445 for (i=0;i<total_faces;i++) faces_adjncy[i]=(MetisInt)(petsc_faces_adjncy[i]/ranks_stretching_ratio); /* cast to MetisInt */ 2446 ierr = PetscFree(faces_displacements);CHKERRQ(ierr); 2447 ierr = PetscFree(number_of_faces);CHKERRQ(ierr); 2448 ierr = PetscFree(petsc_faces_adjncy);CHKERRQ(ierr); 2449 } 2450 2451 if ( rank_prec_comm == master_proc ) { 2452 2453 PetscInt heuristic_for_metis=3; 2454 2455 ncon=1; 2456 faces_nvtxs=n_subdomains; 2457 /* partition graoh induced by face connectivity */ 2458 ierr = PetscMalloc (n_subdomains*sizeof(MetisInt),&metis_coarse_subdivision);CHKERRQ(ierr); 2459 ierr = METIS_SetDefaultOptions(options); 2460 /* we need a contiguous partition of the coarse mesh */ 2461 options[METIS_OPTION_CONTIG]=1; 2462 options[METIS_OPTION_NITER]=30; 2463 if (pcbddc->coarsening_ratio > 1) { 2464 if (n_subdomains>n_parts*heuristic_for_metis) { 2465 options[METIS_OPTION_IPTYPE]=METIS_IPTYPE_EDGE; 2466 options[METIS_OPTION_OBJTYPE]=METIS_OBJTYPE_CUT; 2467 ierr = METIS_PartGraphKway(&faces_nvtxs,&ncon,faces_xadj,faces_adjncy,NULL,NULL,NULL,&n_parts,NULL,NULL,options,&objval,metis_coarse_subdivision); 2468 if (ierr != METIS_OK) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in METIS_PartGraphKway (metis error code %D) called from PCBDDCSetUpCoarseEnvironment\n",ierr); 2469 } else { 2470 ierr = METIS_PartGraphRecursive(&faces_nvtxs,&ncon,faces_xadj,faces_adjncy,NULL,NULL,NULL,&n_parts,NULL,NULL,options,&objval,metis_coarse_subdivision); 2471 if (ierr != METIS_OK) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in METIS_PartGraphRecursive (metis error code %D) called from PCBDDCSetUpCoarseEnvironment\n",ierr); 2472 } 2473 } else { 2474 for (i=0;i<n_subdomains;i++) metis_coarse_subdivision[i]=i; 2475 } 2476 ierr = PetscFree(faces_xadj);CHKERRQ(ierr); 2477 ierr = PetscFree(faces_adjncy);CHKERRQ(ierr); 2478 ierr = PetscMalloc(size_prec_comm*sizeof(PetscMPIInt),&coarse_subdivision);CHKERRQ(ierr); 2479 2480 /* copy/cast values avoiding possible type conflicts between PETSc, MPI and METIS */ 2481 for (i=0;i<size_prec_comm;i++) coarse_subdivision[i]=MPI_PROC_NULL; 2482 for (i=0;i<n_subdomains;i++) coarse_subdivision[ranks_stretching_ratio*i]=(PetscInt)(metis_coarse_subdivision[i]); 2483 ierr = PetscFree(metis_coarse_subdivision);CHKERRQ(ierr); 2484 } 2485 2486 /* Create new communicator for coarse problem splitting the old one */ 2487 if ( !(rank_prec_comm%procs_jumps_coarse_comm) && rank_prec_comm < procs_jumps_coarse_comm*n_parts ){ 2488 coarse_color=0; /* for communicator splitting */ 2489 active_rank=rank_prec_comm; /* for insertion of matrix values */ 2490 } 2491 /* procs with coarse_color = MPI_UNDEFINED will have coarse_comm = MPI_COMM_NULL (from mpi standards) 2492 key = rank_prec_comm -> keep same ordering of ranks from the old to the new communicator */ 2493 ierr = MPI_Comm_split(prec_comm,coarse_color,rank_prec_comm,&coarse_comm);CHKERRQ(ierr); 2494 2495 if ( coarse_color == 0 ) { 2496 ierr = MPI_Comm_size(coarse_comm,&size_coarse_comm);CHKERRQ(ierr); 2497 ierr = MPI_Comm_rank(coarse_comm,&rank_coarse_comm);CHKERRQ(ierr); 2498 } else { 2499 rank_coarse_comm = MPI_PROC_NULL; 2500 } 2501 2502 /* master proc take care of arranging and distributing coarse information */ 2503 if (rank_coarse_comm == master_proc) { 2504 ierr = PetscMalloc (size_coarse_comm*sizeof(PetscMPIInt),&displacements_recv);CHKERRQ(ierr); 2505 ierr = PetscMalloc (size_coarse_comm*sizeof(PetscMPIInt),&total_count_recv);CHKERRQ(ierr); 2506 ierr = PetscMalloc (n_subdomains*sizeof(PetscMPIInt),&total_ranks_recv);CHKERRQ(ierr); 2507 /* some initializations */ 2508 displacements_recv[0]=0; 2509 ierr = PetscMemzero(total_count_recv,size_coarse_comm*sizeof(PetscMPIInt));CHKERRQ(ierr); 2510 /* count from how many processes the j-th process of the coarse decomposition will receive data */ 2511 for (j=0;j<size_coarse_comm;j++) { 2512 for (i=0;i<size_prec_comm;i++) { 2513 if (coarse_subdivision[i]==j) total_count_recv[j]++; 2514 } 2515 } 2516 /* displacements needed for scatterv of total_ranks_recv */ 2517 for (i=1; i<size_coarse_comm; i++) displacements_recv[i]=displacements_recv[i-1]+total_count_recv[i-1]; 2518 2519 /* Now fill properly total_ranks_recv -> each coarse process will receive the ranks (in prec_comm communicator) of its friend (sending) processes */ 2520 ierr = PetscMemzero(total_count_recv,size_coarse_comm*sizeof(PetscMPIInt));CHKERRQ(ierr); 2521 for (j=0;j<size_coarse_comm;j++) { 2522 for (i=0;i<size_prec_comm;i++) { 2523 if (coarse_subdivision[i]==j) { 2524 total_ranks_recv[displacements_recv[j]+total_count_recv[j]]=i; 2525 total_count_recv[j]+=1; 2526 } 2527 } 2528 } 2529 /*for (j=0;j<size_coarse_comm;j++) { 2530 printf("process %d in new rank will receive from %d processes (original ranks follows)\n",j,total_count_recv[j]); 2531 for (i=0;i<total_count_recv[j];i++) { 2532 printf("%d ",total_ranks_recv[displacements_recv[j]+i]); 2533 } 2534 printf("\n"); 2535 }*/ 2536 2537 /* identify new decomposition in terms of ranks in the old communicator */ 2538 for (i=0;i<n_subdomains;i++) { 2539 coarse_subdivision[ranks_stretching_ratio*i]=coarse_subdivision[ranks_stretching_ratio*i]*procs_jumps_coarse_comm; 2540 } 2541 /*printf("coarse_subdivision in old end new ranks\n"); 2542 for (i=0;i<size_prec_comm;i++) 2543 if (coarse_subdivision[i]!=MPI_PROC_NULL) { 2544 printf("%d=(%d %d), ",i,coarse_subdivision[i],coarse_subdivision[i]/procs_jumps_coarse_comm); 2545 } else { 2546 printf("%d=(%d %d), ",i,coarse_subdivision[i],coarse_subdivision[i]); 2547 } 2548 printf("\n");*/ 2549 } 2550 2551 /* Scatter new decomposition for send details */ 2552 ierr = MPI_Scatter(&coarse_subdivision[0],1,MPIU_INT,&rank_coarse_proc_send_to,1,MPIU_INT,master_proc,prec_comm);CHKERRQ(ierr); 2553 /* Scatter receiving details to members of coarse decomposition */ 2554 if ( coarse_color == 0) { 2555 ierr = MPI_Scatter(&total_count_recv[0],1,MPIU_INT,&count_recv,1,MPIU_INT,master_proc,coarse_comm);CHKERRQ(ierr); 2556 ierr = PetscMalloc (count_recv*sizeof(PetscMPIInt),&ranks_recv);CHKERRQ(ierr); 2557 ierr = MPI_Scatterv(&total_ranks_recv[0],total_count_recv,displacements_recv,MPIU_INT,&ranks_recv[0],count_recv,MPIU_INT,master_proc,coarse_comm);CHKERRQ(ierr); 2558 } 2559 2560 /*printf("I will send my matrix data to proc %d\n",rank_coarse_proc_send_to); 2561 if (coarse_color == 0) { 2562 printf("I will receive some matrix data from %d processes (ranks follows)\n",count_recv); 2563 for (i=0;i<count_recv;i++) 2564 printf("%d ",ranks_recv[i]); 2565 printf("\n"); 2566 }*/ 2567 2568 if (rank_prec_comm == master_proc) { 2569 ierr = PetscFree(coarse_subdivision);CHKERRQ(ierr); 2570 ierr = PetscFree(total_count_recv);CHKERRQ(ierr); 2571 ierr = PetscFree(total_ranks_recv);CHKERRQ(ierr); 2572 ierr = PetscFree(displacements_recv);CHKERRQ(ierr); 2573 } 2574 #endif 2575 break; 2576 } 2577 2578 case(REPLICATED_BDDC): 2579 2580 pcbddc->coarse_communications_type = GATHERS_BDDC; 2581 coarse_mat_type = MATSEQAIJ; 2582 coarse_pc_type = PCLU; 2583 coarse_ksp_type = KSPPREONLY; 2584 coarse_comm = PETSC_COMM_SELF; 2585 active_rank = rank_prec_comm; 2586 break; 2587 2588 case(PARALLEL_BDDC): 2589 2590 pcbddc->coarse_communications_type = SCATTERS_BDDC; 2591 coarse_mat_type = MATAIJ; 2592 coarse_pc_type = PCREDUNDANT; 2593 coarse_ksp_type = KSPPREONLY; 2594 coarse_comm = prec_comm; 2595 active_rank = rank_prec_comm; 2596 break; 2597 2598 case(SEQUENTIAL_BDDC): 2599 pcbddc->coarse_communications_type = GATHERS_BDDC; 2600 coarse_mat_type = MATAIJ; 2601 coarse_pc_type = PCLU; 2602 coarse_ksp_type = KSPPREONLY; 2603 coarse_comm = PETSC_COMM_SELF; 2604 active_rank = master_proc; 2605 break; 2606 } 2607 2608 switch(pcbddc->coarse_communications_type){ 2609 2610 case(SCATTERS_BDDC): 2611 { 2612 if (pcbddc->coarse_problem_type==MULTILEVEL_BDDC) { 2613 2614 IS coarse_IS; 2615 2616 if(pcbddc->coarsening_ratio == 1) { 2617 ins_local_primal_size = pcbddc->local_primal_size; 2618 ins_local_primal_indices = pcbddc->local_primal_indices; 2619 if (coarse_color == 0) { ierr = PetscFree(ranks_recv);CHKERRQ(ierr); } 2620 /* nonzeros */ 2621 ierr = PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&dnz);CHKERRQ(ierr); 2622 ierr = PetscMemzero(dnz,ins_local_primal_size*sizeof(PetscInt));CHKERRQ(ierr); 2623 for (i=0;i<ins_local_primal_size;i++) { 2624 dnz[i] = ins_local_primal_size; 2625 } 2626 } else { 2627 PetscMPIInt send_size; 2628 PetscMPIInt *send_buffer; 2629 PetscInt *aux_ins_indices; 2630 PetscInt ii,jj; 2631 MPI_Request *requests; 2632 2633 ierr = PetscMalloc(count_recv*sizeof(PetscMPIInt),&localdispl2);CHKERRQ(ierr); 2634 /* reusing pcbddc->local_primal_displacements and pcbddc->replicated_primal_size */ 2635 ierr = PetscFree(pcbddc->local_primal_displacements);CHKERRQ(ierr); 2636 ierr = PetscMalloc((count_recv+1)*sizeof(PetscMPIInt),&pcbddc->local_primal_displacements);CHKERRQ(ierr); 2637 pcbddc->replicated_primal_size = count_recv; 2638 j = 0; 2639 for (i=0;i<count_recv;i++) { 2640 pcbddc->local_primal_displacements[i] = j; 2641 j += pcbddc->local_primal_sizes[ranks_recv[i]]; 2642 } 2643 pcbddc->local_primal_displacements[count_recv] = j; 2644 ierr = PetscMalloc(j*sizeof(PetscMPIInt),&pcbddc->replicated_local_primal_indices);CHKERRQ(ierr); 2645 /* allocate auxiliary space */ 2646 ierr = PetscMalloc(count_recv*sizeof(PetscMPIInt),&localsizes2);CHKERRQ(ierr); 2647 ierr = PetscMalloc(pcbddc->coarse_size*sizeof(PetscInt),&aux_ins_indices);CHKERRQ(ierr); 2648 ierr = PetscMemzero(aux_ins_indices,pcbddc->coarse_size*sizeof(PetscInt));CHKERRQ(ierr); 2649 /* allocate stuffs for message massing */ 2650 ierr = PetscMalloc((count_recv+1)*sizeof(MPI_Request),&requests);CHKERRQ(ierr); 2651 for (i=0;i<count_recv+1;i++) { requests[i]=MPI_REQUEST_NULL; } 2652 /* send indices to be inserted */ 2653 for (i=0;i<count_recv;i++) { 2654 send_size = pcbddc->local_primal_sizes[ranks_recv[i]]; 2655 ierr = MPI_Irecv(&pcbddc->replicated_local_primal_indices[pcbddc->local_primal_displacements[i]],send_size,MPIU_INT,ranks_recv[i],999,prec_comm,&requests[i]);CHKERRQ(ierr); 2656 } 2657 if (rank_coarse_proc_send_to != MPI_PROC_NULL ) { 2658 send_size = pcbddc->local_primal_size; 2659 ierr = PetscMalloc(send_size*sizeof(PetscMPIInt),&send_buffer);CHKERRQ(ierr); 2660 for (i=0;i<send_size;i++) { 2661 send_buffer[i]=(PetscMPIInt)pcbddc->local_primal_indices[i]; 2662 } 2663 ierr = MPI_Isend(send_buffer,send_size,MPIU_INT,rank_coarse_proc_send_to,999,prec_comm,&requests[count_recv]);CHKERRQ(ierr); 2664 } 2665 ierr = MPI_Waitall(count_recv+1,requests,MPI_STATUSES_IGNORE);CHKERRQ(ierr); 2666 if (rank_coarse_proc_send_to != MPI_PROC_NULL ) { 2667 ierr = PetscFree(send_buffer);CHKERRQ(ierr); 2668 } 2669 j = 0; 2670 for (i=0;i<count_recv;i++) { 2671 ii = pcbddc->local_primal_displacements[i+1]-pcbddc->local_primal_displacements[i]; 2672 localsizes2[i] = ii*ii; 2673 localdispl2[i] = j; 2674 j += localsizes2[i]; 2675 jj = pcbddc->local_primal_displacements[i]; 2676 /* it counts the coarse subdomains sharing the coarse node */ 2677 for (k=0;k<ii;k++) { 2678 aux_ins_indices[pcbddc->replicated_local_primal_indices[jj+k]] += 1; 2679 } 2680 } 2681 /* temp_coarse_mat_vals used to store matrix values to be received */ 2682 ierr = PetscMalloc(j*sizeof(PetscScalar),&temp_coarse_mat_vals);CHKERRQ(ierr); 2683 /* evaluate how many values I will insert in coarse mat */ 2684 ins_local_primal_size = 0; 2685 for (i=0;i<pcbddc->coarse_size;i++) { 2686 if (aux_ins_indices[i]) { 2687 ins_local_primal_size++; 2688 } 2689 } 2690 /* evaluate indices I will insert in coarse mat */ 2691 ierr = PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);CHKERRQ(ierr); 2692 j = 0; 2693 for(i=0;i<pcbddc->coarse_size;i++) { 2694 if(aux_ins_indices[i]) { 2695 ins_local_primal_indices[j] = i; 2696 j++; 2697 } 2698 } 2699 /* processes partecipating in coarse problem receive matrix data from their friends */ 2700 for (i=0;i<count_recv;i++) { 2701 ierr = MPI_Irecv(&temp_coarse_mat_vals[localdispl2[i]],localsizes2[i],MPIU_SCALAR,ranks_recv[i],666,prec_comm,&requests[i]);CHKERRQ(ierr); 2702 } 2703 if (rank_coarse_proc_send_to != MPI_PROC_NULL ) { 2704 send_size = pcbddc->local_primal_size*pcbddc->local_primal_size; 2705 ierr = MPI_Isend(&coarse_submat_vals[0],send_size,MPIU_SCALAR,rank_coarse_proc_send_to,666,prec_comm,&requests[count_recv]);CHKERRQ(ierr); 2706 } 2707 ierr = MPI_Waitall(count_recv+1,requests,MPI_STATUSES_IGNORE);CHKERRQ(ierr); 2708 /* nonzeros */ 2709 ierr = PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&dnz);CHKERRQ(ierr); 2710 ierr = PetscMemzero(dnz,ins_local_primal_size*sizeof(PetscInt));CHKERRQ(ierr); 2711 /* use aux_ins_indices to realize a global to local mapping */ 2712 j=0; 2713 for(i=0;i<pcbddc->coarse_size;i++){ 2714 if(aux_ins_indices[i]==0){ 2715 aux_ins_indices[i]=-1; 2716 } else { 2717 aux_ins_indices[i]=j; 2718 j++; 2719 } 2720 } 2721 for (i=0;i<count_recv;i++) { 2722 j = pcbddc->local_primal_sizes[ranks_recv[i]]; 2723 for (k=0;k<j;k++) { 2724 dnz[aux_ins_indices[pcbddc->replicated_local_primal_indices[pcbddc->local_primal_displacements[i]+k]]] += j; 2725 } 2726 } 2727 /* check */ 2728 for (i=0;i<ins_local_primal_size;i++) { 2729 if (dnz[i] > ins_local_primal_size) { 2730 dnz[i] = ins_local_primal_size; 2731 } 2732 } 2733 ierr = PetscFree(requests);CHKERRQ(ierr); 2734 ierr = PetscFree(aux_ins_indices);CHKERRQ(ierr); 2735 if (coarse_color == 0) { ierr = PetscFree(ranks_recv);CHKERRQ(ierr); } 2736 } 2737 /* create local to global mapping needed by coarse MATIS */ 2738 if (coarse_comm != MPI_COMM_NULL ) {ierr = MPI_Comm_free(&coarse_comm);CHKERRQ(ierr);} 2739 coarse_comm = prec_comm; 2740 active_rank = rank_prec_comm; 2741 ierr = ISCreateGeneral(coarse_comm,ins_local_primal_size,ins_local_primal_indices,PETSC_COPY_VALUES,&coarse_IS);CHKERRQ(ierr); 2742 ierr = ISLocalToGlobalMappingCreateIS(coarse_IS,&coarse_ISLG);CHKERRQ(ierr); 2743 ierr = ISDestroy(&coarse_IS);CHKERRQ(ierr); 2744 } else if (pcbddc->coarse_problem_type==PARALLEL_BDDC) { 2745 /* arrays for values insertion */ 2746 ins_local_primal_size = pcbddc->local_primal_size; 2747 ierr = PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);CHKERRQ(ierr); 2748 ierr = PetscMalloc(ins_local_primal_size*ins_local_primal_size*sizeof(PetscScalar),&ins_coarse_mat_vals);CHKERRQ(ierr); 2749 for (j=0;j<ins_local_primal_size;j++){ 2750 ins_local_primal_indices[j]=pcbddc->local_primal_indices[j]; 2751 for (i=0;i<ins_local_primal_size;i++) { 2752 ins_coarse_mat_vals[j*ins_local_primal_size+i]=coarse_submat_vals[j*ins_local_primal_size+i]; 2753 } 2754 } 2755 } 2756 break; 2757 2758 } 2759 2760 case(GATHERS_BDDC): 2761 { 2762 2763 PetscMPIInt mysize,mysize2; 2764 PetscMPIInt *send_buffer; 2765 2766 if (rank_prec_comm==active_rank) { 2767 ierr = PetscMalloc ( pcbddc->replicated_primal_size*sizeof(PetscMPIInt),&pcbddc->replicated_local_primal_indices);CHKERRQ(ierr); 2768 ierr = PetscMalloc ( pcbddc->replicated_primal_size*sizeof(PetscScalar),&pcbddc->replicated_local_primal_values);CHKERRQ(ierr); 2769 ierr = PetscMalloc ( size_prec_comm*sizeof(PetscMPIInt),&localsizes2);CHKERRQ(ierr); 2770 ierr = PetscMalloc ( size_prec_comm*sizeof(PetscMPIInt),&localdispl2);CHKERRQ(ierr); 2771 /* arrays for values insertion */ 2772 for (i=0;i<size_prec_comm;i++) localsizes2[i]=pcbddc->local_primal_sizes[i]*pcbddc->local_primal_sizes[i]; 2773 localdispl2[0]=0; 2774 for (i=1;i<size_prec_comm;i++) localdispl2[i]=localsizes2[i-1]+localdispl2[i-1]; 2775 j=0; 2776 for (i=0;i<size_prec_comm;i++) j+=localsizes2[i]; 2777 ierr = PetscMalloc ( j*sizeof(PetscScalar),&temp_coarse_mat_vals);CHKERRQ(ierr); 2778 } 2779 2780 mysize=pcbddc->local_primal_size; 2781 mysize2=pcbddc->local_primal_size*pcbddc->local_primal_size; 2782 ierr = PetscMalloc(mysize*sizeof(PetscMPIInt),&send_buffer);CHKERRQ(ierr); 2783 for (i=0; i<mysize; i++) send_buffer[i]=(PetscMPIInt)pcbddc->local_primal_indices[i]; 2784 2785 if (pcbddc->coarse_problem_type == SEQUENTIAL_BDDC){ 2786 ierr = MPI_Gatherv(send_buffer,mysize,MPIU_INT,&pcbddc->replicated_local_primal_indices[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_INT,master_proc,prec_comm);CHKERRQ(ierr); 2787 ierr = MPI_Gatherv(&coarse_submat_vals[0],mysize2,MPIU_SCALAR,&temp_coarse_mat_vals[0],localsizes2,localdispl2,MPIU_SCALAR,master_proc,prec_comm);CHKERRQ(ierr); 2788 } else { 2789 ierr = MPI_Allgatherv(send_buffer,mysize,MPIU_INT,&pcbddc->replicated_local_primal_indices[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_INT,prec_comm);CHKERRQ(ierr); 2790 ierr = MPI_Allgatherv(&coarse_submat_vals[0],mysize2,MPIU_SCALAR,&temp_coarse_mat_vals[0],localsizes2,localdispl2,MPIU_SCALAR,prec_comm);CHKERRQ(ierr); 2791 } 2792 ierr = PetscFree(send_buffer);CHKERRQ(ierr); 2793 break; 2794 }/* switch on coarse problem and communications associated with finished */ 2795 } 2796 2797 /* Now create and fill up coarse matrix */ 2798 if ( rank_prec_comm == active_rank ) { 2799 2800 Mat matis_coarse_local_mat; 2801 2802 if (pcbddc->coarse_problem_type != MULTILEVEL_BDDC) { 2803 ierr = MatCreate(coarse_comm,&pcbddc->coarse_mat);CHKERRQ(ierr); 2804 ierr = MatSetSizes(pcbddc->coarse_mat,PETSC_DECIDE,PETSC_DECIDE,pcbddc->coarse_size,pcbddc->coarse_size);CHKERRQ(ierr); 2805 ierr = MatSetType(pcbddc->coarse_mat,coarse_mat_type);CHKERRQ(ierr); 2806 ierr = MatSetOptionsPrefix(pcbddc->coarse_mat,"coarse_");CHKERRQ(ierr); 2807 ierr = MatSetFromOptions(pcbddc->coarse_mat);CHKERRQ(ierr); 2808 ierr = MatSetUp(pcbddc->coarse_mat);CHKERRQ(ierr); 2809 ierr = MatSetOption(pcbddc->coarse_mat,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); /* local values stored in column major */ 2810 ierr = MatSetOption(pcbddc->coarse_mat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 2811 } else { 2812 ierr = MatCreateIS(coarse_comm,1,PETSC_DECIDE,PETSC_DECIDE,pcbddc->coarse_size,pcbddc->coarse_size,coarse_ISLG,&pcbddc->coarse_mat);CHKERRQ(ierr); 2813 ierr = MatSetUp(pcbddc->coarse_mat);CHKERRQ(ierr); 2814 ierr = MatISGetLocalMat(pcbddc->coarse_mat,&matis_coarse_local_mat);CHKERRQ(ierr); 2815 ierr = MatSetOptionsPrefix(pcbddc->coarse_mat,"coarse_");CHKERRQ(ierr); 2816 ierr = MatSetFromOptions(pcbddc->coarse_mat);CHKERRQ(ierr); 2817 ierr = MatSetUp(matis_coarse_local_mat);CHKERRQ(ierr); 2818 ierr = MatSetOption(matis_coarse_local_mat,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); /* local values stored in column major */ 2819 ierr = MatSetOption(matis_coarse_local_mat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 2820 } 2821 /* preallocation */ 2822 if (pcbddc->coarse_problem_type != MULTILEVEL_BDDC) { 2823 2824 PetscInt lrows,lcols,bs; 2825 2826 ierr = MatGetLocalSize(pcbddc->coarse_mat,&lrows,&lcols);CHKERRQ(ierr); 2827 ierr = MatPreallocateInitialize(coarse_comm,lrows,lcols,dnz,onz);CHKERRQ(ierr); 2828 ierr = MatGetBlockSize(pcbddc->coarse_mat,&bs);CHKERRQ(ierr); 2829 2830 if (pcbddc->coarse_problem_type == PARALLEL_BDDC) { 2831 2832 Vec vec_dnz,vec_onz; 2833 PetscScalar *my_dnz,*my_onz,*array; 2834 PetscInt *mat_ranges,*row_ownership; 2835 PetscInt coarse_index_row,coarse_index_col,owner; 2836 2837 ierr = VecCreate(prec_comm,&vec_dnz);CHKERRQ(ierr); 2838 ierr = VecSetBlockSize(vec_dnz,bs);CHKERRQ(ierr); 2839 ierr = VecSetSizes(vec_dnz,PETSC_DECIDE,pcbddc->coarse_size);CHKERRQ(ierr); 2840 ierr = VecSetType(vec_dnz,VECMPI);CHKERRQ(ierr); 2841 ierr = VecDuplicate(vec_dnz,&vec_onz);CHKERRQ(ierr); 2842 2843 ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(PetscScalar),&my_dnz);CHKERRQ(ierr); 2844 ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(PetscScalar),&my_onz);CHKERRQ(ierr); 2845 ierr = PetscMemzero(my_dnz,pcbddc->local_primal_size*sizeof(PetscScalar));CHKERRQ(ierr); 2846 ierr = PetscMemzero(my_onz,pcbddc->local_primal_size*sizeof(PetscScalar));CHKERRQ(ierr); 2847 2848 ierr = PetscMalloc(pcbddc->coarse_size*sizeof(PetscInt),&row_ownership);CHKERRQ(ierr); 2849 ierr = MatGetOwnershipRanges(pcbddc->coarse_mat,(const PetscInt**)&mat_ranges);CHKERRQ(ierr); 2850 for (i=0;i<size_prec_comm;i++) { 2851 for (j=mat_ranges[i];j<mat_ranges[i+1];j++) { 2852 row_ownership[j]=i; 2853 } 2854 } 2855 2856 for (i=0;i<pcbddc->local_primal_size;i++) { 2857 coarse_index_row = pcbddc->local_primal_indices[i]; 2858 owner = row_ownership[coarse_index_row]; 2859 for (j=i;j<pcbddc->local_primal_size;j++) { 2860 owner = row_ownership[coarse_index_row]; 2861 coarse_index_col = pcbddc->local_primal_indices[j]; 2862 if (coarse_index_col > mat_ranges[owner]-1 && coarse_index_col < mat_ranges[owner+1] ) { 2863 my_dnz[i] += 1.0; 2864 } else { 2865 my_onz[i] += 1.0; 2866 } 2867 if (i != j) { 2868 owner = row_ownership[coarse_index_col]; 2869 if (coarse_index_row > mat_ranges[owner]-1 && coarse_index_row < mat_ranges[owner+1] ) { 2870 my_dnz[j] += 1.0; 2871 } else { 2872 my_onz[j] += 1.0; 2873 } 2874 } 2875 } 2876 } 2877 ierr = VecSet(vec_dnz,0.0);CHKERRQ(ierr); 2878 ierr = VecSet(vec_onz,0.0);CHKERRQ(ierr); 2879 if (pcbddc->local_primal_size) { 2880 ierr = VecSetValues(vec_dnz,pcbddc->local_primal_size,pcbddc->local_primal_indices,my_dnz,ADD_VALUES);CHKERRQ(ierr); 2881 ierr = VecSetValues(vec_onz,pcbddc->local_primal_size,pcbddc->local_primal_indices,my_onz,ADD_VALUES);CHKERRQ(ierr); 2882 } 2883 ierr = VecAssemblyBegin(vec_dnz);CHKERRQ(ierr); 2884 ierr = VecAssemblyBegin(vec_onz);CHKERRQ(ierr); 2885 ierr = VecAssemblyEnd(vec_dnz);CHKERRQ(ierr); 2886 ierr = VecAssemblyEnd(vec_onz);CHKERRQ(ierr); 2887 j = mat_ranges[rank_prec_comm+1]-mat_ranges[rank_prec_comm]; 2888 ierr = VecGetArray(vec_dnz,&array);CHKERRQ(ierr); 2889 for (i=0; i<j; i++) dnz[i] = (PetscInt)PetscRealPart(array[i]); 2890 2891 ierr = VecRestoreArray(vec_dnz,&array);CHKERRQ(ierr); 2892 ierr = VecGetArray(vec_onz,&array);CHKERRQ(ierr); 2893 for (i=0;i<j;i++) onz[i] = (PetscInt)PetscRealPart(array[i]); 2894 2895 ierr = VecRestoreArray(vec_onz,&array);CHKERRQ(ierr); 2896 ierr = PetscFree(my_dnz);CHKERRQ(ierr); 2897 ierr = PetscFree(my_onz);CHKERRQ(ierr); 2898 ierr = PetscFree(row_ownership);CHKERRQ(ierr); 2899 ierr = VecDestroy(&vec_dnz);CHKERRQ(ierr); 2900 ierr = VecDestroy(&vec_onz);CHKERRQ(ierr); 2901 } else { 2902 for (k=0;k<size_prec_comm;k++){ 2903 offset=pcbddc->local_primal_displacements[k]; 2904 offset2=localdispl2[k]; 2905 ins_local_primal_size = pcbddc->local_primal_sizes[k]; 2906 ierr = PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);CHKERRQ(ierr); 2907 for (j=0;j<ins_local_primal_size;j++){ 2908 ins_local_primal_indices[j]=(PetscInt)pcbddc->replicated_local_primal_indices[offset+j]; 2909 } 2910 for (j=0;j<ins_local_primal_size;j++) { 2911 ierr = MatPreallocateSet(ins_local_primal_indices[j],ins_local_primal_size,ins_local_primal_indices,dnz,onz);CHKERRQ(ierr); 2912 } 2913 ierr = PetscFree(ins_local_primal_indices);CHKERRQ(ierr); 2914 } 2915 } 2916 2917 /* check */ 2918 for (i=0;i<lrows;i++) { 2919 if (dnz[i]>lcols) dnz[i]=lcols; 2920 if (onz[i]>pcbddc->coarse_size-lcols) onz[i]=pcbddc->coarse_size-lcols; 2921 } 2922 ierr = MatSeqAIJSetPreallocation(pcbddc->coarse_mat,0,dnz);CHKERRQ(ierr); 2923 ierr = MatMPIAIJSetPreallocation(pcbddc->coarse_mat,0,dnz,0,onz);CHKERRQ(ierr); 2924 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 2925 } else { 2926 ierr = MatSeqAIJSetPreallocation(matis_coarse_local_mat,0,dnz);CHKERRQ(ierr); 2927 ierr = PetscFree(dnz);CHKERRQ(ierr); 2928 } 2929 /* insert values */ 2930 if (pcbddc->coarse_problem_type == PARALLEL_BDDC) { 2931 ierr = MatSetValues(pcbddc->coarse_mat,ins_local_primal_size,ins_local_primal_indices,ins_local_primal_size,ins_local_primal_indices,ins_coarse_mat_vals,ADD_VALUES);CHKERRQ(ierr); 2932 } else if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 2933 if (pcbddc->coarsening_ratio == 1) { 2934 ins_coarse_mat_vals = coarse_submat_vals; 2935 ierr = MatSetValues(pcbddc->coarse_mat,ins_local_primal_size,ins_local_primal_indices,ins_local_primal_size,ins_local_primal_indices,ins_coarse_mat_vals,INSERT_VALUES);CHKERRQ(ierr); 2936 } else { 2937 ierr = PetscFree(ins_local_primal_indices);CHKERRQ(ierr); 2938 for (k=0;k<pcbddc->replicated_primal_size;k++) { 2939 offset = pcbddc->local_primal_displacements[k]; 2940 offset2 = localdispl2[k]; 2941 ins_local_primal_size = pcbddc->local_primal_displacements[k+1]-pcbddc->local_primal_displacements[k]; 2942 ierr = PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);CHKERRQ(ierr); 2943 for (j=0;j<ins_local_primal_size;j++){ 2944 ins_local_primal_indices[j]=(PetscInt)pcbddc->replicated_local_primal_indices[offset+j]; 2945 } 2946 ins_coarse_mat_vals = &temp_coarse_mat_vals[offset2]; 2947 ierr = MatSetValues(pcbddc->coarse_mat,ins_local_primal_size,ins_local_primal_indices,ins_local_primal_size,ins_local_primal_indices,ins_coarse_mat_vals,ADD_VALUES);CHKERRQ(ierr); 2948 ierr = PetscFree(ins_local_primal_indices);CHKERRQ(ierr); 2949 } 2950 } 2951 ins_local_primal_indices = 0; 2952 ins_coarse_mat_vals = 0; 2953 } else { 2954 for (k=0;k<size_prec_comm;k++){ 2955 offset=pcbddc->local_primal_displacements[k]; 2956 offset2=localdispl2[k]; 2957 ins_local_primal_size = pcbddc->local_primal_sizes[k]; 2958 ierr = PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);CHKERRQ(ierr); 2959 for (j=0;j<ins_local_primal_size;j++){ 2960 ins_local_primal_indices[j]=(PetscInt)pcbddc->replicated_local_primal_indices[offset+j]; 2961 } 2962 ins_coarse_mat_vals = &temp_coarse_mat_vals[offset2]; 2963 ierr = MatSetValues(pcbddc->coarse_mat,ins_local_primal_size,ins_local_primal_indices,ins_local_primal_size,ins_local_primal_indices,ins_coarse_mat_vals,ADD_VALUES);CHKERRQ(ierr); 2964 ierr = PetscFree(ins_local_primal_indices);CHKERRQ(ierr); 2965 } 2966 ins_local_primal_indices = 0; 2967 ins_coarse_mat_vals = 0; 2968 } 2969 ierr = MatAssemblyBegin(pcbddc->coarse_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2970 ierr = MatAssemblyEnd(pcbddc->coarse_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2971 /* symmetry of coarse matrix */ 2972 if (issym) { 2973 ierr = MatSetOption(pcbddc->coarse_mat,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 2974 } 2975 ierr = MatGetVecs(pcbddc->coarse_mat,&pcbddc->coarse_vec,&pcbddc->coarse_rhs);CHKERRQ(ierr); 2976 } 2977 2978 /* create loc to glob scatters if needed */ 2979 if (pcbddc->coarse_communications_type == SCATTERS_BDDC) { 2980 IS local_IS,global_IS; 2981 ierr = ISCreateStride(PETSC_COMM_SELF,pcbddc->local_primal_size,0,1,&local_IS);CHKERRQ(ierr); 2982 ierr = ISCreateGeneral(PETSC_COMM_SELF,pcbddc->local_primal_size,pcbddc->local_primal_indices,PETSC_COPY_VALUES,&global_IS);CHKERRQ(ierr); 2983 ierr = VecScatterCreate(pcbddc->vec1_P,local_IS,pcbddc->coarse_vec,global_IS,&pcbddc->coarse_loc_to_glob);CHKERRQ(ierr); 2984 ierr = ISDestroy(&local_IS);CHKERRQ(ierr); 2985 ierr = ISDestroy(&global_IS);CHKERRQ(ierr); 2986 } 2987 2988 /* free memory no longer needed */ 2989 if (coarse_ISLG) { ierr = ISLocalToGlobalMappingDestroy(&coarse_ISLG);CHKERRQ(ierr); } 2990 if (ins_local_primal_indices) { ierr = PetscFree(ins_local_primal_indices);CHKERRQ(ierr); } 2991 if (ins_coarse_mat_vals) { ierr = PetscFree(ins_coarse_mat_vals);CHKERRQ(ierr); } 2992 if (localsizes2) { ierr = PetscFree(localsizes2);CHKERRQ(ierr); } 2993 if (localdispl2) { ierr = PetscFree(localdispl2);CHKERRQ(ierr); } 2994 if (temp_coarse_mat_vals) { ierr = PetscFree(temp_coarse_mat_vals);CHKERRQ(ierr); } 2995 2996 /* Compute coarse null space */ 2997 CoarseNullSpace = 0; 2998 if (pcbddc->NullSpace) { 2999 ierr = PCBDDCNullSpaceAssembleCoarse(pc,&CoarseNullSpace);CHKERRQ(ierr); 3000 } 3001 3002 /* KSP for coarse problem */ 3003 if (rank_prec_comm == active_rank) { 3004 PetscBool isbddc=PETSC_FALSE; 3005 3006 ierr = KSPCreate(coarse_comm,&pcbddc->coarse_ksp);CHKERRQ(ierr); 3007 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->coarse_ksp,(PetscObject)pc,1);CHKERRQ(ierr); 3008 ierr = KSPSetOperators(pcbddc->coarse_ksp,pcbddc->coarse_mat,pcbddc->coarse_mat,SAME_PRECONDITIONER);CHKERRQ(ierr); 3009 ierr = KSPSetTolerances(pcbddc->coarse_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,max_it_coarse_ksp);CHKERRQ(ierr); 3010 ierr = KSPSetType(pcbddc->coarse_ksp,coarse_ksp_type);CHKERRQ(ierr); 3011 ierr = KSPGetPC(pcbddc->coarse_ksp,&pc_temp);CHKERRQ(ierr); 3012 ierr = PCSetType(pc_temp,coarse_pc_type);CHKERRQ(ierr); 3013 /* Allow user's customization */ 3014 ierr = KSPSetOptionsPrefix(pcbddc->coarse_ksp,"coarse_");CHKERRQ(ierr); 3015 /* Set Up PC for coarse problem BDDC */ 3016 if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 3017 i = pcbddc->current_level+1; 3018 ierr = PCBDDCSetLevel(pc_temp,i);CHKERRQ(ierr); 3019 ierr = PCBDDCSetCoarseningRatio(pc_temp,pcbddc->coarsening_ratio);CHKERRQ(ierr); 3020 ierr = PCBDDCSetMaxLevels(pc_temp,pcbddc->max_levels);CHKERRQ(ierr); 3021 ierr = PCBDDCSetCoarseProblemType(pc_temp,MULTILEVEL_BDDC);CHKERRQ(ierr); 3022 if (CoarseNullSpace) { 3023 ierr = PCBDDCSetNullSpace(pc_temp,CoarseNullSpace);CHKERRQ(ierr); 3024 } 3025 if (dbg_flag) { 3026 ierr = PetscViewerASCIIPrintf(viewer,"----------------Level %d: Setting up level %d---------------\n",pcbddc->current_level,i);CHKERRQ(ierr); 3027 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 3028 } 3029 } else { 3030 if (CoarseNullSpace) { 3031 ierr = KSPSetNullSpace(pcbddc->coarse_ksp,CoarseNullSpace);CHKERRQ(ierr); 3032 } 3033 } 3034 ierr = KSPSetFromOptions(pcbddc->coarse_ksp);CHKERRQ(ierr); 3035 ierr = KSPSetUp(pcbddc->coarse_ksp);CHKERRQ(ierr); 3036 3037 ierr = KSPGetTolerances(pcbddc->coarse_ksp,NULL,NULL,NULL,&j);CHKERRQ(ierr); 3038 ierr = KSPGetPC(pcbddc->coarse_ksp,&pc_temp);CHKERRQ(ierr); 3039 ierr = PetscObjectTypeCompare((PetscObject)pc_temp,PCBDDC,&isbddc);CHKERRQ(ierr); 3040 if (j == 1) { 3041 ierr = KSPSetNormType(pcbddc->coarse_ksp,KSP_NORM_NONE);CHKERRQ(ierr); 3042 if (isbddc) { 3043 ierr = PCBDDCSetUseExactDirichlet(pc_temp,PETSC_FALSE);CHKERRQ(ierr); 3044 } 3045 } 3046 } 3047 /* Check coarse problem if requested */ 3048 if ( dbg_flag && rank_prec_comm == active_rank ) { 3049 KSP check_ksp; 3050 PC check_pc; 3051 Vec check_vec; 3052 PetscReal abs_infty_error,infty_error,lambda_min,lambda_max; 3053 KSPType check_ksp_type; 3054 3055 /* Create ksp object suitable for extreme eigenvalues' estimation */ 3056 ierr = KSPCreate(coarse_comm,&check_ksp);CHKERRQ(ierr); 3057 ierr = KSPSetOperators(check_ksp,pcbddc->coarse_mat,pcbddc->coarse_mat,SAME_PRECONDITIONER);CHKERRQ(ierr); 3058 ierr = KSPSetTolerances(check_ksp,1.e-12,1.e-12,PETSC_DEFAULT,pcbddc->coarse_size);CHKERRQ(ierr); 3059 if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 3060 if (issym) check_ksp_type = KSPCG; 3061 else check_ksp_type = KSPGMRES; 3062 ierr = KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);CHKERRQ(ierr); 3063 } else { 3064 check_ksp_type = KSPPREONLY; 3065 } 3066 ierr = KSPSetType(check_ksp,check_ksp_type);CHKERRQ(ierr); 3067 ierr = KSPGetPC(pcbddc->coarse_ksp,&check_pc);CHKERRQ(ierr); 3068 ierr = KSPSetPC(check_ksp,check_pc);CHKERRQ(ierr); 3069 ierr = KSPSetUp(check_ksp);CHKERRQ(ierr); 3070 /* create random vec */ 3071 ierr = VecDuplicate(pcbddc->coarse_vec,&check_vec);CHKERRQ(ierr); 3072 ierr = VecSetRandom(check_vec,NULL);CHKERRQ(ierr); 3073 if (CoarseNullSpace) { 3074 ierr = MatNullSpaceRemove(CoarseNullSpace,check_vec);CHKERRQ(ierr); 3075 } 3076 ierr = MatMult(pcbddc->coarse_mat,check_vec,pcbddc->coarse_rhs);CHKERRQ(ierr); 3077 /* solve coarse problem */ 3078 ierr = KSPSolve(check_ksp,pcbddc->coarse_rhs,pcbddc->coarse_vec);CHKERRQ(ierr); 3079 if (CoarseNullSpace) { 3080 ierr = MatNullSpaceRemove(CoarseNullSpace,pcbddc->coarse_vec);CHKERRQ(ierr); 3081 } 3082 /* check coarse problem residual error */ 3083 ierr = VecAXPY(check_vec,-1.0,pcbddc->coarse_vec);CHKERRQ(ierr); 3084 ierr = VecNorm(check_vec,NORM_INFINITY,&infty_error);CHKERRQ(ierr); 3085 ierr = MatMult(pcbddc->coarse_mat,check_vec,pcbddc->coarse_rhs);CHKERRQ(ierr); 3086 ierr = VecNorm(pcbddc->coarse_rhs,NORM_INFINITY,&abs_infty_error);CHKERRQ(ierr); 3087 ierr = VecDestroy(&check_vec);CHKERRQ(ierr); 3088 /* get eigenvalue estimation if inexact */ 3089 if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 3090 ierr = KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);CHKERRQ(ierr); 3091 ierr = KSPGetIterationNumber(check_ksp,&k);CHKERRQ(ierr); 3092 ierr = PetscViewerASCIIPrintf(viewer,"Coarse problem eigenvalues estimated with %d iterations of %s.\n",k,check_ksp_type);CHKERRQ(ierr); 3093 ierr = PetscViewerASCIIPrintf(viewer,"Coarse problem eigenvalues: % 1.14e %1.14e\n",lambda_min,lambda_max);CHKERRQ(ierr); 3094 } 3095 ierr = PetscViewerASCIIPrintf(viewer,"Coarse problem exact infty_error : %1.14e\n",infty_error);CHKERRQ(ierr); 3096 ierr = PetscViewerASCIIPrintf(viewer,"Coarse problem residual infty_error: %1.14e\n",abs_infty_error);CHKERRQ(ierr); 3097 ierr = KSPDestroy(&check_ksp);CHKERRQ(ierr); 3098 } 3099 if (dbg_flag) { 3100 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 3101 } 3102 ierr = MatNullSpaceDestroy(&CoarseNullSpace);CHKERRQ(ierr); 3103 3104 PetscFunctionReturn(0); 3105 } 3106 3107