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