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