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 if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE; 1204 1205 computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling); 1206 if (pcbddc->faster_deluxe && pcbddc->adaptive_selection && pcbddc->use_change_of_basis) { 1207 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"); 1208 } 1209 /* Get stdout for dbg */ 1210 if (pcbddc->dbg_flag) { 1211 if (!pcbddc->dbg_viewer) { 1212 pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)); 1213 ierr = PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr); 1214 } 1215 ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1216 } 1217 1218 if (pcbddc->user_ChangeOfBasisMatrix) { 1219 /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */ 1220 pcbddc->use_change_of_basis = PETSC_FALSE; 1221 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 1222 } else { 1223 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1224 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1225 pcbddc->local_mat = matis->A; 1226 } 1227 1228 /* workaround for reals */ 1229 #if !defined(PETSC_USE_COMPLEX) 1230 if (matis->A->symmetric_set) { 1231 ierr = MatSetOption(pcbddc->local_mat,MAT_HERMITIAN,matis->A->symmetric);CHKERRQ(ierr); 1232 } 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->symmetric_primal = PETSC_TRUE; 1946 pcbddc->dbg_flag = 0; 1947 /* private */ 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