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