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