Actual source code: bddc.c

  1: #include <petsc/private/pcbddcimpl.h>
  2: #include <petsc/private/pcbddcprivateimpl.h>
  3: #include <petscblaslapack.h>

  5: static PetscBool PCBDDCPackageInitialized = PETSC_FALSE;

  7: static PetscBool  cited      = PETSC_FALSE;
  8: static const char citation[] = "@article{ZampiniPCBDDC,\n"
  9:                                "author = {Stefano Zampini},\n"
 10:                                "title = {{PCBDDC}: A Class of Robust Dual-Primal Methods in {PETS}c},\n"
 11:                                "journal = {SIAM Journal on Scientific Computing},\n"
 12:                                "volume = {38},\n"
 13:                                "number = {5},\n"
 14:                                "pages = {S282-S306},\n"
 15:                                "year = {2016},\n"
 16:                                "doi = {10.1137/15M1025785},\n"
 17:                                "URL = {http://dx.doi.org/10.1137/15M1025785},\n"
 18:                                "eprint = {http://dx.doi.org/10.1137/15M1025785}\n"
 19:                                "}\n";

 21: PetscLogEvent PC_BDDC_Topology[PETSC_PCBDDC_MAXLEVELS];
 22: PetscLogEvent PC_BDDC_LocalSolvers[PETSC_PCBDDC_MAXLEVELS];
 23: PetscLogEvent PC_BDDC_LocalWork[PETSC_PCBDDC_MAXLEVELS];
 24: PetscLogEvent PC_BDDC_CorrectionSetUp[PETSC_PCBDDC_MAXLEVELS];
 25: PetscLogEvent PC_BDDC_ApproxSetUp[PETSC_PCBDDC_MAXLEVELS];
 26: PetscLogEvent PC_BDDC_ApproxApply[PETSC_PCBDDC_MAXLEVELS];
 27: PetscLogEvent PC_BDDC_CoarseSetUp[PETSC_PCBDDC_MAXLEVELS];
 28: PetscLogEvent PC_BDDC_CoarseSolver[PETSC_PCBDDC_MAXLEVELS];
 29: PetscLogEvent PC_BDDC_AdaptiveSetUp[PETSC_PCBDDC_MAXLEVELS];
 30: PetscLogEvent PC_BDDC_Scaling[PETSC_PCBDDC_MAXLEVELS];
 31: PetscLogEvent PC_BDDC_Schurs[PETSC_PCBDDC_MAXLEVELS];
 32: PetscLogEvent PC_BDDC_Solves[PETSC_PCBDDC_MAXLEVELS][3];

 34: const char *const PCBDDCInterfaceExtTypes[] = {"DIRICHLET", "LUMP", "PCBDDCInterfaceExtType", "PC_BDDC_INTERFACE_EXT_", NULL};

 36: static PetscErrorCode PCApply_BDDC(PC, Vec, Vec);

 38: static PetscErrorCode PCSetFromOptions_BDDC(PC pc, PetscOptionItems *PetscOptionsObject)
 39: {
 40:   PC_BDDC  *pcbddc = (PC_BDDC *)pc->data;
 41:   PetscInt  nt, i;
 42:   char      load[PETSC_MAX_PATH_LEN] = {'\0'};
 43:   PetscBool flg;

 45:   PetscFunctionBegin;
 46:   PetscOptionsHeadBegin(PetscOptionsObject, "BDDC options");
 47:   /* Load customization from binary file (debugging) */
 48:   PetscCall(PetscOptionsString("-pc_bddc_load", "Load customization from file (intended for debug)", "none", load, load, sizeof(load), &flg));
 49:   if (flg) {
 50:     size_t len;

 52:     PetscCall(PetscStrlen(load, &len));
 53:     PetscCall(PCBDDCLoadOrViewCustomization(pc, PETSC_TRUE, len ? load : NULL));
 54:   }
 55:   /* Verbose debugging */
 56:   PetscCall(PetscOptionsInt("-pc_bddc_check_level", "Verbose output for PCBDDC (intended for debug)", "none", pcbddc->dbg_flag, &pcbddc->dbg_flag, NULL));
 57:   /* Approximate solvers */
 58:   PetscCall(PetscOptionsEnum("-pc_bddc_interface_ext_type", "Use DIRICHLET or LUMP to extend interface corrections to interior", "PCBDDCSetInterfaceExtType", PCBDDCInterfaceExtTypes, (PetscEnum)pcbddc->interface_extension, (PetscEnum *)&pcbddc->interface_extension, NULL));
 59:   if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_DIRICHLET) {
 60:     PetscCall(PetscOptionsBool("-pc_bddc_dirichlet_approximate", "Inform PCBDDC that we are using approximate Dirichlet solvers", "none", pcbddc->NullSpace_corr[0], &pcbddc->NullSpace_corr[0], NULL));
 61:     PetscCall(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));
 62:   } else {
 63:     /* This flag is needed/implied by lumping */
 64:     pcbddc->switch_static = PETSC_TRUE;
 65:   }
 66:   PetscCall(PetscOptionsBool("-pc_bddc_neumann_approximate", "Inform PCBDDC that we are using approximate Neumann solvers", "none", pcbddc->NullSpace_corr[2], &pcbddc->NullSpace_corr[2], NULL));
 67:   PetscCall(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));
 68:   /* Primal space customization */
 69:   PetscCall(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));
 70:   PetscCall(PetscOptionsInt("-pc_bddc_local_mat_graph_square", "Square adjacency graph of local mat for interface analysis", "none", pcbddc->local_adj_square, &pcbddc->local_adj_square, NULL));
 71:   PetscCall(PetscOptionsInt("-pc_bddc_graph_maxcount", "Maximum number of shared subdomains for a connected component", "none", pcbddc->graphmaxcount, &pcbddc->graphmaxcount, NULL));
 72:   PetscCall(PetscOptionsBool("-pc_bddc_corner_selection", "Activates face-based corner selection", "none", pcbddc->corner_selection, &pcbddc->corner_selection, NULL));
 73:   PetscCall(PetscOptionsBool("-pc_bddc_use_vertices", "Use or not corner dofs in coarse space", "none", pcbddc->use_vertices, &pcbddc->use_vertices, NULL));
 74:   PetscCall(PetscOptionsBool("-pc_bddc_use_edges", "Use or not edge constraints in coarse space", "none", pcbddc->use_edges, &pcbddc->use_edges, NULL));
 75:   PetscCall(PetscOptionsBool("-pc_bddc_use_faces", "Use or not face constraints in coarse space", "none", pcbddc->use_faces, &pcbddc->use_faces, NULL));
 76:   PetscCall(PetscOptionsInt("-pc_bddc_vertex_size", "Connected components smaller or equal to vertex size will be considered as primal vertices", "none", pcbddc->vertex_size, &pcbddc->vertex_size, NULL));
 77:   PetscCall(PetscOptionsBool("-pc_bddc_use_nnsp", "Use near null space attached to the matrix to compute constraints", "none", pcbddc->use_nnsp, &pcbddc->use_nnsp, NULL));
 78:   PetscCall(PetscOptionsBool("-pc_bddc_use_nnsp_true", "Use near null space attached to the matrix to compute constraints as is", "none", pcbddc->use_nnsp_true, &pcbddc->use_nnsp_true, NULL));
 79:   PetscCall(PetscOptionsBool("-pc_bddc_use_qr_single", "Use QR factorization for single constraints on cc (QR is always used when multiple constraints are present)", "none", pcbddc->use_qr_single, &pcbddc->use_qr_single, NULL));
 80:   /* Change of basis */
 81:   PetscCall(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));
 82:   PetscCall(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));
 83:   if (!pcbddc->use_change_of_basis) pcbddc->use_change_on_faces = PETSC_FALSE;
 84:   /* Switch between M_2 (default) and M_3 preconditioners (as defined by C. Dohrmann in the ref. article) */
 85:   PetscCall(PetscOptionsBool("-pc_bddc_switch_static", "Switch on static condensation ops around the interface preconditioner", "none", pcbddc->switch_static, &pcbddc->switch_static, NULL));
 86:   PetscCall(PetscOptionsInt("-pc_bddc_coarse_eqs_per_proc", "Target number of equations per process for coarse problem redistribution (significant only at the coarsest level)", "none", pcbddc->coarse_eqs_per_proc, &pcbddc->coarse_eqs_per_proc, NULL));
 87:   i = pcbddc->coarsening_ratio;
 88:   PetscCall(PetscOptionsInt("-pc_bddc_coarsening_ratio", "Set coarsening ratio used in multilevel coarsening", "PCBDDCSetCoarseningRatio", i, &i, NULL));
 89:   PetscCall(PCBDDCSetCoarseningRatio(pc, i));
 90:   i = pcbddc->max_levels;
 91:   PetscCall(PetscOptionsInt("-pc_bddc_levels", "Set maximum number of levels for multilevel", "PCBDDCSetLevels", i, &i, NULL));
 92:   PetscCall(PCBDDCSetLevels(pc, i));
 93:   PetscCall(PetscOptionsInt("-pc_bddc_coarse_eqs_limit", "Set maximum number of equations on coarsest grid to aim for", "none", pcbddc->coarse_eqs_limit, &pcbddc->coarse_eqs_limit, NULL));
 94:   PetscCall(PetscOptionsBool("-pc_bddc_use_coarse_estimates", "Use estimated eigenvalues for coarse problem", "none", pcbddc->use_coarse_estimates, &pcbddc->use_coarse_estimates, NULL));
 95:   PetscCall(PetscOptionsBool("-pc_bddc_use_deluxe_scaling", "Use deluxe scaling for BDDC", "none", pcbddc->use_deluxe_scaling, &pcbddc->use_deluxe_scaling, NULL));
 96:   PetscCall(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));
 97:   PetscCall(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));
 98:   PetscCall(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));
 99:   PetscCall(PetscOptionsBool("-pc_bddc_schur_exact", "Whether or not to use the exact Schur complement instead of the reduced one (which excludes size 1 cc)", "none", pcbddc->sub_schurs_exact_schur, &pcbddc->sub_schurs_exact_schur, NULL));
100:   PetscCall(PetscOptionsBool("-pc_bddc_deluxe_zerorows", "Zero rows and columns of deluxe operators associated with primal dofs", "none", pcbddc->deluxe_zerorows, &pcbddc->deluxe_zerorows, NULL));
101:   PetscCall(PetscOptionsBool("-pc_bddc_deluxe_singlemat", "Collapse deluxe operators", "none", pcbddc->deluxe_singlemat, &pcbddc->deluxe_singlemat, NULL));
102:   PetscCall(PetscOptionsBool("-pc_bddc_adaptive_userdefined", "Use user-defined constraints (should be attached via MatSetNearNullSpace to pmat) in addition to those adaptively generated", "none", pcbddc->adaptive_userdefined, &pcbddc->adaptive_userdefined, NULL));
103:   nt = 2;
104:   PetscCall(PetscOptionsRealArray("-pc_bddc_adaptive_threshold", "Thresholds to be used for adaptive selection of constraints", "none", pcbddc->adaptive_threshold, &nt, NULL));
105:   if (nt == 1) pcbddc->adaptive_threshold[1] = pcbddc->adaptive_threshold[0];
106:   PetscCall(PetscOptionsInt("-pc_bddc_adaptive_nmin", "Minimum number of constraints per connected components", "none", pcbddc->adaptive_nmin, &pcbddc->adaptive_nmin, NULL));
107:   PetscCall(PetscOptionsInt("-pc_bddc_adaptive_nmax", "Maximum number of constraints per connected components", "none", pcbddc->adaptive_nmax, &pcbddc->adaptive_nmax, NULL));
108:   PetscCall(PetscOptionsBool("-pc_bddc_symmetric", "Symmetric computation of primal basis functions", "none", pcbddc->symmetric_primal, &pcbddc->symmetric_primal, NULL));
109:   PetscCall(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));
110:   PetscCall(PetscOptionsBool("-pc_bddc_benign_trick", "Apply the benign subspace trick to saddle point problems with discontinuous pressures", "none", pcbddc->benign_saddle_point, &pcbddc->benign_saddle_point, NULL));
111:   PetscCall(PetscOptionsBool("-pc_bddc_benign_change", "Compute the pressure change of basis explicitly", "none", pcbddc->benign_change_explicit, &pcbddc->benign_change_explicit, NULL));
112:   PetscCall(PetscOptionsBool("-pc_bddc_benign_compute_correction", "Compute the benign correction during PreSolve", "none", pcbddc->benign_compute_correction, &pcbddc->benign_compute_correction, NULL));
113:   PetscCall(PetscOptionsBool("-pc_bddc_nonetflux", "Automatic computation of no-net-flux quadrature weights", "none", pcbddc->compute_nonetflux, &pcbddc->compute_nonetflux, NULL));
114:   PetscCall(PetscOptionsBool("-pc_bddc_detect_disconnected", "Detects disconnected subdomains", "none", pcbddc->detect_disconnected, &pcbddc->detect_disconnected, NULL));
115:   PetscCall(PetscOptionsBool("-pc_bddc_detect_disconnected_filter", "Filters out small entries in the local matrix when detecting disconnected subdomains", "none", pcbddc->detect_disconnected_filter, &pcbddc->detect_disconnected_filter, NULL));
116:   PetscCall(PetscOptionsBool("-pc_bddc_eliminate_dirichlet", "Whether or not we want to eliminate dirichlet dofs during presolve", "none", pcbddc->eliminate_dirdofs, &pcbddc->eliminate_dirdofs, NULL));
117:   PetscOptionsHeadEnd();
118:   PetscFunctionReturn(PETSC_SUCCESS);
119: }

121: static PetscErrorCode PCView_BDDC(PC pc, PetscViewer viewer)
122: {
123:   PC_BDDC     *pcbddc = (PC_BDDC *)pc->data;
124:   PC_IS       *pcis   = (PC_IS *)pc->data;
125:   PetscBool    isascii;
126:   PetscSubcomm subcomm;
127:   PetscViewer  subviewer;

129:   PetscFunctionBegin;
130:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
131:   /* ASCII viewer */
132:   if (isascii) {
133:     PetscMPIInt color, rank, size;
134:     PetscInt64  loc[7], gsum[6], gmax[6], gmin[6], totbenign;
135:     PetscScalar interface_size;
136:     PetscReal   ratio1 = 0., ratio2 = 0.;
137:     Vec         counter;

139:     if (!pc->setupcalled) PetscCall(PetscViewerASCIIPrintf(viewer, "  Partial information available: preconditioner has not been setup yet\n"));
140:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use verbose output: %" PetscInt_FMT "\n", pcbddc->dbg_flag));
141:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use user-defined CSR: %d\n", !!pcbddc->mat_graph->nvtxs_csr));
142:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use local mat graph: %d\n", pcbddc->use_local_adj && !pcbddc->mat_graph->nvtxs_csr));
143:     if (pcbddc->mat_graph->twodim) {
144:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Connectivity graph topological dimension: 2\n"));
145:     } else {
146:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Connectivity graph topological dimension: 3\n"));
147:     }
148:     if (pcbddc->graphmaxcount != PETSC_INT_MAX) PetscCall(PetscViewerASCIIPrintf(viewer, "  Graph max count: %" PetscInt_FMT "\n", pcbddc->graphmaxcount));
149:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Corner selection: %d (selected %d)\n", pcbddc->corner_selection, pcbddc->corner_selected));
150:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use vertices: %d (vertex size %" PetscInt_FMT ")\n", pcbddc->use_vertices, pcbddc->vertex_size));
151:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use edges: %d\n", pcbddc->use_edges));
152:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use faces: %d\n", pcbddc->use_faces));
153:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use true near null space: %d\n", pcbddc->use_nnsp_true));
154:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use QR for single constraints on cc: %d\n", pcbddc->use_qr_single));
155:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use change of basis on local edge nodes: %d\n", pcbddc->use_change_of_basis));
156:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use change of basis on local face nodes: %d\n", pcbddc->use_change_on_faces));
157:     PetscCall(PetscViewerASCIIPrintf(viewer, "  User defined change of basis matrix: %d\n", !!pcbddc->user_ChangeOfBasisMatrix));
158:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Has change of basis matrix: %d\n", !!pcbddc->ChangeOfBasisMatrix));
159:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Eliminate dirichlet boundary dofs: %d\n", pcbddc->eliminate_dirdofs));
160:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Switch on static condensation ops around the interface preconditioner: %d\n", pcbddc->switch_static));
161:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use exact dirichlet trick: %d\n", pcbddc->use_exact_dirichlet_trick));
162:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Interface extension: %s\n", PCBDDCInterfaceExtTypes[pcbddc->interface_extension]));
163:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Multilevel max levels: %" PetscInt_FMT "\n", pcbddc->max_levels));
164:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Multilevel coarsening ratio: %" PetscInt_FMT "\n", pcbddc->coarsening_ratio));
165:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use estimated eigs for coarse problem: %d\n", pcbddc->use_coarse_estimates));
166:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use deluxe scaling: %d\n", pcbddc->use_deluxe_scaling));
167:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use deluxe zerorows: %d\n", pcbddc->deluxe_zerorows));
168:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use deluxe singlemat: %d\n", pcbddc->deluxe_singlemat));
169:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Rebuild interface graph for Schur principal minors: %d\n", pcbddc->sub_schurs_rebuild));
170:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Number of dofs' layers for the computation of principal minors: %" PetscInt_FMT "\n", pcbddc->sub_schurs_layers));
171:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Use user CSR graph to compute successive layers: %d\n", pcbddc->sub_schurs_use_useradj));
172:     if (pcbddc->adaptive_threshold[1] != pcbddc->adaptive_threshold[0]) {
173:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Adaptive constraint selection thresholds (active %d, userdefined %d): %g,%g\n", pcbddc->adaptive_selection, pcbddc->adaptive_userdefined, (double)pcbddc->adaptive_threshold[0], (double)pcbddc->adaptive_threshold[1]));
174:     } else {
175:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Adaptive constraint selection threshold (active %d, userdefined %d): %g\n", pcbddc->adaptive_selection, pcbddc->adaptive_userdefined, (double)pcbddc->adaptive_threshold[0]));
176:     }
177:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Min constraints / connected component: %" PetscInt_FMT "\n", pcbddc->adaptive_nmin));
178:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Max constraints / connected component: %" PetscInt_FMT "\n", pcbddc->adaptive_nmax));
179:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Invert exact Schur complement for adaptive selection: %d\n", pcbddc->sub_schurs_exact_schur));
180:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Symmetric computation of primal basis functions: %d\n", pcbddc->symmetric_primal));
181:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Num. Procs. to map coarse adjacency list: %" PetscInt_FMT "\n", pcbddc->coarse_adj_red));
182:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Coarse eqs per proc (significant at the coarsest level): %" PetscInt_FMT "\n", pcbddc->coarse_eqs_per_proc));
183:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Detect disconnected: %d (filter %d)\n", pcbddc->detect_disconnected, pcbddc->detect_disconnected_filter));
184:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Benign subspace trick: %d (change explicit %d)\n", pcbddc->benign_saddle_point, pcbddc->benign_change_explicit));
185:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Benign subspace trick is active: %d\n", pcbddc->benign_have_null));
186:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Algebraic computation of no-net-flux: %d\n", pcbddc->compute_nonetflux));
187:     if (!pc->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);

189:     /* compute interface size */
190:     PetscCall(VecSet(pcis->vec1_B, 1.0));
191:     PetscCall(MatCreateVecs(pc->pmat, &counter, NULL));
192:     PetscCall(VecSet(counter, 0.0));
193:     PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, counter, INSERT_VALUES, SCATTER_REVERSE));
194:     PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, counter, INSERT_VALUES, SCATTER_REVERSE));
195:     PetscCall(VecSum(counter, &interface_size));
196:     PetscCall(VecDestroy(&counter));

