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