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