198:     /* compute some statistics on the domain decomposition */
199:     gsum[0] = 1;
200:     gsum[1] = gsum[2] = gsum[3] = gsum[4] = gsum[5] = 0;
201:     loc[0]                                          = !!pcis->n;
202:     loc[1]                                          = pcis->n - pcis->n_B;
203:     loc[2]                                          = pcis->n_B;
204:     loc[3]                                          = pcbddc->local_primal_size;
205:     loc[4]                                          = pcis->n;
206:     loc[5]                                          = pcbddc->n_local_subs > 0 ? pcbddc->n_local_subs : (pcis->n ? 1 : 0);
207:     loc[6]                                          = pcbddc->benign_n;
208:     PetscCallMPI(MPI_Reduce(loc, gsum, 6, MPIU_INT64, MPI_SUM, 0, PetscObjectComm((PetscObject)pc)));
209:     if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = loc[5] = -1;
210:     PetscCallMPI(MPI_Reduce(loc, gmax, 6, MPIU_INT64, MPI_MAX, 0, PetscObjectComm((PetscObject)pc)));
211:     if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = loc[5] = PETSC_INT_MAX;
212:     PetscCallMPI(MPI_Reduce(loc, gmin, 6, MPIU_INT64, MPI_MIN, 0, PetscObjectComm((PetscObject)pc)));
213:     PetscCallMPI(MPI_Reduce(&loc[6], &totbenign, 1, MPIU_INT64, MPI_SUM, 0, PetscObjectComm((PetscObject)pc)));
214:     if (pcbddc->coarse_size) {
215:       ratio1 = pc->pmat->rmap->N / (1. * pcbddc->coarse_size);
216:       ratio2 = PetscRealPart(interface_size) / pcbddc->coarse_size;
217:     }
218:     PetscCall(PetscViewerASCIIPrintf(viewer, "********************************** STATISTICS AT LEVEL %" PetscInt_FMT " **********************************\n", pcbddc->current_level));
219:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Global dofs sizes: all %" PetscInt_FMT " interface %" PetscInt_FMT " coarse %" PetscInt_FMT "\n", pc->pmat->rmap->N, (PetscInt)PetscRealPart(interface_size), pcbddc->coarse_size));
220:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Coarsening ratios: all/coarse %" PetscInt_FMT " interface/coarse %" PetscInt_FMT "\n", (PetscInt)ratio1, (PetscInt)ratio2));
221:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Active processes : %" PetscInt64_FMT "\n", gsum[0]));
222:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Total subdomains : %" PetscInt64_FMT "\n", gsum[5]));
223:     if (pcbddc->benign_have_null) PetscCall(PetscViewerASCIIPrintf(viewer, "  Benign subs      : %" PetscInt64_FMT "\n", totbenign));
224:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Dofs type        :\tMIN\tMAX\tMEAN\n"));
225:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Interior  dofs   :\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\n", gmin[1], gmax[1], gsum[1] / gsum[0]));
226:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Interface dofs   :\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\n", gmin[2], gmax[2], gsum[2] / gsum[0]));
227:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Primal    dofs   :\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\n", gmin[3], gmax[3], gsum[3] / gsum[0]));
228:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Local     dofs   :\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\n", gmin[4], gmax[4], gsum[4] / gsum[0]));
229:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Local     subs   :\t%" PetscInt64_FMT "\t%" PetscInt64_FMT "\n", gmin[5], gmax[5]));
230:     PetscCall(PetscViewerFlush(viewer));

232:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));

234:     /* local solvers */
235:     PetscCall(PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)pcbddc->ksp_D), &subviewer));
236:     if (rank == 0) {
237:       PetscCall(PetscViewerASCIIPrintf(subviewer, "--- Interior solver (rank 0)\n"));
238:       PetscCall(PetscViewerASCIIPushTab(subviewer));
239:       PetscCall(KSPView(pcbddc->ksp_D, subviewer));
240:       PetscCall(PetscViewerASCIIPopTab(subviewer));
241:       PetscCall(PetscViewerASCIIPrintf(subviewer, "--- Correction solver (rank 0)\n"));
242:       PetscCall(PetscViewerASCIIPushTab(subviewer));
243:       PetscCall(KSPView(pcbddc->ksp_R, subviewer));
244:       PetscCall(PetscViewerASCIIPopTab(subviewer));
245:       PetscCall(PetscViewerFlush(subviewer));
246:     }
247:     PetscCall(PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)pcbddc->ksp_D), &subviewer));
248:     /* the coarse problem can be handled by a different communicator */
249:     if (pcbddc->coarse_ksp) color = 1;
250:     else color = 0;
251:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
252:     PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)pc), &subcomm));
253:     PetscCall(PetscSubcommSetNumber(subcomm, PetscMin(size, 2)));
254:     PetscCall(PetscSubcommSetTypeGeneral(subcomm, color, rank));
255:     PetscCall(PetscViewerGetSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer));
256:     if (color == 1) {
257:       PetscCall(PetscViewerASCIIPrintf(subviewer, "--- Coarse solver\n"));
258:       PetscCall(PetscViewerASCIIPushTab(subviewer));
259:       PetscCall(KSPView(pcbddc->coarse_ksp, subviewer));
260:       PetscCall(PetscViewerASCIIPopTab(subviewer));
261:       PetscCall(PetscViewerFlush(subviewer));
262:     }
263:     PetscCall(PetscViewerRestoreSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer));
264:     PetscCall(PetscSubcommDestroy(&subcomm));
265:     PetscCall(PetscViewerFlush(viewer));
266:   }
267:   PetscFunctionReturn(PETSC_SUCCESS);
268: }

270: static PetscErrorCode PCBDDCSetDiscreteGradient_BDDC(PC pc, Mat G, PetscInt order, PetscInt field, PetscBool global, PetscBool conforming)
271: {
272:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

274:   PetscFunctionBegin;
275:   PetscCall(PetscObjectReference((PetscObject)G));
276:   PetscCall(MatDestroy(&pcbddc->discretegradient));
277:   pcbddc->discretegradient = G;
278:   pcbddc->nedorder         = order > 0 ? order : -order;
279:   pcbddc->nedfield         = field;
280:   pcbddc->nedglobal        = global;
281:   pcbddc->conforming       = conforming;
282:   PetscFunctionReturn(PETSC_SUCCESS);
283: }

285: /*@
286:   PCBDDCSetDiscreteGradient - Sets the discrete gradient to be used by the `PCBDDC` preconditioner

288:   Collective

290:   Input Parameters:
291: + pc         - the preconditioning context
292: . G          - the discrete gradient matrix (in `MATAIJ` format)
293: . order      - the order of the Nedelec space (1 for the lowest order)
294: . field      - the field id of the Nedelec dofs (not used if the fields have not been specified)
295: . global     - the type of global ordering for the rows of `G`
296: - conforming - whether the mesh is conforming or not

298:   Level: advanced

300:   Note:
301:   The discrete gradient matrix `G` is used to analyze the subdomain edges, and it should not contain any zero entry.
302:   For variable order spaces, the order should be set to zero.
303:   If `global` is `PETSC_TRUE`, the rows of `G` should be given in global ordering for the whole dofs;
304:   if `PETSC_FALSE`, the ordering should be global for the Nedelec field.
305:   In the latter case, it should hold gid[i] < gid[j] iff geid[i] < geid[j], with gid the global orderding for all the dofs
306:   and geid the one for the Nedelec field.

308: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDofsSplitting()`, `PCBDDCSetDofsSplittingLocal()`, `MATAIJ`, `PCBDDCSetDivergenceMat()`
309: @*/
310: PetscErrorCode PCBDDCSetDiscreteGradient(PC pc, Mat G, PetscInt order, PetscInt field, PetscBool global, PetscBool conforming)
311: {
312:   PetscFunctionBegin;
319:   PetscCheckSameComm(pc, 1, G, 2);
320:   PetscTryMethod(pc, "PCBDDCSetDiscreteGradient_C", (PC, Mat, PetscInt, PetscInt, PetscBool, PetscBool), (pc, G, order, field, global, conforming));
321:   PetscFunctionReturn(PETSC_SUCCESS);
322: }

324: static PetscErrorCode PCBDDCSetDivergenceMat_BDDC(PC pc, Mat divudotp, PetscBool trans, IS vl2l)
325: {
326:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

328:   PetscFunctionBegin;
329:   PetscCall(PetscObjectReference((PetscObject)divudotp));
330:   PetscCall(MatDestroy(&pcbddc->divudotp));
331:   pcbddc->divudotp          = divudotp;
332:   pcbddc->divudotp_trans    = trans;
333:   pcbddc->compute_nonetflux = PETSC_TRUE;
334:   if (vl2l) {
335:     PetscCall(PetscObjectReference((PetscObject)vl2l));
336:     PetscCall(ISDestroy(&pcbddc->divudotp_vl2l));
337:     pcbddc->divudotp_vl2l = vl2l;
338:   }
339:   PetscFunctionReturn(PETSC_SUCCESS);
340: }

342: /*@
343:   PCBDDCSetDivergenceMat - Sets the linear operator representing \int_\Omega \div {\bf u} \cdot p dx for the `PCBDDC` preconditioner

345:   Collective

347:   Input Parameters:
348: + pc       - the preconditioning context
349: . divudotp - the matrix (must be of type `MATIS`)
350: . trans    - if `PETSC_FALSE` (resp. `PETSC_TRUE`), then pressures are in the test (trial) space and velocities are in the trial (test) space.
351: - vl2l     - optional index set describing the local (wrt the local matrix in `divudotp`) to local (wrt the local matrix
352:    in the preconditioning matrix) map for the velocities

354:   Level: advanced

356:   Notes:
357:   This auxiliary matrix is used to compute quadrature weights representing the net-flux across subdomain boundaries

359:   If `vl2l` is `NULL`, the local ordering for velocities in `divudotp` should match that of the preconditioning matrix

361: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDiscreteGradient()`
362: @*/
363: PetscErrorCode PCBDDCSetDivergenceMat(PC pc, Mat divudotp, PetscBool trans, IS vl2l)
364: {
365:   PetscBool ismatis;

367:   PetscFunctionBegin;
370:   PetscCheckSameComm(pc, 1, divudotp, 2);
373:   PetscCall(PetscObjectTypeCompare((PetscObject)divudotp, MATIS, &ismatis));
374:   PetscCheck(ismatis, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Divergence matrix needs to be of type MATIS");
375:   PetscTryMethod(pc, "PCBDDCSetDivergenceMat_C", (PC, Mat, PetscBool, IS), (pc, divudotp, trans, vl2l));
376:   PetscFunctionReturn(PETSC_SUCCESS);
377: }

379: static PetscErrorCode PCBDDCSetChangeOfBasisMat_BDDC(PC pc, Mat change, PetscBool interior)
380: {
381:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

383:   PetscFunctionBegin;
384:   PetscCall(PetscObjectReference((PetscObject)change));
385:   PetscCall(MatDestroy(&pcbddc->user_ChangeOfBasisMatrix));
386:   pcbddc->user_ChangeOfBasisMatrix = change;
387:   pcbddc->change_interior          = interior;
388:   PetscFunctionReturn(PETSC_SUCCESS);
389: }

391: /*@
392:   PCBDDCSetChangeOfBasisMat - Set user defined change of basis for dofs

394:   Collective

396:   Input Parameters:
397: + pc       - the preconditioning context
398: . change   - the change of basis matrix
399: - interior - whether or not the change of basis modifies interior dofs

401:   Level: intermediate

403: .seealso: [](ch_ksp), `PCBDDC`
404: @*/
405: PetscErrorCode PCBDDCSetChangeOfBasisMat(PC pc, Mat change, PetscBool interior)
406: {
407:   PetscFunctionBegin;
410:   PetscCheckSameComm(pc, 1, change, 2);
411:   if (pc->mat) {
412:     PetscInt rows_c, cols_c, rows, cols;
413:     PetscCall(MatGetSize(pc->mat, &rows, &cols));
414:     PetscCall(MatGetSize(change, &rows_c, &cols_c));
415:     PetscCheck(rows_c == rows, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Invalid number of rows for change of basis matrix! %" PetscInt_FMT " != %" PetscInt_FMT, rows_c, rows);
416:     PetscCheck(cols_c == cols, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Invalid number of columns for change of basis matrix! %" PetscInt_FMT " != %" PetscInt_FMT, cols_c, cols);
417:     PetscCall(MatGetLocalSize(pc->mat, &rows, &cols));
418:     PetscCall(MatGetLocalSize(change, &rows_c, &cols_c));
419:     PetscCheck(rows_c == rows, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Invalid number of local rows for change of basis matrix! %" PetscInt_FMT " != %" PetscInt_FMT, rows_c, rows);
420:     PetscCheck(cols_c == cols, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Invalid number of local columns for change of basis matrix! %" PetscInt_FMT " != %" PetscInt_FMT, cols_c, cols);
421:   }
422:   PetscTryMethod(pc, "PCBDDCSetChangeOfBasisMat_C", (PC, Mat, PetscBool), (pc, change, interior));
423:   PetscFunctionReturn(PETSC_SUCCESS);
424: }

426: static PetscErrorCode PCBDDCSetPrimalVerticesIS_BDDC(PC pc, IS PrimalVertices)
427: {
428:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
429:   PetscBool isequal = PETSC_FALSE;

431:   PetscFunctionBegin;
432:   PetscCall(PetscObjectReference((PetscObject)PrimalVertices));
433:   if (pcbddc->user_primal_vertices) PetscCall(ISEqual(PrimalVertices, pcbddc->user_primal_vertices, &isequal));
434:   PetscCall(ISDestroy(&pcbddc->user_primal_vertices));
435:   PetscCall(ISDestroy(&pcbddc->user_primal_vertices_local));
436:   pcbddc->user_primal_vertices = PrimalVertices;
437:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
438:   PetscFunctionReturn(PETSC_SUCCESS);
439: }

441: /*@
442:   PCBDDCSetPrimalVerticesIS - Set additional user defined primal vertices in `PCBDDC`

444:   Collective

446:   Input Parameters:
447: + pc             - the preconditioning context
448: - PrimalVertices - index set of primal vertices in global numbering (can be empty)

450:   Level: intermediate

452:   Note:
453:   Any process can list any global node

455: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
456: @*/
457: PetscErrorCode PCBDDCSetPrimalVerticesIS(PC pc, IS PrimalVertices)
458: {
459:   PetscFunctionBegin;
462:   PetscCheckSameComm(pc, 1, PrimalVertices, 2);
463:   PetscTryMethod(pc, "PCBDDCSetPrimalVerticesIS_C", (PC, IS), (pc, PrimalVertices));
464:   PetscFunctionReturn(PETSC_SUCCESS);
465: }

467: static PetscErrorCode PCBDDCGetPrimalVerticesIS_BDDC(PC pc, IS *is)
468: {
469:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

471:   PetscFunctionBegin;
472:   *is = pcbddc->user_primal_vertices;
473:   PetscFunctionReturn(PETSC_SUCCESS);
474: }

476: /*@
477:   PCBDDCGetPrimalVerticesIS - Get user defined primal vertices set with `PCBDDCSetPrimalVerticesIS()`

479:   Collective

481:   Input Parameter:
482: . pc - the preconditioning context

484:   Output Parameter:
485: . is - index set of primal vertices in global numbering (`NULL` if not set)

487:   Level: intermediate

489: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
490: @*/
491: PetscErrorCode PCBDDCGetPrimalVerticesIS(PC pc, IS *is)
492: {
493:   PetscFunctionBegin;
495:   PetscAssertPointer(is, 2);
496:   PetscUseMethod(pc, "PCBDDCGetPrimalVerticesIS_C", (PC, IS *), (pc, is));
497:   PetscFunctionReturn(PETSC_SUCCESS);
498: }

500: static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices)
501: {
502:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
503:   PetscBool isequal = PETSC_FALSE;

505:   PetscFunctionBegin;
506:   PetscCall(PetscObjectReference((PetscObject)PrimalVertices));
507:   if (pcbddc->user_primal_vertices_local) PetscCall(ISEqual(PrimalVertices, pcbddc->user_primal_vertices_local, &isequal));
508:   PetscCall(ISDestroy(&pcbddc->user_primal_vertices));
509:   PetscCall(ISDestroy(&pcbddc->user_primal_vertices_local));
510:   pcbddc->user_primal_vertices_local = PrimalVertices;
511:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
512:   PetscFunctionReturn(PETSC_SUCCESS);
513: }

515: /*@
516:   PCBDDCSetPrimalVerticesLocalIS - Set additional user defined primal vertices in `PCBDDC`

518:   Collective

520:   Input Parameters:
521: + pc             - the preconditioning context
522: - PrimalVertices - index set of primal vertices in local numbering (can be empty)

524:   Level: intermediate

526: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
527: @*/
528: PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices)
529: {
530:   PetscFunctionBegin;
533:   PetscCheckSameComm(pc, 1, PrimalVertices, 2);
534:   PetscTryMethod(pc, "PCBDDCSetPrimalVerticesLocalIS_C", (PC, IS), (pc, PrimalVertices));
535:   PetscFunctionReturn(PETSC_SUCCESS);
536: }

538: static PetscErrorCode PCBDDCGetPrimalVerticesLocalIS_BDDC(PC pc, IS *is)
539: {
540:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

542:   PetscFunctionBegin;
543:   *is = pcbddc->user_primal_vertices_local;
544:   PetscFunctionReturn(PETSC_SUCCESS);
545: }

547: /*@
548:   PCBDDCGetPrimalVerticesLocalIS - Get user defined primal vertices set with `PCBDDCSetPrimalVerticesLocalIS()`

550:   Collective

552:   Input Parameter:
553: . pc - the preconditioning context

555:   Output Parameter:
556: . is - index set of primal vertices in local numbering (`NULL` if not set)

558:   Level: intermediate

560: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`
561: @*/
562: PetscErrorCode PCBDDCGetPrimalVerticesLocalIS(PC pc, IS *is)
563: {
564:   PetscFunctionBegin;
566:   PetscAssertPointer(is, 2);
567:   PetscUseMethod(pc, "PCBDDCGetPrimalVerticesLocalIS_C", (PC, IS *), (pc, is));
568:   PetscFunctionReturn(PETSC_SUCCESS);
569: }

571: static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc, PetscInt k)
572: {
573:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

575:   PetscFunctionBegin;
576:   pcbddc->coarsening_ratio = k;
577:   PetscFunctionReturn(PETSC_SUCCESS);
578: }

580: /*@
581:   PCBDDCSetCoarseningRatio - Set coarsening ratio used in the multi-level version of `PCBDDC`

583:   Logically Collective

585:   Input Parameters:
586: + pc - the preconditioning context
587: - k  - coarsening ratio (H/h at the coarser level)

589:   Options Database Key:
590: . -pc_bddc_coarsening_ratio <int> - Set the coarsening ratio used in multi-level coarsening

592:   Level: intermediate

594:   Note:
595:   Approximately `k` subdomains at the finer level will be aggregated into a single subdomain at the coarser level

597: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetLevels()`
598: @*/
599: PetscErrorCode PCBDDCSetCoarseningRatio(PC pc, PetscInt k)
600: {
601:   PetscFunctionBegin;
604:   PetscTryMethod(pc, "PCBDDCSetCoarseningRatio_C", (PC, PetscInt), (pc, k));
605:   PetscFunctionReturn(PETSC_SUCCESS);
606: }

608: /* The following functions (PCBDDCSetUseExactDirichlet PCBDDCSetLevel) are not public */
609: static PetscErrorCode PCBDDCSetUseExactDirichlet_BDDC(PC pc, PetscBool flg)
610: {
611:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

613:   PetscFunctionBegin;
614:   pcbddc->use_exact_dirichlet_trick = flg;
615:   PetscFunctionReturn(PETSC_SUCCESS);
616: }

618: PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc, PetscBool flg)
619: {
620:   PetscFunctionBegin;
623:   PetscTryMethod(pc, "PCBDDCSetUseExactDirichlet_C", (PC, PetscBool), (pc, flg));
624:   PetscFunctionReturn(PETSC_SUCCESS);
625: }

627: static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc, PetscInt level)
628: {
629:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

631:   PetscFunctionBegin;
632:   pcbddc->current_level = level;
633:   PetscFunctionReturn(PETSC_SUCCESS);
634: }

