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