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