636: PetscErrorCode PCBDDCSetLevel(PC pc, PetscInt level)
637: {
638:   PetscFunctionBegin;
641:   PetscTryMethod(pc, "PCBDDCSetLevel_C", (PC, PetscInt), (pc, level));
642:   PetscFunctionReturn(PETSC_SUCCESS);
643: }

645: static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc, PetscInt levels)
646: {
647:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

649:   PetscFunctionBegin;
650:   PetscCheck(levels < PETSC_PCBDDC_MAXLEVELS, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Maximum number of additional levels for BDDC is %d", PETSC_PCBDDC_MAXLEVELS - 1);
651:   pcbddc->max_levels = levels;
652:   PetscFunctionReturn(PETSC_SUCCESS);
653: }

655: /*@
656:   PCBDDCSetLevels - Sets the maximum number of additional levels allowed for multilevel `PCBDDC`

658:   Logically Collective

660:   Input Parameters:
661: + pc     - the preconditioning context
662: - levels - the maximum number of levels

664:   Options Database Key:
665: . -pc_bddc_levels <int> - Set maximum number of levels for multilevel

667:   Level: intermediate

669:   Note:
670:   The default value is 0, that gives the classical two-levels BDDC algorithm

672: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetCoarseningRatio()`
673: @*/
674: PetscErrorCode PCBDDCSetLevels(PC pc, PetscInt levels)
675: {
676:   PetscFunctionBegin;
679:   PetscTryMethod(pc, "PCBDDCSetLevels_C", (PC, PetscInt), (pc, levels));
680:   PetscFunctionReturn(PETSC_SUCCESS);
681: }

683: static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc, IS DirichletBoundaries)
684: {
685:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
686:   PetscBool isequal = PETSC_FALSE;

688:   PetscFunctionBegin;
689:   PetscCall(PetscObjectReference((PetscObject)DirichletBoundaries));
690:   if (pcbddc->DirichletBoundaries) PetscCall(ISEqual(DirichletBoundaries, pcbddc->DirichletBoundaries, &isequal));
691:   /* last user setting takes precedence -> destroy any other customization */
692:   PetscCall(ISDestroy(&pcbddc->DirichletBoundariesLocal));
693:   PetscCall(ISDestroy(&pcbddc->DirichletBoundaries));
694:   pcbddc->DirichletBoundaries = DirichletBoundaries;
695:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
696:   PetscFunctionReturn(PETSC_SUCCESS);
697: }

699: /*@
700:   PCBDDCSetDirichletBoundaries - Set the `IS` defining Dirichlet boundaries for the global problem.

702:   Collective

704:   Input Parameters:
705: + pc                  - the preconditioning context
706: - DirichletBoundaries - parallel `IS` defining the Dirichlet boundaries

708:   Level: intermediate

710:   Note:
711:   Provide the information if you used `MatZeroRows()` or `MatZeroRowsColumns()`. Any process can list any global node

713: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDirichletBoundariesLocal()`, `MatZeroRows()`, `MatZeroRowsColumns()`
714: @*/
715: PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc, IS DirichletBoundaries)
716: {
717:   PetscFunctionBegin;
720:   PetscCheckSameComm(pc, 1, DirichletBoundaries, 2);
721:   PetscTryMethod(pc, "PCBDDCSetDirichletBoundaries_C", (PC, IS), (pc, DirichletBoundaries));
722:   PetscFunctionReturn(PETSC_SUCCESS);
723: }

725: static PetscErrorCode PCBDDCSetDirichletBoundariesLocal_BDDC(PC pc, IS DirichletBoundaries)
726: {
727:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
728:   PetscBool isequal = PETSC_FALSE;

730:   PetscFunctionBegin;
731:   PetscCall(PetscObjectReference((PetscObject)DirichletBoundaries));
732:   if (pcbddc->DirichletBoundariesLocal) PetscCall(ISEqual(DirichletBoundaries, pcbddc->DirichletBoundariesLocal, &isequal));
733:   /* last user setting takes precedence -> destroy any other customization */
734:   PetscCall(ISDestroy(&pcbddc->DirichletBoundariesLocal));
735:   PetscCall(ISDestroy(&pcbddc->DirichletBoundaries));
736:   pcbddc->DirichletBoundariesLocal = DirichletBoundaries;
737:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
738:   PetscFunctionReturn(PETSC_SUCCESS);
739: }

741: /*@
742:   PCBDDCSetDirichletBoundariesLocal - Set the `IS` defining Dirichlet boundaries for the global problem in local ordering.

744:   Collective

746:   Input Parameters:
747: + pc                  - the preconditioning context
748: - DirichletBoundaries - parallel `IS` defining the Dirichlet boundaries (in local ordering)

750:   Level: intermediate

752: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDirichletBoundaries()`, `MatZeroRows()`, `MatZeroRowsColumns()`
753: @*/
754: PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc, IS DirichletBoundaries)
755: {
756:   PetscFunctionBegin;
759:   PetscCheckSameComm(pc, 1, DirichletBoundaries, 2);
760:   PetscTryMethod(pc, "PCBDDCSetDirichletBoundariesLocal_C", (PC, IS), (pc, DirichletBoundaries));
761:   PetscFunctionReturn(PETSC_SUCCESS);
762: }

764: static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc, IS NeumannBoundaries)
765: {
766:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
767:   PetscBool isequal = PETSC_FALSE;

769:   PetscFunctionBegin;
770:   PetscCall(PetscObjectReference((PetscObject)NeumannBoundaries));
771:   if (pcbddc->NeumannBoundaries) PetscCall(ISEqual(NeumannBoundaries, pcbddc->NeumannBoundaries, &isequal));
772:   /* last user setting takes precedence -> destroy any other customization */
773:   PetscCall(ISDestroy(&pcbddc->NeumannBoundariesLocal));
774:   PetscCall(ISDestroy(&pcbddc->NeumannBoundaries));
775:   pcbddc->NeumannBoundaries = NeumannBoundaries;
776:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
777:   PetscFunctionReturn(PETSC_SUCCESS);
778: }

780: /*@
781:   PCBDDCSetNeumannBoundaries - Set the `IS` defining Neumann boundaries for the global problem.

783:   Collective

785:   Input Parameters:
786: + pc                - the preconditioning context
787: - NeumannBoundaries - parallel `IS` defining the Neumann boundaries

789:   Level: intermediate

791:   Note:
792:   Any process can list any global node

794: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetNeumannBoundariesLocal()`
795: @*/
796: PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc, IS NeumannBoundaries)
797: {
798:   PetscFunctionBegin;
801:   PetscCheckSameComm(pc, 1, NeumannBoundaries, 2);
802:   PetscTryMethod(pc, "PCBDDCSetNeumannBoundaries_C", (PC, IS), (pc, NeumannBoundaries));
803:   PetscFunctionReturn(PETSC_SUCCESS);
804: }

806: static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc, IS NeumannBoundaries)
807: {
808:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
809:   PetscBool isequal = PETSC_FALSE;

811:   PetscFunctionBegin;
812:   PetscCall(PetscObjectReference((PetscObject)NeumannBoundaries));
813:   if (pcbddc->NeumannBoundariesLocal) PetscCall(ISEqual(NeumannBoundaries, pcbddc->NeumannBoundariesLocal, &isequal));
814:   /* last user setting takes precedence -> destroy any other customization */
815:   PetscCall(ISDestroy(&pcbddc->NeumannBoundariesLocal));
816:   PetscCall(ISDestroy(&pcbddc->NeumannBoundaries));
817:   pcbddc->NeumannBoundariesLocal = NeumannBoundaries;
818:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
819:   PetscFunctionReturn(PETSC_SUCCESS);
820: }

822: /*@
823:   PCBDDCSetNeumannBoundariesLocal - Set the `IS` defining Neumann boundaries for the global problem in local ordering.

825:   Collective

827:   Input Parameters:
828: + pc                - the preconditioning context
829: - NeumannBoundaries - parallel `IS` defining the subdomain part of Neumann boundaries (in local ordering)

831:   Level: intermediate

833: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCGetDirichletBoundaries()`
834: @*/
835: PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc, IS NeumannBoundaries)
836: {
837:   PetscFunctionBegin;
840:   PetscCheckSameComm(pc, 1, NeumannBoundaries, 2);
841:   PetscTryMethod(pc, "PCBDDCSetNeumannBoundariesLocal_C", (PC, IS), (pc, NeumannBoundaries));
842:   PetscFunctionReturn(PETSC_SUCCESS);
843: }

845: static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc, IS *DirichletBoundaries)
846: {
847:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

849:   PetscFunctionBegin;
850:   *DirichletBoundaries = pcbddc->DirichletBoundaries;
851:   PetscFunctionReturn(PETSC_SUCCESS);
852: }

854: /*@
855:   PCBDDCGetDirichletBoundaries - Get parallel `IS` for Dirichlet boundaries

857:   Collective

859:   Input Parameter:
860: . pc - the preconditioning context

862:   Output Parameter:
863: . DirichletBoundaries - index set defining the Dirichlet boundaries

865:   Level: intermediate

867:   Note:
868:   The `IS` returned (if any) is the same passed in earlier by the user with `PCBDDCSetDirichletBoundaries()`

870: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDirichletBoundaries()`
871: @*/
872: PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc, IS *DirichletBoundaries)
873: {
874:   PetscFunctionBegin;
876:   PetscUseMethod(pc, "PCBDDCGetDirichletBoundaries_C", (PC, IS *), (pc, DirichletBoundaries));
877:   PetscFunctionReturn(PETSC_SUCCESS);
878: }

880: static PetscErrorCode PCBDDCGetDirichletBoundariesLocal_BDDC(PC pc, IS *DirichletBoundaries)
881: {
882:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

884:   PetscFunctionBegin;
885:   *DirichletBoundaries = pcbddc->DirichletBoundariesLocal;
886:   PetscFunctionReturn(PETSC_SUCCESS);
887: }

889: /*@
890:   PCBDDCGetDirichletBoundariesLocal - Get parallel `IS` for Dirichlet boundaries (in local ordering)

892:   Collective

894:   Input Parameter:
895: . pc - the preconditioning context

897:   Output Parameter:
898: . DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries

900:   Level: intermediate

902:   Note:
903:   The `IS` returned could be the same passed in earlier by the user (if provided with `PCBDDCSetDirichletBoundariesLocal()`)
904:   or a global-to-local map of the global `IS` (if provided with `PCBDDCSetDirichletBoundaries()`).
905:   In the latter case, the `IS` will be available only after `PCSetUp()`.

907: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCGetDirichletBoundaries()`, `PCBDDCSetDirichletBoundaries()`
908: @*/
909: PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc, IS *DirichletBoundaries)
910: {
911:   PetscFunctionBegin;
913:   PetscUseMethod(pc, "PCBDDCGetDirichletBoundariesLocal_C", (PC, IS *), (pc, DirichletBoundaries));
914:   PetscFunctionReturn(PETSC_SUCCESS);
915: }

917: static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc, IS *NeumannBoundaries)
918: {
919:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

921:   PetscFunctionBegin;
922:   *NeumannBoundaries = pcbddc->NeumannBoundaries;
923:   PetscFunctionReturn(PETSC_SUCCESS);
924: }

926: /*@
927:   PCBDDCGetNeumannBoundaries - Get parallel `IS` for Neumann boundaries

929:   Not Collective

931:   Input Parameter:
932: . pc - the preconditioning context

934:   Output Parameter:
935: . NeumannBoundaries - index set defining the Neumann boundaries

937:   Level: intermediate

939:   Note:
940:   The `IS` returned (if any) is the same passed in earlier by the user with `PCBDDCSetNeumannBoundaries()`

942: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCGetDirichletBoundaries()`, `PCBDDCSetDirichletBoundaries()`
943: @*/
944: PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc, IS *NeumannBoundaries)
945: {
946:   PetscFunctionBegin;
948:   PetscUseMethod(pc, "PCBDDCGetNeumannBoundaries_C", (PC, IS *), (pc, NeumannBoundaries));
949:   PetscFunctionReturn(PETSC_SUCCESS);
950: }

952: static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc, IS *NeumannBoundaries)
953: {
954:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

956:   PetscFunctionBegin;
957:   *NeumannBoundaries = pcbddc->NeumannBoundariesLocal;
958:   PetscFunctionReturn(PETSC_SUCCESS);
959: }

961: /*@
962:   PCBDDCGetNeumannBoundariesLocal - Get parallel `IS` for Neumann boundaries (in local ordering)

964:   Not Collective

966:   Input Parameter:
967: . pc - the preconditioning context

969:   Output Parameter:
970: . NeumannBoundaries - index set defining the subdomain part of Neumann boundaries

972:   Level: intermediate

974:   Note:
975:   The `IS` returned could be the same passed in earlier by the user (if provided with `PCBDDCSetNeumannBoundariesLocal()`
976:   or a global-to-local map of the global `IS` (if provided with `PCBDDCSetNeumannBoundaries()`).
977:   In the latter case, the `IS` will be available after `PCSetUp()`.

979: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCSetNeumannBoundariesLocal)`, `PCBDDCGetNeumannBoundaries()`
980: @*/
981: PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc, IS *NeumannBoundaries)
982: {
983:   PetscFunctionBegin;
985:   PetscUseMethod(pc, "PCBDDCGetNeumannBoundariesLocal_C", (PC, IS *), (pc, NeumannBoundaries));
986:   PetscFunctionReturn(PETSC_SUCCESS);
987: }

989: static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs, const PetscInt xadj[], const PetscInt adjncy[], PetscCopyMode copymode)
990: {
991:   PC_BDDC    *pcbddc    = (PC_BDDC *)pc->data;
992:   PCBDDCGraph mat_graph = pcbddc->mat_graph;
993:   PetscBool   same_data = PETSC_FALSE;

995:   PetscFunctionBegin;
996:   if (!nvtxs) {
997:     if (copymode == PETSC_OWN_POINTER) {
998:       PetscCall(PetscFree(xadj));
999:       PetscCall(PetscFree(adjncy));
1000:     }
1001:     PetscCall(PCBDDCGraphResetCSR(mat_graph));
1002:     PetscFunctionReturn(PETSC_SUCCESS);
1003:   }
1004:   if (mat_graph->nvtxs == nvtxs && mat_graph->freecsr) { /* we own the data */
1005:     if (mat_graph->xadj == xadj && mat_graph->adjncy == adjncy) same_data = PETSC_TRUE;
1006:     if (!same_data && mat_graph->xadj[nvtxs] == xadj[nvtxs]) {
1007:       PetscCall(PetscArraycmp(xadj, mat_graph->xadj, nvtxs + 1, &same_data));
1008:       if (same_data) PetscCall(PetscArraycmp(adjncy, mat_graph->adjncy, xadj[nvtxs], &same_data));
1009:     }
1010:   }
1011:   if (!same_data) {
1012:     /* free old CSR */
1013:     PetscCall(PCBDDCGraphResetCSR(mat_graph));
1014:     /* get CSR into graph structure */
1015:     if (copymode == PETSC_COPY_VALUES) {
1016:       PetscCall(PetscMalloc1(nvtxs + 1, &mat_graph->xadj));
1017:       PetscCall(PetscMalloc1(xadj[nvtxs], &mat_graph->adjncy));
1018:       PetscCall(PetscArraycpy(mat_graph->xadj, xadj, nvtxs + 1));
1019:       PetscCall(PetscArraycpy(mat_graph->adjncy, adjncy, xadj[nvtxs]));
1020:       mat_graph->freecsr = PETSC_TRUE;
1021:     } else if (copymode == PETSC_OWN_POINTER) {
1022:       mat_graph->xadj    = (PetscInt *)xadj;
1023:       mat_graph->adjncy  = (PetscInt *)adjncy;
1024:       mat_graph->freecsr = PETSC_TRUE;
1025:     } else if (copymode == PETSC_USE_POINTER) {
1026:       mat_graph->xadj    = (PetscInt *)xadj;
1027:       mat_graph->adjncy  = (PetscInt *)adjncy;
1028:       mat_graph->freecsr = PETSC_FALSE;
1029:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported copy mode %d", copymode);
1030:     mat_graph->nvtxs_csr         = nvtxs;
1031:     pcbddc->recompute_topography = PETSC_TRUE;
1032:   }
1033:   PetscFunctionReturn(PETSC_SUCCESS);
1034: }

1036: /*@
1037:   PCBDDCSetLocalAdjacencyGraph - Set adjacency structure (CSR graph) of the local degrees of freedom.

1039:   Not collective

1041:   Input Parameters:
1042: + pc       - the preconditioning context.
1043: . nvtxs    - number of local vertices of the graph (i.e., the number of local dofs).
1044: . xadj     - CSR format row pointers for the connectivity of the dofs
1045: . adjncy   - CSR format column pointers for the connectivity of the dofs
1046: - copymode - supported modes are `PETSC_COPY_VALUES`, `PETSC_USE_POINTER` or `PETSC_OWN_POINTER`.

1048:   Level: intermediate

1050:   Note:
1051:   A dof is considered connected with all local dofs if xadj[dof+1]-xadj[dof] == 1 and adjncy[xadj[dof]] is negative.

1053: .seealso: [](ch_ksp), `PCBDDC`, `PetscCopyMode`
1054: @*/
1055: PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc, PetscInt nvtxs, const PetscInt xadj[], const PetscInt adjncy[], PetscCopyMode copymode)
1056: {
1057:   PetscBool f = PETSC_FALSE;

1059:   PetscFunctionBegin;
1061:   if (nvtxs) {
1062:     PetscAssertPointer(xadj, 3);
1063:     if (xadj[nvtxs]) PetscAssertPointer(adjncy, 4);
1064:   }
1065:   PetscTryMethod(pc, "PCBDDCSetLocalAdjacencyGraph_C", (PC, PetscInt, const PetscInt[], const PetscInt[], PetscCopyMode), (pc, nvtxs, xadj, adjncy, copymode));
1066:   /* free arrays if PCBDDC is not the PC type */
1067:   PetscCall(PetscObjectHasFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", &f));
1068:   if (!f && copymode == PETSC_OWN_POINTER) {
1069:     PetscCall(PetscFree(xadj));
1070:     PetscCall(PetscFree(adjncy));
1071:   }
1072:   PetscFunctionReturn(PETSC_SUCCESS);
1073: }

1075: static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc, PetscInt n_is, IS ISForDofs[])
1076: {
1077:   PC_BDDC  *pcbddc = (PC_BDDC *)pc->data;
1078:   PetscInt  i;
1079:   PetscBool isequal = PETSC_FALSE;

1081:   PetscFunctionBegin;
1082:   if (pcbddc->n_ISForDofsLocal == n_is) {
1083:     for (i = 0; i < n_is; i++) {
1084:       PetscBool isequalt;
1085:       PetscCall(ISEqual(ISForDofs[i], pcbddc->ISForDofsLocal[i], &isequalt));
1086:       if (!isequalt) break;
1087:     }
1088:     if (i == n_is) isequal = PETSC_TRUE;
1089:   }
1090:   for (i = 0; i < n_is; i++) PetscCall(PetscObjectReference((PetscObject)ISForDofs[i]));
1091:   /* Destroy ISes if they were already set */
1092:   for (i = 0; i < pcbddc->n_ISForDofsLocal; i++) PetscCall(ISDestroy(&pcbddc->ISForDofsLocal[i]));
1093:   PetscCall(PetscFree(pcbddc->ISForDofsLocal));
1094:   /* last user setting takes precedence -> destroy any other customization */
1095:   for (i = 0; i < pcbddc->n_ISForDofs; i++) PetscCall(ISDestroy(&pcbddc->ISForDofs[i]));
1096:   PetscCall(PetscFree(pcbddc->ISForDofs));
1097:   pcbddc->n_ISForDofs = 0;
1098:   /* allocate space then set */
1099:   if (n_is) PetscCall(PetscMalloc1(n_is, &pcbddc->ISForDofsLocal));
1100:   for (i = 0; i < n_is; i++) pcbddc->ISForDofsLocal[i] = ISForDofs[i];
1101:   pcbddc->n_ISForDofsLocal = n_is;
1102:   if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
1103:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1104:   PetscFunctionReturn(PETSC_SUCCESS);
1105: }

