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