1 /* TODOLIST 2 3 ConstraintsSetup 4 - tolerances for constraints as an option (take care of single precision!) 5 - Can MAT_IGNORE_ZERO_ENTRIES be used for Constraints Matrix? 6 7 Solvers 8 - Add support for reuse fill and cholecky factor for coarse solver (similar to local solvers) 9 - Propagate ksp prefixes for solvers to mat objects? 10 - Propagate nearnullspace info among levels 11 12 User interface 13 - Negative indices in dirichlet and Neumann is should be skipped (now they cause out-of-bounds access) 14 - DofSplitting and DM attached to pc? 15 16 Debugging output 17 - Better management of verbosity levels of debugging output 18 19 Build 20 - make runexe59 21 22 Extra 23 - Is it possible to work with PCBDDCGraph on boundary indices only (less memory consumed)? 24 - add support for computing h,H and related using coordinates? 25 - Change of basis approach does not work with my nonlinear mechanics example. why? (seems not an issue with l2gmap) 26 - Better management in PCIS code 27 - BDDC with MG framework? 28 29 FETIDP 30 - Move FETIDP code to its own classes 31 32 MATIS related operations contained in BDDC code 33 - Provide general case for subassembling 34 - Preallocation routines in MatISGetMPIAXAIJ 35 36 */ 37 38 /* ---------------------------------------------------------------------------------------------------------------------------------------------- 39 Implementation of BDDC preconditioner based on: 40 C. Dohrmann "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 41 ---------------------------------------------------------------------------------------------------------------------------------------------- */ 42 43 #include <../src/ksp/pc/impls/bddc/bddc.h> /*I "petscpc.h" I*/ /* includes for fortran wrappers */ 44 #include <../src/ksp/pc/impls/bddc/bddcprivate.h> 45 #include <petscblaslapack.h> 46 47 /* -------------------------------------------------------------------------- */ 48 #undef __FUNCT__ 49 #define __FUNCT__ "PCSetFromOptions_BDDC" 50 PetscErrorCode PCSetFromOptions_BDDC(PC pc) 51 { 52 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 53 PetscErrorCode ierr; 54 55 PetscFunctionBegin; 56 ierr = PetscOptionsHead("BDDC options");CHKERRQ(ierr); 57 /* Verbose debugging */ 58 ierr = PetscOptionsInt("-pc_bddc_check_level","Verbose output for PCBDDC (intended for debug)","none",pcbddc->dbg_flag,&pcbddc->dbg_flag,NULL);CHKERRQ(ierr); 59 /* Primal space cumstomization */ 60 ierr = PetscOptionsBool("-pc_bddc_use_vertices","Use or not corner dofs in coarse space","none",pcbddc->use_vertices,&pcbddc->use_vertices,NULL);CHKERRQ(ierr); 61 ierr = PetscOptionsBool("-pc_bddc_use_edges","Use or not edge constraints in coarse space","none",pcbddc->use_edges,&pcbddc->use_edges,NULL);CHKERRQ(ierr); 62 ierr = PetscOptionsBool("-pc_bddc_use_faces","Use or not face constraints in coarse space","none",pcbddc->use_faces,&pcbddc->use_faces,NULL);CHKERRQ(ierr); 63 /* Change of basis */ 64 ierr = PetscOptionsBool("-pc_bddc_use_change_of_basis","Use or not change of basis on local edge nodes","none",pcbddc->use_change_of_basis,&pcbddc->use_change_of_basis,NULL);CHKERRQ(ierr); 65 ierr = PetscOptionsBool("-pc_bddc_use_change_on_faces","Use or not change of basis on local face nodes","none",pcbddc->use_change_on_faces,&pcbddc->use_change_on_faces,NULL);CHKERRQ(ierr); 66 if (!pcbddc->use_change_of_basis) { 67 pcbddc->use_change_on_faces = PETSC_FALSE; 68 } 69 /* Switch between M_2 (default) and M_3 preconditioners (as defined by C. Dohrmann in the ref. article) */ 70 ierr = PetscOptionsBool("-pc_bddc_switch_static","Switch on static condensation ops around the interface preconditioner","none",pcbddc->switch_static,&pcbddc->switch_static,NULL);CHKERRQ(ierr); 71 ierr = PetscOptionsInt("-pc_bddc_coarsening_ratio","Set coarsening ratio used in multilevel coarsening","none",pcbddc->coarsening_ratio,&pcbddc->coarsening_ratio,NULL);CHKERRQ(ierr); 72 ierr = PetscOptionsInt("-pc_bddc_levels","Set maximum number of levels for multilevel","none",pcbddc->max_levels,&pcbddc->max_levels,NULL);CHKERRQ(ierr); 73 ierr = PetscOptionsBool("-pc_bddc_use_coarse_estimates","Use estimated eigenvalues for coarse problem","none",pcbddc->use_coarse_estimates,&pcbddc->use_coarse_estimates,NULL);CHKERRQ(ierr); 74 ierr = PetscOptionsBool("-pc_bddc_use_deluxe_scaling","Use deluxe scaling for BDDC","none",pcbddc->use_deluxe_scaling,&pcbddc->use_deluxe_scaling,NULL);CHKERRQ(ierr); 75 ierr = PetscOptionsTail();CHKERRQ(ierr); 76 PetscFunctionReturn(0); 77 } 78 /* -------------------------------------------------------------------------- */ 79 #undef __FUNCT__ 80 #define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS_BDDC" 81 static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices) 82 { 83 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 84 PetscErrorCode ierr; 85 86 PetscFunctionBegin; 87 ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr); 88 ierr = PetscObjectReference((PetscObject)PrimalVertices);CHKERRQ(ierr); 89 pcbddc->user_primal_vertices = PrimalVertices; 90 PetscFunctionReturn(0); 91 } 92 #undef __FUNCT__ 93 #define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS" 94 /*@ 95 PCBDDCSetPrimalVerticesLocalIS - Set additional user defined primal vertices in PCBDDC 96 97 Not collective 98 99 Input Parameters: 100 + pc - the preconditioning context 101 - PrimalVertices - index set of primal vertices in local numbering 102 103 Level: intermediate 104 105 Notes: 106 107 .seealso: PCBDDC 108 @*/ 109 PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices) 110 { 111 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 115 PetscValidHeaderSpecific(PrimalVertices,IS_CLASSID,2); 116 ierr = PetscTryMethod(pc,"PCBDDCSetPrimalVerticesLocalIS_C",(PC,IS),(pc,PrimalVertices));CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 /* -------------------------------------------------------------------------- */ 120 #undef __FUNCT__ 121 #define __FUNCT__ "PCBDDCSetCoarseningRatio_BDDC" 122 static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc,PetscInt k) 123 { 124 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 125 126 PetscFunctionBegin; 127 pcbddc->coarsening_ratio = k; 128 PetscFunctionReturn(0); 129 } 130 131 #undef __FUNCT__ 132 #define __FUNCT__ "PCBDDCSetCoarseningRatio" 133 /*@ 134 PCBDDCSetCoarseningRatio - Set coarsening ratio used in multilevel 135 136 Logically collective on PC 137 138 Input Parameters: 139 + pc - the preconditioning context 140 - k - coarsening ratio (H/h at the coarser level) 141 142 Approximatively k subdomains at the finer level will be aggregated into a single subdomain at the coarser level 143 144 Level: intermediate 145 146 Notes: 147 148 .seealso: PCBDDC 149 @*/ 150 PetscErrorCode PCBDDCSetCoarseningRatio(PC pc,PetscInt k) 151 { 152 PetscErrorCode ierr; 153 154 PetscFunctionBegin; 155 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 156 PetscValidLogicalCollectiveInt(pc,k,2); 157 ierr = PetscTryMethod(pc,"PCBDDCSetCoarseningRatio_C",(PC,PetscInt),(pc,k));CHKERRQ(ierr); 158 PetscFunctionReturn(0); 159 } 160 161 /* The following functions (PCBDDCSetUseExactDirichlet PCBDDCSetLevel) are not public */ 162 #undef __FUNCT__ 163 #define __FUNCT__ "PCBDDCSetUseExactDirichlet_BDDC" 164 static PetscErrorCode PCBDDCSetUseExactDirichlet_BDDC(PC pc,PetscBool flg) 165 { 166 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 167 168 PetscFunctionBegin; 169 pcbddc->use_exact_dirichlet_trick = flg; 170 PetscFunctionReturn(0); 171 } 172 173 #undef __FUNCT__ 174 #define __FUNCT__ "PCBDDCSetUseExactDirichlet" 175 PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool flg) 176 { 177 PetscErrorCode ierr; 178 179 PetscFunctionBegin; 180 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 181 PetscValidLogicalCollectiveBool(pc,flg,2); 182 ierr = PetscTryMethod(pc,"PCBDDCSetUseExactDirichlet_C",(PC,PetscBool),(pc,flg));CHKERRQ(ierr); 183 PetscFunctionReturn(0); 184 } 185 186 #undef __FUNCT__ 187 #define __FUNCT__ "PCBDDCSetLevel_BDDC" 188 static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc,PetscInt level) 189 { 190 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 191 192 PetscFunctionBegin; 193 pcbddc->current_level = level; 194 PetscFunctionReturn(0); 195 } 196 197 #undef __FUNCT__ 198 #define __FUNCT__ "PCBDDCSetLevel" 199 PetscErrorCode PCBDDCSetLevel(PC pc,PetscInt level) 200 { 201 PetscErrorCode ierr; 202 203 PetscFunctionBegin; 204 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 205 PetscValidLogicalCollectiveInt(pc,level,2); 206 ierr = PetscTryMethod(pc,"PCBDDCSetLevel_C",(PC,PetscInt),(pc,level));CHKERRQ(ierr); 207 PetscFunctionReturn(0); 208 } 209 210 #undef __FUNCT__ 211 #define __FUNCT__ "PCBDDCSetLevels_BDDC" 212 static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc,PetscInt levels) 213 { 214 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 215 216 PetscFunctionBegin; 217 pcbddc->max_levels = levels; 218 PetscFunctionReturn(0); 219 } 220 221 #undef __FUNCT__ 222 #define __FUNCT__ "PCBDDCSetLevels" 223 /*@ 224 PCBDDCSetLevels - Sets the maximum number of levels for multilevel 225 226 Logically collective on PC 227 228 Input Parameters: 229 + pc - the preconditioning context 230 - levels - the maximum number of levels (max 9) 231 232 Default value is 0, i.e. traditional one-level BDDC 233 234 Level: intermediate 235 236 Notes: 237 238 .seealso: PCBDDC 239 @*/ 240 PetscErrorCode PCBDDCSetLevels(PC pc,PetscInt levels) 241 { 242 PetscErrorCode ierr; 243 244 PetscFunctionBegin; 245 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 246 PetscValidLogicalCollectiveInt(pc,levels,2); 247 if (levels > 99) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Maximum number of levels for bddc is 99\n"); 248 ierr = PetscTryMethod(pc,"PCBDDCSetLevels_C",(PC,PetscInt),(pc,levels));CHKERRQ(ierr); 249 PetscFunctionReturn(0); 250 } 251 /* -------------------------------------------------------------------------- */ 252 253 #undef __FUNCT__ 254 #define __FUNCT__ "PCBDDCSetNullSpace_BDDC" 255 static PetscErrorCode PCBDDCSetNullSpace_BDDC(PC pc,MatNullSpace NullSpace) 256 { 257 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 258 PetscErrorCode ierr; 259 260 PetscFunctionBegin; 261 ierr = PetscObjectReference((PetscObject)NullSpace);CHKERRQ(ierr); 262 ierr = MatNullSpaceDestroy(&pcbddc->NullSpace);CHKERRQ(ierr); 263 pcbddc->NullSpace = NullSpace; 264 PetscFunctionReturn(0); 265 } 266 267 #undef __FUNCT__ 268 #define __FUNCT__ "PCBDDCSetNullSpace" 269 /*@ 270 PCBDDCSetNullSpace - Set nullspace for BDDC operator 271 272 Logically collective on PC and MatNullSpace 273 274 Input Parameters: 275 + pc - the preconditioning context 276 - NullSpace - Null space of the linear operator to be preconditioned (Pmat) 277 278 Level: intermediate 279 280 Notes: 281 282 .seealso: PCBDDC 283 @*/ 284 PetscErrorCode PCBDDCSetNullSpace(PC pc,MatNullSpace NullSpace) 285 { 286 PetscErrorCode ierr; 287 288 PetscFunctionBegin; 289 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 290 PetscValidHeaderSpecific(NullSpace,MAT_NULLSPACE_CLASSID,2); 291 PetscCheckSameComm(pc,1,NullSpace,2); 292 ierr = PetscTryMethod(pc,"PCBDDCSetNullSpace_C",(PC,MatNullSpace),(pc,NullSpace));CHKERRQ(ierr); 293 PetscFunctionReturn(0); 294 } 295 /* -------------------------------------------------------------------------- */ 296 297 #undef __FUNCT__ 298 #define __FUNCT__ "PCBDDCSetDirichletBoundaries_BDDC" 299 static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc,IS DirichletBoundaries) 300 { 301 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 302 PetscErrorCode ierr; 303 304 PetscFunctionBegin; 305 /* last user setting takes precendence -> destroy any other customization */ 306 ierr = ISDestroy(&pcbddc->DirichletBoundariesLocal);CHKERRQ(ierr); 307 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 308 ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr); 309 pcbddc->DirichletBoundaries = DirichletBoundaries; 310 pcbddc->recompute_topography = PETSC_TRUE; 311 PetscFunctionReturn(0); 312 } 313 314 #undef __FUNCT__ 315 #define __FUNCT__ "PCBDDCSetDirichletBoundaries" 316 /*@ 317 PCBDDCSetDirichletBoundaries - Set IS defining Dirichlet boundaries for the global problem. 318 319 Collective 320 321 Input Parameters: 322 + pc - the preconditioning context 323 - DirichletBoundaries - parallel IS defining the Dirichlet boundaries 324 325 Level: intermediate 326 327 Notes: Any process can list any global node 328 329 .seealso: PCBDDC 330 @*/ 331 PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc,IS DirichletBoundaries) 332 { 333 PetscErrorCode ierr; 334 335 PetscFunctionBegin; 336 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 337 PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2); 338 PetscCheckSameComm(pc,1,DirichletBoundaries,2); 339 ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundaries_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr); 340 PetscFunctionReturn(0); 341 } 342 /* -------------------------------------------------------------------------- */ 343 344 #undef __FUNCT__ 345 #define __FUNCT__ "PCBDDCSetDirichletBoundariesLocal_BDDC" 346 static PetscErrorCode PCBDDCSetDirichletBoundariesLocal_BDDC(PC pc,IS DirichletBoundaries) 347 { 348 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 349 PetscErrorCode ierr; 350 351 PetscFunctionBegin; 352 /* last user setting takes precendence -> destroy any other customization */ 353 ierr = ISDestroy(&pcbddc->DirichletBoundariesLocal);CHKERRQ(ierr); 354 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 355 ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr); 356 pcbddc->DirichletBoundariesLocal = DirichletBoundaries; 357 pcbddc->recompute_topography = PETSC_TRUE; 358 PetscFunctionReturn(0); 359 } 360 361 #undef __FUNCT__ 362 #define __FUNCT__ "PCBDDCSetDirichletBoundariesLocal" 363 /*@ 364 PCBDDCSetDirichletBoundariesLocal - Set IS defining Dirichlet boundaries for the global problem in local ordering. 365 366 Collective 367 368 Input Parameters: 369 + pc - the preconditioning context 370 - DirichletBoundaries - parallel IS defining the Dirichlet boundaries (in local ordering) 371 372 Level: intermediate 373 374 Notes: 375 376 .seealso: PCBDDC 377 @*/ 378 PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc,IS DirichletBoundaries) 379 { 380 PetscErrorCode ierr; 381 382 PetscFunctionBegin; 383 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 384 PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2); 385 PetscCheckSameComm(pc,1,DirichletBoundaries,2); 386 ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundariesLocal_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr); 387 PetscFunctionReturn(0); 388 } 389 /* -------------------------------------------------------------------------- */ 390 391 #undef __FUNCT__ 392 #define __FUNCT__ "PCBDDCSetNeumannBoundaries_BDDC" 393 static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc,IS NeumannBoundaries) 394 { 395 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 396 PetscErrorCode ierr; 397 398 PetscFunctionBegin; 399 /* last user setting takes precendence -> destroy any other customization */ 400 ierr = ISDestroy(&pcbddc->NeumannBoundariesLocal);CHKERRQ(ierr); 401 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 402 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 403 pcbddc->NeumannBoundaries = NeumannBoundaries; 404 pcbddc->recompute_topography = PETSC_TRUE; 405 PetscFunctionReturn(0); 406 } 407 408 #undef __FUNCT__ 409 #define __FUNCT__ "PCBDDCSetNeumannBoundaries" 410 /*@ 411 PCBDDCSetNeumannBoundaries - Set IS defining Neumann boundaries for the global problem. 412 413 Collective 414 415 Input Parameters: 416 + pc - the preconditioning context 417 - NeumannBoundaries - parallel IS defining the Neumann boundaries 418 419 Level: intermediate 420 421 Notes: Any process can list any global node 422 423 .seealso: PCBDDC 424 @*/ 425 PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc,IS NeumannBoundaries) 426 { 427 PetscErrorCode ierr; 428 429 PetscFunctionBegin; 430 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 431 PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2); 432 PetscCheckSameComm(pc,1,NeumannBoundaries,2); 433 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundaries_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 434 PetscFunctionReturn(0); 435 } 436 /* -------------------------------------------------------------------------- */ 437 438 #undef __FUNCT__ 439 #define __FUNCT__ "PCBDDCSetNeumannBoundariesLocal_BDDC" 440 static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc,IS NeumannBoundaries) 441 { 442 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 443 PetscErrorCode ierr; 444 445 PetscFunctionBegin; 446 /* last user setting takes precendence -> destroy any other customization */ 447 ierr = ISDestroy(&pcbddc->NeumannBoundariesLocal);CHKERRQ(ierr); 448 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 449 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 450 pcbddc->NeumannBoundariesLocal = NeumannBoundaries; 451 pcbddc->recompute_topography = PETSC_TRUE; 452 PetscFunctionReturn(0); 453 } 454 455 #undef __FUNCT__ 456 #define __FUNCT__ "PCBDDCSetNeumannBoundariesLocal" 457 /*@ 458 PCBDDCSetNeumannBoundariesLocal - Set IS defining Neumann boundaries for the global problem in local ordering. 459 460 Collective 461 462 Input Parameters: 463 + pc - the preconditioning context 464 - NeumannBoundaries - parallel IS defining the subdomain part of Neumann boundaries (in local ordering) 465 466 Level: intermediate 467 468 Notes: 469 470 .seealso: PCBDDC 471 @*/ 472 PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc,IS NeumannBoundaries) 473 { 474 PetscErrorCode ierr; 475 476 PetscFunctionBegin; 477 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 478 PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2); 479 PetscCheckSameComm(pc,1,NeumannBoundaries,2); 480 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundariesLocal_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 481 PetscFunctionReturn(0); 482 } 483 /* -------------------------------------------------------------------------- */ 484 485 #undef __FUNCT__ 486 #define __FUNCT__ "PCBDDCGetDirichletBoundaries_BDDC" 487 static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc,IS *DirichletBoundaries) 488 { 489 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 490 491 PetscFunctionBegin; 492 *DirichletBoundaries = pcbddc->DirichletBoundaries; 493 PetscFunctionReturn(0); 494 } 495 496 #undef __FUNCT__ 497 #define __FUNCT__ "PCBDDCGetDirichletBoundaries" 498 /*@ 499 PCBDDCGetDirichletBoundaries - Get parallel IS for Dirichlet boundaries 500 501 Collective 502 503 Input Parameters: 504 . pc - the preconditioning context 505 506 Output Parameters: 507 . DirichletBoundaries - index set defining the Dirichlet boundaries 508 509 Level: intermediate 510 511 Notes: The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetDirichletBoundaries 512 513 .seealso: PCBDDC 514 @*/ 515 PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc,IS *DirichletBoundaries) 516 { 517 PetscErrorCode ierr; 518 519 PetscFunctionBegin; 520 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 521 ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundaries_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr); 522 PetscFunctionReturn(0); 523 } 524 /* -------------------------------------------------------------------------- */ 525 526 #undef __FUNCT__ 527 #define __FUNCT__ "PCBDDCGetDirichletBoundariesLocal_BDDC" 528 static PetscErrorCode PCBDDCGetDirichletBoundariesLocal_BDDC(PC pc,IS *DirichletBoundaries) 529 { 530 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 531 532 PetscFunctionBegin; 533 *DirichletBoundaries = pcbddc->DirichletBoundariesLocal; 534 PetscFunctionReturn(0); 535 } 536 537 #undef __FUNCT__ 538 #define __FUNCT__ "PCBDDCGetDirichletBoundariesLocal" 539 /*@ 540 PCBDDCGetDirichletBoundariesLocal - Get parallel IS for Dirichlet boundaries (in local ordering) 541 542 Collective 543 544 Input Parameters: 545 . pc - the preconditioning context 546 547 Output Parameters: 548 . DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries 549 550 Level: intermediate 551 552 Notes: 553 554 .seealso: PCBDDC 555 @*/ 556 PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc,IS *DirichletBoundaries) 557 { 558 PetscErrorCode ierr; 559 560 PetscFunctionBegin; 561 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 562 ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundariesLocal_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr); 563 PetscFunctionReturn(0); 564 } 565 /* -------------------------------------------------------------------------- */ 566 567 #undef __FUNCT__ 568 #define __FUNCT__ "PCBDDCGetNeumannBoundaries_BDDC" 569 static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc,IS *NeumannBoundaries) 570 { 571 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 572 573 PetscFunctionBegin; 574 *NeumannBoundaries = pcbddc->NeumannBoundaries; 575 PetscFunctionReturn(0); 576 } 577 578 #undef __FUNCT__ 579 #define __FUNCT__ "PCBDDCGetNeumannBoundaries" 580 /*@ 581 PCBDDCGetNeumannBoundaries - Get parallel IS for Neumann boundaries 582 583 Collective 584 585 Input Parameters: 586 . pc - the preconditioning context 587 588 Output Parameters: 589 . NeumannBoundaries - index set defining the Neumann boundaries 590 591 Level: intermediate 592 593 Notes: The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetNeumannBoundaries 594 595 .seealso: PCBDDC 596 @*/ 597 PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc,IS *NeumannBoundaries) 598 { 599 PetscErrorCode ierr; 600 601 PetscFunctionBegin; 602 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 603 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundaries_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 604 PetscFunctionReturn(0); 605 } 606 /* -------------------------------------------------------------------------- */ 607 608 #undef __FUNCT__ 609 #define __FUNCT__ "PCBDDCGetNeumannBoundariesLocal_BDDC" 610 static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc,IS *NeumannBoundaries) 611 { 612 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 613 614 PetscFunctionBegin; 615 *NeumannBoundaries = pcbddc->NeumannBoundariesLocal; 616 PetscFunctionReturn(0); 617 } 618 619 #undef __FUNCT__ 620 #define __FUNCT__ "PCBDDCGetNeumannBoundariesLocal" 621 /*@ 622 PCBDDCGetNeumannBoundariesLocal - Get parallel IS for Neumann boundaries (in local ordering) 623 624 Collective 625 626 Input Parameters: 627 . pc - the preconditioning context 628 629 Output Parameters: 630 . NeumannBoundaries - index set defining the subdomain part of Neumann boundaries 631 632 Level: intermediate 633 634 Notes: 635 636 .seealso: PCBDDC 637 @*/ 638 PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc,IS *NeumannBoundaries) 639 { 640 PetscErrorCode ierr; 641 642 PetscFunctionBegin; 643 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 644 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundariesLocal_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 645 PetscFunctionReturn(0); 646 } 647 /* -------------------------------------------------------------------------- */ 648 649 #undef __FUNCT__ 650 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph_BDDC" 651 static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 652 { 653 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 654 PCBDDCGraph mat_graph = pcbddc->mat_graph; 655 PetscErrorCode ierr; 656 657 PetscFunctionBegin; 658 /* free old CSR */ 659 ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr); 660 /* TODO: PCBDDCGraphSetAdjacency */ 661 /* get CSR into graph structure */ 662 if (copymode == PETSC_COPY_VALUES) { 663 ierr = PetscMalloc1((nvtxs+1),&mat_graph->xadj);CHKERRQ(ierr); 664 ierr = PetscMalloc1(xadj[nvtxs],&mat_graph->adjncy);CHKERRQ(ierr); 665 ierr = PetscMemcpy(mat_graph->xadj,xadj,(nvtxs+1)*sizeof(PetscInt));CHKERRQ(ierr); 666 ierr = PetscMemcpy(mat_graph->adjncy,adjncy,xadj[nvtxs]*sizeof(PetscInt));CHKERRQ(ierr); 667 } else if (copymode == PETSC_OWN_POINTER) { 668 mat_graph->xadj = (PetscInt*)xadj; 669 mat_graph->adjncy = (PetscInt*)adjncy; 670 } 671 mat_graph->nvtxs_csr = nvtxs; 672 PetscFunctionReturn(0); 673 } 674 675 #undef __FUNCT__ 676 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph" 677 /*@ 678 PCBDDCSetLocalAdjacencyGraph - Set adjacency structure (CSR graph) of the local Neumann matrix 679 680 Not collective 681 682 Input Parameters: 683 + pc - the preconditioning context 684 . nvtxs - number of local vertices of the graph (i.e., the local size of your problem) 685 . xadj, adjncy - the CSR graph 686 - copymode - either PETSC_COPY_VALUES or PETSC_OWN_POINTER. 687 688 Level: intermediate 689 690 Notes: 691 692 .seealso: PCBDDC,PetscCopyMode 693 @*/ 694 PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc,PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 695 { 696 void (*f)(void) = 0; 697 PetscErrorCode ierr; 698 699 PetscFunctionBegin; 700 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 701 PetscValidIntPointer(xadj,3); 702 PetscValidIntPointer(xadj,4); 703 if (copymode != PETSC_COPY_VALUES && copymode != PETSC_OWN_POINTER) { 704 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported copy mode %d in %s\n",copymode,__FUNCT__); 705 } 706 ierr = PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));CHKERRQ(ierr); 707 /* free arrays if PCBDDC is not the PC type */ 708 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);CHKERRQ(ierr); 709 if (!f && copymode == PETSC_OWN_POINTER) { 710 ierr = PetscFree(xadj);CHKERRQ(ierr); 711 ierr = PetscFree(adjncy);CHKERRQ(ierr); 712 } 713 PetscFunctionReturn(0); 714 } 715 /* -------------------------------------------------------------------------- */ 716 717 #undef __FUNCT__ 718 #define __FUNCT__ "PCBDDCSetDofsSplittingLocal_BDDC" 719 static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 720 { 721 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 722 PetscInt i; 723 PetscErrorCode ierr; 724 725 PetscFunctionBegin; 726 /* Destroy ISes if they were already set */ 727 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 728 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 729 } 730 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 731 /* last user setting takes precendence -> destroy any other customization */ 732 for (i=0;i<pcbddc->n_ISForDofs;i++) { 733 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 734 } 735 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 736 pcbddc->n_ISForDofs = 0; 737 /* allocate space then set */ 738 ierr = PetscMalloc(n_is*sizeof(IS),&pcbddc->ISForDofsLocal);CHKERRQ(ierr); 739 for (i=0;i<n_is;i++) { 740 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 741 pcbddc->ISForDofsLocal[i]=ISForDofs[i]; 742 } 743 pcbddc->n_ISForDofsLocal=n_is; 744 if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE; 745 pcbddc->recompute_topography = PETSC_TRUE; 746 PetscFunctionReturn(0); 747 } 748 749 #undef __FUNCT__ 750 #define __FUNCT__ "PCBDDCSetDofsSplittingLocal" 751 /*@ 752 PCBDDCSetDofsSplittingLocal - Set index sets defining fields of the local subdomain matrix 753 754 Collective 755 756 Input Parameters: 757 + pc - the preconditioning context 758 - n_is - number of index sets defining the fields 759 . ISForDofs - array of IS describing the fields in local ordering 760 761 Level: intermediate 762 763 Notes: n_is should be the same among processes. Not all nodes need to be listed: unlisted nodes will belong to a different field. 764 765 .seealso: PCBDDC 766 @*/ 767 PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc,PetscInt n_is, IS ISForDofs[]) 768 { 769 PetscInt i; 770 PetscErrorCode ierr; 771 772 PetscFunctionBegin; 773 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 774 PetscValidLogicalCollectiveInt(pc,n_is,2); 775 for (i=0;i<n_is;i++) { 776 PetscCheckSameComm(pc,1,ISForDofs[i],3); 777 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 778 } 779 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 780 PetscFunctionReturn(0); 781 } 782 /* -------------------------------------------------------------------------- */ 783 784 #undef __FUNCT__ 785 #define __FUNCT__ "PCBDDCSetDofsSplitting_BDDC" 786 static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 787 { 788 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 789 PetscInt i; 790 PetscErrorCode ierr; 791 792 PetscFunctionBegin; 793 /* Destroy ISes if they were already set */ 794 for (i=0;i<pcbddc->n_ISForDofs;i++) { 795 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 796 } 797 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 798 /* last user setting takes precendence -> destroy any other customization */ 799 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 800 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 801 } 802 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 803 pcbddc->n_ISForDofsLocal = 0; 804 /* allocate space then set */ 805 ierr = PetscMalloc1(n_is,&pcbddc->ISForDofs);CHKERRQ(ierr); 806 for (i=0;i<n_is;i++) { 807 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 808 pcbddc->ISForDofs[i]=ISForDofs[i]; 809 } 810 pcbddc->n_ISForDofs=n_is; 811 if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE; 812 pcbddc->recompute_topography = PETSC_TRUE; 813 PetscFunctionReturn(0); 814 } 815 816 #undef __FUNCT__ 817 #define __FUNCT__ "PCBDDCSetDofsSplitting" 818 /*@ 819 PCBDDCSetDofsSplitting - Set index sets defining fields of the global matrix 820 821 Collective 822 823 Input Parameters: 824 + pc - the preconditioning context 825 - n_is - number of index sets defining the fields 826 . ISForDofs - array of IS describing the fields in global ordering 827 828 Level: intermediate 829 830 Notes: Any process can list any global node. Not all nodes need to be listed: unlisted nodes will belong to a different field. 831 832 .seealso: PCBDDC 833 @*/ 834 PetscErrorCode PCBDDCSetDofsSplitting(PC pc,PetscInt n_is, IS ISForDofs[]) 835 { 836 PetscInt i; 837 PetscErrorCode ierr; 838 839 PetscFunctionBegin; 840 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 841 PetscValidLogicalCollectiveInt(pc,n_is,2); 842 for (i=0;i<n_is;i++) { 843 PetscCheckSameComm(pc,1,ISForDofs[i],3); 844 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 845 } 846 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 847 PetscFunctionReturn(0); 848 } 849 /* -------------------------------------------------------------------------- */ 850 #undef __FUNCT__ 851 #define __FUNCT__ "PCPreSolve_BDDC" 852 /* -------------------------------------------------------------------------- */ 853 /* 854 PCPreSolve_BDDC - Changes the right hand side and (if necessary) the initial 855 guess if a transformation of basis approach has been selected. 856 857 Input Parameter: 858 + pc - the preconditioner contex 859 860 Application Interface Routine: PCPreSolve() 861 862 Notes: 863 The interface routine PCPreSolve() is not usually called directly by 864 the user, but instead is called by KSPSolve(). 865 */ 866 static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 867 { 868 PetscErrorCode ierr; 869 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 870 PC_IS *pcis = (PC_IS*)(pc->data); 871 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 872 Mat temp_mat; 873 IS dirIS; 874 Vec used_vec; 875 PetscBool guess_nonzero; 876 877 PetscFunctionBegin; 878 /* if we are working with cg, one dirichlet solve can be avoided during Krylov iterations */ 879 if (ksp) { 880 PetscBool iscg; 881 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPCG,&iscg);CHKERRQ(ierr); 882 if (!iscg) { 883 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 884 } 885 } 886 /* Creates parallel work vectors used in presolve */ 887 if (!pcbddc->original_rhs) { 888 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 889 } 890 if (!pcbddc->temp_solution) { 891 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);CHKERRQ(ierr); 892 } 893 if (x) { 894 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr); 895 used_vec = x; 896 } else { 897 ierr = PetscObjectReference((PetscObject)pcbddc->temp_solution);CHKERRQ(ierr); 898 used_vec = pcbddc->temp_solution; 899 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 900 } 901 /* hack into ksp data structure PCPreSolve comes earlier in src/ksp/ksp/interface/itfunc.c */ 902 if (ksp) { 903 ierr = KSPGetInitialGuessNonzero(ksp,&guess_nonzero);CHKERRQ(ierr); 904 if (!guess_nonzero) { 905 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 906 } 907 } 908 909 /* store the original rhs */ 910 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 911 912 /* Take into account zeroed rows -> change rhs and store solution removed */ 913 /* note that Dirichlet boundaries in global ordering (if any) has already been translated into local ordering in PCBDDCAnalyzeInterface */ 914 ierr = PCBDDCGetDirichletBoundariesLocal(pc,&dirIS);CHKERRQ(ierr); 915 if (rhs && dirIS) { 916 PetscInt dirsize,i,*is_indices; 917 PetscScalar *array_x,*array_diagonal; 918 919 ierr = MatGetDiagonal(pc->pmat,pcis->vec1_global);CHKERRQ(ierr); 920 ierr = VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);CHKERRQ(ierr); 921 ierr = VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 922 ierr = VecScatterEnd(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 923 ierr = VecScatterBegin(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 924 ierr = VecScatterEnd(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 925 ierr = ISGetLocalSize(dirIS,&dirsize);CHKERRQ(ierr); 926 ierr = VecGetArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 927 ierr = VecGetArray(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 928 ierr = ISGetIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 929 for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]]; 930 ierr = ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 931 ierr = VecRestoreArray(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 932 ierr = VecRestoreArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 933 ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 934 ierr = VecScatterEnd(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 935 936 /* remove the computed solution from the rhs */ 937 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 938 ierr = MatMultAdd(pc->pmat,used_vec,rhs,rhs);CHKERRQ(ierr); 939 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 940 } 941 942 /* store partially computed solution and set initial guess */ 943 if (x) { 944 ierr = VecCopy(used_vec,pcbddc->temp_solution);CHKERRQ(ierr); 945 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 946 if (pcbddc->use_exact_dirichlet_trick) { 947 ierr = VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 948 ierr = VecScatterEnd (pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 949 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 950 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 951 ierr = VecScatterEnd (pcis->global_to_D,pcis->vec2_D,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 952 if (ksp) { 953 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 954 } 955 } 956 } 957 958 /* prepare MatMult and rhs for solver */ 959 if (pcbddc->use_change_of_basis) { 960 /* swap pointers for local matrices */ 961 temp_mat = matis->A; 962 matis->A = pcbddc->local_mat; 963 pcbddc->local_mat = temp_mat; 964 if (rhs) { 965 /* Get local rhs and apply transformation of basis */ 966 ierr = VecScatterBegin(pcis->global_to_B,rhs,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 967 ierr = VecScatterEnd (pcis->global_to_B,rhs,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 968 /* from original basis to modified basis */ 969 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,pcis->vec1_B,pcis->vec2_B);CHKERRQ(ierr); 970 /* put back modified values into the global vec using INSERT_VALUES copy mode */ 971 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec2_B,rhs,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 972 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec2_B,rhs,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 973 } 974 } 975 976 /* remove nullspace if present */ 977 if (ksp && pcbddc->NullSpace) { 978 ierr = MatNullSpaceRemove(pcbddc->NullSpace,used_vec);CHKERRQ(ierr); 979 ierr = MatNullSpaceRemove(pcbddc->NullSpace,rhs);CHKERRQ(ierr); 980 } 981 ierr = VecDestroy(&used_vec);CHKERRQ(ierr); 982 PetscFunctionReturn(0); 983 } 984 /* -------------------------------------------------------------------------- */ 985 #undef __FUNCT__ 986 #define __FUNCT__ "PCPostSolve_BDDC" 987 /* -------------------------------------------------------------------------- */ 988 /* 989 PCPostSolve_BDDC - Changes the computed solution if a transformation of basis 990 approach has been selected. Also, restores rhs to its original state. 991 992 Input Parameter: 993 + pc - the preconditioner contex 994 995 Application Interface Routine: PCPostSolve() 996 997 Notes: 998 The interface routine PCPostSolve() is not usually called directly by 999 the user, but instead is called by KSPSolve(). 1000 */ 1001 static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 1002 { 1003 PetscErrorCode ierr; 1004 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1005 PC_IS *pcis = (PC_IS*)(pc->data); 1006 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 1007 Mat temp_mat; 1008 1009 PetscFunctionBegin; 1010 if (pcbddc->use_change_of_basis) { 1011 /* swap pointers for local matrices */ 1012 temp_mat = matis->A; 1013 matis->A = pcbddc->local_mat; 1014 pcbddc->local_mat = temp_mat; 1015 } 1016 if (pcbddc->use_change_of_basis && x) { 1017 /* Get Local boundary and apply transformation of basis to solution vector */ 1018 ierr = VecScatterBegin(pcis->global_to_B,x,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1019 ierr = VecScatterEnd (pcis->global_to_B,x,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1020 /* from modified basis to original basis */ 1021 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcis->vec1_B,pcis->vec2_B);CHKERRQ(ierr); 1022 /* put back modified values into the global vec using INSERT_VALUES copy mode */ 1023 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec2_B,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1024 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec2_B,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1025 } 1026 /* add solution removed in presolve */ 1027 if (x) { 1028 ierr = VecAXPY(x,1.0,pcbddc->temp_solution);CHKERRQ(ierr); 1029 } 1030 /* restore rhs to its original state */ 1031 if (rhs) { 1032 ierr = VecCopy(pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1033 } 1034 PetscFunctionReturn(0); 1035 } 1036 /* -------------------------------------------------------------------------- */ 1037 #undef __FUNCT__ 1038 #define __FUNCT__ "PCSetUp_BDDC" 1039 /* -------------------------------------------------------------------------- */ 1040 /* 1041 PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner 1042 by setting data structures and options. 1043 1044 Input Parameter: 1045 + pc - the preconditioner context 1046 1047 Application Interface Routine: PCSetUp() 1048 1049 Notes: 1050 The interface routine PCSetUp() is not usually called directly by 1051 the user, but instead is called by PCApply() if necessary. 1052 */ 1053 PetscErrorCode PCSetUp_BDDC(PC pc) 1054 { 1055 PetscErrorCode ierr; 1056 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 1057 MatNullSpace nearnullspace; 1058 PetscBool computeis,computetopography,computesolvers; 1059 PetscBool new_nearnullspace_provided; 1060 1061 PetscFunctionBegin; 1062 /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */ 1063 /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup 1064 Also, BDDC directly build the Dirichlet problem */ 1065 1066 /* split work */ 1067 if (pc->setupcalled) { 1068 computeis = PETSC_FALSE; 1069 if (pc->flag == SAME_NONZERO_PATTERN) { 1070 computetopography = PETSC_FALSE; 1071 computesolvers = PETSC_TRUE; 1072 } else { /* DIFFERENT_NONZERO_PATTERN */ 1073 computetopography = PETSC_TRUE; 1074 computesolvers = PETSC_TRUE; 1075 } 1076 } else { 1077 computeis = PETSC_TRUE; 1078 computetopography = PETSC_TRUE; 1079 computesolvers = PETSC_TRUE; 1080 } 1081 if (pcbddc->recompute_topography) { 1082 computetopography = PETSC_TRUE; 1083 } 1084 1085 /* Get stdout for dbg */ 1086 if (pcbddc->dbg_flag) { 1087 if (!pcbddc->dbg_viewer) { 1088 pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)); 1089 ierr = PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr); 1090 } 1091 ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1092 } 1093 1094 /* Set up all the "iterative substructuring" common block without computing solvers */ 1095 if (computeis) { 1096 /* HACK INTO PCIS */ 1097 PC_IS* pcis = (PC_IS*)pc->data; 1098 pcis->computesolvers = PETSC_FALSE; 1099 ierr = PCISSetUp(pc);CHKERRQ(ierr); 1100 ierr = ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&pcbddc->BtoNmap);CHKERRQ(ierr); 1101 } 1102 1103 /* Analyze interface */ 1104 if (computetopography) { 1105 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 1106 } 1107 1108 /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */ 1109 new_nearnullspace_provided = PETSC_FALSE; 1110 ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr); 1111 if (pcbddc->onearnullspace) { /* already used nearnullspace */ 1112 if (!nearnullspace) { /* near null space attached to mat has been destroyed */ 1113 new_nearnullspace_provided = PETSC_TRUE; 1114 } else { 1115 /* determine if the two nullspaces are different (should be lightweight) */ 1116 if (nearnullspace != pcbddc->onearnullspace) { 1117 new_nearnullspace_provided = PETSC_TRUE; 1118 } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */ 1119 PetscInt i; 1120 const Vec *nearnullvecs; 1121 PetscObjectState state; 1122 PetscInt nnsp_size; 1123 ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 1124 for (i=0;i<nnsp_size;i++) { 1125 ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr); 1126 if (pcbddc->onearnullvecs_state[i] != state) { 1127 new_nearnullspace_provided = PETSC_TRUE; 1128 break; 1129 } 1130 } 1131 } 1132 } 1133 } else { 1134 if (!nearnullspace) { /* both nearnullspaces are null */ 1135 new_nearnullspace_provided = PETSC_FALSE; 1136 } else { /* nearnullspace attached later */ 1137 new_nearnullspace_provided = PETSC_TRUE; 1138 } 1139 } 1140 1141 /* Setup constraints and related work vectors */ 1142 /* reset primal space flags */ 1143 pcbddc->new_primal_space = PETSC_FALSE; 1144 pcbddc->new_primal_space_local = PETSC_FALSE; 1145 if (computetopography || new_nearnullspace_provided) { 1146 /* It also sets the primal space flags */ 1147 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 1148 /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */ 1149 ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr); 1150 } 1151 1152 if (computesolvers || pcbddc->new_primal_space) { 1153 /* reset data */ 1154 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 1155 /* Create coarse and local stuffs */ 1156 ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr); 1157 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1158 } 1159 1160 if (pcbddc->dbg_flag) { 1161 ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1162 } 1163 PetscFunctionReturn(0); 1164 } 1165 1166 /* -------------------------------------------------------------------------- */ 1167 /* 1168 PCApply_BDDC - Applies the BDDC operator to a vector. 1169 1170 Input Parameters: 1171 . pc - the preconditioner context 1172 . r - input vector (global) 1173 1174 Output Parameter: 1175 . z - output vector (global) 1176 1177 Application Interface Routine: PCApply() 1178 */ 1179 #undef __FUNCT__ 1180 #define __FUNCT__ "PCApply_BDDC" 1181 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 1182 { 1183 PC_IS *pcis = (PC_IS*)(pc->data); 1184 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1185 PetscErrorCode ierr; 1186 const PetscScalar one = 1.0; 1187 const PetscScalar m_one = -1.0; 1188 const PetscScalar zero = 0.0; 1189 1190 /* This code is similar to that provided in nn.c for PCNN 1191 NN interface preconditioner changed to BDDC 1192 Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static = PETSC_TRUE) */ 1193 1194 PetscFunctionBegin; 1195 if (!pcbddc->use_exact_dirichlet_trick) { 1196 /* First Dirichlet solve */ 1197 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1198 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1199 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1200 /* 1201 Assembling right hand side for BDDC operator 1202 - pcis->vec1_D for the Dirichlet part (if needed, i.e. prec_flag=PETSC_TRUE) 1203 - pcis->vec1_B the interface part of the global vector z 1204 */ 1205 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1206 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1207 if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 1208 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1209 ierr = VecCopy(r,z);CHKERRQ(ierr); 1210 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1211 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1212 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1213 } else { 1214 ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr); 1215 ierr = VecSet(pcis->vec2_D,zero);CHKERRQ(ierr); 1216 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1217 } 1218 1219 /* Apply interface preconditioner 1220 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1221 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr); 1222 1223 /* Apply transpose of partition of unity operator */ 1224 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1225 1226 /* Second Dirichlet solve and assembling of output */ 1227 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1228 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1229 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1230 if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 1231 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1232 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1233 if (pcbddc->switch_static) { ierr = VecAXPY(pcis->vec4_D,one,pcis->vec1_D);CHKERRQ(ierr); } 1234 ierr = VecAXPY(pcis->vec2_D,one,pcis->vec4_D);CHKERRQ(ierr); 1235 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1236 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1237 PetscFunctionReturn(0); 1238 } 1239 1240 /* -------------------------------------------------------------------------- */ 1241 /* 1242 PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector. 1243 1244 Input Parameters: 1245 . pc - the preconditioner context 1246 . r - input vector (global) 1247 1248 Output Parameter: 1249 . z - output vector (global) 1250 1251 Application Interface Routine: PCApplyTranspose() 1252 */ 1253 #undef __FUNCT__ 1254 #define __FUNCT__ "PCApplyTranspose_BDDC" 1255 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z) 1256 { 1257 PC_IS *pcis = (PC_IS*)(pc->data); 1258 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1259 PetscErrorCode ierr; 1260 const PetscScalar one = 1.0; 1261 const PetscScalar m_one = -1.0; 1262 const PetscScalar zero = 0.0; 1263 1264 PetscFunctionBegin; 1265 if (!pcbddc->use_exact_dirichlet_trick) { 1266 /* First Dirichlet solve */ 1267 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1268 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1269 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1270 /* 1271 Assembling right hand side for BDDC operator 1272 - pcis->vec1_D for the Dirichlet part (if needed, i.e. prec_flag=PETSC_TRUE) 1273 - pcis->vec1_B the interface part of the global vector z 1274 */ 1275 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1276 ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1277 if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 1278 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1279 ierr = VecCopy(r,z);CHKERRQ(ierr); 1280 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1281 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1282 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1283 } else { 1284 ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr); 1285 ierr = VecSet(pcis->vec2_D,zero);CHKERRQ(ierr); 1286 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1287 } 1288 1289 /* Apply interface preconditioner 1290 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1291 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr); 1292 1293 /* Apply transpose of partition of unity operator */ 1294 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1295 1296 /* Second Dirichlet solve and assembling of output */ 1297 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1298 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1299 ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1300 if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 1301 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1302 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1303 if (pcbddc->switch_static) { ierr = VecAXPY(pcis->vec4_D,one,pcis->vec1_D);CHKERRQ(ierr); } 1304 ierr = VecAXPY(pcis->vec2_D,one,pcis->vec4_D);CHKERRQ(ierr); 1305 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1306 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1307 PetscFunctionReturn(0); 1308 } 1309 /* -------------------------------------------------------------------------- */ 1310 1311 #undef __FUNCT__ 1312 #define __FUNCT__ "PCDestroy_BDDC" 1313 PetscErrorCode PCDestroy_BDDC(PC pc) 1314 { 1315 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1316 PetscErrorCode ierr; 1317 1318 PetscFunctionBegin; 1319 /* free data created by PCIS */ 1320 ierr = PCISDestroy(pc);CHKERRQ(ierr); 1321 /* free BDDC custom data */ 1322 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 1323 /* destroy objects related to topography */ 1324 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 1325 /* free allocated graph structure */ 1326 ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr); 1327 /* free data for scaling operator */ 1328 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 1329 /* free solvers stuff */ 1330 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 1331 /* free global vectors needed in presolve */ 1332 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 1333 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 1334 ierr = ISLocalToGlobalMappingDestroy(&pcbddc->BtoNmap);CHKERRQ(ierr); 1335 /* remove functions */ 1336 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 1337 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 1338 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr); 1339 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr); 1340 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr); 1341 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",NULL);CHKERRQ(ierr); 1342 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 1343 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 1344 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 1345 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 1346 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 1347 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 1348 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 1349 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 1350 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 1351 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr); 1352 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 1353 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 1354 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 1355 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 1356 /* Free the private data structure */ 1357 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1358 PetscFunctionReturn(0); 1359 } 1360 /* -------------------------------------------------------------------------- */ 1361 1362 #undef __FUNCT__ 1363 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC" 1364 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 1365 { 1366 FETIDPMat_ctx mat_ctx; 1367 PC_IS* pcis; 1368 PC_BDDC* pcbddc; 1369 PetscErrorCode ierr; 1370 1371 PetscFunctionBegin; 1372 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1373 pcis = (PC_IS*)mat_ctx->pc->data; 1374 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 1375 1376 /* change of basis for physical rhs if needed 1377 It also changes the rhs in case of dirichlet boundaries */ 1378 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,standard_rhs,NULL);CHKERRQ(ierr); 1379 /* store vectors for computation of fetidp final solution */ 1380 ierr = VecScatterBegin(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1381 ierr = VecScatterEnd(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1382 /* scale rhs since it should be unassembled */ 1383 /* TODO use counter scaling? (also below) */ 1384 ierr = VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1385 ierr = VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1386 /* Apply partition of unity */ 1387 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1388 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 1389 if (!pcbddc->switch_static) { 1390 /* compute partially subassembled Schur complement right-hand side */ 1391 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1392 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 1393 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 1394 ierr = VecSet(standard_rhs,0.0);CHKERRQ(ierr); 1395 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1396 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1397 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 1398 ierr = VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1399 ierr = VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1400 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1401 } 1402 /* BDDC rhs */ 1403 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 1404 if (pcbddc->switch_static) { 1405 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1406 } 1407 /* apply BDDC */ 1408 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 1409 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 1410 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 1411 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 1412 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1413 ierr = VecScatterEnd (mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1414 /* restore original rhs */ 1415 ierr = VecCopy(pcbddc->original_rhs,standard_rhs);CHKERRQ(ierr); 1416 PetscFunctionReturn(0); 1417 } 1418 1419 #undef __FUNCT__ 1420 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS" 1421 /*@ 1422 PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETIDP linear system 1423 1424 Collective 1425 1426 Input Parameters: 1427 + fetidp_mat - the FETIDP matrix object obtained by calling PCBDDCCreateFETIDPOperators 1428 . standard_rhs - the right-hand side for your linear system 1429 1430 Output Parameters: 1431 - fetidp_flux_rhs - the right-hand side for the FETIDP linear system 1432 1433 Level: developer 1434 1435 Notes: 1436 1437 .seealso: PCBDDC,PCBDDCCreateFETIDPOperators 1438 @*/ 1439 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 1440 { 1441 FETIDPMat_ctx mat_ctx; 1442 PetscErrorCode ierr; 1443 1444 PetscFunctionBegin; 1445 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1446 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 1447 PetscFunctionReturn(0); 1448 } 1449 /* -------------------------------------------------------------------------- */ 1450 1451 #undef __FUNCT__ 1452 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC" 1453 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 1454 { 1455 FETIDPMat_ctx mat_ctx; 1456 PC_IS* pcis; 1457 PC_BDDC* pcbddc; 1458 PetscErrorCode ierr; 1459 1460 PetscFunctionBegin; 1461 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1462 pcis = (PC_IS*)mat_ctx->pc->data; 1463 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 1464 1465 /* apply B_delta^T */ 1466 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1467 ierr = VecScatterEnd (mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1468 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 1469 /* compute rhs for BDDC application */ 1470 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1471 if (pcbddc->switch_static) { 1472 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1473 } 1474 /* apply BDDC */ 1475 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 1476 /* put values into standard global vector */ 1477 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1478 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1479 if (!pcbddc->switch_static) { 1480 /* compute values into the interior if solved for the partially subassembled Schur complement */ 1481 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 1482 ierr = VecAXPY(mat_ctx->temp_solution_D,-1.0,pcis->vec1_D);CHKERRQ(ierr); 1483 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1484 } 1485 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1486 ierr = VecScatterEnd (pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1487 /* final change of basis if needed 1488 Is also sums the dirichlet part removed during RHS assembling */ 1489 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 1490 PetscFunctionReturn(0); 1491 } 1492 1493 #undef __FUNCT__ 1494 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution" 1495 /*@ 1496 PCBDDCMatFETIDPGetSolution - Compute the physical solution from the solution of the FETIDP linear system 1497 1498 Collective 1499 1500 Input Parameters: 1501 + fetidp_mat - the FETIDP matrix obtained by calling PCBDDCCreateFETIDPOperators 1502 . fetidp_flux_sol - the solution of the FETIDP linear system 1503 1504 Output Parameters: 1505 - standard_sol - the solution defined on the physical domain 1506 1507 Level: developer 1508 1509 Notes: 1510 1511 .seealso: PCBDDC,PCBDDCCreateFETIDPOperators 1512 @*/ 1513 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 1514 { 1515 FETIDPMat_ctx mat_ctx; 1516 PetscErrorCode ierr; 1517 1518 PetscFunctionBegin; 1519 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1520 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 1521 PetscFunctionReturn(0); 1522 } 1523 /* -------------------------------------------------------------------------- */ 1524 1525 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec); 1526 extern PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec); 1527 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat); 1528 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec); 1529 extern PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec); 1530 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC); 1531 1532 #undef __FUNCT__ 1533 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC" 1534 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 1535 { 1536 1537 FETIDPMat_ctx fetidpmat_ctx; 1538 Mat newmat; 1539 FETIDPPC_ctx fetidppc_ctx; 1540 PC newpc; 1541 MPI_Comm comm; 1542 PetscErrorCode ierr; 1543 1544 PetscFunctionBegin; 1545 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 1546 /* FETIDP linear matrix */ 1547 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 1548 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 1549 ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 1550 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 1551 ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr); 1552 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 1553 ierr = MatSetUp(newmat);CHKERRQ(ierr); 1554 /* FETIDP preconditioner */ 1555 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 1556 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 1557 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 1558 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 1559 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 1560 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 1561 ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 1562 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 1563 ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr); 1564 ierr = PCSetUp(newpc);CHKERRQ(ierr); 1565 /* return pointers for objects created */ 1566 *fetidp_mat=newmat; 1567 *fetidp_pc=newpc; 1568 PetscFunctionReturn(0); 1569 } 1570 1571 #undef __FUNCT__ 1572 #define __FUNCT__ "PCBDDCCreateFETIDPOperators" 1573 /*@ 1574 PCBDDCCreateFETIDPOperators - Create operators for FETIDP 1575 1576 Collective 1577 1578 Input Parameters: 1579 + pc - the BDDC preconditioning context already setup 1580 1581 Output Parameters: 1582 . fetidp_mat - shell FETIDP matrix object 1583 . fetidp_pc - shell Dirichlet preconditioner for FETIDP matrix 1584 1585 Options Database Keys: 1586 - -fetidp_fullyredundant: use or not a fully redundant set of Lagrange multipliers 1587 1588 Level: developer 1589 1590 Notes: 1591 Currently the only operation provided for FETIDP matrix is MatMult 1592 1593 .seealso: PCBDDC 1594 @*/ 1595 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 1596 { 1597 PetscErrorCode ierr; 1598 1599 PetscFunctionBegin; 1600 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1601 if (pc->setupcalled) { 1602 ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,Mat*,PC*),(pc,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 1603 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n"); 1604 PetscFunctionReturn(0); 1605 } 1606 /* -------------------------------------------------------------------------- */ 1607 /*MC 1608 PCBDDC - Balancing Domain Decomposition by Constraints. 1609 1610 An implementation of the BDDC preconditioner based on 1611 1612 .vb 1613 [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 1614 [2] A. Klawonn and O. B. Widlund. "Dual-Primal FETI Methods for Linear Elasticity", http://cs.nyu.edu/csweb/Research/TechReports/TR2004-855/TR2004-855.pdf 1615 [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977 1616 .ve 1617 1618 The matrix to be preconditioned (Pmat) must be of type MATIS. 1619 1620 Currently works with MATIS matrices with local Neumann matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers. 1621 1622 It also works with unsymmetric and indefinite problems. 1623 1624 Unlike 'conventional' interface preconditioners, PCBDDC iterates over all degrees of freedom, not just those on the interface. This allows the use of approximate solvers on the subdomains. 1625 1626 Approximate local solvers are automatically adapted for singular linear problems (see [1]) if the user has provided the nullspace using PCBDDCSetNullSpace 1627 1628 Boundary nodes are split in vertices, edges and faces using information from the local to global mapping of dofs and the local connectivity graph of nodes. The latter can be customized by using PCBDDCSetLocalAdjacencyGraph 1629 1630 Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace. 1631 1632 Change of basis is performed similarly to [2] when requested. When more the one constraint is present on a single connected component (i.e. an edge or a face), a robust method based on local QR factorizations is used. 1633 1634 The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using MatPartitioning object. 1635 1636 Options Database Keys: 1637 1638 . -pc_bddc_use_vertices <1> - use or not vertices in primal space 1639 . -pc_bddc_use_edges <1> - use or not edges in primal space 1640 . -pc_bddc_use_faces <0> - use or not faces in primal space 1641 . -pc_bddc_use_change_of_basis <0> - use change of basis approach (on edges only) 1642 . -pc_bddc_use_change_on_faces <0> - use change of basis approach on faces if change of basis has been requested 1643 . -pc_bddc_switch_static <0> - switches from M_2 to M_3 operator (see reference article [1]) 1644 . -pc_bddc_levels <0> - maximum number of levels for multilevel 1645 . -pc_bddc_coarsening_ratio - H/h ratio at the coarser level 1646 - -pc_bddc_check_level <0> - set verbosity level of debugging output 1647 1648 Options for Dirichlet, Neumann or coarse solver can be set with 1649 .vb 1650 -pc_bddc_dirichlet_ 1651 -pc_bddc_neumann_ 1652 -pc_bddc_coarse_ 1653 .ve 1654 e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg 1655 1656 When using a multilevel approach, solvers' options at the N-th level can be specified as 1657 .vb 1658 -pc_bddc_dirichlet_lN_ 1659 -pc_bddc_neumann_lN_ 1660 -pc_bddc_coarse_lN_ 1661 .ve 1662 Note that level number ranges from the finest 0 to the coarsest N. 1663 1664 Level: intermediate 1665 1666 Developer notes: 1667 1668 New deluxe scaling operator will be available soon. 1669 1670 Contributed by Stefano Zampini 1671 1672 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 1673 M*/ 1674 1675 #undef __FUNCT__ 1676 #define __FUNCT__ "PCCreate_BDDC" 1677 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 1678 { 1679 PetscErrorCode ierr; 1680 PC_BDDC *pcbddc; 1681 1682 PetscFunctionBegin; 1683 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 1684 ierr = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr); 1685 pc->data = (void*)pcbddc; 1686 1687 /* create PCIS data structure */ 1688 ierr = PCISCreate(pc);CHKERRQ(ierr); 1689 1690 /* BDDC customization */ 1691 pcbddc->use_vertices = PETSC_TRUE; 1692 pcbddc->use_edges = PETSC_TRUE; 1693 pcbddc->use_faces = PETSC_FALSE; 1694 pcbddc->use_change_of_basis = PETSC_FALSE; 1695 pcbddc->use_change_on_faces = PETSC_FALSE; 1696 pcbddc->switch_static = PETSC_FALSE; 1697 pcbddc->use_nnsp_true = PETSC_FALSE; /* not yet exposed */ 1698 pcbddc->dbg_flag = 0; 1699 1700 pcbddc->BtoNmap = 0; 1701 pcbddc->local_primal_size = 0; 1702 pcbddc->n_vertices = 0; 1703 pcbddc->n_actual_vertices = 0; 1704 pcbddc->n_constraints = 0; 1705 pcbddc->primal_indices_local_idxs = 0; 1706 pcbddc->recompute_topography = PETSC_FALSE; 1707 pcbddc->coarse_size = -1; 1708 pcbddc->new_primal_space = PETSC_FALSE; 1709 pcbddc->new_primal_space_local = PETSC_FALSE; 1710 pcbddc->global_primal_indices = 0; 1711 pcbddc->onearnullspace = 0; 1712 pcbddc->onearnullvecs_state = 0; 1713 pcbddc->user_primal_vertices = 0; 1714 pcbddc->NullSpace = 0; 1715 pcbddc->temp_solution = 0; 1716 pcbddc->original_rhs = 0; 1717 pcbddc->local_mat = 0; 1718 pcbddc->ChangeOfBasisMatrix = 0; 1719 pcbddc->coarse_vec = 0; 1720 pcbddc->coarse_rhs = 0; 1721 pcbddc->coarse_ksp = 0; 1722 pcbddc->coarse_phi_B = 0; 1723 pcbddc->coarse_phi_D = 0; 1724 pcbddc->coarse_psi_B = 0; 1725 pcbddc->coarse_psi_D = 0; 1726 pcbddc->vec1_P = 0; 1727 pcbddc->vec1_R = 0; 1728 pcbddc->vec2_R = 0; 1729 pcbddc->local_auxmat1 = 0; 1730 pcbddc->local_auxmat2 = 0; 1731 pcbddc->R_to_B = 0; 1732 pcbddc->R_to_D = 0; 1733 pcbddc->ksp_D = 0; 1734 pcbddc->ksp_R = 0; 1735 pcbddc->NeumannBoundaries = 0; 1736 pcbddc->NeumannBoundariesLocal = 0; 1737 pcbddc->DirichletBoundaries = 0; 1738 pcbddc->DirichletBoundariesLocal = 0; 1739 pcbddc->user_provided_isfordofs = PETSC_FALSE; 1740 pcbddc->n_ISForDofs = 0; 1741 pcbddc->n_ISForDofsLocal = 0; 1742 pcbddc->ISForDofs = 0; 1743 pcbddc->ISForDofsLocal = 0; 1744 pcbddc->ConstraintMatrix = 0; 1745 pcbddc->use_exact_dirichlet_trick = PETSC_TRUE; 1746 pcbddc->coarse_loc_to_glob = 0; 1747 pcbddc->coarsening_ratio = 8; 1748 pcbddc->current_level = 0; 1749 pcbddc->max_levels = 0; 1750 pcbddc->use_coarse_estimates = PETSC_FALSE; 1751 pcbddc->coarse_subassembling = 0; 1752 pcbddc->coarse_subassembling_init = 0; 1753 1754 /* create local graph structure */ 1755 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 1756 1757 /* scaling */ 1758 pcbddc->use_deluxe_scaling = PETSC_FALSE; 1759 pcbddc->work_scaling = 0; 1760 1761 /* function pointers */ 1762 pc->ops->apply = PCApply_BDDC; 1763 pc->ops->applytranspose = PCApplyTranspose_BDDC; 1764 pc->ops->setup = PCSetUp_BDDC; 1765 pc->ops->destroy = PCDestroy_BDDC; 1766 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 1767 pc->ops->view = 0; 1768 pc->ops->applyrichardson = 0; 1769 pc->ops->applysymmetricleft = 0; 1770 pc->ops->applysymmetricright = 0; 1771 pc->ops->presolve = PCPreSolve_BDDC; 1772 pc->ops->postsolve = PCPostSolve_BDDC; 1773 1774 /* composing function */ 1775 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 1776 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 1777 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr); 1778 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr); 1779 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr); 1780 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",PCBDDCSetNullSpace_BDDC);CHKERRQ(ierr); 1781 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 1782 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 1783 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 1784 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 1785 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 1786 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 1787 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 1788 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 1789 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 1790 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr); 1791 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 1792 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 1793 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 1794 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 1795 PetscFunctionReturn(0); 1796 } 1797 1798