1107: /*@
1108:   PCBDDCSetDofsSplittingLocal - Set the `IS` defining fields of the local subdomain matrix

1110:   Collective

1112:   Input Parameters:
1113: + pc        - the preconditioning context
1114: . n_is      - number of index sets defining the fields, must be the same on all MPI processes
1115: - ISForDofs - array of `IS` describing the fields in local ordering

1117:   Level: intermediate

1119:   Note:
1120:   Not all nodes need to be listed, unlisted nodes will belong to the complement field.

1122: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDofsSplitting()`
1123: @*/
1124: PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc, PetscInt n_is, IS ISForDofs[])
1125: {
1126:   PetscInt i;

1128:   PetscFunctionBegin;
1131:   for (i = 0; i < n_is; i++) {
1132:     PetscCheckSameComm(pc, 1, ISForDofs[i], 3);
1134:   }
1135:   PetscTryMethod(pc, "PCBDDCSetDofsSplittingLocal_C", (PC, PetscInt, IS[]), (pc, n_is, ISForDofs));
1136:   PetscFunctionReturn(PETSC_SUCCESS);
1137: }

1139: static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc, PetscInt n_is, IS ISForDofs[])
1140: {
1141:   PC_BDDC  *pcbddc = (PC_BDDC *)pc->data;
1142:   PetscInt  i;
1143:   PetscBool isequal = PETSC_FALSE;

1145:   PetscFunctionBegin;
1146:   if (pcbddc->n_ISForDofs == n_is) {
1147:     for (i = 0; i < n_is; i++) {
1148:       PetscBool isequalt;
1149:       PetscCall(ISEqual(ISForDofs[i], pcbddc->ISForDofs[i], &isequalt));
1150:       if (!isequalt) break;
1151:     }
1152:     if (i == n_is) isequal = PETSC_TRUE;
1153:   }
1154:   for (i = 0; i < n_is; i++) PetscCall(PetscObjectReference((PetscObject)ISForDofs[i]));
1155:   /* Destroy ISes if they were already set */
1156:   for (i = 0; i < pcbddc->n_ISForDofs; i++) PetscCall(ISDestroy(&pcbddc->ISForDofs[i]));
1157:   PetscCall(PetscFree(pcbddc->ISForDofs));
1158:   /* last user setting takes precedence -> destroy any other customization */
1159:   for (i = 0; i < pcbddc->n_ISForDofsLocal; i++) PetscCall(ISDestroy(&pcbddc->ISForDofsLocal[i]));
1160:   PetscCall(PetscFree(pcbddc->ISForDofsLocal));
1161:   pcbddc->n_ISForDofsLocal = 0;
1162:   /* allocate space then set */
1163:   if (n_is) PetscCall(PetscMalloc1(n_is, &pcbddc->ISForDofs));
1164:   for (i = 0; i < n_is; i++) pcbddc->ISForDofs[i] = ISForDofs[i];
1165:   pcbddc->n_ISForDofs = n_is;
1166:   if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
1167:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1168:   PetscFunctionReturn(PETSC_SUCCESS);
1169: }

1171: /*@
1172:   PCBDDCSetDofsSplitting - Set the `IS` defining fields of the global matrix

1174:   Collective

1176:   Input Parameters:
1177: + pc        - the preconditioning context
1178: . n_is      - number of index sets defining the fields
1179: - ISForDofs - array of `IS` describing the fields in global ordering

1181:   Level: intermediate

1183:   Note:
1184:   Any process can list any global node. Not all nodes need to be listed: unlisted nodes will belong to the complement field.

1186: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDofsSplittingLocal()`
1187: @*/
1188: PetscErrorCode PCBDDCSetDofsSplitting(PC pc, PetscInt n_is, IS ISForDofs[])
1189: {
1190:   PetscInt i;

1192:   PetscFunctionBegin;
1195:   for (i = 0; i < n_is; i++) {
1197:     PetscCheckSameComm(pc, 1, ISForDofs[i], 3);
1198:   }
1199:   PetscTryMethod(pc, "PCBDDCSetDofsSplitting_C", (PC, PetscInt, IS[]), (pc, n_is, ISForDofs));
1200:   PetscFunctionReturn(PETSC_SUCCESS);
1201: }

1203: static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
1204: {
1205:   PC_BDDC  *pcbddc = (PC_BDDC *)pc->data;
1206:   PC_IS    *pcis   = (PC_IS *)pc->data;
1207:   Vec       used_vec;
1208:   PetscBool iscg, save_rhs = PETSC_TRUE, benign_correction_computed;

1210:   PetscFunctionBegin;
1211:   /* if we are working with CG, one dirichlet solve can be avoided during Krylov iterations */
1212:   if (ksp) {
1213:     PetscCall(PetscObjectTypeCompareAny((PetscObject)ksp, &iscg, KSPCG, KSPGROPPCG, KSPPIPECG, KSPPIPELCG, KSPPIPECGRR, ""));
1214:     if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || !iscg || pc->mat != pc->pmat) PetscCall(PCBDDCSetUseExactDirichlet(pc, PETSC_FALSE));
1215:   }
1216:   if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || pc->mat != pc->pmat) PetscCall(PCBDDCSetUseExactDirichlet(pc, PETSC_FALSE));

1218:   /* Creates parallel work vectors used in presolve */
1219:   if (!pcbddc->original_rhs) PetscCall(VecDuplicate(pcis->vec1_global, &pcbddc->original_rhs));
1220:   if (!pcbddc->temp_solution) PetscCall(VecDuplicate(pcis->vec1_global, &pcbddc->temp_solution));

1222:   pcbddc->temp_solution_used = PETSC_FALSE;
1223:   if (x) {
1224:     PetscCall(PetscObjectReference((PetscObject)x));
1225:     used_vec = x;
1226:   } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */
1227:     PetscCall(PetscObjectReference((PetscObject)pcbddc->temp_solution));
1228:     used_vec = pcbddc->temp_solution;
1229:     PetscCall(VecSet(used_vec, 0.0));
1230:     pcbddc->temp_solution_used = PETSC_TRUE;
1231:     PetscCall(VecCopy(rhs, pcbddc->original_rhs));
1232:     save_rhs                  = PETSC_FALSE;
1233:     pcbddc->eliminate_dirdofs = PETSC_TRUE;
1234:   }

1236:   /* hack into ksp data structure since PCPreSolve comes earlier than setting to zero the guess in src/ksp/ksp/interface/itfunc.c */
1237:   if (ksp) {
1238:     /* store the flag for the initial guess since it will be restored back during PCPostSolve_BDDC */
1239:     PetscCall(KSPGetInitialGuessNonzero(ksp, &pcbddc->ksp_guess_nonzero));
1240:     if (!pcbddc->ksp_guess_nonzero) PetscCall(VecSet(used_vec, 0.0));
1241:   }

1243:   pcbddc->rhs_change = PETSC_FALSE;
1244:   /* Take into account zeroed rows -> change rhs and store solution removed */
1245:   if (rhs && pcbddc->eliminate_dirdofs) {
1246:     IS dirIS = NULL;

1248:     /* DirichletBoundariesLocal may not be consistent among neighbours; gets a dirichlet dofs IS from graph (may be cached) */
1249:     PetscCall(PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph, &dirIS));
1250:     if (dirIS) {
1251:       Mat_IS            *matis = (Mat_IS *)pc->pmat->data;
1252:       PetscInt           dirsize, i, *is_indices;
1253:       PetscScalar       *array_x;
1254:       const PetscScalar *array_diagonal;

1256:       PetscCall(MatGetDiagonal(pc->pmat, pcis->vec1_global));
1257:       PetscCall(VecPointwiseDivide(pcis->vec1_global, rhs, pcis->vec1_global));
1258:       PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_global, pcis->vec2_N, INSERT_VALUES, SCATTER_FORWARD));
1259:       PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_global, pcis->vec2_N, INSERT_VALUES, SCATTER_FORWARD));
1260:       PetscCall(VecScatterBegin(matis->rctx, used_vec, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD));
1261:       PetscCall(VecScatterEnd(matis->rctx, used_vec, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD));
1262:       PetscCall(ISGetLocalSize(dirIS, &dirsize));
1263:       PetscCall(VecGetArray(pcis->vec1_N, &array_x));
1264:       PetscCall(VecGetArrayRead(pcis->vec2_N, &array_diagonal));
1265:       PetscCall(ISGetIndices(dirIS, (const PetscInt **)&is_indices));
1266:       for (i = 0; i < dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]];
1267:       PetscCall(ISRestoreIndices(dirIS, (const PetscInt **)&is_indices));
1268:       PetscCall(VecRestoreArrayRead(pcis->vec2_N, &array_diagonal));
1269:       PetscCall(VecRestoreArray(pcis->vec1_N, &array_x));
1270:       PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_N, used_vec, INSERT_VALUES, SCATTER_REVERSE));
1271:       PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_N, used_vec, INSERT_VALUES, SCATTER_REVERSE));
1272:       pcbddc->rhs_change = PETSC_TRUE;
1273:       PetscCall(ISDestroy(&dirIS));
1274:     }
1275:   }

1277:   /* remove the computed solution or the initial guess from the rhs */
1278:   if (pcbddc->rhs_change || (ksp && pcbddc->ksp_guess_nonzero)) {
1279:     /* save the original rhs */
1280:     if (save_rhs) {
1281:       PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1282:       save_rhs = PETSC_FALSE;
1283:     }
1284:     pcbddc->rhs_change = PETSC_TRUE;
1285:     PetscCall(VecScale(used_vec, -1.0));
1286:     PetscCall(MatMultAdd(pc->mat, used_vec, pcbddc->original_rhs, rhs));
1287:     PetscCall(VecScale(used_vec, -1.0));
1288:     PetscCall(VecCopy(used_vec, pcbddc->temp_solution));
1289:     pcbddc->temp_solution_used = PETSC_TRUE;
1290:     if (ksp) PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_FALSE));
1291:   }
1292:   PetscCall(VecDestroy(&used_vec));

1294:   /* compute initial vector in benign space if needed
1295:      and remove non-benign solution from the rhs */
1296:   benign_correction_computed = PETSC_FALSE;
1297:   if (rhs && pcbddc->benign_compute_correction && (pcbddc->benign_have_null || pcbddc->benign_apply_coarse_only)) {
1298:     /* compute u^*_h using ideas similar to those in Xuemin Tu's PhD thesis (see Section 4.8.1)
1299:        Recursively apply BDDC in the multilevel case */
1300:     if (!pcbddc->benign_vec) PetscCall(VecDuplicate(rhs, &pcbddc->benign_vec));
1301:     /* keep applying coarse solver unless we no longer have benign subdomains */
1302:     pcbddc->benign_apply_coarse_only = pcbddc->benign_have_null ? PETSC_TRUE : PETSC_FALSE;
1303:     if (!pcbddc->benign_skip_correction) {
1304:       PetscCall(PCApply_BDDC(pc, rhs, pcbddc->benign_vec));
1305:       benign_correction_computed = PETSC_TRUE;
1306:       if (pcbddc->temp_solution_used) PetscCall(VecAXPY(pcbddc->temp_solution, 1.0, pcbddc->benign_vec));
1307:       PetscCall(VecScale(pcbddc->benign_vec, -1.0));
1308:       /* store the original rhs if not done earlier */
1309:       if (save_rhs) PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1310:       if (pcbddc->rhs_change) {
1311:         PetscCall(MatMultAdd(pc->mat, pcbddc->benign_vec, rhs, rhs));
1312:       } else {
1313:         PetscCall(MatMultAdd(pc->mat, pcbddc->benign_vec, pcbddc->original_rhs, rhs));
1314:       }
1315:       pcbddc->rhs_change = PETSC_TRUE;
1316:     }
1317:     pcbddc->benign_apply_coarse_only = PETSC_FALSE;
1318:   } else {
1319:     PetscCall(VecDestroy(&pcbddc->benign_vec));
1320:   }

1322:   /* dbg output */
1323:   if (pcbddc->dbg_flag && benign_correction_computed) {
1324:     Vec v;

1326:     PetscCall(VecDuplicate(pcis->vec1_global, &v));
1327:     if (pcbddc->ChangeOfBasisMatrix) {
1328:       PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, rhs, v));
1329:     } else {
1330:       PetscCall(VecCopy(rhs, v));
1331:     }
1332:     PetscCall(PCBDDCBenignGetOrSetP0(pc, v, PETSC_TRUE));
1333:     PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "LEVEL %" PetscInt_FMT ": is the correction benign?\n", pcbddc->current_level));
1334:     PetscCall(PetscScalarView(pcbddc->benign_n, pcbddc->benign_p0, pcbddc->dbg_viewer));
1335:     PetscCall(PetscViewerFlush(pcbddc->dbg_viewer));
1336:     PetscCall(VecDestroy(&v));
1337:   }

1339:   /* set initial guess if using PCG */
1340:   pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1341:   if (x && pcbddc->use_exact_dirichlet_trick) {
1342:     PetscCall(VecSet(x, 0.0));
1343:     if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) {
1344:       if (benign_correction_computed) { /* we have already saved the changed rhs */
1345:         PetscCall(VecLockReadPop(pcis->vec1_global));
1346:       } else {
1347:         PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, rhs, pcis->vec1_global));
1348:       }
1349:       PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_global, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1350:       PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_global, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1351:     } else {
1352:       PetscCall(VecScatterBegin(pcis->global_to_D, rhs, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1353:       PetscCall(VecScatterEnd(pcis->global_to_D, rhs, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1354:     }
1355:     PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1356:     PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
1357:     PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1358:     PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
1359:     if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) {
1360:       PetscCall(VecSet(pcis->vec1_global, 0.));
1361:       PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, pcis->vec1_global, INSERT_VALUES, SCATTER_REVERSE));
1362:       PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, pcis->vec1_global, INSERT_VALUES, SCATTER_REVERSE));
1363:       PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcis->vec1_global, x));
1364:     } else {
1365:       PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, x, INSERT_VALUES, SCATTER_REVERSE));
1366:       PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, x, INSERT_VALUES, SCATTER_REVERSE));
1367:     }
1368:     if (ksp) PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
1369:     pcbddc->exact_dirichlet_trick_app = PETSC_TRUE;
1370:   } else if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior && benign_correction_computed && pcbddc->use_exact_dirichlet_trick) {
1371:     PetscCall(VecLockReadPop(pcis->vec1_global));
1372:   }
1373:   PetscFunctionReturn(PETSC_SUCCESS);
1374: }

1376: static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
1377: {
1378:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

1380:   PetscFunctionBegin;
1381:   /* add solution removed in presolve */
1382:   if (x && pcbddc->rhs_change) {
1383:     if (pcbddc->temp_solution_used) {
1384:       PetscCall(VecAXPY(x, 1.0, pcbddc->temp_solution));
1385:     } else if (pcbddc->benign_compute_correction && pcbddc->benign_vec) {
1386:       PetscCall(VecAXPY(x, -1.0, pcbddc->benign_vec));
1387:     }
1388:     /* restore to original state (not for FETI-DP) */
1389:     if (ksp) pcbddc->temp_solution_used = PETSC_FALSE;
1390:   }

1392:   /* restore rhs to its original state (not needed for FETI-DP) */
1393:   if (rhs && pcbddc->rhs_change) {
1394:     PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1395:     pcbddc->rhs_change = PETSC_FALSE;
1396:   }
1397:   /* restore ksp guess state */
1398:   if (ksp) {
1399:     PetscCall(KSPSetInitialGuessNonzero(ksp, pcbddc->ksp_guess_nonzero));
1400:     /* reset flag for exact dirichlet trick */
1401:     pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1402:   }
1403:   PetscFunctionReturn(PETSC_SUCCESS);
1404: }

1406: static PetscErrorCode PCSetUp_BDDC(PC pc)
1407: {
1408:   PC_BDDC        *pcbddc = (PC_BDDC *)pc->data;
1409:   PCBDDCSubSchurs sub_schurs;
1410:   Mat_IS         *matis;
1411:   MatNullSpace    nearnullspace;
1412:   Mat             lA;
1413:   IS              lP, zerodiag = NULL;
1414:   PetscInt        nrows, ncols;
1415:   PetscMPIInt     size;
1416:   PetscBool       computesubschurs;
1417:   PetscBool       computeconstraintsmatrix;
1418:   PetscBool       new_nearnullspace_provided, ismatis, rl;
1419:   PetscBool       isset, issym, isspd;

1421:   PetscFunctionBegin;
1422:   PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATIS, &ismatis));
1423:   PetscCheck(ismatis, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PCBDDC preconditioner requires matrix of type MATIS");
1424:   PetscCall(MatGetSize(pc->pmat, &nrows, &ncols));
1425:   PetscCheck(nrows == ncols, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCBDDC preconditioner requires a square preconditioning matrix");
1426:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));

1428:   matis = (Mat_IS *)pc->pmat->data;
1429:   /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */
1430:   /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup
1431:      Also, BDDC builds its own KSP for the Dirichlet problem */
1432:   rl = pcbddc->recompute_topography;
1433:   if (!pc->setupcalled || pc->flag == DIFFERENT_NONZERO_PATTERN) rl = PETSC_TRUE;
1434:   PetscCallMPI(MPIU_Allreduce(&rl, &pcbddc->recompute_topography, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)pc)));
1435:   if (pcbddc->recompute_topography) {
1436:     pcbddc->graphanalyzed    = PETSC_FALSE;
1437:     computeconstraintsmatrix = PETSC_TRUE;
1438:   } else {
1439:     computeconstraintsmatrix = PETSC_FALSE;
1440:   }

1442:   /* check parameters' compatibility */
1443:   if (!pcbddc->use_deluxe_scaling) pcbddc->deluxe_zerorows = PETSC_FALSE;
1444:   pcbddc->adaptive_selection   = (PetscBool)(pcbddc->adaptive_threshold[0] != 0.0 || pcbddc->adaptive_threshold[1] != 0.0);
1445:   pcbddc->use_deluxe_scaling   = (PetscBool)(pcbddc->use_deluxe_scaling && size > 1);
1446:   pcbddc->adaptive_selection   = (PetscBool)(pcbddc->adaptive_selection && size > 1);
1447:   pcbddc->adaptive_userdefined = (PetscBool)(pcbddc->adaptive_selection && pcbddc->adaptive_userdefined);
1448:   if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE;

1450:   computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling);

1452:   /* activate all connected components if the netflux has been requested */
1453:   if (pcbddc->compute_nonetflux) {
1454:     pcbddc->use_vertices = PETSC_TRUE;
1455:     pcbddc->use_edges    = PETSC_TRUE;
1456:     pcbddc->use_faces    = PETSC_TRUE;
1457:   }

