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