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