1459:   /* Get stdout for dbg */
1460:   if (pcbddc->dbg_flag) {
1461:     if (!pcbddc->dbg_viewer) pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc));
1462:     PetscCall(PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer));
1463:     PetscCall(PetscViewerASCIIAddTab(pcbddc->dbg_viewer, 2 * pcbddc->current_level));
1464:   }

1466:   /* process topology information */
1467:   PetscCall(PetscLogEventBegin(PC_BDDC_Topology[pcbddc->current_level], pc, 0, 0, 0));
1468:   if (pcbddc->recompute_topography) {
1469:     PetscCall(PCBDDCComputeLocalTopologyInfo(pc));
1470:     if (pcbddc->discretegradient) PetscCall(PCBDDCNedelecSupport(pc));
1471:   }
1472:   if (pcbddc->corner_selected) pcbddc->use_vertices = PETSC_TRUE;

1474:   /* change basis if requested by the user */
1475:   if (pcbddc->user_ChangeOfBasisMatrix) {
1476:     /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */
1477:     pcbddc->use_change_of_basis = PETSC_FALSE;
1478:     PetscCall(PCBDDCComputeLocalMatrix(pc, pcbddc->user_ChangeOfBasisMatrix));
1479:   } else {
1480:     PetscCall(MatDestroy(&pcbddc->local_mat));
1481:     PetscCall(PetscObjectReference((PetscObject)matis->A));
1482:     pcbddc->local_mat = matis->A;
1483:   }

1485:   /*
1486:      Compute change of basis on local pressures (aka zerodiag dofs) with the benign trick
1487:      This should come earlier then PCISSetUp for extracting the correct subdomain matrices
1488:   */
1489:   PetscCall(PCBDDCBenignShellMat(pc, PETSC_TRUE));
1490:   if (pcbddc->benign_saddle_point) {
1491:     PC_IS *pcis = (PC_IS *)pc->data;

1493:     if (pcbddc->user_ChangeOfBasisMatrix || pcbddc->use_change_of_basis || !computesubschurs) pcbddc->benign_change_explicit = PETSC_TRUE;
1494:     /* detect local saddle point and change the basis in pcbddc->local_mat */
1495:     PetscCall(PCBDDCBenignDetectSaddlePoint(pc, (PetscBool)(!pcbddc->recompute_topography), &zerodiag));
1496:     /* pop B0 mat from local mat */
1497:     PetscCall(PCBDDCBenignPopOrPushB0(pc, PETSC_TRUE));
1498:     /* give pcis a hint to not reuse submatrices during PCISCreate */
1499:     if (pc->flag == SAME_NONZERO_PATTERN && pcis->reusesubmatrices == PETSC_TRUE) {
1500:       if (pcbddc->benign_n && (pcbddc->benign_change_explicit || pcbddc->dbg_flag)) {
1501:         pcis->reusesubmatrices = PETSC_FALSE;
1502:       } else {
1503:         pcis->reusesubmatrices = PETSC_TRUE;
1504:       }
1505:     } else {
1506:       pcis->reusesubmatrices = PETSC_FALSE;
1507:     }
1508:   }

1510:   /* propagate relevant information */
1511:   PetscCall(MatIsSymmetricKnown(matis->A, &isset, &issym));
1512:   if (isset) PetscCall(MatSetOption(pcbddc->local_mat, MAT_SYMMETRIC, issym));
1513:   PetscCall(MatIsSPDKnown(matis->A, &isset, &isspd));
1514:   if (isset) PetscCall(MatSetOption(pcbddc->local_mat, MAT_SPD, isspd));

1516:   /* Set up all the "iterative substructuring" common block without computing solvers */
1517:   {
1518:     Mat temp_mat;

1520:     temp_mat = matis->A;
1521:     matis->A = pcbddc->local_mat;
1522:     PetscCall(PCISSetUp(pc, PETSC_TRUE, PETSC_FALSE));
1523:     pcbddc->local_mat = matis->A;
1524:     matis->A          = temp_mat;
1525:   }

1527:   /* Analyze interface */
1528:   if (!pcbddc->graphanalyzed) {
1529:     PetscCall(PCBDDCAnalyzeInterface(pc));
1530:     computeconstraintsmatrix = PETSC_TRUE;
1531:     PetscCheck(!(pcbddc->adaptive_selection && !pcbddc->use_deluxe_scaling && !pcbddc->mat_graph->twodim), PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot compute the adaptive primal space for a problem with 3D edges without deluxe scaling");
1532:     if (pcbddc->compute_nonetflux) {
1533:       MatNullSpace nnfnnsp;

1535:       PetscCheck(pcbddc->divudotp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Missing divudotp operator");
1536:       PetscCall(PCBDDCComputeNoNetFlux(pc->pmat, pcbddc->divudotp, pcbddc->divudotp_trans, pcbddc->divudotp_vl2l, pcbddc->mat_graph, &nnfnnsp));
1537:       /* TODO what if a nearnullspace is already attached? */
1538:       if (nnfnnsp) {
1539:         PetscCall(MatSetNearNullSpace(pc->pmat, nnfnnsp));
1540:         PetscCall(MatNullSpaceDestroy(&nnfnnsp));
1541:       }
1542:     }
1543:   }
1544:   PetscCall(PetscLogEventEnd(PC_BDDC_Topology[pcbddc->current_level], pc, 0, 0, 0));

1546:   /* check existence of a divergence free extension, i.e.
1547:      b(v_I,p_0) = 0 for all v_I (raise error if not).
1548:      Also, check that PCBDDCBenignGetOrSetP0 works */
1549:   if (pcbddc->benign_saddle_point && pcbddc->dbg_flag > 1) PetscCall(PCBDDCBenignCheck(pc, zerodiag));
1550:   PetscCall(ISDestroy(&zerodiag));

1552:   /* Setup local dirichlet solver ksp_D and sub_schurs solvers */
1553:   if (computesubschurs && pcbddc->recompute_topography) PetscCall(PCBDDCInitSubSchurs(pc));
1554:   /* SetUp Scaling operator (scaling matrices could be needed in SubSchursSetUp)*/
1555:   if (!pcbddc->use_deluxe_scaling) PetscCall(PCBDDCScalingSetUp(pc));

1557:   /* finish setup solvers and do adaptive selection of constraints */
1558:   sub_schurs = pcbddc->sub_schurs;
1559:   if (sub_schurs && sub_schurs->schur_explicit) {
1560:     if (computesubschurs) PetscCall(PCBDDCSetUpSubSchurs(pc));
1561:     PetscCall(PCBDDCSetUpLocalSolvers(pc, PETSC_TRUE, PETSC_FALSE));
1562:   } else {
1563:     PetscCall(PCBDDCSetUpLocalSolvers(pc, PETSC_TRUE, PETSC_FALSE));
1564:     if (computesubschurs) PetscCall(PCBDDCSetUpSubSchurs(pc));
1565:   }
1566:   if (pcbddc->adaptive_selection) {
1567:     PetscCall(PCBDDCAdaptiveSelection(pc));
1568:     computeconstraintsmatrix = PETSC_TRUE;
1569:   }

1571:   /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */
1572:   new_nearnullspace_provided = PETSC_FALSE;
1573:   PetscCall(MatGetNearNullSpace(pc->pmat, &nearnullspace));
1574:   if (pcbddc->onearnullspace) { /* already used nearnullspace */
1575:     if (!nearnullspace) {       /* near null space attached to mat has been destroyed */
1576:       new_nearnullspace_provided = PETSC_TRUE;
1577:     } else {
1578:       /* determine if the two nullspaces are different (should be lightweight) */
1579:       if (nearnullspace != pcbddc->onearnullspace) {
1580:         new_nearnullspace_provided = PETSC_TRUE;
1581:       } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */
1582:         PetscInt         i;
1583:         const Vec       *nearnullvecs;
1584:         PetscObjectState state;
1585:         PetscInt         nnsp_size;
1586:         PetscCall(MatNullSpaceGetVecs(nearnullspace, NULL, &nnsp_size, &nearnullvecs));
1587:         for (i = 0; i < nnsp_size; i++) {
1588:           PetscCall(PetscObjectStateGet((PetscObject)nearnullvecs[i], &state));
1589:           if (pcbddc->onearnullvecs_state[i] != state) {
1590:             new_nearnullspace_provided = PETSC_TRUE;
1591:             break;
1592:           }
1593:         }
1594:       }
1595:     }
1596:   } else {
1597:     if (!nearnullspace) { /* both nearnullspaces are null */
1598:       new_nearnullspace_provided = PETSC_FALSE;
1599:     } else { /* nearnullspace attached later */
1600:       new_nearnullspace_provided = PETSC_TRUE;
1601:     }
1602:   }

1604:   /* Setup constraints and related work vectors */
1605:   /* reset primal space flags */
1606:   PetscCall(PetscLogEventBegin(PC_BDDC_LocalWork[pcbddc->current_level], pc, 0, 0, 0));
1607:   pcbddc->new_primal_space       = PETSC_FALSE;
1608:   pcbddc->new_primal_space_local = PETSC_FALSE;
1609:   if (computeconstraintsmatrix || new_nearnullspace_provided) {
1610:     /* It also sets the primal space flags */
1611:     PetscCall(PCBDDCConstraintsSetUp(pc));
1612:   }
1613:   /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */
1614:   PetscCall(PCBDDCSetUpLocalWorkVectors(pc));

1616:   if (pcbddc->use_change_of_basis) {
1617:     PC_IS *pcis = (PC_IS *)pc->data;

1619:     PetscCall(PCBDDCComputeLocalMatrix(pc, pcbddc->ChangeOfBasisMatrix));
1620:     if (pcbddc->benign_change) {
1621:       PetscCall(MatDestroy(&pcbddc->benign_B0));
1622:       /* pop B0 from pcbddc->local_mat */
1623:       PetscCall(PCBDDCBenignPopOrPushB0(pc, PETSC_TRUE));
1624:     }
1625:     /* get submatrices */
1626:     PetscCall(MatDestroy(&pcis->A_IB));
1627:     PetscCall(MatDestroy(&pcis->A_BI));
1628:     PetscCall(MatDestroy(&pcis->A_BB));
1629:     PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_B_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &pcis->A_BB));
1630:     PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_I_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &pcis->A_IB));
1631:     PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_B_local, pcis->is_I_local, MAT_INITIAL_MATRIX, &pcis->A_BI));
1632:     /* set flag in pcis to not reuse submatrices during PCISCreate */
1633:     pcis->reusesubmatrices = PETSC_FALSE;
1634:   } else if (!pcbddc->user_ChangeOfBasisMatrix && !pcbddc->benign_change) {
1635:     PetscCall(MatDestroy(&pcbddc->local_mat));
1636:     PetscCall(PetscObjectReference((PetscObject)matis->A));
1637:     pcbddc->local_mat = matis->A;
1638:   }

1640:   /* interface pressure block row for B_C */
1641:   PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_lP", (PetscObject *)&lP));
1642:   PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_lA", (PetscObject *)&lA));
1643:   if (lA && lP) {
1644:     PC_IS    *pcis = (PC_IS *)pc->data;
1645:     Mat       B_BI, B_BB, Bt_BI, Bt_BB;
1646:     PetscBool issym;

1648:     PetscCall(MatIsSymmetric(lA, PETSC_SMALL, &issym));
1649:     if (issym) {
1650:       PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_I_local, MAT_INITIAL_MATRIX, &B_BI));
1651:       PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_B_local, MAT_INITIAL_MATRIX, &B_BB));
1652:       PetscCall(MatCreateTranspose(B_BI, &Bt_BI));
1653:       PetscCall(MatCreateTranspose(B_BB, &Bt_BB));
1654:     } else {
1655:       PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_I_local, MAT_INITIAL_MATRIX, &B_BI));
1656:       PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_B_local, MAT_INITIAL_MATRIX, &B_BB));
1657:       PetscCall(MatCreateSubMatrix(lA, pcis->is_I_local, lP, MAT_INITIAL_MATRIX, &Bt_BI));
1658:       PetscCall(MatCreateSubMatrix(lA, pcis->is_B_local, lP, MAT_INITIAL_MATRIX, &Bt_BB));
1659:     }
1660:     PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_B_BI", (PetscObject)B_BI));
1661:     PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_B_BB", (PetscObject)B_BB));
1662:     PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_Bt_BI", (PetscObject)Bt_BI));
1663:     PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_Bt_BB", (PetscObject)Bt_BB));
1664:     PetscCall(MatDestroy(&B_BI));
1665:     PetscCall(MatDestroy(&B_BB));
1666:     PetscCall(MatDestroy(&Bt_BI));
1667:     PetscCall(MatDestroy(&Bt_BB));
1668:   }
1669:   PetscCall(PetscLogEventEnd(PC_BDDC_LocalWork[pcbddc->current_level], pc, 0, 0, 0));

1671:   /* SetUp coarse and local Neumann solvers */
1672:   PetscCall(PCBDDCSetUpSolvers(pc));
1673:   /* SetUp Scaling operator */
1674:   if (pcbddc->use_deluxe_scaling) PetscCall(PCBDDCScalingSetUp(pc));

1676:   /* mark topography as done */
1677:   pcbddc->recompute_topography = PETSC_FALSE;

1679:   /* wrap pcis->A_IB and pcis->A_BI if we did not change explicitly the variables on the pressures */
1680:   PetscCall(PCBDDCBenignShellMat(pc, PETSC_FALSE));

1682:   if (pcbddc->dbg_flag) {
1683:     PetscCall(PetscViewerASCIISubtractTab(pcbddc->dbg_viewer, 2 * pcbddc->current_level));
1684:     PetscCall(PetscViewerASCIIPopSynchronized(pcbddc->dbg_viewer));
1685:   }

1687:   { /* Dump customization */
1688:     PetscBool flg;
1689:     char      save[PETSC_MAX_PATH_LEN] = {'\0'};

1691:     PetscCall(PetscOptionsGetString(NULL, ((PetscObject)pc)->prefix, "-pc_bddc_save", save, sizeof(save), &flg));
1692:     if (flg) {
1693:       size_t len;

1695:       PetscCall(PetscStrlen(save, &len));
1696:       PetscCall(PCBDDCLoadOrViewCustomization(pc, PETSC_FALSE, len ? save : NULL));
1697:     }
1698:   }
1699:   PetscFunctionReturn(PETSC_SUCCESS);
1700: }

1702: static PetscErrorCode PCApply_BDDC(PC pc, Vec r, Vec z)
1703: {
1704:   PC_IS            *pcis   = (PC_IS *)pc->data;
1705:   PC_BDDC          *pcbddc = (PC_BDDC *)pc->data;
1706:   Mat               lA     = NULL;
1707:   PetscInt          n_B = pcis->n_B, n_D = pcis->n - n_B;
1708:   const PetscScalar one   = 1.0;
1709:   const PetscScalar m_one = -1.0;
1710:   const PetscScalar zero  = 0.0;
1711:   /* This code is similar to that provided in nn.c for PCNN
1712:    NN interface preconditioner changed to BDDC
1713:    Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */

1715:   PetscFunctionBegin;
1716:   PetscCall(PetscCitationsRegister(citation, &cited));
1717:   if (pcbddc->switch_static) PetscCall(MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat, &lA));

1719:   if (pcbddc->ChangeOfBasisMatrix) {
1720:     Vec swap;

1722:     PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, r, pcbddc->work_change));
1723:     swap                = pcbddc->work_change;
1724:     pcbddc->work_change = r;
1725:     r                   = swap;
1726:     /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */
1727:     if (pcbddc->benign_apply_coarse_only && pcbddc->use_exact_dirichlet_trick && pcbddc->change_interior) {
1728:       PetscCall(VecCopy(r, pcis->vec1_global));
1729:       PetscCall(VecLockReadPush(pcis->vec1_global));
1730:     }
1731:   }
1732:   if (pcbddc->benign_have_null) { /* get p0 from r */
1733:     PetscCall(PCBDDCBenignGetOrSetP0(pc, r, PETSC_TRUE));
1734:   }
1735:   if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_DIRICHLET && !pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
1736:     PetscCall(VecCopy(r, z));
1737:     /* First Dirichlet solve */
1738:     PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1739:     PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1740:     /*
1741:       Assembling right-hand side for BDDC operator
1742:       - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
1743:       - pcis->vec1_B the interface part of the global vector z
1744:     */
1745:     PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1746:     if (n_D) {
1747:       PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
1748:       PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1749:       PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
1750:       PetscCall(VecScale(pcis->vec2_D, m_one));
1751:       if (pcbddc->switch_static) {
1752:         PetscCall(VecSet(pcis->vec1_N, 0.));
1753:         PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1754:         PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1755:         if (!pcbddc->switch_static_change) {
1756:           PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
1757:         } else {
1758:           PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1759:           PetscCall(MatMult(lA, pcis->vec2_N, pcis->vec1_N));
1760:           PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1761:         }
1762:         PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
1763:         PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
1764:         PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1765:         PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1766:       } else {
1767:         PetscCall(MatMult(pcis->A_BI, pcis->vec2_D, pcis->vec1_B));
1768:       }
1769:     } else {
1770:       PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1771:       PetscCall(VecSet(pcis->vec1_B, zero));
1772:     }
1773:     PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
1774:     PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
1775:     PetscCall(PCBDDCScalingRestriction(pc, z, pcis->vec1_B));
1776:   } else {
1777:     if (!pcbddc->benign_apply_coarse_only) PetscCall(PCBDDCScalingRestriction(pc, r, pcis->vec1_B));
1778:   }
1779:   if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_LUMP) {
1780:     PetscCheck(pcbddc->switch_static, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "You forgot to pass -pc_bddc_switch_static");
1781:     PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1782:     PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1783:   }

1785:   /* Apply interface preconditioner
1786:      input/output vecs: pcis->vec1_B and pcis->vec1_D */
1787:   PetscCall(PCBDDCApplyInterfacePreconditioner(pc, PETSC_FALSE));

1789:   /* Apply transpose of partition of unity operator */
1790:   PetscCall(PCBDDCScalingExtension(pc, pcis->vec1_B, z));
1791:   if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_LUMP) {
1792:     PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_D, z, INSERT_VALUES, SCATTER_REVERSE));
1793:     PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_D, z, INSERT_VALUES, SCATTER_REVERSE));
1794:     PetscFunctionReturn(PETSC_SUCCESS);
1795:   }
1796:   /* Second Dirichlet solve and assembling of output */
1797:   PetscCall(VecScatterBegin(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1798:   PetscCall(VecScatterEnd(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1799:   if (n_B) {
1800:     if (pcbddc->switch_static) {
1801:       PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1802:       PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1803:       PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1804:       PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1805:       if (!pcbddc->switch_static_change) {
1806:         PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
1807:       } else {
1808:         PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1809:         PetscCall(MatMult(lA, pcis->vec2_N, pcis->vec1_N));
1810:         PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1811:       }
1812:       PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
1813:       PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
1814:     } else {
1815:       PetscCall(MatMult(pcis->A_IB, pcis->vec1_B, pcis->vec3_D));
1816:     }
1817:   } else if (pcbddc->switch_static) { /* n_B is zero */
1818:     if (!pcbddc->switch_static_change) {
1819:       PetscCall(MatMult(lA, pcis->vec1_D, pcis->vec3_D));
1820:     } else {
1821:       PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_D, pcis->vec1_N));
1822:       PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
1823:       PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec2_N, pcis->vec3_D));
1824:     }
1825:   }
1826:   PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1827:   PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec3_D, pcis->vec4_D));
1828:   PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1829:   PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec4_D));

