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