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 && pcbddc->DirichletBoundariesLocal) { 982 IS dirIS = NULL; 983 984 /* DirichletBoundariesLocal may not be consistent among neighbours */ 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 } 1010 ierr = ISDestroy(&dirIS);CHKERRQ(ierr); 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 1204 /* Get stdout for dbg */ 1205 if (pcbddc->dbg_flag) { 1206 if (!pcbddc->dbg_viewer) { 1207 pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)); 1208 ierr = PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr); 1209 } 1210 ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1211 } 1212 1213 /*ierr = MatIsSymmetric(pc->pmat,0.,&pcbddc->issym);CHKERRQ(ierr);*/ 1214 { /* this is a temporary workaround since seqbaij matrices does not have support for symmetry checking */ 1215 PetscBool setsym; 1216 ierr = MatIsSymmetricKnown(pc->pmat,&setsym,&pcbddc->issym);CHKERRQ(ierr); 1217 if (!setsym) pcbddc->issym = PETSC_FALSE; 1218 } 1219 1220 if (pcbddc->user_ChangeOfBasisMatrix) { 1221 /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */ 1222 pcbddc->use_change_of_basis = PETSC_FALSE; 1223 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 1224 } else { 1225 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1226 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1227 pcbddc->local_mat = matis->A; 1228 } 1229 1230 /* workaround for reals */ 1231 #if !defined(PETSC_USE_COMPLEX) 1232 ierr = MatSetOption(pcbddc->local_mat,MAT_HERMITIAN,pcbddc->issym);CHKERRQ(ierr); 1233 #endif 1234 1235 /* Set up all the "iterative substructuring" common block without computing solvers */ 1236 { 1237 Mat temp_mat; 1238 1239 temp_mat = matis->A; 1240 matis->A = pcbddc->local_mat; 1241 ierr = PCISSetUp(pc,PETSC_FALSE);CHKERRQ(ierr); 1242 pcbddc->local_mat = matis->A; 1243 matis->A = temp_mat; 1244 } 1245 1246 /* Analyze interface and setup sub_schurs data */ 1247 if (computetopography) { 1248 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 1249 computeconstraintsmatrix = PETSC_TRUE; 1250 } 1251 1252 /* Setup local dirichlet solver ksp_D and sub_schurs solvers */ 1253 if (computesolvers) { 1254 PCBDDCSubSchurs sub_schurs=pcbddc->sub_schurs; 1255 1256 if (computesubschurs && computetopography) { 1257 ierr = PCBDDCInitSubSchurs(pc);CHKERRQ(ierr); 1258 } 1259 if (sub_schurs->use_mumps) { 1260 if (computesubschurs) { 1261 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1262 } 1263 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1264 } else { 1265 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1266 if (computesubschurs) { 1267 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1268 } 1269 } 1270 if (pcbddc->adaptive_selection) { 1271 ierr = PCBDDCAdaptiveSelection(pc);CHKERRQ(ierr); 1272 computeconstraintsmatrix = PETSC_TRUE; 1273 } 1274 } 1275 1276 /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */ 1277 new_nearnullspace_provided = PETSC_FALSE; 1278 ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr); 1279 if (pcbddc->onearnullspace) { /* already used nearnullspace */ 1280 if (!nearnullspace) { /* near null space attached to mat has been destroyed */ 1281 new_nearnullspace_provided = PETSC_TRUE; 1282 } else { 1283 /* determine if the two nullspaces are different (should be lightweight) */ 1284 if (nearnullspace != pcbddc->onearnullspace) { 1285 new_nearnullspace_provided = PETSC_TRUE; 1286 } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */ 1287 PetscInt i; 1288 const Vec *nearnullvecs; 1289 PetscObjectState state; 1290 PetscInt nnsp_size; 1291 ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 1292 for (i=0;i<nnsp_size;i++) { 1293 ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr); 1294 if (pcbddc->onearnullvecs_state[i] != state) { 1295 new_nearnullspace_provided = PETSC_TRUE; 1296 break; 1297 } 1298 } 1299 } 1300 } 1301 } else { 1302 if (!nearnullspace) { /* both nearnullspaces are null */ 1303 new_nearnullspace_provided = PETSC_FALSE; 1304 } else { /* nearnullspace attached later */ 1305 new_nearnullspace_provided = PETSC_TRUE; 1306 } 1307 } 1308 1309 /* Setup constraints and related work vectors */ 1310 /* reset primal space flags */ 1311 pcbddc->new_primal_space = PETSC_FALSE; 1312 pcbddc->new_primal_space_local = PETSC_FALSE; 1313 if (computeconstraintsmatrix || new_nearnullspace_provided) { 1314 /* It also sets the primal space flags */ 1315 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 1316 /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */ 1317 ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr); 1318 } 1319 1320 if (computesolvers || pcbddc->new_primal_space) { 1321 if (pcbddc->use_change_of_basis) { 1322 PC_IS *pcis = (PC_IS*)(pc->data); 1323 1324 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 1325 /* get submatrices */ 1326 ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr); 1327 ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr); 1328 ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr); 1329 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);CHKERRQ(ierr); 1330 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);CHKERRQ(ierr); 1331 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);CHKERRQ(ierr); 1332 /* set flag in pcis to not reuse submatrices during PCISCreate */ 1333 pcis->reusesubmatrices = PETSC_FALSE; 1334 } else if (!pcbddc->user_ChangeOfBasisMatrix) { 1335 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1336 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1337 pcbddc->local_mat = matis->A; 1338 } 1339 /* SetUp coarse and local Neumann solvers */ 1340 ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr); 1341 /* SetUp Scaling operator */ 1342 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1343 } 1344 1345 if (pcbddc->dbg_flag) { 1346 ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1347 } 1348 PetscFunctionReturn(0); 1349 } 1350 1351 /* -------------------------------------------------------------------------- */ 1352 /* 1353 PCApply_BDDC - Applies the BDDC operator to a vector. 1354 1355 Input Parameters: 1356 . pc - the preconditioner context 1357 . r - input vector (global) 1358 1359 Output Parameter: 1360 . z - output vector (global) 1361 1362 Application Interface Routine: PCApply() 1363 */ 1364 #undef __FUNCT__ 1365 #define __FUNCT__ "PCApply_BDDC" 1366 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 1367 { 1368 PC_IS *pcis = (PC_IS*)(pc->data); 1369 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1370 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1371 PetscErrorCode ierr; 1372 const PetscScalar one = 1.0; 1373 const PetscScalar m_one = -1.0; 1374 const PetscScalar zero = 0.0; 1375 1376 /* This code is similar to that provided in nn.c for PCNN 1377 NN interface preconditioner changed to BDDC 1378 Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */ 1379 1380 PetscFunctionBegin; 1381 if (!pcbddc->use_exact_dirichlet_trick) { 1382 ierr = VecCopy(r,z);CHKERRQ(ierr); 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 /* 1387 Assembling right hand side for BDDC operator 1388 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1389 - pcis->vec1_B the interface part of the global vector z 1390 */ 1391 if (n_D) { 1392 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1393 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1394 if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 1395 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1396 } else { 1397 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1398 } 1399 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1400 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1401 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1402 } else { 1403 if (pcbddc->switch_static) { 1404 ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr); 1405 } 1406 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1407 } 1408 1409 /* Apply interface preconditioner 1410 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1411 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr); 1412 1413 /* Apply transpose of partition of unity operator */ 1414 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1415 1416 /* Second Dirichlet solve and assembling of output */ 1417 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1418 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1419 if (n_B) { 1420 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1421 if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 1422 } else if (pcbddc->switch_static) { 1423 ierr = MatMult(pcis->A_II,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 1424 } 1425 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1426 if (!pcbddc->use_exact_dirichlet_trick) { 1427 if (pcbddc->switch_static) { 1428 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 1429 } else { 1430 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 1431 } 1432 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1433 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1434 } else { 1435 if (pcbddc->switch_static) { 1436 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 1437 } else { 1438 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1439 } 1440 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1441 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1442 } 1443 PetscFunctionReturn(0); 1444 } 1445 1446 /* -------------------------------------------------------------------------- */ 1447 /* 1448 PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector. 1449 1450 Input Parameters: 1451 . pc - the preconditioner context 1452 . r - input vector (global) 1453 1454 Output Parameter: 1455 . z - output vector (global) 1456 1457 Application Interface Routine: PCApplyTranspose() 1458 */ 1459 #undef __FUNCT__ 1460 #define __FUNCT__ "PCApplyTranspose_BDDC" 1461 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z) 1462 { 1463 PC_IS *pcis = (PC_IS*)(pc->data); 1464 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1465 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1466 PetscErrorCode ierr; 1467 const PetscScalar one = 1.0; 1468 const PetscScalar m_one = -1.0; 1469 const PetscScalar zero = 0.0; 1470 1471 PetscFunctionBegin; 1472 if (!pcbddc->use_exact_dirichlet_trick) { 1473 ierr = VecCopy(r,z);CHKERRQ(ierr); 1474 /* First Dirichlet solve */ 1475 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1476 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1477 /* 1478 Assembling right hand side for BDDC operator 1479 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1480 - pcis->vec1_B the interface part of the global vector z 1481 */ 1482 if (n_D) { 1483 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1484 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1485 if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 1486 ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1487 } else { 1488 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1489 } 1490 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1491 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1492 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1493 } else { 1494 if (pcbddc->switch_static) { 1495 ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr); 1496 } 1497 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1498 } 1499 1500 /* Apply interface preconditioner 1501 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1502 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr); 1503 1504 /* Apply transpose of partition of unity operator */ 1505 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1506 1507 /* Second Dirichlet solve and assembling of output */ 1508 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1509 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1510 if (n_B) { 1511 ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1512 if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 1513 } else if (pcbddc->switch_static) { 1514 ierr = MatMultTranspose(pcis->A_II,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 1515 } 1516 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1517 if (!pcbddc->use_exact_dirichlet_trick) { 1518 if (pcbddc->switch_static) { 1519 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 1520 } else { 1521 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 1522 } 1523 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1524 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1525 } else { 1526 if (pcbddc->switch_static) { 1527 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 1528 } else { 1529 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1530 } 1531 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1532 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1533 } 1534 PetscFunctionReturn(0); 1535 } 1536 /* -------------------------------------------------------------------------- */ 1537 1538 #undef __FUNCT__ 1539 #define __FUNCT__ "PCDestroy_BDDC" 1540 PetscErrorCode PCDestroy_BDDC(PC pc) 1541 { 1542 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1543 PetscErrorCode ierr; 1544 1545 PetscFunctionBegin; 1546 /* free data created by PCIS */ 1547 ierr = PCISDestroy(pc);CHKERRQ(ierr); 1548 /* free BDDC custom data */ 1549 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 1550 /* destroy objects related to topography */ 1551 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 1552 /* free allocated graph structure */ 1553 ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr); 1554 /* free allocated sub schurs structure */ 1555 ierr = PetscFree(pcbddc->sub_schurs);CHKERRQ(ierr); 1556 /* destroy objects for scaling operator */ 1557 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 1558 ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr); 1559 /* free solvers stuff */ 1560 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 1561 /* free global vectors needed in presolve */ 1562 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 1563 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 1564 /* free stuff for change of basis hooks */ 1565 if (pcbddc->new_global_mat) { 1566 PCBDDCChange_ctx change_ctx; 1567 ierr = MatShellGetContext(pcbddc->new_global_mat,&change_ctx);CHKERRQ(ierr); 1568 ierr = MatDestroy(&change_ctx->original_mat);CHKERRQ(ierr); 1569 ierr = MatDestroy(&change_ctx->global_change);CHKERRQ(ierr); 1570 ierr = VecDestroyVecs(2,&change_ctx->work);CHKERRQ(ierr); 1571 ierr = PetscFree(change_ctx);CHKERRQ(ierr); 1572 } 1573 ierr = MatDestroy(&pcbddc->new_global_mat);CHKERRQ(ierr); 1574 /* remove functions */ 1575 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);CHKERRQ(ierr); 1576 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 1577 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 1578 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr); 1579 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr); 1580 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr); 1581 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",NULL);CHKERRQ(ierr); 1582 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 1583 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 1584 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 1585 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 1586 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 1587 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 1588 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 1589 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 1590 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 1591 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr); 1592 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 1593 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 1594 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 1595 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 1596 /* Free the private data structure */ 1597 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1598 PetscFunctionReturn(0); 1599 } 1600 /* -------------------------------------------------------------------------- */ 1601 1602 #undef __FUNCT__ 1603 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC" 1604 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 1605 { 1606 FETIDPMat_ctx mat_ctx; 1607 Vec copy_standard_rhs; 1608 PC_IS* pcis; 1609 PC_BDDC* pcbddc; 1610 PetscErrorCode ierr; 1611 1612 PetscFunctionBegin; 1613 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1614 pcis = (PC_IS*)mat_ctx->pc->data; 1615 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 1616 1617 /* 1618 change of basis for physical rhs if needed 1619 It also changes the rhs in case of dirichlet boundaries 1620 TODO: better management when FETIDP will have its own class 1621 */ 1622 ierr = VecDuplicate(standard_rhs,©_standard_rhs);CHKERRQ(ierr); 1623 ierr = VecCopy(standard_rhs,copy_standard_rhs);CHKERRQ(ierr); 1624 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,copy_standard_rhs,NULL);CHKERRQ(ierr); 1625 /* store vectors for computation of fetidp final solution */ 1626 ierr = VecScatterBegin(pcis->global_to_D,copy_standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1627 ierr = VecScatterEnd(pcis->global_to_D,copy_standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1628 /* scale rhs since it should be unassembled */ 1629 /* TODO use counter scaling? (also below) */ 1630 ierr = VecScatterBegin(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1631 ierr = VecScatterEnd(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1632 /* Apply partition of unity */ 1633 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1634 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,copy_standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 1635 if (!pcbddc->switch_static) { 1636 /* compute partially subassembled Schur complement right-hand side */ 1637 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1638 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 1639 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 1640 ierr = VecSet(copy_standard_rhs,0.0);CHKERRQ(ierr); 1641 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,copy_standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1642 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,copy_standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1643 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,copy_standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 1644 ierr = VecScatterBegin(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1645 ierr = VecScatterEnd(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1646 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1647 } 1648 ierr = VecDestroy(©_standard_rhs);CHKERRQ(ierr); 1649 /* BDDC rhs */ 1650 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 1651 if (pcbddc->switch_static) { 1652 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1653 } 1654 /* apply BDDC */ 1655 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 1656 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 1657 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 1658 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 1659 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1660 ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1661 PetscFunctionReturn(0); 1662 } 1663 1664 #undef __FUNCT__ 1665 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS" 1666 /*@ 1667 PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETIDP linear system 1668 1669 Collective 1670 1671 Input Parameters: 1672 + fetidp_mat - the FETIDP matrix object obtained by calling PCBDDCCreateFETIDPOperators 1673 . standard_rhs - the right-hand side for your linear system 1674 1675 Output Parameters: 1676 - fetidp_flux_rhs - the right-hand side for the FETIDP linear system 1677 1678 Level: developer 1679 1680 Notes: 1681 1682 .seealso: PCBDDC,PCBDDCCreateFETIDPOperators 1683 @*/ 1684 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 1685 { 1686 FETIDPMat_ctx mat_ctx; 1687 PetscErrorCode ierr; 1688 1689 PetscFunctionBegin; 1690 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1691 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 1692 PetscFunctionReturn(0); 1693 } 1694 /* -------------------------------------------------------------------------- */ 1695 1696 #undef __FUNCT__ 1697 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC" 1698 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 1699 { 1700 FETIDPMat_ctx mat_ctx; 1701 PC_IS* pcis; 1702 PC_BDDC* pcbddc; 1703 PetscErrorCode ierr; 1704 1705 PetscFunctionBegin; 1706 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1707 pcis = (PC_IS*)mat_ctx->pc->data; 1708 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 1709 1710 /* apply B_delta^T */ 1711 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1712 ierr = VecScatterEnd (mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1713 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 1714 /* compute rhs for BDDC application */ 1715 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1716 if (pcbddc->switch_static) { 1717 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1718 } 1719 /* apply BDDC */ 1720 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 1721 /* put values into standard global vector */ 1722 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1723 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1724 if (!pcbddc->switch_static) { 1725 /* compute values into the interior if solved for the partially subassembled Schur complement */ 1726 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 1727 ierr = VecAXPY(mat_ctx->temp_solution_D,-1.0,pcis->vec1_D);CHKERRQ(ierr); 1728 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1729 } 1730 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1731 ierr = VecScatterEnd (pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1732 /* final change of basis if needed 1733 Is also sums the dirichlet part removed during RHS assembling */ 1734 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 1735 PetscFunctionReturn(0); 1736 } 1737 1738 #undef __FUNCT__ 1739 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution" 1740 /*@ 1741 PCBDDCMatFETIDPGetSolution - Compute the physical solution from the solution of the FETIDP linear system 1742 1743 Collective 1744 1745 Input Parameters: 1746 + fetidp_mat - the FETIDP matrix obtained by calling PCBDDCCreateFETIDPOperators 1747 . fetidp_flux_sol - the solution of the FETIDP linear system 1748 1749 Output Parameters: 1750 - standard_sol - the solution defined on the physical domain 1751 1752 Level: developer 1753 1754 Notes: 1755 1756 .seealso: PCBDDC,PCBDDCCreateFETIDPOperators 1757 @*/ 1758 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 1759 { 1760 FETIDPMat_ctx mat_ctx; 1761 PetscErrorCode ierr; 1762 1763 PetscFunctionBegin; 1764 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1765 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 1766 PetscFunctionReturn(0); 1767 } 1768 /* -------------------------------------------------------------------------- */ 1769 1770 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec); 1771 extern PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec); 1772 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat); 1773 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec); 1774 extern PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec); 1775 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC); 1776 1777 #undef __FUNCT__ 1778 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC" 1779 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 1780 { 1781 1782 FETIDPMat_ctx fetidpmat_ctx; 1783 Mat newmat; 1784 FETIDPPC_ctx fetidppc_ctx; 1785 PC newpc; 1786 MPI_Comm comm; 1787 PetscErrorCode ierr; 1788 1789 PetscFunctionBegin; 1790 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 1791 /* FETIDP linear matrix */ 1792 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 1793 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 1794 ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 1795 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 1796 ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr); 1797 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 1798 ierr = MatSetUp(newmat);CHKERRQ(ierr); 1799 /* FETIDP preconditioner */ 1800 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 1801 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 1802 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 1803 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 1804 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 1805 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 1806 ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 1807 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 1808 ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr); 1809 ierr = PCSetUp(newpc);CHKERRQ(ierr); 1810 /* return pointers for objects created */ 1811 *fetidp_mat=newmat; 1812 *fetidp_pc=newpc; 1813 PetscFunctionReturn(0); 1814 } 1815 1816 #undef __FUNCT__ 1817 #define __FUNCT__ "PCBDDCCreateFETIDPOperators" 1818 /*@ 1819 PCBDDCCreateFETIDPOperators - Create operators for FETIDP 1820 1821 Collective 1822 1823 Input Parameters: 1824 + pc - the BDDC preconditioning context already setup 1825 1826 Output Parameters: 1827 . fetidp_mat - shell FETIDP matrix object 1828 . fetidp_pc - shell Dirichlet preconditioner for FETIDP matrix 1829 1830 Options Database Keys: 1831 - -fetidp_fullyredundant: use or not a fully redundant set of Lagrange multipliers 1832 1833 Level: developer 1834 1835 Notes: 1836 Currently the only operation provided for FETIDP matrix is MatMult 1837 1838 .seealso: PCBDDC 1839 @*/ 1840 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 1841 { 1842 PetscErrorCode ierr; 1843 1844 PetscFunctionBegin; 1845 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1846 if (pc->setupcalled) { 1847 ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,Mat*,PC*),(pc,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 1848 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n"); 1849 PetscFunctionReturn(0); 1850 } 1851 /* -------------------------------------------------------------------------- */ 1852 /*MC 1853 PCBDDC - Balancing Domain Decomposition by Constraints. 1854 1855 An implementation of the BDDC preconditioner based on 1856 1857 .vb 1858 [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 1859 [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 1860 [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977 1861 .ve 1862 1863 The matrix to be preconditioned (Pmat) must be of type MATIS. 1864 1865 Currently works with MATIS matrices with local Neumann matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers. 1866 1867 It also works with unsymmetric and indefinite problems. 1868 1869 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. 1870 1871 Approximate local solvers are automatically adapted for singular linear problems (see [1]) if the user has provided the nullspace using PCBDDCSetNullSpace 1872 1873 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() 1874 1875 Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). 1876 1877 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. 1878 1879 The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using MatPartitioning object. 1880 1881 Options Database Keys: 1882 1883 . -pc_bddc_use_vertices <1> - use or not vertices in primal space 1884 . -pc_bddc_use_edges <1> - use or not edges in primal space 1885 . -pc_bddc_use_faces <0> - use or not faces in primal space 1886 . -pc_bddc_use_change_of_basis <0> - use change of basis approach (on edges only) 1887 . -pc_bddc_use_change_on_faces <0> - use change of basis approach on faces if change of basis has been requested 1888 . -pc_bddc_switch_static <0> - switches from M_2 to M_3 operator (see reference article [1]) 1889 . -pc_bddc_levels <0> - maximum number of levels for multilevel 1890 . -pc_bddc_coarsening_ratio - H/h ratio at the coarser level 1891 - -pc_bddc_check_level <0> - set verbosity level of debugging output 1892 1893 Options for Dirichlet, Neumann or coarse solver can be set with 1894 .vb 1895 -pc_bddc_dirichlet_ 1896 -pc_bddc_neumann_ 1897 -pc_bddc_coarse_ 1898 .ve 1899 e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg 1900 1901 When using a multilevel approach, solvers' options at the N-th level can be specified as 1902 .vb 1903 -pc_bddc_dirichlet_lN_ 1904 -pc_bddc_neumann_lN_ 1905 -pc_bddc_coarse_lN_ 1906 .ve 1907 Note that level number ranges from the finest 0 to the coarsest N. 1908 1909 Level: intermediate 1910 1911 Developer notes: 1912 1913 New deluxe scaling operator will be available soon. 1914 1915 Contributed by Stefano Zampini 1916 1917 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 1918 M*/ 1919 1920 #undef __FUNCT__ 1921 #define __FUNCT__ "PCCreate_BDDC" 1922 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 1923 { 1924 PetscErrorCode ierr; 1925 PC_BDDC *pcbddc; 1926 1927 PetscFunctionBegin; 1928 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 1929 ierr = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr); 1930 pc->data = (void*)pcbddc; 1931 1932 /* create PCIS data structure */ 1933 ierr = PCISCreate(pc);CHKERRQ(ierr); 1934 1935 /* BDDC customization */ 1936 pcbddc->use_local_adj = PETSC_TRUE; 1937 pcbddc->use_vertices = PETSC_TRUE; 1938 pcbddc->use_edges = PETSC_TRUE; 1939 pcbddc->use_faces = PETSC_FALSE; 1940 pcbddc->use_change_of_basis = PETSC_FALSE; 1941 pcbddc->use_change_on_faces = PETSC_FALSE; 1942 pcbddc->switch_static = PETSC_FALSE; 1943 pcbddc->use_nnsp_true = PETSC_FALSE; 1944 pcbddc->use_qr_single = PETSC_FALSE; 1945 pcbddc->dbg_flag = 0; 1946 /* private */ 1947 pcbddc->issym = PETSC_FALSE; 1948 pcbddc->local_primal_size = 0; 1949 pcbddc->local_primal_size_cc = 0; 1950 pcbddc->local_primal_ref_node = 0; 1951 pcbddc->local_primal_ref_mult = 0; 1952 pcbddc->n_vertices = 0; 1953 pcbddc->primal_indices_local_idxs = 0; 1954 pcbddc->recompute_topography = PETSC_FALSE; 1955 pcbddc->coarse_size = -1; 1956 pcbddc->new_primal_space = PETSC_FALSE; 1957 pcbddc->new_primal_space_local = PETSC_FALSE; 1958 pcbddc->global_primal_indices = 0; 1959 pcbddc->onearnullspace = 0; 1960 pcbddc->onearnullvecs_state = 0; 1961 pcbddc->user_primal_vertices = 0; 1962 pcbddc->NullSpace = 0; 1963 pcbddc->temp_solution = 0; 1964 pcbddc->original_rhs = 0; 1965 pcbddc->local_mat = 0; 1966 pcbddc->ChangeOfBasisMatrix = 0; 1967 pcbddc->user_ChangeOfBasisMatrix = 0; 1968 pcbddc->new_global_mat = 0; 1969 pcbddc->coarse_vec = 0; 1970 pcbddc->coarse_ksp = 0; 1971 pcbddc->coarse_phi_B = 0; 1972 pcbddc->coarse_phi_D = 0; 1973 pcbddc->coarse_psi_B = 0; 1974 pcbddc->coarse_psi_D = 0; 1975 pcbddc->vec1_P = 0; 1976 pcbddc->vec1_R = 0; 1977 pcbddc->vec2_R = 0; 1978 pcbddc->local_auxmat1 = 0; 1979 pcbddc->local_auxmat2 = 0; 1980 pcbddc->R_to_B = 0; 1981 pcbddc->R_to_D = 0; 1982 pcbddc->ksp_D = 0; 1983 pcbddc->ksp_R = 0; 1984 pcbddc->NeumannBoundaries = 0; 1985 pcbddc->NeumannBoundariesLocal = 0; 1986 pcbddc->DirichletBoundaries = 0; 1987 pcbddc->DirichletBoundariesLocal = 0; 1988 pcbddc->user_provided_isfordofs = PETSC_FALSE; 1989 pcbddc->n_ISForDofs = 0; 1990 pcbddc->n_ISForDofsLocal = 0; 1991 pcbddc->ISForDofs = 0; 1992 pcbddc->ISForDofsLocal = 0; 1993 pcbddc->ConstraintMatrix = 0; 1994 pcbddc->use_exact_dirichlet_trick = PETSC_TRUE; 1995 pcbddc->coarse_loc_to_glob = 0; 1996 pcbddc->coarsening_ratio = 8; 1997 pcbddc->current_level = 0; 1998 pcbddc->max_levels = 0; 1999 pcbddc->use_coarse_estimates = PETSC_FALSE; 2000 pcbddc->redistribute_coarse = 0; 2001 pcbddc->coarse_subassembling = 0; 2002 pcbddc->coarse_subassembling_init = 0; 2003 2004 /* create local graph structure */ 2005 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 2006 2007 /* scaling */ 2008 pcbddc->work_scaling = 0; 2009 pcbddc->use_deluxe_scaling = PETSC_FALSE; 2010 pcbddc->faster_deluxe = PETSC_FALSE; 2011 2012 /* create sub schurs structure */ 2013 ierr = PCBDDCSubSchursCreate(&pcbddc->sub_schurs);CHKERRQ(ierr); 2014 pcbddc->sub_schurs_rebuild = PETSC_FALSE; 2015 pcbddc->sub_schurs_layers = -1; 2016 pcbddc->sub_schurs_use_useradj = PETSC_FALSE; 2017 2018 pcbddc->computed_rowadj = PETSC_FALSE; 2019 2020 /* adaptivity */ 2021 pcbddc->adaptive_threshold = 0.0; 2022 pcbddc->adaptive_nmax = 0; 2023 pcbddc->adaptive_nmin = 0; 2024 2025 /* function pointers */ 2026 pc->ops->apply = PCApply_BDDC; 2027 pc->ops->applytranspose = PCApplyTranspose_BDDC; 2028 pc->ops->setup = PCSetUp_BDDC; 2029 pc->ops->destroy = PCDestroy_BDDC; 2030 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 2031 pc->ops->view = 0; 2032 pc->ops->applyrichardson = 0; 2033 pc->ops->applysymmetricleft = 0; 2034 pc->ops->applysymmetricright = 0; 2035 pc->ops->presolve = PCPreSolve_BDDC; 2036 pc->ops->postsolve = PCPostSolve_BDDC; 2037 2038 /* composing function */ 2039 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);CHKERRQ(ierr); 2040 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 2041 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 2042 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr); 2043 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr); 2044 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr); 2045 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",PCBDDCSetNullSpace_BDDC);CHKERRQ(ierr); 2046 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2047 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2048 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2049 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2050 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2051 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2052 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2053 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2054 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 2055 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr); 2056 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 2057 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 2058 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 2059 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 2060 PetscFunctionReturn(0); 2061 } 2062 2063