1831:   if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
1832:     if (pcbddc->switch_static) {
1833:       PetscCall(VecAXPBYPCZ(pcis->vec2_D, m_one, one, m_one, pcis->vec4_D, pcis->vec1_D));
1834:     } else {
1835:       PetscCall(VecAXPBY(pcis->vec2_D, m_one, m_one, pcis->vec4_D));
1836:     }
1837:     PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
1838:     PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
1839:   } else {
1840:     if (pcbddc->switch_static) {
1841:       PetscCall(VecAXPBY(pcis->vec4_D, one, m_one, pcis->vec1_D));
1842:     } else {
1843:       PetscCall(VecScale(pcis->vec4_D, m_one));
1844:     }
1845:     PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
1846:     PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
1847:   }
1848:   if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */
1849:     if (pcbddc->benign_apply_coarse_only) PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
1850:     PetscCall(PCBDDCBenignGetOrSetP0(pc, z, PETSC_FALSE));
1851:   }

1853:   if (pcbddc->ChangeOfBasisMatrix) {
1854:     pcbddc->work_change = r;
1855:     PetscCall(VecCopy(z, pcbddc->work_change));
1856:     PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcbddc->work_change, z));
1857:   }
1858:   PetscFunctionReturn(PETSC_SUCCESS);
1859: }

1861: static PetscErrorCode PCApplyTranspose_BDDC(PC pc, Vec r, Vec z)
1862: {
1863:   PC_IS            *pcis   = (PC_IS *)pc->data;
1864:   PC_BDDC          *pcbddc = (PC_BDDC *)pc->data;
1865:   Mat               lA     = NULL;
1866:   PetscInt          n_B = pcis->n_B, n_D = pcis->n - n_B;
1867:   const PetscScalar one   = 1.0;
1868:   const PetscScalar m_one = -1.0;
1869:   const PetscScalar zero  = 0.0;

1871:   PetscFunctionBegin;
1872:   PetscCall(PetscCitationsRegister(citation, &cited));
1873:   if (pcbddc->switch_static) PetscCall(MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat, &lA));
1874:   if (pcbddc->ChangeOfBasisMatrix) {
1875:     Vec swap;

1877:     PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, r, pcbddc->work_change));
1878:     swap                = pcbddc->work_change;
1879:     pcbddc->work_change = r;
1880:     r                   = swap;
1881:     /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */
1882:     if (pcbddc->benign_apply_coarse_only && pcbddc->exact_dirichlet_trick_app && pcbddc->change_interior) {
1883:       PetscCall(VecCopy(r, pcis->vec1_global));
1884:       PetscCall(VecLockReadPush(pcis->vec1_global));
1885:     }
1886:   }
1887:   if (pcbddc->benign_have_null) { /* get p0 from r */
1888:     PetscCall(PCBDDCBenignGetOrSetP0(pc, r, PETSC_TRUE));
1889:   }
1890:   if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
1891:     PetscCall(VecCopy(r, z));
1892:     /* First Dirichlet solve */
1893:     PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1894:     PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1895:     /*
1896:       Assembling right-hand side for BDDC operator
1897:       - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
1898:       - pcis->vec1_B the interface part of the global vector z
1899:     */
1900:     if (n_D) {
1901:       PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1902:       PetscCall(KSPSolveTranspose(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
1903:       PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1904:       PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
1905:       PetscCall(VecScale(pcis->vec2_D, m_one));
1906:       if (pcbddc->switch_static) {
1907:         PetscCall(VecSet(pcis->vec1_N, 0.));
1908:         PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1909:         PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1910:         if (!pcbddc->switch_static_change) {
1911:           PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
1912:         } else {
1913:           PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1914:           PetscCall(MatMultTranspose(lA, pcis->vec2_N, pcis->vec1_N));
1915:           PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1916:         }
1917:         PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
1918:         PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
1919:         PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1920:         PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1921:       } else {
1922:         PetscCall(MatMultTranspose(pcis->A_IB, pcis->vec2_D, pcis->vec1_B));
1923:       }
1924:     } else {
1925:       PetscCall(VecSet(pcis->vec1_B, zero));
1926:     }
1927:     PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
1928:     PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
1929:     PetscCall(PCBDDCScalingRestriction(pc, z, pcis->vec1_B));
1930:   } else {
1931:     PetscCall(PCBDDCScalingRestriction(pc, r, pcis->vec1_B));
1932:   }

1934:   /* Apply interface preconditioner
1935:      input/output vecs: pcis->vec1_B and pcis->vec1_D */
1936:   PetscCall(PCBDDCApplyInterfacePreconditioner(pc, PETSC_TRUE));

1938:   /* Apply transpose of partition of unity operator */
1939:   PetscCall(PCBDDCScalingExtension(pc, pcis->vec1_B, z));

1941:   /* Second Dirichlet solve and assembling of output */
1942:   PetscCall(VecScatterBegin(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1943:   PetscCall(VecScatterEnd(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1944:   if (n_B) {
1945:     if (pcbddc->switch_static) {
1946:       PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1947:       PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1948:       PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1949:       PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1950:       if (!pcbddc->switch_static_change) {
1951:         PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
1952:       } else {
1953:         PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1954:         PetscCall(MatMultTranspose(lA, pcis->vec2_N, pcis->vec1_N));
1955:         PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1956:       }
1957:       PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
1958:       PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
1959:     } else {
1960:       PetscCall(MatMultTranspose(pcis->A_BI, pcis->vec1_B, pcis->vec3_D));
1961:     }
1962:   } else if (pcbddc->switch_static) { /* n_B is zero */
1963:     if (!pcbddc->switch_static_change) {
1964:       PetscCall(MatMultTranspose(lA, pcis->vec1_D, pcis->vec3_D));
1965:     } else {
1966:       PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_D, pcis->vec1_N));
1967:       PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
1968:       PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec2_N, pcis->vec3_D));
1969:     }
1970:   }
1971:   PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1972:   PetscCall(KSPSolveTranspose(pcbddc->ksp_D, pcis->vec3_D, pcis->vec4_D));
1973:   PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1974:   PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec4_D));
1975:   if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
1976:     if (pcbddc->switch_static) {
1977:       PetscCall(VecAXPBYPCZ(pcis->vec2_D, m_one, one, m_one, pcis->vec4_D, pcis->vec1_D));
1978:     } else {
1979:       PetscCall(VecAXPBY(pcis->vec2_D, m_one, m_one, pcis->vec4_D));
1980:     }
1981:     PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
1982:     PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
1983:   } else {
1984:     if (pcbddc->switch_static) {
1985:       PetscCall(VecAXPBY(pcis->vec4_D, one, m_one, pcis->vec1_D));
1986:     } else {
1987:       PetscCall(VecScale(pcis->vec4_D, m_one));
1988:     }
1989:     PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
1990:     PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
1991:   }
1992:   if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */
1993:     PetscCall(PCBDDCBenignGetOrSetP0(pc, z, PETSC_FALSE));
1994:   }
1995:   if (pcbddc->ChangeOfBasisMatrix) {
1996:     pcbddc->work_change = r;
1997:     PetscCall(VecCopy(z, pcbddc->work_change));
1998:     PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcbddc->work_change, z));
1999:   }
2000:   PetscFunctionReturn(PETSC_SUCCESS);
2001: }

2003: static PetscErrorCode PCReset_BDDC(PC pc)
2004: {
2005:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
2006:   PC_IS   *pcis   = (PC_IS *)pc->data;
2007:   KSP      kspD, kspR, kspC;

2009:   PetscFunctionBegin;
2010:   /* free BDDC custom data  */
2011:   PetscCall(PCBDDCResetCustomization(pc));
2012:   /* destroy objects related to topography */
2013:   PetscCall(PCBDDCResetTopography(pc));
2014:   /* destroy objects for scaling operator */
2015:   PetscCall(PCBDDCScalingDestroy(pc));
2016:   /* free solvers stuff */
2017:   PetscCall(PCBDDCResetSolvers(pc));
2018:   /* free global vectors needed in presolve */
2019:   PetscCall(VecDestroy(&pcbddc->temp_solution));
2020:   PetscCall(VecDestroy(&pcbddc->original_rhs));
2021:   /* free data created by PCIS */
2022:   PetscCall(PCISReset(pc));

2024:   /* restore defaults */
2025:   kspD = pcbddc->ksp_D;
2026:   kspR = pcbddc->ksp_R;
2027:   kspC = pcbddc->coarse_ksp;
2028:   PetscCall(PetscMemzero(pc->data, sizeof(*pcbddc)));
2029:   pcis->n_neigh                     = -1;
2030:   pcis->scaling_factor              = 1.0;
2031:   pcis->reusesubmatrices            = PETSC_TRUE;
2032:   pcbddc->use_local_adj             = PETSC_TRUE;
2033:   pcbddc->use_vertices              = PETSC_TRUE;
2034:   pcbddc->use_edges                 = PETSC_TRUE;
2035:   pcbddc->symmetric_primal          = PETSC_TRUE;
2036:   pcbddc->vertex_size               = 1;
2037:   pcbddc->recompute_topography      = PETSC_TRUE;
2038:   pcbddc->coarse_size               = -1;
2039:   pcbddc->use_exact_dirichlet_trick = PETSC_TRUE;
2040:   pcbddc->coarsening_ratio          = 8;
2041:   pcbddc->coarse_eqs_per_proc       = 1;
2042:   pcbddc->benign_compute_correction = PETSC_TRUE;
2043:   pcbddc->nedfield                  = -1;
2044:   pcbddc->nedglobal                 = PETSC_TRUE;
2045:   pcbddc->graphmaxcount             = PETSC_INT_MAX;
2046:   pcbddc->sub_schurs_layers         = -1;
2047:   pcbddc->ksp_D                     = kspD;
2048:   pcbddc->ksp_R                     = kspR;
2049:   pcbddc->coarse_ksp                = kspC;
2050:   PetscFunctionReturn(PETSC_SUCCESS);
2051: }

2053: static PetscErrorCode PCDestroy_BDDC(PC pc)
2054: {
2055:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

2057:   PetscFunctionBegin;
2058:   PetscCall(PCReset_BDDC(pc));
2059:   PetscCall(KSPDestroy(&pcbddc->ksp_D));
2060:   PetscCall(KSPDestroy(&pcbddc->ksp_R));
2061:   PetscCall(KSPDestroy(&pcbddc->coarse_ksp));
2062:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDiscreteGradient_C", NULL));
2063:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDivergenceMat_C", NULL));
2064:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetChangeOfBasisMat_C", NULL));
2065:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesLocalIS_C", NULL));
2066:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesIS_C", NULL));
2067:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesLocalIS_C", NULL));
2068:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesIS_C", NULL));
2069:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetCoarseningRatio_C", NULL));
2070:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevel_C", NULL));
2071:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetUseExactDirichlet_C", NULL));
2072:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevels_C", NULL));
2073:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundaries_C", NULL));
2074:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundariesLocal_C", NULL));
2075:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundaries_C", NULL));
2076:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundariesLocal_C", NULL));
2077:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundaries_C", NULL));
2078:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundariesLocal_C", NULL));
2079:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundaries_C", NULL));
2080:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundariesLocal_C", NULL));
2081:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplitting_C", NULL));
2082:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplittingLocal_C", NULL));
2083:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", NULL));
2084:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCCreateFETIDPOperators_C", NULL));
2085:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetRHS_C", NULL));
2086:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetSolution_C", NULL));
2087:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", NULL));
2088:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
2089:   PetscCall(PetscFree(pc->data));
2090:   PetscFunctionReturn(PETSC_SUCCESS);
2091: }

2093: static PetscErrorCode PCSetCoordinates_BDDC(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
2094: {
2095:   PC_BDDC    *pcbddc    = (PC_BDDC *)pc->data;
2096:   PCBDDCGraph mat_graph = pcbddc->mat_graph;

2098:   PetscFunctionBegin;
2099:   PetscCall(PetscFree(mat_graph->coords));
2100:   PetscCall(PetscMalloc1(nloc * dim, &mat_graph->coords));
2101:   PetscCall(PetscArraycpy(mat_graph->coords, coords, nloc * dim));
2102:   mat_graph->cnloc = nloc;
2103:   mat_graph->cdim  = dim;
2104:   mat_graph->cloc  = PETSC_FALSE;
2105:   /* flg setup */
2106:   pcbddc->recompute_topography = PETSC_TRUE;
2107:   pcbddc->corner_selected      = PETSC_FALSE;
2108:   PetscFunctionReturn(PETSC_SUCCESS);
2109: }

2111: static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool *change)
2112: {
2113:   PetscFunctionBegin;
2114:   *change = PETSC_TRUE;
2115:   PetscFunctionReturn(PETSC_SUCCESS);
2116: }

2118: static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
2119: {
2120:   FETIDPMat_ctx mat_ctx;
2121:   Vec           work;
2122:   PC_IS        *pcis;
2123:   PC_BDDC      *pcbddc;

2125:   PetscFunctionBegin;
2126:   PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2127:   pcis   = (PC_IS *)mat_ctx->pc->data;
2128:   pcbddc = (PC_BDDC *)mat_ctx->pc->data;

2130:   PetscCall(VecSet(fetidp_flux_rhs, 0.0));
2131:   /* copy rhs since we may change it during PCPreSolve_BDDC */
2132:   if (!pcbddc->original_rhs) PetscCall(VecDuplicate(pcis->vec1_global, &pcbddc->original_rhs));
2133:   if (mat_ctx->rhs_flip) {
2134:     PetscCall(VecPointwiseMult(pcbddc->original_rhs, standard_rhs, mat_ctx->rhs_flip));
2135:   } else {
2136:     PetscCall(VecCopy(standard_rhs, pcbddc->original_rhs));
2137:   }
2138:   if (mat_ctx->g2g_p) {
2139:     /* interface pressure rhs */
2140:     PetscCall(VecScatterBegin(mat_ctx->g2g_p, fetidp_flux_rhs, pcbddc->original_rhs, INSERT_VALUES, SCATTER_REVERSE));
2141:     PetscCall(VecScatterEnd(mat_ctx->g2g_p, fetidp_flux_rhs, pcbddc->original_rhs, INSERT_VALUES, SCATTER_REVERSE));
2142:     PetscCall(VecScatterBegin(mat_ctx->g2g_p, standard_rhs, fetidp_flux_rhs, INSERT_VALUES, SCATTER_FORWARD));
2143:     PetscCall(VecScatterEnd(mat_ctx->g2g_p, standard_rhs, fetidp_flux_rhs, INSERT_VALUES, SCATTER_FORWARD));
2144:     if (!mat_ctx->rhs_flip) PetscCall(VecScale(fetidp_flux_rhs, -1.));
2145:   }
2146:   /*
2147:      change of basis for physical rhs if needed
2148:      It also changes the rhs in case of dirichlet boundaries
2149:   */
2150:   PetscCall(PCPreSolve_BDDC(mat_ctx->pc, NULL, pcbddc->original_rhs, NULL));
2151:   if (pcbddc->ChangeOfBasisMatrix) {
2152:     PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, pcbddc->original_rhs, pcbddc->work_change));
2153:     work = pcbddc->work_change;
2154:   } else {
2155:     work = pcbddc->original_rhs;
2156:   }
2157:   /* store vectors for computation of fetidp final solution */
2158:   PetscCall(VecScatterBegin(pcis->global_to_D, work, mat_ctx->temp_solution_D, INSERT_VALUES, SCATTER_FORWARD));
2159:   PetscCall(VecScatterEnd(pcis->global_to_D, work, mat_ctx->temp_solution_D, INSERT_VALUES, SCATTER_FORWARD));
2160:   /* scale rhs since it should be unassembled */
2161:   /* TODO use counter scaling? (also below) */
2162:   PetscCall(VecScatterBegin(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2163:   PetscCall(VecScatterEnd(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2164:   /* Apply partition of unity */
2165:   PetscCall(VecPointwiseMult(mat_ctx->temp_solution_B, pcis->D, mat_ctx->temp_solution_B));
2166:   /* PetscCall(PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B)); */
2167:   if (!pcbddc->switch_static) {
2168:     /* compute partially subassembled Schur complement right-hand side */
2169:     PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2170:     PetscCall(KSPSolve(pcbddc->ksp_D, mat_ctx->temp_solution_D, pcis->vec1_D));
2171:     PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2172:     /* Cannot propagate up error in KSPSolve() because there is no access to the PC */
2173:     PetscCall(MatMult(pcis->A_BI, pcis->vec1_D, pcis->vec1_B));
2174:     PetscCall(VecAXPY(mat_ctx->temp_solution_B, -1.0, pcis->vec1_B));
2175:     PetscCall(VecSet(work, 0.0));
2176:     PetscCall(VecScatterBegin(pcis->global_to_B, mat_ctx->temp_solution_B, work, ADD_VALUES, SCATTER_REVERSE));
2177:     PetscCall(VecScatterEnd(pcis->global_to_B, mat_ctx->temp_solution_B, work, ADD_VALUES, SCATTER_REVERSE));
2178:     /* PetscCall(PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B)); */
2179:     PetscCall(VecScatterBegin(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2180:     PetscCall(VecScatterEnd(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2181:     PetscCall(VecPointwiseMult(mat_ctx->temp_solution_B, pcis->D, mat_ctx->temp_solution_B));
2182:   }
2183:   /* BDDC rhs */
2184:   PetscCall(VecCopy(mat_ctx->temp_solution_B, pcis->vec1_B));
2185:   if (pcbddc->switch_static) PetscCall(VecCopy(mat_ctx->temp_solution_D, pcis->vec1_D));
2186:   /* apply BDDC */
2187:   PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2188:   PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, PETSC_FALSE));
2189:   PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));

2191:   /* Application of B_delta and assembling of rhs for fetidp fluxes */
2192:   PetscCall(MatMult(mat_ctx->B_delta, pcis->vec1_B, mat_ctx->lambda_local));
2193:   PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, mat_ctx->lambda_local, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2194:   PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, mat_ctx->lambda_local, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2195:   /* Add contribution to interface pressures */
2196:   if (mat_ctx->l2g_p) {
2197:     PetscCall(MatMult(mat_ctx->B_BB, pcis->vec1_B, mat_ctx->vP));
2198:     if (pcbddc->switch_static) PetscCall(MatMultAdd(mat_ctx->B_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP));
2199:     PetscCall(VecScatterBegin(mat_ctx->l2g_p, mat_ctx->vP, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2200:     PetscCall(VecScatterEnd(mat_ctx->l2g_p, mat_ctx->vP, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2201:   }
2202:   PetscFunctionReturn(PETSC_SUCCESS);
2203: }

2205: /*@
2206:   PCBDDCMatFETIDPGetRHS - Compute the right-hand side for a FETI-DP linear system using the physical right-hand side

2208:   Collective

2210:   Input Parameters:
2211: + fetidp_mat   - the FETI-DP matrix object obtained by a call to `PCBDDCCreateFETIDPOperators()`
2212: - standard_rhs - the right-hand side of the original linear system

2214:   Output Parameter:
2215: . fetidp_flux_rhs - the right-hand side for the FETI-DP linear system

2217:   Level: developer

2219:   Note:
2220:   Most users should employ the `KSP` interface for linear solvers and create a solver of type `KSPFETIDP`.

2222: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCCreateFETIDPOperators()`, `PCBDDCMatFETIDPGetSolution()`
2223: @*/
2224: PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
2225: {
2226:   FETIDPMat_ctx mat_ctx;

2228:   PetscFunctionBegin;
2232:   PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2233:   PetscUseMethod(mat_ctx->pc, "PCBDDCMatFETIDPGetRHS_C", (Mat, Vec, Vec), (fetidp_mat, standard_rhs, fetidp_flux_rhs));
2234:   PetscFunctionReturn(PETSC_SUCCESS);
2235: }

2237: static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
2238: {
2239:   FETIDPMat_ctx mat_ctx;
2240:   PC_IS        *pcis;
2241:   PC_BDDC      *pcbddc;
2242:   Vec           work;

2244:   PetscFunctionBegin;
2245:   PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2246:   pcis   = (PC_IS *)mat_ctx->pc->data;
2247:   pcbddc = (PC_BDDC *)mat_ctx->pc->data;

2249:   /* apply B_delta^T */
2250:   PetscCall(VecSet(pcis->vec1_B, 0.));
2251:   PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, fetidp_flux_sol, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
2252:   PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, fetidp_flux_sol, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
2253:   PetscCall(MatMultTranspose(mat_ctx->B_delta, mat_ctx->lambda_local, pcis->vec1_B));
2254:   if (mat_ctx->l2g_p) {
2255:     PetscCall(VecScatterBegin(mat_ctx->l2g_p, fetidp_flux_sol, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
2256:     PetscCall(VecScatterEnd(mat_ctx->l2g_p, fetidp_flux_sol, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
2257:     PetscCall(MatMultAdd(mat_ctx->Bt_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B));
2258:   }

2260:   /* compute rhs for BDDC application */
2261:   PetscCall(VecAYPX(pcis->vec1_B, -1.0, mat_ctx->temp_solution_B));
2262:   if (pcbddc->switch_static) {
2263:     PetscCall(VecCopy(mat_ctx->temp_solution_D, pcis->vec1_D));
2264:     if (mat_ctx->l2g_p) {
2265:       PetscCall(VecScale(mat_ctx->vP, -1.));
2266:       PetscCall(MatMultAdd(mat_ctx->Bt_BI, mat_ctx->vP, pcis->vec1_D, pcis->vec1_D));
2267:     }
2268:   }

2270:   /* apply BDDC */
2271:   PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2272:   PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, PETSC_FALSE));

2274:   /* put values into global vector */
2275:   if (pcbddc->ChangeOfBasisMatrix) work = pcbddc->work_change;
2276:   else work = standard_sol;
2277:   PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, work, INSERT_VALUES, SCATTER_REVERSE));
2278:   PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, work, INSERT_VALUES, SCATTER_REVERSE));
2279:   if (!pcbddc->switch_static) {
2280:     /* compute values into the interior if solved for the partially subassembled Schur complement */
2281:     PetscCall(MatMult(pcis->A_IB, pcis->vec1_B, pcis->vec1_D));
2282:     PetscCall(VecAYPX(pcis->vec1_D, -1.0, mat_ctx->temp_solution_D));
2283:     PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2284:     PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec1_D));
2285:     PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2286:     /* Cannot propagate up error in KSPSolve() because there is no access to the PC */
2287:   }

2289:   PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_D, work, INSERT_VALUES, SCATTER_REVERSE));
2290:   PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_D, work, INSERT_VALUES, SCATTER_REVERSE));
2291:   /* add p0 solution to final solution */
2292:   PetscCall(PCBDDCBenignGetOrSetP0(mat_ctx->pc, work, PETSC_FALSE));
2293:   if (pcbddc->ChangeOfBasisMatrix) PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, work, standard_sol));
2294:   PetscCall(PCPostSolve_BDDC(mat_ctx->pc, NULL, NULL, standard_sol));
2295:   if (mat_ctx->g2g_p) {
2296:     PetscCall(VecScatterBegin(mat_ctx->g2g_p, fetidp_flux_sol, standard_sol, INSERT_VALUES, SCATTER_REVERSE));
2297:     PetscCall(VecScatterEnd(mat_ctx->g2g_p, fetidp_flux_sol, standard_sol, INSERT_VALUES, SCATTER_REVERSE));
2298:   }
2299:   PetscFunctionReturn(PETSC_SUCCESS);
2300: }

