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