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