2302: static PetscErrorCode PCView_BDDCIPC(PC pc, PetscViewer viewer)
2303: {
2304:   BDDCIPC_ctx bddcipc_ctx;
2305:   PetscBool   isascii;

2307:   PetscFunctionBegin;
2308:   PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2309:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2310:   if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "BDDC interface preconditioner\n"));
2311:   PetscCall(PetscViewerASCIIPushTab(viewer));
2312:   PetscCall(PCView(bddcipc_ctx->bddc, viewer));
2313:   PetscCall(PetscViewerASCIIPopTab(viewer));
2314:   PetscFunctionReturn(PETSC_SUCCESS);
2315: }

2317: static PetscErrorCode PCSetUp_BDDCIPC(PC pc)
2318: {
2319:   BDDCIPC_ctx bddcipc_ctx;
2320:   PetscBool   isbddc;
2321:   Vec         vv;
2322:   IS          is;
2323:   PC_IS      *pcis;

2325:   PetscFunctionBegin;
2326:   PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2327:   PetscCall(PetscObjectTypeCompare((PetscObject)bddcipc_ctx->bddc, PCBDDC, &isbddc));
2328:   PetscCheck(isbddc, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Invalid type %s. Must be of type bddc", ((PetscObject)bddcipc_ctx->bddc)->type_name);
2329:   PetscCall(PCSetUp(bddcipc_ctx->bddc));

2331:   /* create interface scatter */
2332:   pcis = (PC_IS *)bddcipc_ctx->bddc->data;
2333:   PetscCall(VecScatterDestroy(&bddcipc_ctx->g2l));
2334:   PetscCall(MatCreateVecs(pc->pmat, &vv, NULL));
2335:   PetscCall(ISRenumber(pcis->is_B_global, NULL, NULL, &is));
2336:   PetscCall(VecScatterCreate(vv, is, pcis->vec1_B, NULL, &bddcipc_ctx->g2l));
2337:   PetscCall(ISDestroy(&is));
2338:   PetscCall(VecDestroy(&vv));
2339:   PetscFunctionReturn(PETSC_SUCCESS);
2340: }

2342: static PetscErrorCode PCApply_BDDCIPC(PC pc, Vec r, Vec x)
2343: {
2344:   BDDCIPC_ctx bddcipc_ctx;
2345:   PC_IS      *pcis;
2346:   VecScatter  tmps;

2348:   PetscFunctionBegin;
2349:   PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2350:   pcis              = (PC_IS *)bddcipc_ctx->bddc->data;
2351:   tmps              = pcis->global_to_B;
2352:   pcis->global_to_B = bddcipc_ctx->g2l;
2353:   PetscCall(PCBDDCScalingRestriction(bddcipc_ctx->bddc, r, pcis->vec1_B));
2354:   PetscCall(PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc, PETSC_FALSE));
2355:   PetscCall(PCBDDCScalingExtension(bddcipc_ctx->bddc, pcis->vec1_B, x));
2356:   pcis->global_to_B = tmps;
2357:   PetscFunctionReturn(PETSC_SUCCESS);
2358: }

2360: static PetscErrorCode PCApplyTranspose_BDDCIPC(PC pc, Vec r, Vec x)
2361: {
2362:   BDDCIPC_ctx bddcipc_ctx;
2363:   PC_IS      *pcis;
2364:   VecScatter  tmps;

2366:   PetscFunctionBegin;
2367:   PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2368:   pcis              = (PC_IS *)bddcipc_ctx->bddc->data;
2369:   tmps              = pcis->global_to_B;
2370:   pcis->global_to_B = bddcipc_ctx->g2l;
2371:   PetscCall(PCBDDCScalingRestriction(bddcipc_ctx->bddc, r, pcis->vec1_B));
2372:   PetscCall(PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc, PETSC_TRUE));
2373:   PetscCall(PCBDDCScalingExtension(bddcipc_ctx->bddc, pcis->vec1_B, x));
2374:   pcis->global_to_B = tmps;
2375:   PetscFunctionReturn(PETSC_SUCCESS);
2376: }

2378: static PetscErrorCode PCDestroy_BDDCIPC(PC pc)
2379: {
2380:   BDDCIPC_ctx bddcipc_ctx;

2382:   PetscFunctionBegin;
2383:   PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2384:   PetscCall(PCDestroy(&bddcipc_ctx->bddc));
2385:   PetscCall(VecScatterDestroy(&bddcipc_ctx->g2l));
2386:   PetscCall(PetscFree(bddcipc_ctx));
2387:   PetscFunctionReturn(PETSC_SUCCESS);
2388: }

2390: /*@
2391:   PCBDDCMatFETIDPGetSolution - Compute the physical solution using the solution of the FETI-DP linear system

2393:   Collective

2395:   Input Parameters:
2396: + fetidp_mat      - the FETI-DP matrix obtained by a call to `PCBDDCCreateFETIDPOperators()`
2397: - fetidp_flux_sol - the solution of the FETI-DP linear system`

2399:   Output Parameter:
2400: . standard_sol - the solution defined on the physical domain

2402:   Level: developer

2404:   Note:
2405:   Most users should employ the `KSP` interface for linear solvers and create a solver of type `KSPFETIDP`.

2407: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCCreateFETIDPOperators()`, `PCBDDCMatFETIDPGetRHS()`
2408: @*/
2409: PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
2410: {
2411:   FETIDPMat_ctx mat_ctx;

2413:   PetscFunctionBegin;
2417:   PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2418:   PetscUseMethod(mat_ctx->pc, "PCBDDCMatFETIDPGetSolution_C", (Mat, Vec, Vec), (fetidp_mat, fetidp_flux_sol, standard_sol));
2419:   PetscFunctionReturn(PETSC_SUCCESS);
2420: }

2422: static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc)
2423: {
2424:   FETIDPMat_ctx fetidpmat_ctx;
2425:   Mat           newmat;
2426:   FETIDPPC_ctx  fetidppc_ctx;
2427:   PC            newpc;
2428:   MPI_Comm      comm;

2430:   PetscFunctionBegin;
2431:   PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
2432:   /* FETI-DP matrix */
2433:   PetscCall(PCBDDCCreateFETIDPMatContext(pc, &fetidpmat_ctx));
2434:   fetidpmat_ctx->fully_redundant = fully_redundant;
2435:   PetscCall(PCBDDCSetupFETIDPMatContext(fetidpmat_ctx));
2436:   PetscCall(MatCreateShell(comm, fetidpmat_ctx->n, fetidpmat_ctx->n, fetidpmat_ctx->N, fetidpmat_ctx->N, fetidpmat_ctx, &newmat));
2437:   PetscCall(PetscObjectSetName((PetscObject)newmat, !fetidpmat_ctx->l2g_lambda_only ? "F" : "G"));
2438:   PetscCall(MatShellSetOperation(newmat, MATOP_MULT, (void (*)(void))FETIDPMatMult));
2439:   PetscCall(MatShellSetOperation(newmat, MATOP_MULT_TRANSPOSE, (void (*)(void))FETIDPMatMultTranspose));
2440:   PetscCall(MatShellSetOperation(newmat, MATOP_DESTROY, (void (*)(void))PCBDDCDestroyFETIDPMat));
2441:   /* propagate MatOptions */
2442:   {
2443:     PC_BDDC  *pcbddc = (PC_BDDC *)fetidpmat_ctx->pc->data;
2444:     PetscBool isset, issym;

2446:     PetscCall(MatIsSymmetricKnown(pc->mat, &isset, &issym));
2447:     if ((isset && issym) || pcbddc->symmetric_primal) PetscCall(MatSetOption(newmat, MAT_SYMMETRIC, PETSC_TRUE));
2448:   }
2449:   PetscCall(MatSetOptionsPrefix(newmat, prefix));
2450:   PetscCall(MatAppendOptionsPrefix(newmat, "fetidp_"));
2451:   PetscCall(MatSetUp(newmat));
2452:   /* FETI-DP preconditioner */
2453:   PetscCall(PCBDDCCreateFETIDPPCContext(pc, &fetidppc_ctx));
2454:   PetscCall(PCBDDCSetupFETIDPPCContext(newmat, fetidppc_ctx));
2455:   PetscCall(PCCreate(comm, &newpc));
2456:   PetscCall(PCSetOperators(newpc, newmat, newmat));
2457:   PetscCall(PCSetOptionsPrefix(newpc, prefix));
2458:   PetscCall(PCAppendOptionsPrefix(newpc, "fetidp_"));
2459:   PetscCall(PCSetErrorIfFailure(newpc, pc->erroriffailure));
2460:   if (!fetidpmat_ctx->l2g_lambda_only) { /* standard FETI-DP */
2461:     PetscCall(PCSetType(newpc, PCSHELL));
2462:     PetscCall(PCShellSetName(newpc, "FETI-DP multipliers"));
2463:     PetscCall(PCShellSetContext(newpc, fetidppc_ctx));
2464:     PetscCall(PCShellSetApply(newpc, FETIDPPCApply));
2465:     PetscCall(PCShellSetApplyTranspose(newpc, FETIDPPCApplyTranspose));
2466:     PetscCall(PCShellSetView(newpc, FETIDPPCView));
2467:     PetscCall(PCShellSetDestroy(newpc, PCBDDCDestroyFETIDPPC));
2468:   } else { /* saddle-point FETI-DP */
2469:     Mat       M;
2470:     PetscInt  psize;
2471:     PetscBool fake = PETSC_FALSE, isfieldsplit;

2473:     PetscCall(ISViewFromOptions(fetidpmat_ctx->lagrange, NULL, "-lag_view"));
2474:     PetscCall(ISViewFromOptions(fetidpmat_ctx->pressure, NULL, "-press_view"));
2475:     PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_PPmat", (PetscObject *)&M));
2476:     PetscCall(PCSetType(newpc, PCFIELDSPLIT));
2477:     PetscCall(PCFieldSplitSetIS(newpc, "lag", fetidpmat_ctx->lagrange));
2478:     PetscCall(PCFieldSplitSetIS(newpc, "p", fetidpmat_ctx->pressure));
2479:     PetscCall(PCFieldSplitSetType(newpc, PC_COMPOSITE_SCHUR));
2480:     PetscCall(PCFieldSplitSetSchurFactType(newpc, PC_FIELDSPLIT_SCHUR_FACT_DIAG));
2481:     PetscCall(ISGetSize(fetidpmat_ctx->pressure, &psize));
2482:     if (psize != M->rmap->N) {
2483:       Mat      M2;
2484:       PetscInt lpsize;

2486:       fake = PETSC_TRUE;
2487:       PetscCall(ISGetLocalSize(fetidpmat_ctx->pressure, &lpsize));
2488:       PetscCall(MatCreate(comm, &M2));
2489:       PetscCall(MatSetType(M2, MATAIJ));
2490:       PetscCall(MatSetSizes(M2, lpsize, lpsize, psize, psize));
2491:       PetscCall(MatSetUp(M2));
2492:       PetscCall(MatAssemblyBegin(M2, MAT_FINAL_ASSEMBLY));
2493:       PetscCall(MatAssemblyEnd(M2, MAT_FINAL_ASSEMBLY));
2494:       PetscCall(PCFieldSplitSetSchurPre(newpc, PC_FIELDSPLIT_SCHUR_PRE_USER, M2));
2495:       PetscCall(MatDestroy(&M2));
2496:     } else {
2497:       PetscCall(PCFieldSplitSetSchurPre(newpc, PC_FIELDSPLIT_SCHUR_PRE_USER, M));
2498:     }
2499:     PetscCall(PCFieldSplitSetSchurScale(newpc, 1.0));

2501:     /* we need to setfromoptions and setup here to access the blocks */
2502:     PetscCall(PCSetFromOptions(newpc));
2503:     PetscCall(PCSetUp(newpc));

2505:     /* user may have changed the type (e.g. -fetidp_pc_type none) */
2506:     PetscCall(PetscObjectTypeCompare((PetscObject)newpc, PCFIELDSPLIT, &isfieldsplit));
2507:     if (isfieldsplit) {
2508:       KSP      *ksps;
2509:       PC        ppc, lagpc;
2510:       PetscInt  nn;
2511:       PetscBool ismatis, matisok = PETSC_FALSE, check = PETSC_FALSE;

2513:       /* set the solver for the (0,0) block */
2514:       PetscCall(PCFieldSplitSchurGetSubKSP(newpc, &nn, &ksps));
2515:       if (!nn) { /* not of type PC_COMPOSITE_SCHUR */
2516:         PetscCall(PCFieldSplitGetSubKSP(newpc, &nn, &ksps));
2517:         if (!fake) { /* pass pmat to the pressure solver */
2518:           Mat F;

2520:           PetscCall(KSPGetOperators(ksps[1], &F, NULL));
2521:           PetscCall(KSPSetOperators(ksps[1], F, M));
2522:         }
2523:       } else {
2524:         PetscBool issym, isset;
2525:         Mat       S;

2527:         PetscCall(PCFieldSplitSchurGetS(newpc, &S));
2528:         PetscCall(MatIsSymmetricKnown(newmat, &isset, &issym));
2529:         if (isset) PetscCall(MatSetOption(S, MAT_SYMMETRIC, issym));
2530:       }
2531:       PetscCall(KSPGetPC(ksps[0], &lagpc));
2532:       PetscCall(PCSetType(lagpc, PCSHELL));
2533:       PetscCall(PCShellSetName(lagpc, "FETI-DP multipliers"));
2534:       PetscCall(PCShellSetContext(lagpc, fetidppc_ctx));
2535:       PetscCall(PCShellSetApply(lagpc, FETIDPPCApply));
2536:       PetscCall(PCShellSetApplyTranspose(lagpc, FETIDPPCApplyTranspose));
2537:       PetscCall(PCShellSetView(lagpc, FETIDPPCView));
2538:       PetscCall(PCShellSetDestroy(lagpc, PCBDDCDestroyFETIDPPC));

2540:       /* Olof's idea: interface Schur complement preconditioner for the mass matrix */
2541:       PetscCall(KSPGetPC(ksps[1], &ppc));
2542:       if (fake) {
2543:         BDDCIPC_ctx    bddcipc_ctx;
2544:         PetscContainer c;

2546:         matisok = PETSC_TRUE;

2548:         /* create inner BDDC solver */
2549:         PetscCall(PetscNew(&bddcipc_ctx));
2550:         PetscCall(PCCreate(comm, &bddcipc_ctx->bddc));
2551:         PetscCall(PCSetType(bddcipc_ctx->bddc, PCBDDC));
2552:         PetscCall(PCSetOperators(bddcipc_ctx->bddc, M, M));
2553:         PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_pCSR", (PetscObject *)&c));
2554:         PetscCall(PetscObjectTypeCompare((PetscObject)M, MATIS, &ismatis));
2555:         if (c && ismatis) {
2556:           Mat       lM;
2557:           PetscInt *csr, n;

2559:           PetscCall(MatISGetLocalMat(M, &lM));
2560:           PetscCall(MatGetSize(lM, &n, NULL));
2561:           PetscCall(PetscContainerGetPointer(c, (void **)&csr));
2562:           PetscCall(PCBDDCSetLocalAdjacencyGraph(bddcipc_ctx->bddc, n, csr, csr + (n + 1), PETSC_COPY_VALUES));
2563:           PetscCall(MatISRestoreLocalMat(M, &lM));
2564:         }
2565:         PetscCall(PCSetOptionsPrefix(bddcipc_ctx->bddc, ((PetscObject)ksps[1])->prefix));
2566:         PetscCall(PCSetErrorIfFailure(bddcipc_ctx->bddc, pc->erroriffailure));
2567:         PetscCall(PCSetFromOptions(bddcipc_ctx->bddc));

2569:         /* wrap the interface application */
2570:         PetscCall(PCSetType(ppc, PCSHELL));
2571:         PetscCall(PCShellSetName(ppc, "FETI-DP pressure"));
2572:         PetscCall(PCShellSetContext(ppc, bddcipc_ctx));
2573:         PetscCall(PCShellSetSetUp(ppc, PCSetUp_BDDCIPC));
2574:         PetscCall(PCShellSetApply(ppc, PCApply_BDDCIPC));
2575:         PetscCall(PCShellSetApplyTranspose(ppc, PCApplyTranspose_BDDCIPC));
2576:         PetscCall(PCShellSetView(ppc, PCView_BDDCIPC));
2577:         PetscCall(PCShellSetDestroy(ppc, PCDestroy_BDDCIPC));
2578:       }

2580:       /* determine if we need to assemble M to construct a preconditioner */
2581:       if (!matisok) {
2582:         PetscCall(PetscObjectTypeCompare((PetscObject)M, MATIS, &ismatis));
2583:         PetscCall(PetscObjectTypeCompareAny((PetscObject)ppc, &matisok, PCBDDC, PCJACOBI, PCNONE, PCMG, ""));
2584:         if (ismatis && !matisok) PetscCall(MatConvert(M, MATAIJ, MAT_INPLACE_MATRIX, &M));
2585:       }

2587:       /* run the subproblems to check convergence */
2588:       PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)newmat)->prefix, "-check_saddlepoint", &check, NULL));
2589:       if (check) {
2590:         PetscInt i;

2592:         for (i = 0; i < nn; i++) {
2593:           KSP       kspC;
2594:           PC        npc;
2595:           Mat       F, pF;
2596:           Vec       x, y;
2597:           PetscBool isschur, prec = PETSC_TRUE;

2599:           PetscCall(KSPCreate(PetscObjectComm((PetscObject)ksps[i]), &kspC));
2600:           PetscCall(KSPSetNestLevel(kspC, pc->kspnestlevel));
2601:           PetscCall(KSPSetOptionsPrefix(kspC, ((PetscObject)ksps[i])->prefix));
2602:           PetscCall(KSPAppendOptionsPrefix(kspC, "check_"));
2603:           PetscCall(KSPGetOperators(ksps[i], &F, &pF));
2604:           PetscCall(PetscObjectTypeCompare((PetscObject)F, MATSCHURCOMPLEMENT, &isschur));
2605:           if (isschur) {
2606:             KSP  kspS, kspS2;
2607:             Mat  A00, pA00, A10, A01, A11;
2608:             char prefix[256];

2610:             PetscCall(MatSchurComplementGetKSP(F, &kspS));
2611:             PetscCall(MatSchurComplementGetSubMatrices(F, &A00, &pA00, &A01, &A10, &A11));
2612:             PetscCall(MatCreateSchurComplement(A00, pA00, A01, A10, A11, &F));
2613:             PetscCall(MatSchurComplementGetKSP(F, &kspS2));
2614:             PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sschur_", ((PetscObject)kspC)->prefix));
2615:             PetscCall(KSPSetOptionsPrefix(kspS2, prefix));
2616:             PetscCall(KSPGetPC(kspS2, &npc));
2617:             PetscCall(PCSetType(npc, PCKSP));
2618:             PetscCall(PCKSPSetKSP(npc, kspS));
2619:             PetscCall(KSPSetFromOptions(kspS2));
2620:             PetscCall(KSPGetPC(kspS2, &npc));
2621:             PetscCall(PCSetUseAmat(npc, PETSC_TRUE));
2622:           } else {
2623:             PetscCall(PetscObjectReference((PetscObject)F));
2624:           }
2625:           PetscCall(KSPSetFromOptions(kspC));
2626:           PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)kspC)->prefix, "-preconditioned", &prec, NULL));
2627:           if (prec) {
2628:             PetscCall(KSPGetPC(ksps[i], &npc));
2629:             PetscCall(KSPSetPC(kspC, npc));
2630:           }
2631:           PetscCall(KSPSetOperators(kspC, F, pF));
2632:           PetscCall(MatCreateVecs(F, &x, &y));
2633:           PetscCall(VecSetRandom(x, NULL));
2634:           PetscCall(MatMult(F, x, y));
2635:           PetscCall(KSPSolve(kspC, y, x));
2636:           PetscCall(KSPCheckSolve(kspC, npc, x));
2637:           PetscCall(KSPDestroy(&kspC));
2638:           PetscCall(MatDestroy(&F));
2639:           PetscCall(VecDestroy(&x));
2640:           PetscCall(VecDestroy(&y));
2641:         }
2642:       }
2643:       PetscCall(PetscFree(ksps));
2644:     }
2645:   }
2646:   /* return pointers for objects created */
2647:   *fetidp_mat = newmat;
2648:   *fetidp_pc  = newpc;
2649:   PetscFunctionReturn(PETSC_SUCCESS);
2650: }

