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