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