2652: /*@C
2653:   PCBDDCCreateFETIDPOperators - Create FETI-DP operators

2655:   Collective

2657:   Input Parameters:
2658: + pc              - the `PCBDDC` preconditioning context (setup should have been called before)
2659: . fully_redundant - true for a fully redundant set of Lagrange multipliers
2660: - prefix          - optional options database prefix for the objects to be created (can be `NULL`)

2662:   Output Parameters:
2663: + fetidp_mat - shell FETI-DP matrix object
2664: - fetidp_pc  - shell Dirichlet preconditioner for FETI-DP matrix

2666:   Level: developer

2668:   Notes:
2669:   Most users should employ the `KSP` interface for linear solvers and create a solver of type `KSPFETIDP`.
2670:   Currently the only operations provided for the FETI-DP matrix are `MatMult()` and `MatMultTranspose()`

2672: .seealso: [](ch_ksp), `KSPFETIDP`, `PCBDDC`, `PCBDDCMatFETIDPGetRHS()`, `PCBDDCMatFETIDPGetSolution()`
2673: @*/
2674: PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc)
2675: {
2676:   PetscFunctionBegin;
2678:   PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "You must call PCSetup_BDDC() first");
2679:   PetscUseMethod(pc, "PCBDDCCreateFETIDPOperators_C", (PC, PetscBool, const char *, Mat *, PC *), (pc, fully_redundant, prefix, fetidp_mat, fetidp_pc));
2680:   PetscFunctionReturn(PETSC_SUCCESS);
2681: }

2683: /*MC
2684:    PCBDDC - Balancing Domain Decomposition by Constraints preconditioners, {cite}`dohrmann2007approximate`, {cite}`klawonn2006dual`, {cite}`mandel2008multispace`

2686:    Requires `MATIS` matrices (Pmat) with local matrices (inside the `MATIS`) of type `MATSEQAIJ`, `MATSEQBAIJ` or `MATSEQSBAIJ`

2688:    It also works with unsymmetric and indefinite problems.

2690:    Unlike 'conventional' interface preconditioners, `PCBDDC` iterates over all degrees of freedom, not just those on the interface. This allows the use
2691:    of approximate solvers on the subdomains.

2693:    Approximate local solvers are automatically adapted (see {cite}`dohrmann2007approximate`,) if the user has attached a nullspace object to the subdomain matrices, and informed
2694:    `PCBDDC` of using approximate solvers (via the command line).

2696:    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.
2697:    The latter can be customized by using `PCBDDCSetLocalAdjacencyGraph()`

2699:    Additional information on dofs can be provided by using `PCBDDCSetDofsSplitting()`, `PCBDDCSetDirichletBoundaries()`, `PCBDDCSetNeumannBoundaries()`, and
2700:    `PCBDDCSetPrimalVerticesIS()` and their local counterparts.

2702:    Constraints can be customized by attaching a `MatNullSpace` object to the `MATIS` matrix via `MatSetNearNullSpace()`. Non-singular modes are retained via SVD.

2704:    Change of basis is performed similarly to {cite}`klawonn2006dual` when requested. When more than one constraint is present on a single connected component
2705:    (i.e. an edge or a face), a robust method based on local QR factorizations is used.
2706:    User defined change of basis can be passed to `PCBDDC` with `PCBDDCSetChangeOfBasisMat()`

2708:    The PETSc implementation also supports multilevel `PCBDDC` {cite}`mandel2008multispace`. Coarse grids are partitioned using a `MatPartitioning` object.

2710:    Adaptive selection of primal constraints is supported for SPD systems with high-contrast in the coefficients if MUMPS or MKL_PARDISO are present.

2712:    Options Database Keys:
2713: +    -pc_bddc_use_vertices <true>         - use or not vertices in primal space
2714: .    -pc_bddc_use_edges <true>            - use or not edges in primal space
2715: .    -pc_bddc_use_faces <false>           - use or not faces in primal space
2716: .    -pc_bddc_symmetric <true>            - symmetric computation of primal basis functions. Specify false for unsymmetric problems
2717: .    -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only)
2718: .    -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested
2719: .    -pc_bddc_switch_static <false>       - switches from M_2 (default) to M_3 operator (see reference article [1])
2720: .    -pc_bddc_levels <0>                  - maximum number of levels for multilevel
2721: .    -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)
2722: .    -pc_bddc_coarse_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)
2723: .    -pc_bddc_use_deluxe_scaling <false>  - use deluxe scaling
2724: .    -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)
2725: .    -pc_bddc_adaptive_threshold <0.0>    - when a value different than zero is specified, adaptive selection of constraints is performed on edges and faces (requires deluxe scaling and MUMPS or MKL_PARDISO installed)
2726: -    -pc_bddc_check_level <0>             - set verbosity level of debugging output

2728:    Options for Dirichlet, Neumann or coarse solver can be set using the appropriate options prefix
2729: .vb
2730:       -pc_bddc_dirichlet_
2731:       -pc_bddc_neumann_
2732:       -pc_bddc_coarse_
2733: .ve
2734:    e.g. -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. `PCBDDC` uses by default `KSPPREONLY` and `PCLU`.

2736:    When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified using the options prefix
2737: .vb
2738:       -pc_bddc_dirichlet_lN_
2739:       -pc_bddc_neumann_lN_
2740:       -pc_bddc_coarse_lN_
2741: .ve
2742:    Note that level number ranges from the finest (0) to the coarsest (N).
2743:    In order to specify options for the `PCBDDC` operators at the coarser levels (and not for the solvers), prepend -pc_bddc_coarse_ or -pc_bddc_coarse_l
2744:    to the option, e.g.
2745: .vb
2746:      -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3
2747: .ve
2748:    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

2750:    Level: intermediate

2752: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `MATIS`, `KSPFETIDP`, `PCLU`, `PCGAMG`, `PC`, `PCBDDCSetLocalAdjacencyGraph()`, `PCBDDCSetDofsSplitting()`,
2753:           `PCBDDCSetDirichletBoundaries()`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCSetPrimalVerticesIS()`, `MatNullSpace`, `MatSetNearNullSpace()`,
2754:           `PCBDDCSetChangeOfBasisMat()`, `PCBDDCCreateFETIDPOperators()`, `PCNN`
2755: M*/

2757: PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc)
2758: {
2759:   PC_BDDC *pcbddc;

2761:   PetscFunctionBegin;
2762:   PetscCall(PetscNew(&pcbddc));
2763:   pc->data = pcbddc;

2765:   PetscCall(PCISInitialize(pc));

2767:   /* create local graph structure */
2768:   PetscCall(PCBDDCGraphCreate(&pcbddc->mat_graph));

2770:   /* BDDC nonzero defaults */
2771:   pcbddc->use_nnsp                  = PETSC_TRUE;
2772:   pcbddc->use_local_adj             = PETSC_TRUE;
2773:   pcbddc->use_vertices              = PETSC_TRUE;
2774:   pcbddc->use_edges                 = PETSC_TRUE;
2775:   pcbddc->symmetric_primal          = PETSC_TRUE;
2776:   pcbddc->vertex_size               = 1;
2777:   pcbddc->recompute_topography      = PETSC_TRUE;
2778:   pcbddc->coarse_size               = -1;
2779:   pcbddc->use_exact_dirichlet_trick = PETSC_TRUE;
2780:   pcbddc->coarsening_ratio          = 8;
2781:   pcbddc->coarse_eqs_per_proc       = 1;
2782:   pcbddc->benign_compute_correction = PETSC_TRUE;
2783:   pcbddc->nedfield                  = -1;
2784:   pcbddc->nedglobal                 = PETSC_TRUE;
2785:   pcbddc->graphmaxcount             = PETSC_INT_MAX;
2786:   pcbddc->sub_schurs_layers         = -1;
2787:   pcbddc->adaptive_threshold[0]     = 0.0;
2788:   pcbddc->adaptive_threshold[1]     = 0.0;

2790:   /* function pointers */
2791:   pc->ops->apply               = PCApply_BDDC;
2792:   pc->ops->applytranspose      = PCApplyTranspose_BDDC;
2793:   pc->ops->setup               = PCSetUp_BDDC;
2794:   pc->ops->destroy             = PCDestroy_BDDC;
2795:   pc->ops->setfromoptions      = PCSetFromOptions_BDDC;
2796:   pc->ops->view                = PCView_BDDC;
2797:   pc->ops->applyrichardson     = NULL;
2798:   pc->ops->applysymmetricleft  = NULL;
2799:   pc->ops->applysymmetricright = NULL;
2800:   pc->ops->presolve            = PCPreSolve_BDDC;
2801:   pc->ops->postsolve           = PCPostSolve_BDDC;
2802:   pc->ops->reset               = PCReset_BDDC;

2804:   /* composing function */
2805:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDiscreteGradient_C", PCBDDCSetDiscreteGradient_BDDC));
2806:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDivergenceMat_C", PCBDDCSetDivergenceMat_BDDC));
2807:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetChangeOfBasisMat_C", PCBDDCSetChangeOfBasisMat_BDDC));
2808:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesLocalIS_C", PCBDDCSetPrimalVerticesLocalIS_BDDC));
2809:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesIS_C", PCBDDCSetPrimalVerticesIS_BDDC));
2810:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesLocalIS_C", PCBDDCGetPrimalVerticesLocalIS_BDDC));
2811:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesIS_C", PCBDDCGetPrimalVerticesIS_BDDC));
2812:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetCoarseningRatio_C", PCBDDCSetCoarseningRatio_BDDC));
2813:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevel_C", PCBDDCSetLevel_BDDC));
2814:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetUseExactDirichlet_C", PCBDDCSetUseExactDirichlet_BDDC));
2815:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevels_C", PCBDDCSetLevels_BDDC));
2816:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundaries_C", PCBDDCSetDirichletBoundaries_BDDC));
2817:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundariesLocal_C", PCBDDCSetDirichletBoundariesLocal_BDDC));
2818:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundaries_C", PCBDDCSetNeumannBoundaries_BDDC));
2819:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundariesLocal_C", PCBDDCSetNeumannBoundariesLocal_BDDC));
2820:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundaries_C", PCBDDCGetDirichletBoundaries_BDDC));
2821:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundariesLocal_C", PCBDDCGetDirichletBoundariesLocal_BDDC));
2822:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundaries_C", PCBDDCGetNeumannBoundaries_BDDC));
2823:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundariesLocal_C", PCBDDCGetNeumannBoundariesLocal_BDDC));
2824:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplitting_C", PCBDDCSetDofsSplitting_BDDC));
2825:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplittingLocal_C", PCBDDCSetDofsSplittingLocal_BDDC));
2826:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", PCBDDCSetLocalAdjacencyGraph_BDDC));
2827:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCCreateFETIDPOperators_C", PCBDDCCreateFETIDPOperators_BDDC));
2828:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetRHS_C", PCBDDCMatFETIDPGetRHS_BDDC));
2829:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetSolution_C", PCBDDCMatFETIDPGetSolution_BDDC));
2830:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", PCPreSolveChangeRHS_BDDC));
2831:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_BDDC));
2832:   PetscFunctionReturn(PETSC_SUCCESS);
2833: }

2835: /*@C
2836:   PCBDDCInitializePackage - This function initializes everything in the `PCBDDC` package. It is called
2837:   from `PCInitializePackage()`.

2839:   Level: developer

2841: .seealso: [](ch_ksp), `PetscInitialize()`, `PCBDDCFinalizePackage()`
2842: @*/
2843: PetscErrorCode PCBDDCInitializePackage(void)
2844: {
2845:   int i;

2847:   PetscFunctionBegin;
2848:   if (PCBDDCPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
2849:   PCBDDCPackageInitialized = PETSC_TRUE;
2850:   PetscCall(PetscRegisterFinalize(PCBDDCFinalizePackage));

2852:   /* general events */
2853:   PetscCall(PetscLogEventRegister("PCBDDCTopo", PC_CLASSID, &PC_BDDC_Topology[0]));
2854:   PetscCall(PetscLogEventRegister("PCBDDCLKSP", PC_CLASSID, &PC_BDDC_LocalSolvers[0]));
2855:   PetscCall(PetscLogEventRegister("PCBDDCLWor", PC_CLASSID, &PC_BDDC_LocalWork[0]));
2856:   PetscCall(PetscLogEventRegister("PCBDDCCorr", PC_CLASSID, &PC_BDDC_CorrectionSetUp[0]));
2857:   PetscCall(PetscLogEventRegister("PCBDDCASet", PC_CLASSID, &PC_BDDC_ApproxSetUp[0]));
2858:   PetscCall(PetscLogEventRegister("PCBDDCAApp", PC_CLASSID, &PC_BDDC_ApproxApply[0]));
2859:   PetscCall(PetscLogEventRegister("PCBDDCCSet", PC_CLASSID, &PC_BDDC_CoarseSetUp[0]));
2860:   PetscCall(PetscLogEventRegister("PCBDDCCKSP", PC_CLASSID, &PC_BDDC_CoarseSolver[0]));
2861:   PetscCall(PetscLogEventRegister("PCBDDCAdap", PC_CLASSID, &PC_BDDC_AdaptiveSetUp[0]));
2862:   PetscCall(PetscLogEventRegister("PCBDDCScal", PC_CLASSID, &PC_BDDC_Scaling[0]));
2863:   PetscCall(PetscLogEventRegister("PCBDDCSchr", PC_CLASSID, &PC_BDDC_Schurs[0]));
2864:   PetscCall(PetscLogEventRegister("PCBDDCDirS", PC_CLASSID, &PC_BDDC_Solves[0][0]));
2865:   PetscCall(PetscLogEventRegister("PCBDDCNeuS", PC_CLASSID, &PC_BDDC_Solves[0][1]));
2866:   PetscCall(PetscLogEventRegister("PCBDDCCoaS", PC_CLASSID, &PC_BDDC_Solves[0][2]));
2867:   for (i = 1; i < PETSC_PCBDDC_MAXLEVELS; i++) {
2868:     char ename[32];

2870:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCTopo l%02d", i));
2871:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Topology[i]));
2872:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCLKSP l%02d", i));
2873:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_LocalSolvers[i]));
2874:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCLWor l%02d", i));
2875:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_LocalWork[i]));
2876:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCorr l%02d", i));
2877:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CorrectionSetUp[i]));
2878:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCASet l%02d", i));
2879:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_ApproxSetUp[i]));
2880:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCAApp l%02d", i));
2881:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_ApproxApply[i]));
2882:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCSet l%02d", i));
2883:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CoarseSetUp[i]));
2884:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCKSP l%02d", i));
2885:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CoarseSolver[i]));
2886:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCAdap l%02d", i));
2887:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_AdaptiveSetUp[i]));
2888:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCScal l%02d", i));
2889:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Scaling[i]));
2890:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCSchr l%02d", i));
2891:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Schurs[i]));
2892:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCDirS l%02d", i));
2893:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][0]));
2894:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCNeuS l%02d", i));
2895:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][1]));
2896:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCoaS l%02d", i));
2897:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][2]));
2898:   }
2899:   PetscFunctionReturn(PETSC_SUCCESS);
2900: }

2902: /*@C
2903:   PCBDDCFinalizePackage - This function frees everything from the `PCBDDC` package. It is
2904:   called from `PetscFinalize()` automatically.

2906:   Level: developer

2908: .seealso: [](ch_ksp), `PetscFinalize()`, `PCBDDCInitializePackage()`
2909: @*/
2910: PetscErrorCode PCBDDCFinalizePackage(void)
2911: {
2912:   PetscFunctionBegin;
2913:   PCBDDCPackageInitialized = PETSC_FALSE;
2914:   PetscFunctionReturn(PETSC_SUCCESS);
2915: }