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