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 } 1218 } 1219 1220 /* change rhs and iteration matrix if using the change of basis */ 1221 if (pcbddc->ChangeOfBasisMatrix) { 1222 PCBDDCChange_ctx change_ctx; 1223 1224 /* get change ctx */ 1225 ierr = MatShellGetContext(pcbddc->new_global_mat,&change_ctx);CHKERRQ(ierr); 1226 1227 /* set current iteration matrix inside change context (change of basis has been already set into the ctx during PCSetUp) */ 1228 if (!pcbddc->benign_original_mat) { 1229 ierr = MatDestroy(&change_ctx->original_mat);CHKERRQ(ierr); 1230 ierr = PetscObjectReference((PetscObject)pc->mat);CHKERRQ(ierr); 1231 change_ctx->original_mat = pc->mat; 1232 1233 /* change iteration matrix */ 1234 ierr = MatDestroy(&pc->mat);CHKERRQ(ierr); 1235 ierr = PetscObjectReference((PetscObject)pcbddc->new_global_mat);CHKERRQ(ierr); 1236 pc->mat = pcbddc->new_global_mat; 1237 } 1238 /* store the original rhs */ 1239 if (copy_rhs) { 1240 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1241 copy_rhs = PETSC_FALSE; 1242 } 1243 1244 /* change rhs */ 1245 ierr = MatMultTranspose(change_ctx->global_change,rhs,pcis->vec1_global);CHKERRQ(ierr); 1246 ierr = VecCopy(pcis->vec1_global,rhs);CHKERRQ(ierr); 1247 pcbddc->rhs_change = PETSC_TRUE; 1248 } 1249 1250 /* remove non-benign solution from the rhs */ 1251 if (pcbddc->benign_saddle_point) { 1252 /* store the original rhs */ 1253 if (copy_rhs) { 1254 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1255 copy_rhs = PETSC_FALSE; 1256 } 1257 if (benign_correction_is_zero) { /* still need to understand why it works great */ 1258 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 1259 ierr = PCApply_BDDC(pc,rhs,pcbddc->benign_vec);CHKERRQ(ierr); 1260 } 1261 ierr = VecScale(pcbddc->benign_vec,-1.0);CHKERRQ(ierr); 1262 ierr = MatMultAdd(pc->mat,pcbddc->benign_vec,rhs,rhs);CHKERRQ(ierr); 1263 pcbddc->rhs_change = PETSC_TRUE; 1264 } 1265 1266 /* set initial guess if using PCG */ 1267 if (x && pcbddc->use_exact_dirichlet_trick) { 1268 ierr = VecSet(x,0.0);CHKERRQ(ierr); 1269 ierr = VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1270 ierr = VecScatterEnd(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1271 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1272 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1273 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1274 if (ksp) { 1275 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 1276 } 1277 } 1278 if (pcbddc->benign_null) { 1279 MatNullSpace null_space; 1280 Vec nullv; 1281 PetscBool isnull; 1282 PetscInt i; 1283 1284 for (i=0;i<pcbddc->benign_n;i++) pcbddc->benign_p0[i] = 1.; 1285 ierr = VecDuplicate(pcis->vec1_global,&nullv);CHKERRQ(ierr); 1286 ierr = VecSet(nullv,0.);CHKERRQ(ierr); 1287 ierr = PCBDDCBenignGetOrSetP0(pc,nullv,PETSC_FALSE);CHKERRQ(ierr); 1288 ierr = VecNormalize(nullv,NULL);CHKERRQ(ierr); 1289 ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)pc),PETSC_FALSE,1,&nullv,&null_space);CHKERRQ(ierr); 1290 ierr = MatNullSpaceTest(null_space,pc->mat,&isnull);CHKERRQ(ierr); 1291 if (isnull) { 1292 ierr = MatSetNullSpace(pc->mat,null_space);CHKERRQ(ierr); 1293 } 1294 ierr = MatNullSpaceDestroy(&null_space);CHKERRQ(ierr); 1295 ierr = VecDestroy(&nullv);CHKERRQ(ierr); 1296 } 1297 1298 /* remove nullspace if present */ 1299 if (ksp && x && pcbddc->NullSpace) { 1300 ierr = MatNullSpaceRemove(pcbddc->NullSpace,x);CHKERRQ(ierr); 1301 /* store the original rhs */ 1302 if (copy_rhs) { 1303 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1304 copy_rhs = PETSC_FALSE; 1305 } 1306 pcbddc->rhs_change = PETSC_TRUE; 1307 ierr = MatNullSpaceRemove(pcbddc->NullSpace,rhs);CHKERRQ(ierr); 1308 } 1309 PetscFunctionReturn(0); 1310 } 1311 1312 /* -------------------------------------------------------------------------- */ 1313 #undef __FUNCT__ 1314 #define __FUNCT__ "PCPostSolve_BDDC" 1315 /* -------------------------------------------------------------------------- */ 1316 /* 1317 PCPostSolve_BDDC - Changes the computed solution if a transformation of basis 1318 approach has been selected. Also, restores rhs to its original state. 1319 1320 Input Parameter: 1321 + pc - the preconditioner contex 1322 1323 Application Interface Routine: PCPostSolve() 1324 1325 Notes: 1326 The interface routine PCPostSolve() is not usually called directly by 1327 the user, but instead is called by KSPSolve(). 1328 */ 1329 static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 1330 { 1331 PetscErrorCode ierr; 1332 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1333 1334 PetscFunctionBegin; 1335 if (pcbddc->ChangeOfBasisMatrix && !pcbddc->benign_original_mat) { 1336 PCBDDCChange_ctx change_ctx; 1337 1338 /* get change ctx */ 1339 ierr = MatShellGetContext(pcbddc->new_global_mat,&change_ctx);CHKERRQ(ierr); 1340 1341 /* restore iteration matrix */ 1342 ierr = MatDestroy(&pc->mat);CHKERRQ(ierr); 1343 ierr = PetscObjectReference((PetscObject)change_ctx->original_mat);CHKERRQ(ierr); 1344 pc->mat = change_ctx->original_mat; 1345 } 1346 1347 /* need to restore the local matrices */ 1348 if (pcbddc->benign_original_mat) { 1349 Mat_IS *matis = (Mat_IS*)pc->mat->data; 1350 1351 if (pcbddc->benign_n && !pcbddc->benign_change_explicit) { 1352 ierr = PCBDDCBenignShellMat(pc,PETSC_TRUE);CHKERRQ(ierr); 1353 } else { 1354 pcbddc->local_mat = matis->A; 1355 matis->A = pcbddc->benign_original_mat; 1356 pcbddc->benign_original_mat = NULL; 1357 } 1358 } 1359 1360 /* get solution in original basis */ 1361 if (x) { 1362 PC_IS *pcis = (PC_IS*)(pc->data); 1363 1364 /* restore solution on pressures */ 1365 if (pcbddc->benign_saddle_point) { 1366 Mat_IS *matis = (Mat_IS*)pc->mat->data; 1367 1368 /* add non-benign solution */ 1369 ierr = VecAXPY(x,-1.0,pcbddc->benign_vec);CHKERRQ(ierr); 1370 1371 /* change basis on pressures for x */ 1372 ierr = VecScatterBegin(matis->rctx,x,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1373 ierr = VecScatterEnd(matis->rctx,x,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1374 if (pcbddc->benign_change) { 1375 ierr = MatMult(pcbddc->benign_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1376 ierr = VecScatterBegin(matis->rctx,pcis->vec2_N,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1377 ierr = VecScatterEnd(matis->rctx,pcis->vec2_N,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1378 } else { 1379 ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1380 ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1381 } 1382 } 1383 1384 /* change basis on x */ 1385 if (pcbddc->ChangeOfBasisMatrix) { 1386 PCBDDCChange_ctx change_ctx; 1387 1388 ierr = MatShellGetContext(pcbddc->new_global_mat,&change_ctx);CHKERRQ(ierr); 1389 ierr = MatMult(change_ctx->global_change,x,pcis->vec1_global);CHKERRQ(ierr); 1390 ierr = VecCopy(pcis->vec1_global,x);CHKERRQ(ierr); 1391 } 1392 } 1393 1394 /* add solution removed in presolve */ 1395 if (x && pcbddc->rhs_change) { 1396 ierr = VecAXPY(x,1.0,pcbddc->temp_solution);CHKERRQ(ierr); 1397 } 1398 1399 /* restore rhs to its original state */ 1400 if (rhs && pcbddc->rhs_change) { 1401 ierr = VecCopy(pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1402 } 1403 pcbddc->rhs_change = PETSC_FALSE; 1404 1405 /* restore ksp guess state */ 1406 if (ksp) { 1407 ierr = KSPSetInitialGuessNonzero(ksp,pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 1408 } 1409 PetscFunctionReturn(0); 1410 } 1411 /* -------------------------------------------------------------------------- */ 1412 #undef __FUNCT__ 1413 #define __FUNCT__ "PCSetUp_BDDC" 1414 /* -------------------------------------------------------------------------- */ 1415 /* 1416 PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner 1417 by setting data structures and options. 1418 1419 Input Parameter: 1420 + pc - the preconditioner context 1421 1422 Application Interface Routine: PCSetUp() 1423 1424 Notes: 1425 The interface routine PCSetUp() is not usually called directly by 1426 the user, but instead is called by PCApply() if necessary. 1427 */ 1428 PetscErrorCode PCSetUp_BDDC(PC pc) 1429 { 1430 PetscErrorCode ierr; 1431 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 1432 Mat_IS* matis; 1433 MatNullSpace nearnullspace; 1434 IS zerodiag = NULL; 1435 PetscInt nrows,ncols; 1436 PetscBool computetopography,computesolvers,computesubschurs; 1437 PetscBool computeconstraintsmatrix; 1438 PetscBool new_nearnullspace_provided,ismatis; 1439 1440 PetscFunctionBegin; 1441 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&ismatis);CHKERRQ(ierr); 1442 if (!ismatis) { 1443 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PCBDDC preconditioner requires matrix of type MATIS"); 1444 } 1445 ierr = MatGetSize(pc->pmat,&nrows,&ncols);CHKERRQ(ierr); 1446 if (nrows != ncols) { 1447 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCBDDC preconditioner requires a square preconditioning matrix"); 1448 } 1449 matis = (Mat_IS*)pc->pmat->data; 1450 /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */ 1451 /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup 1452 Also, BDDC directly build the Dirichlet problem */ 1453 /* split work */ 1454 if (pc->setupcalled) { 1455 if (pc->flag == SAME_NONZERO_PATTERN) { 1456 computetopography = PETSC_FALSE; 1457 computesolvers = PETSC_TRUE; 1458 } else { /* DIFFERENT_NONZERO_PATTERN */ 1459 computetopography = PETSC_TRUE; 1460 computesolvers = PETSC_TRUE; 1461 } 1462 } else { 1463 computetopography = PETSC_TRUE; 1464 computesolvers = PETSC_TRUE; 1465 } 1466 if (pcbddc->recompute_topography) { 1467 computetopography = PETSC_TRUE; 1468 } 1469 pcbddc->recompute_topography = computetopography; 1470 computeconstraintsmatrix = PETSC_FALSE; 1471 1472 /* check parameters' compatibility */ 1473 if (pcbddc->adaptive_threshold > 0.0 && !pcbddc->use_deluxe_scaling) { 1474 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot compute adaptive constraints without deluxe scaling. Rerun with -pc_bddc_use_deluxe_scaling"); 1475 } 1476 pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_threshold > 0.0 && pcbddc->use_deluxe_scaling); 1477 pcbddc->adaptive_userdefined = (PetscBool)(pcbddc->adaptive_selection && pcbddc->adaptive_userdefined); 1478 if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE; 1479 1480 computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling); 1481 if (pcbddc->faster_deluxe && pcbddc->adaptive_selection && pcbddc->use_change_of_basis) { 1482 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"); 1483 } 1484 if (pcbddc->benign_saddle_point && pcbddc->switch_static) { 1485 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Benign subspace trick cannot be used with M_3"); 1486 } 1487 1488 /* check if the iteration matrix is of type MATIS in case the benign trick has been requested */ 1489 ierr = PetscObjectTypeCompare((PetscObject)pc->mat,MATIS,&ismatis);CHKERRQ(ierr); 1490 if (pcbddc->benign_saddle_point && !ismatis) { 1491 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PCBDDC preconditioner with benign subspace trick requires the iteration matrix to be of type MATIS"); 1492 } 1493 if (pcbddc->benign_saddle_point && pc->mat != pc->pmat) { 1494 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PCBDDC preconditioner with benign subspace trick requires Amat == Pmat"); 1495 } 1496 1497 /* Get stdout for dbg */ 1498 if (pcbddc->dbg_flag) { 1499 if (!pcbddc->dbg_viewer) { 1500 pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)); 1501 ierr = PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer);CHKERRQ(ierr); 1502 } 1503 ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1504 } 1505 1506 if (pcbddc->user_ChangeOfBasisMatrix) { 1507 /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */ 1508 pcbddc->use_change_of_basis = PETSC_FALSE; 1509 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 1510 } else { 1511 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1512 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1513 pcbddc->local_mat = matis->A; 1514 } 1515 1516 /* detect local disconnected subdomains if requested and not done before */ 1517 if (pcbddc->detect_disconnected && !pcbddc->n_local_subs) { 1518 ierr = MatDetectDisconnectedComponents(pcbddc->local_mat,PETSC_FALSE,&pcbddc->n_local_subs,&pcbddc->local_subs);CHKERRQ(ierr); 1519 } 1520 1521 /* 1522 change basis on local pressures (aka zerodiag dofs) 1523 This should come earlier then PCISSetUp for extracting the correct subdomain matrices 1524 */ 1525 if (pcbddc->benign_saddle_point) { 1526 PC_IS* pcis = (PC_IS*)pc->data; 1527 1528 if (pcbddc->user_ChangeOfBasisMatrix || pcbddc->use_change_of_basis || !computesubschurs) pcbddc->benign_change_explicit = PETSC_TRUE; 1529 /* detect local saddle point and change the basis in pcbddc->local_mat */ 1530 ierr = PCBDDCBenignDetectSaddlePoint(pc,&zerodiag);CHKERRQ(ierr); 1531 /* pop B0 mat from local mat */ 1532 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1533 /* set flag in pcis to not reuse submatrices during PCISCreate */ 1534 if (pc->flag == SAME_NONZERO_PATTERN && pcbddc->benign_n && !pcbddc->benign_change_explicit && !pcbddc->dbg_flag) { 1535 pcis->reusesubmatrices = PETSC_TRUE; 1536 } else { 1537 pcis->reusesubmatrices = PETSC_FALSE; 1538 } 1539 } 1540 /* propagate relevant information */ 1541 #if !defined(PETSC_USE_COMPLEX) /* workaround for reals */ 1542 if (matis->A->symmetric_set) { 1543 ierr = MatSetOption(pcbddc->local_mat,MAT_HERMITIAN,matis->A->symmetric);CHKERRQ(ierr); 1544 } 1545 #endif 1546 if (matis->A->symmetric_set) { 1547 ierr = MatSetOption(pcbddc->local_mat,MAT_SYMMETRIC,matis->A->symmetric);CHKERRQ(ierr); 1548 } 1549 if (matis->A->spd_set) { 1550 ierr = MatSetOption(pcbddc->local_mat,MAT_SPD,matis->A->spd);CHKERRQ(ierr); 1551 } 1552 1553 /* Set up all the "iterative substructuring" common block without computing solvers */ 1554 { 1555 Mat temp_mat; 1556 1557 temp_mat = matis->A; 1558 matis->A = pcbddc->local_mat; 1559 ierr = PCISSetUp(pc,PETSC_FALSE);CHKERRQ(ierr); 1560 pcbddc->local_mat = matis->A; 1561 matis->A = temp_mat; 1562 } 1563 1564 /* Analyze interface */ 1565 if (computetopography) { 1566 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 1567 computeconstraintsmatrix = PETSC_TRUE; 1568 } 1569 1570 /* check existence of a divergence free extension, i.e. 1571 b(v_I,p_0) = 0 for all v_I (raise error if not). 1572 Also, check that PCBDDCBenignGetOrSetP0 works */ 1573 #if defined(PETSC_USE_DEBUG) 1574 if (pcbddc->benign_saddle_point) { 1575 ierr = PCBDDCBenignCheck(pc,zerodiag);CHKERRQ(ierr); 1576 } 1577 #endif 1578 ierr = ISDestroy(&zerodiag);CHKERRQ(ierr); 1579 1580 /* Setup local dirichlet solver ksp_D and sub_schurs solvers */ 1581 if (computesolvers) { 1582 PCBDDCSubSchurs sub_schurs=pcbddc->sub_schurs; 1583 1584 if (computesubschurs && computetopography) { 1585 ierr = PCBDDCInitSubSchurs(pc);CHKERRQ(ierr); 1586 } 1587 if (sub_schurs->schur_explicit) { 1588 if (computesubschurs) { 1589 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1590 } 1591 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1592 } else { 1593 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1594 if (computesubschurs) { 1595 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1596 } 1597 } 1598 if (pcbddc->adaptive_selection) { 1599 ierr = PCBDDCAdaptiveSelection(pc);CHKERRQ(ierr); 1600 computeconstraintsmatrix = PETSC_TRUE; 1601 } 1602 } 1603 1604 /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */ 1605 new_nearnullspace_provided = PETSC_FALSE; 1606 ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr); 1607 if (pcbddc->onearnullspace) { /* already used nearnullspace */ 1608 if (!nearnullspace) { /* near null space attached to mat has been destroyed */ 1609 new_nearnullspace_provided = PETSC_TRUE; 1610 } else { 1611 /* determine if the two nullspaces are different (should be lightweight) */ 1612 if (nearnullspace != pcbddc->onearnullspace) { 1613 new_nearnullspace_provided = PETSC_TRUE; 1614 } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */ 1615 PetscInt i; 1616 const Vec *nearnullvecs; 1617 PetscObjectState state; 1618 PetscInt nnsp_size; 1619 ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 1620 for (i=0;i<nnsp_size;i++) { 1621 ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr); 1622 if (pcbddc->onearnullvecs_state[i] != state) { 1623 new_nearnullspace_provided = PETSC_TRUE; 1624 break; 1625 } 1626 } 1627 } 1628 } 1629 } else { 1630 if (!nearnullspace) { /* both nearnullspaces are null */ 1631 new_nearnullspace_provided = PETSC_FALSE; 1632 } else { /* nearnullspace attached later */ 1633 new_nearnullspace_provided = PETSC_TRUE; 1634 } 1635 } 1636 1637 /* Setup constraints and related work vectors */ 1638 /* reset primal space flags */ 1639 pcbddc->new_primal_space = PETSC_FALSE; 1640 pcbddc->new_primal_space_local = PETSC_FALSE; 1641 if (computeconstraintsmatrix || new_nearnullspace_provided) { 1642 /* It also sets the primal space flags */ 1643 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 1644 /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */ 1645 ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr); 1646 } 1647 1648 if (computesolvers || pcbddc->new_primal_space) { 1649 if (pcbddc->use_change_of_basis) { 1650 PC_IS *pcis = (PC_IS*)(pc->data); 1651 Mat temp_mat = NULL; 1652 1653 if (pcbddc->benign_change) { 1654 /* insert B0 in pcbddc->local_mat */ 1655 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_FALSE);CHKERRQ(ierr); 1656 /* swap local matrices */ 1657 ierr = MatISGetLocalMat(pc->pmat,&temp_mat);CHKERRQ(ierr); 1658 ierr = PetscObjectReference((PetscObject)temp_mat);CHKERRQ(ierr); 1659 ierr = MatISSetLocalMat(pc->pmat,pcbddc->local_mat);CHKERRQ(ierr); 1660 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1661 } 1662 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 1663 if (pcbddc->benign_change) { 1664 /* restore original matrix */ 1665 ierr = MatISSetLocalMat(pc->pmat,temp_mat);CHKERRQ(ierr); 1666 ierr = PetscObjectDereference((PetscObject)temp_mat);CHKERRQ(ierr); 1667 /* pop B0 from pcbddc->local_mat */ 1668 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1669 } 1670 /* get submatrices */ 1671 ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr); 1672 ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr); 1673 ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr); 1674 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);CHKERRQ(ierr); 1675 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);CHKERRQ(ierr); 1676 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);CHKERRQ(ierr); 1677 /* set flag in pcis to not reuse submatrices during PCISCreate */ 1678 pcis->reusesubmatrices = PETSC_FALSE; 1679 } else if (!pcbddc->user_ChangeOfBasisMatrix && !pcbddc->benign_change) { 1680 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1681 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1682 pcbddc->local_mat = matis->A; 1683 } 1684 /* SetUp coarse and local Neumann solvers */ 1685 ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr); 1686 /* SetUp Scaling operator */ 1687 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1688 } 1689 pcbddc->recompute_topography = PETSC_FALSE; 1690 1691 if (pcbddc->dbg_flag) { 1692 ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1693 } 1694 PetscFunctionReturn(0); 1695 } 1696 1697 /* -------------------------------------------------------------------------- */ 1698 /* 1699 PCApply_BDDC - Applies the BDDC operator to a vector. 1700 1701 Input Parameters: 1702 + pc - the preconditioner context 1703 - r - input vector (global) 1704 1705 Output Parameter: 1706 . z - output vector (global) 1707 1708 Application Interface Routine: PCApply() 1709 */ 1710 #undef __FUNCT__ 1711 #define __FUNCT__ "PCApply_BDDC" 1712 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 1713 { 1714 PC_IS *pcis = (PC_IS*)(pc->data); 1715 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1716 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1717 PetscErrorCode ierr; 1718 const PetscScalar one = 1.0; 1719 const PetscScalar m_one = -1.0; 1720 const PetscScalar zero = 0.0; 1721 1722 /* This code is similar to that provided in nn.c for PCNN 1723 NN interface preconditioner changed to BDDC 1724 Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */ 1725 1726 PetscFunctionBegin; 1727 if (pcbddc->benign_saddle_point) { /* get p0 from r */ 1728 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 1729 } 1730 if (!pcbddc->use_exact_dirichlet_trick) { 1731 ierr = VecCopy(r,z);CHKERRQ(ierr); 1732 /* First Dirichlet solve */ 1733 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1734 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1735 /* 1736 Assembling right hand side for BDDC operator 1737 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1738 - pcis->vec1_B the interface part of the global vector z 1739 */ 1740 if (n_D) { 1741 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1742 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1743 if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 1744 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1745 } else { 1746 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1747 } 1748 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1749 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1750 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1751 } else { 1752 if (!pcbddc->benign_apply_coarse_only) { 1753 if (pcbddc->switch_static) { 1754 ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr); 1755 } 1756 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1757 } 1758 } 1759 1760 /* Apply interface preconditioner 1761 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1762 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr); 1763 1764 /* Apply transpose of partition of unity operator */ 1765 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1766 1767 /* Second Dirichlet solve and assembling of output */ 1768 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1769 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1770 if (n_B) { 1771 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1772 if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 1773 } else if (pcbddc->switch_static) { 1774 ierr = MatMult(pcis->A_II,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 1775 } 1776 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1777 1778 if (!pcbddc->use_exact_dirichlet_trick) { 1779 if (pcbddc->switch_static) { 1780 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 1781 } else { 1782 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 1783 } 1784 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1785 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1786 } else { 1787 if (pcbddc->switch_static) { 1788 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 1789 } else { 1790 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1791 } 1792 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1793 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1794 } 1795 1796 if (pcbddc->benign_saddle_point) { /* set p0 (computed in PCBDDCApplyInterface) */ 1797 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 1798 } 1799 PetscFunctionReturn(0); 1800 } 1801 1802 /* -------------------------------------------------------------------------- */ 1803 /* 1804 PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector. 1805 1806 Input Parameters: 1807 + pc - the preconditioner context 1808 - r - input vector (global) 1809 1810 Output Parameter: 1811 . z - output vector (global) 1812 1813 Application Interface Routine: PCApplyTranspose() 1814 */ 1815 #undef __FUNCT__ 1816 #define __FUNCT__ "PCApplyTranspose_BDDC" 1817 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z) 1818 { 1819 PC_IS *pcis = (PC_IS*)(pc->data); 1820 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1821 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1822 PetscErrorCode ierr; 1823 const PetscScalar one = 1.0; 1824 const PetscScalar m_one = -1.0; 1825 const PetscScalar zero = 0.0; 1826 1827 PetscFunctionBegin; 1828 if (pcbddc->benign_saddle_point) { /* get p0 from r */ 1829 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 1830 } 1831 if (!pcbddc->use_exact_dirichlet_trick) { 1832 ierr = VecCopy(r,z);CHKERRQ(ierr); 1833 /* First Dirichlet solve */ 1834 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1835 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1836 /* 1837 Assembling right hand side for BDDC operator 1838 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1839 - pcis->vec1_B the interface part of the global vector z 1840 */ 1841 if (n_D) { 1842 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1843 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1844 if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 1845 ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1846 } else { 1847 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1848 } 1849 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1850 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1851 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1852 } else { 1853 if (pcbddc->switch_static) { 1854 ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr); 1855 } 1856 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1857 } 1858 1859 /* Apply interface preconditioner 1860 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1861 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr); 1862 1863 /* Apply transpose of partition of unity operator */ 1864 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1865 1866 /* Second Dirichlet solve and assembling of output */ 1867 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1868 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1869 if (n_B) { 1870 ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1871 if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 1872 } else if (pcbddc->switch_static) { 1873 ierr = MatMultTranspose(pcis->A_II,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 1874 } 1875 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1876 if (!pcbddc->use_exact_dirichlet_trick) { 1877 if (pcbddc->switch_static) { 1878 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 1879 } else { 1880 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 1881 } 1882 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1883 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1884 } else { 1885 if (pcbddc->switch_static) { 1886 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 1887 } else { 1888 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1889 } 1890 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1891 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1892 } 1893 if (pcbddc->benign_saddle_point) { /* set p0 (computed in PCBDDCApplyInterface) */ 1894 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 1895 } 1896 PetscFunctionReturn(0); 1897 } 1898 /* -------------------------------------------------------------------------- */ 1899 1900 #undef __FUNCT__ 1901 #define __FUNCT__ "PCDestroy_BDDC" 1902 PetscErrorCode PCDestroy_BDDC(PC pc) 1903 { 1904 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1905 PetscErrorCode ierr; 1906 1907 PetscFunctionBegin; 1908 /* free data created by PCIS */ 1909 ierr = PCISDestroy(pc);CHKERRQ(ierr); 1910 /* free BDDC custom data */ 1911 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 1912 /* destroy objects related to topography */ 1913 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 1914 /* free allocated graph structure */ 1915 ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr); 1916 /* free allocated sub schurs structure */ 1917 ierr = PetscFree(pcbddc->sub_schurs);CHKERRQ(ierr); 1918 /* destroy objects for scaling operator */ 1919 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 1920 ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr); 1921 /* free solvers stuff */ 1922 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 1923 /* free global vectors needed in presolve */ 1924 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 1925 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 1926 /* free stuff for change of basis hooks */ 1927 if (pcbddc->new_global_mat) { 1928 PCBDDCChange_ctx change_ctx; 1929 ierr = MatShellGetContext(pcbddc->new_global_mat,&change_ctx);CHKERRQ(ierr); 1930 ierr = MatDestroy(&change_ctx->original_mat);CHKERRQ(ierr); 1931 ierr = MatDestroy(&change_ctx->global_change);CHKERRQ(ierr); 1932 ierr = VecDestroyVecs(2,&change_ctx->work);CHKERRQ(ierr); 1933 ierr = PetscFree(change_ctx);CHKERRQ(ierr); 1934 } 1935 ierr = MatDestroy(&pcbddc->new_global_mat);CHKERRQ(ierr); 1936 /* remove functions */ 1937 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);CHKERRQ(ierr); 1938 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 1939 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",NULL);CHKERRQ(ierr); 1940 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 1941 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr); 1942 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr); 1943 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr); 1944 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",NULL);CHKERRQ(ierr); 1945 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 1946 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 1947 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 1948 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 1949 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 1950 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 1951 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 1952 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 1953 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 1954 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr); 1955 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 1956 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 1957 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 1958 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 1959 /* Free the private data structure */ 1960 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1961 PetscFunctionReturn(0); 1962 } 1963 /* -------------------------------------------------------------------------- */ 1964 1965 #undef __FUNCT__ 1966 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC" 1967 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 1968 { 1969 FETIDPMat_ctx mat_ctx; 1970 Vec copy_standard_rhs; 1971 PC_IS* pcis; 1972 PC_BDDC* pcbddc; 1973 PetscErrorCode ierr; 1974 1975 PetscFunctionBegin; 1976 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1977 pcis = (PC_IS*)mat_ctx->pc->data; 1978 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 1979 1980 /* 1981 change of basis for physical rhs if needed 1982 It also changes the rhs in case of dirichlet boundaries 1983 TODO: better management when FETIDP will have its own class 1984 */ 1985 ierr = VecDuplicate(standard_rhs,©_standard_rhs);CHKERRQ(ierr); 1986 ierr = VecCopy(standard_rhs,copy_standard_rhs);CHKERRQ(ierr); 1987 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,copy_standard_rhs,NULL);CHKERRQ(ierr); 1988 /* store vectors for computation of fetidp final solution */ 1989 ierr = VecScatterBegin(pcis->global_to_D,copy_standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1990 ierr = VecScatterEnd(pcis->global_to_D,copy_standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1991 /* scale rhs since it should be unassembled */ 1992 /* TODO use counter scaling? (also below) */ 1993 ierr = VecScatterBegin(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1994 ierr = VecScatterEnd(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1995 /* Apply partition of unity */ 1996 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1997 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,copy_standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 1998 if (!pcbddc->switch_static) { 1999 /* compute partially subassembled Schur complement right-hand side */ 2000 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2001 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 2002 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 2003 ierr = VecSet(copy_standard_rhs,0.0);CHKERRQ(ierr); 2004 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,copy_standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2005 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,copy_standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2006 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,copy_standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2007 ierr = VecScatterBegin(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2008 ierr = VecScatterEnd(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2009 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2010 } 2011 ierr = VecDestroy(©_standard_rhs);CHKERRQ(ierr); 2012 /* BDDC rhs */ 2013 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 2014 if (pcbddc->switch_static) { 2015 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2016 } 2017 /* apply BDDC */ 2018 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2019 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 2020 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 2021 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 2022 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2023 ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2024 PetscFunctionReturn(0); 2025 } 2026 2027 #undef __FUNCT__ 2028 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS" 2029 /*@ 2030 PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETI-DP linear system using the physical right-hand side 2031 2032 Collective 2033 2034 Input Parameters: 2035 + fetidp_mat - the FETI-DP matrix object obtained by a call to PCBDDCCreateFETIDPOperators 2036 - standard_rhs - the right-hand side of the original linear system 2037 2038 Output Parameters: 2039 . fetidp_flux_rhs - the right-hand side for the FETI-DP linear system 2040 2041 Level: developer 2042 2043 Notes: 2044 2045 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetSolution 2046 @*/ 2047 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2048 { 2049 FETIDPMat_ctx mat_ctx; 2050 PetscErrorCode ierr; 2051 2052 PetscFunctionBegin; 2053 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2054 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 2055 PetscFunctionReturn(0); 2056 } 2057 /* -------------------------------------------------------------------------- */ 2058 2059 #undef __FUNCT__ 2060 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC" 2061 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2062 { 2063 FETIDPMat_ctx mat_ctx; 2064 PC_IS* pcis; 2065 PC_BDDC* pcbddc; 2066 PetscErrorCode ierr; 2067 2068 PetscFunctionBegin; 2069 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2070 pcis = (PC_IS*)mat_ctx->pc->data; 2071 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2072 2073 /* apply B_delta^T */ 2074 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2075 ierr = VecScatterEnd (mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2076 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 2077 /* compute rhs for BDDC application */ 2078 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2079 if (pcbddc->switch_static) { 2080 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2081 } 2082 /* apply BDDC */ 2083 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2084 /* put values into standard global vector */ 2085 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2086 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2087 if (!pcbddc->switch_static) { 2088 /* compute values into the interior if solved for the partially subassembled Schur complement */ 2089 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 2090 ierr = VecAXPY(mat_ctx->temp_solution_D,-1.0,pcis->vec1_D);CHKERRQ(ierr); 2091 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2092 } 2093 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2094 ierr = VecScatterEnd (pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2095 /* final change of basis if needed 2096 Is also sums the dirichlet part removed during RHS assembling */ 2097 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 2098 PetscFunctionReturn(0); 2099 } 2100 2101 #undef __FUNCT__ 2102 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution" 2103 /*@ 2104 PCBDDCMatFETIDPGetSolution - Compute the physical solution using the solution of the FETI-DP linear system 2105 2106 Collective 2107 2108 Input Parameters: 2109 + fetidp_mat - the FETI-DP matrix obtained by a call to PCBDDCCreateFETIDPOperators 2110 - fetidp_flux_sol - the solution of the FETI-DP linear system 2111 2112 Output Parameters: 2113 . standard_sol - the solution defined on the physical domain 2114 2115 Level: developer 2116 2117 Notes: 2118 2119 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetRHS 2120 @*/ 2121 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2122 { 2123 FETIDPMat_ctx mat_ctx; 2124 PetscErrorCode ierr; 2125 2126 PetscFunctionBegin; 2127 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2128 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 2129 PetscFunctionReturn(0); 2130 } 2131 /* -------------------------------------------------------------------------- */ 2132 2133 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec); 2134 extern PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec); 2135 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat); 2136 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec); 2137 extern PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec); 2138 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC); 2139 2140 #undef __FUNCT__ 2141 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC" 2142 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 2143 { 2144 2145 FETIDPMat_ctx fetidpmat_ctx; 2146 Mat newmat; 2147 FETIDPPC_ctx fetidppc_ctx; 2148 PC newpc; 2149 MPI_Comm comm; 2150 PetscErrorCode ierr; 2151 2152 PetscFunctionBegin; 2153 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 2154 /* FETIDP linear matrix */ 2155 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 2156 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 2157 ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 2158 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 2159 ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr); 2160 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 2161 ierr = MatSetUp(newmat);CHKERRQ(ierr); 2162 /* FETIDP preconditioner */ 2163 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 2164 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 2165 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 2166 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 2167 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 2168 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 2169 ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 2170 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 2171 ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr); 2172 ierr = PCSetUp(newpc);CHKERRQ(ierr); 2173 /* return pointers for objects created */ 2174 *fetidp_mat=newmat; 2175 *fetidp_pc=newpc; 2176 PetscFunctionReturn(0); 2177 } 2178 2179 #undef __FUNCT__ 2180 #define __FUNCT__ "PCBDDCCreateFETIDPOperators" 2181 /*@ 2182 PCBDDCCreateFETIDPOperators - Create FETI-DP operators 2183 2184 Collective 2185 2186 Input Parameters: 2187 . pc - the BDDC preconditioning context (setup should have been called before) 2188 2189 Output Parameters: 2190 + fetidp_mat - shell FETI-DP matrix object 2191 - fetidp_pc - shell Dirichlet preconditioner for FETI-DP matrix 2192 2193 Options Database Keys: 2194 . -fetidp_fullyredundant <false> - use or not a fully redundant set of Lagrange multipliers 2195 2196 Level: developer 2197 2198 Notes: 2199 Currently the only operations provided for FETI-DP matrix are MatMult and MatMultTranspose 2200 2201 .seealso: PCBDDC, PCBDDCMatFETIDPGetRHS, PCBDDCMatFETIDPGetSolution 2202 @*/ 2203 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 2204 { 2205 PetscErrorCode ierr; 2206 2207 PetscFunctionBegin; 2208 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2209 if (pc->setupcalled) { 2210 ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,Mat*,PC*),(pc,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 2211 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n"); 2212 PetscFunctionReturn(0); 2213 } 2214 /* -------------------------------------------------------------------------- */ 2215 /*MC 2216 PCBDDC - Balancing Domain Decomposition by Constraints. 2217 2218 An implementation of the BDDC preconditioner based on 2219 2220 .vb 2221 [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 2222 [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 2223 [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977 2224 [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 2225 .ve 2226 2227 The matrix to be preconditioned (Pmat) must be of type MATIS. 2228 2229 Currently works with MATIS matrices with local matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers. 2230 2231 It also works with unsymmetric and indefinite problems. 2232 2233 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. 2234 2235 Approximate local solvers are automatically adapted for singular linear problems (see [1]) if the user has provided the nullspace using PCBDDCSetNullSpace() 2236 2237 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() 2238 Additional information on dofs can be provided by using PCBDDCSetDofsSplitting(), PCBDDCSetDirichletBoundaries(), PCBDDCSetNeumannBoundaries(), and PCBDDCSetPrimalVerticesIS() and their local counterparts. 2239 2240 Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). Non-singular modes are retained via SVD. 2241 2242 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. 2243 User defined change of basis can be passed to PCBDDC by using PCBDDCSetChangeOfBasisMat() 2244 2245 The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using a MatPartitioning object. 2246 2247 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. 2248 2249 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. 2250 Deluxe scaling is not supported yet for FETI-DP. 2251 2252 Options Database Keys (some of them, run with -h for a complete list): 2253 2254 . -pc_bddc_use_vertices <true> - use or not vertices in primal space 2255 . -pc_bddc_use_edges <true> - use or not edges in primal space 2256 . -pc_bddc_use_faces <false> - use or not faces in primal space 2257 . -pc_bddc_symmetric <true> - symmetric computation of primal basis functions. Specify false for unsymmetric problems 2258 . -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only) 2259 . -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested 2260 . -pc_bddc_switch_static <false> - switches from M_2 (default) to M_3 operator (see reference article [1]) 2261 . -pc_bddc_levels <0> - maximum number of levels for multilevel 2262 . -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) 2263 . -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) 2264 . -pc_bddc_use_deluxe_scaling <false> - use deluxe scaling 2265 . -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) 2266 . -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) 2267 - -pc_bddc_check_level <0> - set verbosity level of debugging output 2268 2269 Options for Dirichlet, Neumann or coarse solver can be set with 2270 .vb 2271 -pc_bddc_dirichlet_ 2272 -pc_bddc_neumann_ 2273 -pc_bddc_coarse_ 2274 .ve 2275 e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. PCBDDC uses by default KPSPREONLY and PCLU. 2276 2277 When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified as 2278 .vb 2279 -pc_bddc_dirichlet_lN_ 2280 -pc_bddc_neumann_lN_ 2281 -pc_bddc_coarse_lN_ 2282 .ve 2283 Note that level number ranges from the finest (0) to the coarsest (N). 2284 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. 2285 .vb 2286 -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3 2287 .ve 2288 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 2289 2290 Level: intermediate 2291 2292 Developer notes: 2293 2294 Contributed by Stefano Zampini 2295 2296 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 2297 M*/ 2298 2299 #undef __FUNCT__ 2300 #define __FUNCT__ "PCCreate_BDDC" 2301 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 2302 { 2303 PetscErrorCode ierr; 2304 PC_BDDC *pcbddc; 2305 2306 PetscFunctionBegin; 2307 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 2308 ierr = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr); 2309 pc->data = (void*)pcbddc; 2310 2311 /* create PCIS data structure */ 2312 ierr = PCISCreate(pc);CHKERRQ(ierr); 2313 2314 /* BDDC customization */ 2315 pcbddc->use_local_adj = PETSC_TRUE; 2316 pcbddc->use_vertices = PETSC_TRUE; 2317 pcbddc->use_edges = PETSC_TRUE; 2318 pcbddc->use_faces = PETSC_FALSE; 2319 pcbddc->use_change_of_basis = PETSC_FALSE; 2320 pcbddc->use_change_on_faces = PETSC_FALSE; 2321 pcbddc->switch_static = PETSC_FALSE; 2322 pcbddc->use_nnsp_true = PETSC_FALSE; 2323 pcbddc->use_qr_single = PETSC_FALSE; 2324 pcbddc->symmetric_primal = PETSC_TRUE; 2325 pcbddc->benign_saddle_point = PETSC_FALSE; 2326 pcbddc->vertex_size = 1; 2327 pcbddc->dbg_flag = 0; 2328 /* private */ 2329 pcbddc->local_primal_size = 0; 2330 pcbddc->local_primal_size_cc = 0; 2331 pcbddc->local_primal_ref_node = 0; 2332 pcbddc->local_primal_ref_mult = 0; 2333 pcbddc->n_vertices = 0; 2334 pcbddc->primal_indices_local_idxs = 0; 2335 pcbddc->recompute_topography = PETSC_FALSE; 2336 pcbddc->coarse_size = -1; 2337 pcbddc->new_primal_space = PETSC_FALSE; 2338 pcbddc->new_primal_space_local = PETSC_FALSE; 2339 pcbddc->global_primal_indices = 0; 2340 pcbddc->onearnullspace = 0; 2341 pcbddc->onearnullvecs_state = 0; 2342 pcbddc->user_primal_vertices = 0; 2343 pcbddc->user_primal_vertices_local = 0; 2344 pcbddc->NullSpace = 0; 2345 pcbddc->temp_solution = 0; 2346 pcbddc->original_rhs = 0; 2347 pcbddc->local_mat = 0; 2348 pcbddc->ChangeOfBasisMatrix = 0; 2349 pcbddc->user_ChangeOfBasisMatrix = 0; 2350 pcbddc->new_global_mat = 0; 2351 pcbddc->coarse_vec = 0; 2352 pcbddc->coarse_ksp = 0; 2353 pcbddc->coarse_phi_B = 0; 2354 pcbddc->coarse_phi_D = 0; 2355 pcbddc->coarse_psi_B = 0; 2356 pcbddc->coarse_psi_D = 0; 2357 pcbddc->vec1_P = 0; 2358 pcbddc->vec1_R = 0; 2359 pcbddc->vec2_R = 0; 2360 pcbddc->local_auxmat1 = 0; 2361 pcbddc->local_auxmat2 = 0; 2362 pcbddc->R_to_B = 0; 2363 pcbddc->R_to_D = 0; 2364 pcbddc->ksp_D = 0; 2365 pcbddc->ksp_R = 0; 2366 pcbddc->NeumannBoundaries = 0; 2367 pcbddc->NeumannBoundariesLocal = 0; 2368 pcbddc->DirichletBoundaries = 0; 2369 pcbddc->DirichletBoundariesLocal = 0; 2370 pcbddc->user_provided_isfordofs = PETSC_FALSE; 2371 pcbddc->n_ISForDofs = 0; 2372 pcbddc->n_ISForDofsLocal = 0; 2373 pcbddc->ISForDofs = 0; 2374 pcbddc->ISForDofsLocal = 0; 2375 pcbddc->ConstraintMatrix = 0; 2376 pcbddc->use_exact_dirichlet_trick = PETSC_TRUE; 2377 pcbddc->coarse_loc_to_glob = 0; 2378 pcbddc->coarsening_ratio = 8; 2379 pcbddc->coarse_adj_red = 0; 2380 pcbddc->current_level = 0; 2381 pcbddc->max_levels = 0; 2382 pcbddc->use_coarse_estimates = PETSC_FALSE; 2383 pcbddc->redistribute_coarse = 0; 2384 pcbddc->coarse_subassembling = 0; 2385 pcbddc->coarse_subassembling_init = 0; 2386 pcbddc->detect_disconnected = PETSC_FALSE; 2387 pcbddc->n_local_subs = 0; 2388 pcbddc->local_subs = NULL; 2389 2390 /* benign subspace trick */ 2391 pcbddc->benign_change = 0; 2392 pcbddc->benign_vec = 0; 2393 pcbddc->benign_original_mat = 0; 2394 pcbddc->benign_sf = 0; 2395 pcbddc->benign_B0 = 0; 2396 pcbddc->benign_n = 0; 2397 pcbddc->benign_p0 = NULL; 2398 pcbddc->benign_p0_lidx = NULL; 2399 pcbddc->benign_p0_gidx = NULL; 2400 pcbddc->benign_null = PETSC_FALSE; 2401 2402 /* create local graph structure */ 2403 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 2404 2405 /* scaling */ 2406 pcbddc->work_scaling = 0; 2407 pcbddc->use_deluxe_scaling = PETSC_FALSE; 2408 pcbddc->faster_deluxe = PETSC_FALSE; 2409 2410 /* create sub schurs structure */ 2411 ierr = PCBDDCSubSchursCreate(&pcbddc->sub_schurs);CHKERRQ(ierr); 2412 pcbddc->sub_schurs_rebuild = PETSC_FALSE; 2413 pcbddc->sub_schurs_layers = -1; 2414 pcbddc->sub_schurs_use_useradj = PETSC_FALSE; 2415 2416 pcbddc->computed_rowadj = PETSC_FALSE; 2417 2418 /* adaptivity */ 2419 pcbddc->adaptive_threshold = 0.0; 2420 pcbddc->adaptive_nmax = 0; 2421 pcbddc->adaptive_nmin = 0; 2422 2423 /* function pointers */ 2424 pc->ops->apply = PCApply_BDDC; 2425 pc->ops->applytranspose = PCApplyTranspose_BDDC; 2426 pc->ops->setup = PCSetUp_BDDC; 2427 pc->ops->destroy = PCDestroy_BDDC; 2428 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 2429 pc->ops->view = 0; 2430 pc->ops->applyrichardson = 0; 2431 pc->ops->applysymmetricleft = 0; 2432 pc->ops->applysymmetricright = 0; 2433 pc->ops->presolve = PCPreSolve_BDDC; 2434 pc->ops->postsolve = PCPostSolve_BDDC; 2435 2436 /* composing function */ 2437 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);CHKERRQ(ierr); 2438 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 2439 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",PCBDDCSetPrimalVerticesIS_BDDC);CHKERRQ(ierr); 2440 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 2441 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr); 2442 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr); 2443 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr); 2444 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",PCBDDCSetNullSpace_BDDC);CHKERRQ(ierr); 2445 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2446 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2447 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2448 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2449 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2450 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2451 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2452 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2453 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 2454 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr); 2455 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 2456 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 2457 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 2458 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 2459 PetscFunctionReturn(0); 2460 } 2461 2462