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