Actual source code: bddc.c

  1: /* TODOLIST

  3:    Solvers
  4:    - Add support for cholesky for coarse solver (similar to local solvers)
  5:    - Propagate ksp prefixes for solvers to mat objects?

  7:    User interface
  8:    - ** DM attached to pc?

 10:    Debugging output
 11:    - * Better management of verbosity levels of debugging output

 13:    Extra
 14:    - *** Is it possible to work with PCBDDCGraph on boundary indices only (less memory consumed)?
 15:    - BDDC with MG framework?

 17:    MATIS related operations contained in BDDC code
 18:    - Provide general case for subassembling

 20: */

 22: #include <petsc/private/pcbddcimpl.h>
 23: #include <petsc/private/pcbddcprivateimpl.h>
 24: #include <petscblaslapack.h>

 26: static PetscBool PCBDDCPackageInitialized = PETSC_FALSE;

 28: static PetscBool  cited      = PETSC_FALSE;
 29: static const char citation[] = "@article{ZampiniPCBDDC,\n"
 30:                                "author = {Stefano Zampini},\n"
 31:                                "title = {{PCBDDC}: A Class of Robust Dual-Primal Methods in {PETS}c},\n"
 32:                                "journal = {SIAM Journal on Scientific Computing},\n"
 33:                                "volume = {38},\n"
 34:                                "number = {5},\n"
 35:                                "pages = {S282-S306},\n"
 36:                                "year = {2016},\n"
 37:                                "doi = {10.1137/15M1025785},\n"
 38:                                "URL = {http://dx.doi.org/10.1137/15M1025785},\n"
 39:                                "eprint = {http://dx.doi.org/10.1137/15M1025785}\n"
 40:                                "}\n";

 42: PetscLogEvent PC_BDDC_Topology[PETSC_PCBDDC_MAXLEVELS];
 43: PetscLogEvent PC_BDDC_LocalSolvers[PETSC_PCBDDC_MAXLEVELS];
 44: PetscLogEvent PC_BDDC_LocalWork[PETSC_PCBDDC_MAXLEVELS];
 45: PetscLogEvent PC_BDDC_CorrectionSetUp[PETSC_PCBDDC_MAXLEVELS];
 46: PetscLogEvent PC_BDDC_ApproxSetUp[PETSC_PCBDDC_MAXLEVELS];
 47: PetscLogEvent PC_BDDC_ApproxApply[PETSC_PCBDDC_MAXLEVELS];
 48: PetscLogEvent PC_BDDC_CoarseSetUp[PETSC_PCBDDC_MAXLEVELS];
 49: PetscLogEvent PC_BDDC_CoarseSolver[PETSC_PCBDDC_MAXLEVELS];
 50: PetscLogEvent PC_BDDC_AdaptiveSetUp[PETSC_PCBDDC_MAXLEVELS];
 51: PetscLogEvent PC_BDDC_Scaling[PETSC_PCBDDC_MAXLEVELS];
 52: PetscLogEvent PC_BDDC_Schurs[PETSC_PCBDDC_MAXLEVELS];
 53: PetscLogEvent PC_BDDC_Solves[PETSC_PCBDDC_MAXLEVELS][3];

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

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

 59: static PetscErrorCode PCSetFromOptions_BDDC(PC pc, PetscOptionItems *PetscOptionsObject)
 60: {
 61:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
 62:   PetscInt nt, i;

 64:   PetscFunctionBegin;
 65:   PetscOptionsHeadBegin(PetscOptionsObject, "BDDC options");
 66:   /* Verbose debugging */
 67:   PetscCall(PetscOptionsInt("-pc_bddc_check_level", "Verbose output for PCBDDC (intended for debug)", "none", pcbddc->dbg_flag, &pcbddc->dbg_flag, NULL));
 68:   /* Approximate solvers */
 69:   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));
 70:   if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_DIRICHLET) {
 71:     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));
 72:     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));
 73:   } else {
 74:     /* This flag is needed/implied by lumping */
 75:     pcbddc->switch_static = PETSC_TRUE;
 76:   }
 77:   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));
 78:   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));
 79:   /* Primal space customization */
 80:   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));
 81:   PetscCall(PetscOptionsInt("-pc_bddc_graph_maxcount", "Maximum number of shared subdomains for a connected component", "none", pcbddc->graphmaxcount, &pcbddc->graphmaxcount, NULL));
 82:   PetscCall(PetscOptionsBool("-pc_bddc_corner_selection", "Activates face-based corner selection", "none", pcbddc->corner_selection, &pcbddc->corner_selection, NULL));
 83:   PetscCall(PetscOptionsBool("-pc_bddc_use_vertices", "Use or not corner dofs in coarse space", "none", pcbddc->use_vertices, &pcbddc->use_vertices, NULL));
 84:   PetscCall(PetscOptionsBool("-pc_bddc_use_edges", "Use or not edge constraints in coarse space", "none", pcbddc->use_edges, &pcbddc->use_edges, NULL));
 85:   PetscCall(PetscOptionsBool("-pc_bddc_use_faces", "Use or not face constraints in coarse space", "none", pcbddc->use_faces, &pcbddc->use_faces, NULL));
 86:   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));
 87:   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));
 88:   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));
 89:   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));
 90:   /* Change of basis */
 91:   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));
 92:   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));
 93:   if (!pcbddc->use_change_of_basis) pcbddc->use_change_on_faces = PETSC_FALSE;
 94:   /* Switch between M_2 (default) and M_3 preconditioners (as defined by C. Dohrmann in the ref. article) */
 95:   PetscCall(PetscOptionsBool("-pc_bddc_switch_static", "Switch on static condensation ops around the interface preconditioner", "none", pcbddc->switch_static, &pcbddc->switch_static, NULL));
 96:   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));
 97:   i = pcbddc->coarsening_ratio;
 98:   PetscCall(PetscOptionsInt("-pc_bddc_coarsening_ratio", "Set coarsening ratio used in multilevel coarsening", "PCBDDCSetCoarseningRatio", i, &i, NULL));
 99:   PetscCall(PCBDDCSetCoarseningRatio(pc, i));
100:   i = pcbddc->max_levels;
101:   PetscCall(PetscOptionsInt("-pc_bddc_levels", "Set maximum number of levels for multilevel", "PCBDDCSetLevels", i, &i, NULL));
102:   PetscCall(PCBDDCSetLevels(pc, i));
103:   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));
104:   PetscCall(PetscOptionsBool("-pc_bddc_use_coarse_estimates", "Use estimated eigenvalues for coarse problem", "none", pcbddc->use_coarse_estimates, &pcbddc->use_coarse_estimates, NULL));
105:   PetscCall(PetscOptionsBool("-pc_bddc_use_deluxe_scaling", "Use deluxe scaling for BDDC", "none", pcbddc->use_deluxe_scaling, &pcbddc->use_deluxe_scaling, NULL));
106:   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));
107:   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));
108:   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));
109:   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));
110:   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));
111:   PetscCall(PetscOptionsBool("-pc_bddc_deluxe_singlemat", "Collapse deluxe operators", "none", pcbddc->deluxe_singlemat, &pcbddc->deluxe_singlemat, NULL));
112:   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));
113:   nt = 2;
114:   PetscCall(PetscOptionsRealArray("-pc_bddc_adaptive_threshold", "Thresholds to be used for adaptive selection of constraints", "none", pcbddc->adaptive_threshold, &nt, NULL));
115:   if (nt == 1) pcbddc->adaptive_threshold[1] = pcbddc->adaptive_threshold[0];
116:   PetscCall(PetscOptionsInt("-pc_bddc_adaptive_nmin", "Minimum number of constraints per connected components", "none", pcbddc->adaptive_nmin, &pcbddc->adaptive_nmin, NULL));
117:   PetscCall(PetscOptionsInt("-pc_bddc_adaptive_nmax", "Maximum number of constraints per connected components", "none", pcbddc->adaptive_nmax, &pcbddc->adaptive_nmax, NULL));
118:   PetscCall(PetscOptionsBool("-pc_bddc_symmetric", "Symmetric computation of primal basis functions", "none", pcbddc->symmetric_primal, &pcbddc->symmetric_primal, NULL));
119:   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));
120:   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));
121:   PetscCall(PetscOptionsBool("-pc_bddc_benign_change", "Compute the pressure change of basis explicitly", "none", pcbddc->benign_change_explicit, &pcbddc->benign_change_explicit, NULL));
122:   PetscCall(PetscOptionsBool("-pc_bddc_benign_compute_correction", "Compute the benign correction during PreSolve", "none", pcbddc->benign_compute_correction, &pcbddc->benign_compute_correction, NULL));
123:   PetscCall(PetscOptionsBool("-pc_bddc_nonetflux", "Automatic computation of no-net-flux quadrature weights", "none", pcbddc->compute_nonetflux, &pcbddc->compute_nonetflux, NULL));
124:   PetscCall(PetscOptionsBool("-pc_bddc_detect_disconnected", "Detects disconnected subdomains", "none", pcbddc->detect_disconnected, &pcbddc->detect_disconnected, NULL));
125:   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));
126:   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));
127:   PetscOptionsHeadEnd();
128:   PetscFunctionReturn(PETSC_SUCCESS);
129: }

131: static PetscErrorCode PCView_BDDC(PC pc, PetscViewer viewer)
132: {
133:   PC_BDDC     *pcbddc = (PC_BDDC *)pc->data;
134:   PC_IS       *pcis   = (PC_IS *)pc->data;
135:   PetscBool    isascii;
136:   PetscSubcomm subcomm;
137:   PetscViewer  subviewer;

139:   PetscFunctionBegin;
140:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
141:   /* ASCII viewer */
142:   if (isascii) {
143:     PetscMPIInt color, rank, size;
144:     PetscInt64  loc[7], gsum[6], gmax[6], gmin[6], totbenign;
145:     PetscScalar interface_size;
146:     PetscReal   ratio1 = 0., ratio2 = 0.;
147:     Vec         counter;

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

199:     /* compute interface size */
200:     PetscCall(VecSet(pcis->vec1_B, 1.0));
201:     PetscCall(MatCreateVecs(pc->pmat, &counter, NULL));
202:     PetscCall(VecSet(counter, 0.0));
203:     PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, counter, INSERT_VALUES, SCATTER_REVERSE));
204:     PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, counter, INSERT_VALUES, SCATTER_REVERSE));
205:     PetscCall(VecSum(counter, &interface_size));
206:     PetscCall(VecDestroy(&counter));

208:     /* compute some statistics on the domain decomposition */
209:     gsum[0] = 1;
210:     gsum[1] = gsum[2] = gsum[3] = gsum[4] = gsum[5] = 0;
211:     loc[0]                                          = !!pcis->n;
212:     loc[1]                                          = pcis->n - pcis->n_B;
213:     loc[2]                                          = pcis->n_B;
214:     loc[3]                                          = pcbddc->local_primal_size;
215:     loc[4]                                          = pcis->n;
216:     loc[5]                                          = pcbddc->n_local_subs > 0 ? pcbddc->n_local_subs : (pcis->n ? 1 : 0);
217:     loc[6]                                          = pcbddc->benign_n;
218:     PetscCallMPI(MPI_Reduce(loc, gsum, 6, MPIU_INT64, MPI_SUM, 0, PetscObjectComm((PetscObject)pc)));
219:     if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = loc[5] = -1;
220:     PetscCallMPI(MPI_Reduce(loc, gmax, 6, MPIU_INT64, MPI_MAX, 0, PetscObjectComm((PetscObject)pc)));
221:     if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = loc[5] = PETSC_MAX_INT;
222:     PetscCallMPI(MPI_Reduce(loc, gmin, 6, MPIU_INT64, MPI_MIN, 0, PetscObjectComm((PetscObject)pc)));
223:     PetscCallMPI(MPI_Reduce(&loc[6], &totbenign, 1, MPIU_INT64, MPI_SUM, 0, PetscObjectComm((PetscObject)pc)));
224:     if (pcbddc->coarse_size) {
225:       ratio1 = pc->pmat->rmap->N / (1. * pcbddc->coarse_size);
226:       ratio2 = PetscRealPart(interface_size) / pcbddc->coarse_size;
227:     }
228:     PetscCall(PetscViewerASCIIPrintf(viewer, "********************************** STATISTICS AT LEVEL %" PetscInt_FMT " **********************************\n", pcbddc->current_level));
229:     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));
230:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Coarsening ratios: all/coarse %" PetscInt_FMT " interface/coarse %" PetscInt_FMT "\n", (PetscInt)ratio1, (PetscInt)ratio2));
231:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Active processes : %" PetscInt_FMT "\n", (PetscInt)gsum[0]));
232:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Total subdomains : %" PetscInt_FMT "\n", (PetscInt)gsum[5]));
233:     if (pcbddc->benign_have_null) PetscCall(PetscViewerASCIIPrintf(viewer, "  Benign subs      : %" PetscInt_FMT "\n", (PetscInt)totbenign));
234:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Dofs type        :\tMIN\tMAX\tMEAN\n"));
235:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Interior  dofs   :\t%" PetscInt_FMT "\t%" PetscInt_FMT "\t%" PetscInt_FMT "\n", (PetscInt)gmin[1], (PetscInt)gmax[1], (PetscInt)(gsum[1] / gsum[0])));
236:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Interface dofs   :\t%" PetscInt_FMT "\t%" PetscInt_FMT "\t%" PetscInt_FMT "\n", (PetscInt)gmin[2], (PetscInt)gmax[2], (PetscInt)(gsum[2] / gsum[0])));
237:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Primal    dofs   :\t%" PetscInt_FMT "\t%" PetscInt_FMT "\t%" PetscInt_FMT "\n", (PetscInt)gmin[3], (PetscInt)gmax[3], (PetscInt)(gsum[3] / gsum[0])));
238:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Local     dofs   :\t%" PetscInt_FMT "\t%" PetscInt_FMT "\t%" PetscInt_FMT "\n", (PetscInt)gmin[4], (PetscInt)gmax[4], (PetscInt)(gsum[4] / gsum[0])));
239:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Local     subs   :\t%" PetscInt_FMT "\t%" PetscInt_FMT "\n", (PetscInt)gmin[5], (PetscInt)gmax[5]));
240:     PetscCall(PetscViewerFlush(viewer));

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

244:     /* local solvers */
245:     PetscCall(PetscViewerGetSubViewer(viewer, PetscObjectComm((PetscObject)pcbddc->ksp_D), &subviewer));
246:     if (rank == 0) {
247:       PetscCall(PetscViewerASCIIPrintf(subviewer, "--- Interior solver (rank 0)\n"));
248:       PetscCall(PetscViewerASCIIPushTab(subviewer));
249:       PetscCall(KSPView(pcbddc->ksp_D, subviewer));
250:       PetscCall(PetscViewerASCIIPopTab(subviewer));
251:       PetscCall(PetscViewerASCIIPrintf(subviewer, "--- Correction solver (rank 0)\n"));
252:       PetscCall(PetscViewerASCIIPushTab(subviewer));
253:       PetscCall(KSPView(pcbddc->ksp_R, subviewer));
254:       PetscCall(PetscViewerASCIIPopTab(subviewer));
255:       PetscCall(PetscViewerFlush(subviewer));
256:     }
257:     PetscCall(PetscViewerRestoreSubViewer(viewer, PetscObjectComm((PetscObject)pcbddc->ksp_D), &subviewer));
258:     PetscCall(PetscViewerFlush(viewer));

260:     /* the coarse problem can be handled by a different communicator */
261:     if (pcbddc->coarse_ksp) color = 1;
262:     else color = 0;
263:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
264:     PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)pc), &subcomm));
265:     PetscCall(PetscSubcommSetNumber(subcomm, PetscMin(size, 2)));
266:     PetscCall(PetscSubcommSetTypeGeneral(subcomm, color, rank));
267:     PetscCall(PetscViewerGetSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer));
268:     if (color == 1) {
269:       PetscCall(PetscViewerASCIIPrintf(subviewer, "--- Coarse solver\n"));
270:       PetscCall(PetscViewerASCIIPushTab(subviewer));
271:       PetscCall(KSPView(pcbddc->coarse_ksp, subviewer));
272:       PetscCall(PetscViewerASCIIPopTab(subviewer));
273:       PetscCall(PetscViewerFlush(subviewer));
274:     }
275:     PetscCall(PetscViewerRestoreSubViewer(viewer, PetscSubcommChild(subcomm), &subviewer));
276:     PetscCall(PetscSubcommDestroy(&subcomm));
277:     PetscCall(PetscViewerFlush(viewer));
278:   }
279:   PetscFunctionReturn(PETSC_SUCCESS);
280: }

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

286:   PetscFunctionBegin;
287:   PetscCall(PetscObjectReference((PetscObject)G));
288:   PetscCall(MatDestroy(&pcbddc->discretegradient));
289:   pcbddc->discretegradient = G;
290:   pcbddc->nedorder         = order > 0 ? order : -order;
291:   pcbddc->nedfield         = field;
292:   pcbddc->nedglobal        = global;
293:   pcbddc->conforming       = conforming;
294:   PetscFunctionReturn(PETSC_SUCCESS);
295: }

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

300:   Collective

302:   Input Parameters:
303: + pc         - the preconditioning context
304: . G          - the discrete gradient matrix (in `MATAIJ` format)
305: . order      - the order of the Nedelec space (1 for the lowest order)
306: . field      - the field id of the Nedelec dofs (not used if the fields have not been specified)
307: . global     - the type of global ordering for the rows of `G`
308: - conforming - whether the mesh is conforming or not

310:   Level: advanced

312:   Note:
313:   The discrete gradient matrix `G` is used to analyze the subdomain edges, and it should not contain any zero entry.
314:   For variable order spaces, the order should be set to zero.
315:   If `global` is `PETSC_TRUE`, the rows of `G` should be given in global ordering for the whole dofs;
316:   if `PETSC_FALSE`, the ordering should be global for the Nedelec field.
317:   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
318:   and geid the one for the Nedelec field.

320: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDofsSplitting()`, `PCBDDCSetDofsSplittingLocal()`, `MATAIJ`, `PCBDDCSetDivergenceMat()`
321: @*/
322: PetscErrorCode PCBDDCSetDiscreteGradient(PC pc, Mat G, PetscInt order, PetscInt field, PetscBool global, PetscBool conforming)
323: {
324:   PetscFunctionBegin;
331:   PetscCheckSameComm(pc, 1, G, 2);
332:   PetscTryMethod(pc, "PCBDDCSetDiscreteGradient_C", (PC, Mat, PetscInt, PetscInt, PetscBool, PetscBool), (pc, G, order, field, global, conforming));
333:   PetscFunctionReturn(PETSC_SUCCESS);
334: }

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

340:   PetscFunctionBegin;
341:   PetscCall(PetscObjectReference((PetscObject)divudotp));
342:   PetscCall(MatDestroy(&pcbddc->divudotp));
343:   pcbddc->divudotp          = divudotp;
344:   pcbddc->divudotp_trans    = trans;
345:   pcbddc->compute_nonetflux = PETSC_TRUE;
346:   if (vl2l) {
347:     PetscCall(PetscObjectReference((PetscObject)vl2l));
348:     PetscCall(ISDestroy(&pcbddc->divudotp_vl2l));
349:     pcbddc->divudotp_vl2l = vl2l;
350:   }
351:   PetscFunctionReturn(PETSC_SUCCESS);
352: }

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

357:   Collective

359:   Input Parameters:
360: + pc       - the preconditioning context
361: . divudotp - the matrix (must be of type `MATIS`)
362: . trans    - if `PETSC_FALSE` (resp. `PETSC_TRUE`), then pressures are in the test (trial) space and velocities are in the trial (test) space.
363: - vl2l     - optional index set describing the local (wrt the local matrix in `divudotp`) to local (wrt the local matrix
364:    in the preconditioning matrix) map for the velocities

366:   Level: advanced

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

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

373: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDiscreteGradient()`
374: @*/
375: PetscErrorCode PCBDDCSetDivergenceMat(PC pc, Mat divudotp, PetscBool trans, IS vl2l)
376: {
377:   PetscBool ismatis;

379:   PetscFunctionBegin;
382:   PetscCheckSameComm(pc, 1, divudotp, 2);
385:   PetscCall(PetscObjectTypeCompare((PetscObject)divudotp, MATIS, &ismatis));
386:   PetscCheck(ismatis, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Divergence matrix needs to be of type MATIS");
387:   PetscTryMethod(pc, "PCBDDCSetDivergenceMat_C", (PC, Mat, PetscBool, IS), (pc, divudotp, trans, vl2l));
388:   PetscFunctionReturn(PETSC_SUCCESS);
389: }

391: static PetscErrorCode PCBDDCSetChangeOfBasisMat_BDDC(PC pc, Mat change, PetscBool interior)
392: {
393:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

395:   PetscFunctionBegin;
396:   PetscCall(PetscObjectReference((PetscObject)change));
397:   PetscCall(MatDestroy(&pcbddc->user_ChangeOfBasisMatrix));
398:   pcbddc->user_ChangeOfBasisMatrix = change;
399:   pcbddc->change_interior          = interior;
400:   PetscFunctionReturn(PETSC_SUCCESS);
401: }

403: /*@
404:   PCBDDCSetChangeOfBasisMat - Set user defined change of basis for dofs

406:   Collective

408:   Input Parameters:
409: + pc       - the preconditioning context
410: . change   - the change of basis matrix
411: - interior - whether or not the change of basis modifies interior dofs

413:   Level: intermediate

415: .seealso: [](ch_ksp), `PCBDDC`
416: @*/
417: PetscErrorCode PCBDDCSetChangeOfBasisMat(PC pc, Mat change, PetscBool interior)
418: {
419:   PetscFunctionBegin;
422:   PetscCheckSameComm(pc, 1, change, 2);
423:   if (pc->mat) {
424:     PetscInt rows_c, cols_c, rows, cols;
425:     PetscCall(MatGetSize(pc->mat, &rows, &cols));
426:     PetscCall(MatGetSize(change, &rows_c, &cols_c));
427:     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);
428:     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);
429:     PetscCall(MatGetLocalSize(pc->mat, &rows, &cols));
430:     PetscCall(MatGetLocalSize(change, &rows_c, &cols_c));
431:     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);
432:     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);
433:   }
434:   PetscTryMethod(pc, "PCBDDCSetChangeOfBasisMat_C", (PC, Mat, PetscBool), (pc, change, interior));
435:   PetscFunctionReturn(PETSC_SUCCESS);
436: }

438: static PetscErrorCode PCBDDCSetPrimalVerticesIS_BDDC(PC pc, IS PrimalVertices)
439: {
440:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
441:   PetscBool isequal = PETSC_FALSE;

443:   PetscFunctionBegin;
444:   PetscCall(PetscObjectReference((PetscObject)PrimalVertices));
445:   if (pcbddc->user_primal_vertices) PetscCall(ISEqual(PrimalVertices, pcbddc->user_primal_vertices, &isequal));
446:   PetscCall(ISDestroy(&pcbddc->user_primal_vertices));
447:   PetscCall(ISDestroy(&pcbddc->user_primal_vertices_local));
448:   pcbddc->user_primal_vertices = PrimalVertices;
449:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
450:   PetscFunctionReturn(PETSC_SUCCESS);
451: }

453: /*@
454:   PCBDDCSetPrimalVerticesIS - Set additional user defined primal vertices in `PCBDDC`

456:   Collective

458:   Input Parameters:
459: + pc             - the preconditioning context
460: - PrimalVertices - index set of primal vertices in global numbering (can be empty)

462:   Level: intermediate

464:   Note:
465:   Any process can list any global node

467: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
468: @*/
469: PetscErrorCode PCBDDCSetPrimalVerticesIS(PC pc, IS PrimalVertices)
470: {
471:   PetscFunctionBegin;
474:   PetscCheckSameComm(pc, 1, PrimalVertices, 2);
475:   PetscTryMethod(pc, "PCBDDCSetPrimalVerticesIS_C", (PC, IS), (pc, PrimalVertices));
476:   PetscFunctionReturn(PETSC_SUCCESS);
477: }

479: static PetscErrorCode PCBDDCGetPrimalVerticesIS_BDDC(PC pc, IS *is)
480: {
481:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

483:   PetscFunctionBegin;
484:   *is = pcbddc->user_primal_vertices;
485:   PetscFunctionReturn(PETSC_SUCCESS);
486: }

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

491:   Collective

493:   Input Parameter:
494: . pc - the preconditioning context

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

499:   Level: intermediate

501: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
502: @*/
503: PetscErrorCode PCBDDCGetPrimalVerticesIS(PC pc, IS *is)
504: {
505:   PetscFunctionBegin;
507:   PetscAssertPointer(is, 2);
508:   PetscUseMethod(pc, "PCBDDCGetPrimalVerticesIS_C", (PC, IS *), (pc, is));
509:   PetscFunctionReturn(PETSC_SUCCESS);
510: }

512: static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices)
513: {
514:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
515:   PetscBool isequal = PETSC_FALSE;

517:   PetscFunctionBegin;
518:   PetscCall(PetscObjectReference((PetscObject)PrimalVertices));
519:   if (pcbddc->user_primal_vertices_local) PetscCall(ISEqual(PrimalVertices, pcbddc->user_primal_vertices_local, &isequal));
520:   PetscCall(ISDestroy(&pcbddc->user_primal_vertices));
521:   PetscCall(ISDestroy(&pcbddc->user_primal_vertices_local));
522:   pcbddc->user_primal_vertices_local = PrimalVertices;
523:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
524:   PetscFunctionReturn(PETSC_SUCCESS);
525: }

527: /*@
528:   PCBDDCSetPrimalVerticesLocalIS - Set additional user defined primal vertices in `PCBDDC`

530:   Collective

532:   Input Parameters:
533: + pc             - the preconditioning context
534: - PrimalVertices - index set of primal vertices in local numbering (can be empty)

536:   Level: intermediate

538: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
539: @*/
540: PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices)
541: {
542:   PetscFunctionBegin;
545:   PetscCheckSameComm(pc, 1, PrimalVertices, 2);
546:   PetscTryMethod(pc, "PCBDDCSetPrimalVerticesLocalIS_C", (PC, IS), (pc, PrimalVertices));
547:   PetscFunctionReturn(PETSC_SUCCESS);
548: }

550: static PetscErrorCode PCBDDCGetPrimalVerticesLocalIS_BDDC(PC pc, IS *is)
551: {
552:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

554:   PetscFunctionBegin;
555:   *is = pcbddc->user_primal_vertices_local;
556:   PetscFunctionReturn(PETSC_SUCCESS);
557: }

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

562:   Collective

564:   Input Parameter:
565: . pc - the preconditioning context

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

570:   Level: intermediate

572: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`
573: @*/
574: PetscErrorCode PCBDDCGetPrimalVerticesLocalIS(PC pc, IS *is)
575: {
576:   PetscFunctionBegin;
578:   PetscAssertPointer(is, 2);
579:   PetscUseMethod(pc, "PCBDDCGetPrimalVerticesLocalIS_C", (PC, IS *), (pc, is));
580:   PetscFunctionReturn(PETSC_SUCCESS);
581: }

583: static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc, PetscInt k)
584: {
585:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

587:   PetscFunctionBegin;
588:   pcbddc->coarsening_ratio = k;
589:   PetscFunctionReturn(PETSC_SUCCESS);
590: }

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

595:   Logically Collective

597:   Input Parameters:
598: + pc - the preconditioning context
599: - k  - coarsening ratio (H/h at the coarser level)

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

604:   Level: intermediate

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

609: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetLevels()`
610: @*/
611: PetscErrorCode PCBDDCSetCoarseningRatio(PC pc, PetscInt k)
612: {
613:   PetscFunctionBegin;
616:   PetscTryMethod(pc, "PCBDDCSetCoarseningRatio_C", (PC, PetscInt), (pc, k));
617:   PetscFunctionReturn(PETSC_SUCCESS);
618: }

620: /* The following functions (PCBDDCSetUseExactDirichlet PCBDDCSetLevel) are not public */
621: static PetscErrorCode PCBDDCSetUseExactDirichlet_BDDC(PC pc, PetscBool flg)
622: {
623:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

625:   PetscFunctionBegin;
626:   pcbddc->use_exact_dirichlet_trick = flg;
627:   PetscFunctionReturn(PETSC_SUCCESS);
628: }

630: PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc, PetscBool flg)
631: {
632:   PetscFunctionBegin;
635:   PetscTryMethod(pc, "PCBDDCSetUseExactDirichlet_C", (PC, PetscBool), (pc, flg));
636:   PetscFunctionReturn(PETSC_SUCCESS);
637: }

639: static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc, PetscInt level)
640: {
641:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

643:   PetscFunctionBegin;
644:   pcbddc->current_level = level;
645:   PetscFunctionReturn(PETSC_SUCCESS);
646: }

648: PetscErrorCode PCBDDCSetLevel(PC pc, PetscInt level)
649: {
650:   PetscFunctionBegin;
653:   PetscTryMethod(pc, "PCBDDCSetLevel_C", (PC, PetscInt), (pc, level));
654:   PetscFunctionReturn(PETSC_SUCCESS);
655: }

657: static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc, PetscInt levels)
658: {
659:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

661:   PetscFunctionBegin;
662:   PetscCheck(levels < PETSC_PCBDDC_MAXLEVELS, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Maximum number of additional levels for BDDC is %d", PETSC_PCBDDC_MAXLEVELS - 1);
663:   pcbddc->max_levels = levels;
664:   PetscFunctionReturn(PETSC_SUCCESS);
665: }

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

670:   Logically Collective

672:   Input Parameters:
673: + pc     - the preconditioning context
674: - levels - the maximum number of levels

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

679:   Level: intermediate

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

684: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetCoarseningRatio()`
685: @*/
686: PetscErrorCode PCBDDCSetLevels(PC pc, PetscInt levels)
687: {
688:   PetscFunctionBegin;
691:   PetscTryMethod(pc, "PCBDDCSetLevels_C", (PC, PetscInt), (pc, levels));
692:   PetscFunctionReturn(PETSC_SUCCESS);
693: }

695: static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc, IS DirichletBoundaries)
696: {
697:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
698:   PetscBool isequal = PETSC_FALSE;

700:   PetscFunctionBegin;
701:   PetscCall(PetscObjectReference((PetscObject)DirichletBoundaries));
702:   if (pcbddc->DirichletBoundaries) PetscCall(ISEqual(DirichletBoundaries, pcbddc->DirichletBoundaries, &isequal));
703:   /* last user setting takes precedence -> destroy any other customization */
704:   PetscCall(ISDestroy(&pcbddc->DirichletBoundariesLocal));
705:   PetscCall(ISDestroy(&pcbddc->DirichletBoundaries));
706:   pcbddc->DirichletBoundaries = DirichletBoundaries;
707:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
708:   PetscFunctionReturn(PETSC_SUCCESS);
709: }

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

714:   Collective

716:   Input Parameters:
717: + pc                  - the preconditioning context
718: - DirichletBoundaries - parallel `IS` defining the Dirichlet boundaries

720:   Level: intermediate

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

725: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDirichletBoundariesLocal()`, `MatZeroRows()`, `MatZeroRowsColumns()`
726: @*/
727: PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc, IS DirichletBoundaries)
728: {
729:   PetscFunctionBegin;
732:   PetscCheckSameComm(pc, 1, DirichletBoundaries, 2);
733:   PetscTryMethod(pc, "PCBDDCSetDirichletBoundaries_C", (PC, IS), (pc, DirichletBoundaries));
734:   PetscFunctionReturn(PETSC_SUCCESS);
735: }

737: static PetscErrorCode PCBDDCSetDirichletBoundariesLocal_BDDC(PC pc, IS DirichletBoundaries)
738: {
739:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
740:   PetscBool isequal = PETSC_FALSE;

742:   PetscFunctionBegin;
743:   PetscCall(PetscObjectReference((PetscObject)DirichletBoundaries));
744:   if (pcbddc->DirichletBoundariesLocal) PetscCall(ISEqual(DirichletBoundaries, pcbddc->DirichletBoundariesLocal, &isequal));
745:   /* last user setting takes precedence -> destroy any other customization */
746:   PetscCall(ISDestroy(&pcbddc->DirichletBoundariesLocal));
747:   PetscCall(ISDestroy(&pcbddc->DirichletBoundaries));
748:   pcbddc->DirichletBoundariesLocal = DirichletBoundaries;
749:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
750:   PetscFunctionReturn(PETSC_SUCCESS);
751: }

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

756:   Collective

758:   Input Parameters:
759: + pc                  - the preconditioning context
760: - DirichletBoundaries - parallel `IS` defining the Dirichlet boundaries (in local ordering)

762:   Level: intermediate

764: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDirichletBoundaries()`, `MatZeroRows()`, `MatZeroRowsColumns()`
765: @*/
766: PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc, IS DirichletBoundaries)
767: {
768:   PetscFunctionBegin;
771:   PetscCheckSameComm(pc, 1, DirichletBoundaries, 2);
772:   PetscTryMethod(pc, "PCBDDCSetDirichletBoundariesLocal_C", (PC, IS), (pc, DirichletBoundaries));
773:   PetscFunctionReturn(PETSC_SUCCESS);
774: }

776: static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc, IS NeumannBoundaries)
777: {
778:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
779:   PetscBool isequal = PETSC_FALSE;

781:   PetscFunctionBegin;
782:   PetscCall(PetscObjectReference((PetscObject)NeumannBoundaries));
783:   if (pcbddc->NeumannBoundaries) PetscCall(ISEqual(NeumannBoundaries, pcbddc->NeumannBoundaries, &isequal));
784:   /* last user setting takes precedence -> destroy any other customization */
785:   PetscCall(ISDestroy(&pcbddc->NeumannBoundariesLocal));
786:   PetscCall(ISDestroy(&pcbddc->NeumannBoundaries));
787:   pcbddc->NeumannBoundaries = NeumannBoundaries;
788:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
789:   PetscFunctionReturn(PETSC_SUCCESS);
790: }

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

795:   Collective

797:   Input Parameters:
798: + pc                - the preconditioning context
799: - NeumannBoundaries - parallel `IS` defining the Neumann boundaries

801:   Level: intermediate

803:   Note:
804:   Any process can list any global node

806: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetNeumannBoundariesLocal()`
807: @*/
808: PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc, IS NeumannBoundaries)
809: {
810:   PetscFunctionBegin;
813:   PetscCheckSameComm(pc, 1, NeumannBoundaries, 2);
814:   PetscTryMethod(pc, "PCBDDCSetNeumannBoundaries_C", (PC, IS), (pc, NeumannBoundaries));
815:   PetscFunctionReturn(PETSC_SUCCESS);
816: }

818: static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc, IS NeumannBoundaries)
819: {
820:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
821:   PetscBool isequal = PETSC_FALSE;

823:   PetscFunctionBegin;
824:   PetscCall(PetscObjectReference((PetscObject)NeumannBoundaries));
825:   if (pcbddc->NeumannBoundariesLocal) PetscCall(ISEqual(NeumannBoundaries, pcbddc->NeumannBoundariesLocal, &isequal));
826:   /* last user setting takes precedence -> destroy any other customization */
827:   PetscCall(ISDestroy(&pcbddc->NeumannBoundariesLocal));
828:   PetscCall(ISDestroy(&pcbddc->NeumannBoundaries));
829:   pcbddc->NeumannBoundariesLocal = NeumannBoundaries;
830:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
831:   PetscFunctionReturn(PETSC_SUCCESS);
832: }

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

837:   Collective

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

843:   Level: intermediate

845: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCGetDirichletBoundaries()`
846: @*/
847: PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc, IS NeumannBoundaries)
848: {
849:   PetscFunctionBegin;
852:   PetscCheckSameComm(pc, 1, NeumannBoundaries, 2);
853:   PetscTryMethod(pc, "PCBDDCSetNeumannBoundariesLocal_C", (PC, IS), (pc, NeumannBoundaries));
854:   PetscFunctionReturn(PETSC_SUCCESS);
855: }

857: static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc, IS *DirichletBoundaries)
858: {
859:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

861:   PetscFunctionBegin;
862:   *DirichletBoundaries = pcbddc->DirichletBoundaries;
863:   PetscFunctionReturn(PETSC_SUCCESS);
864: }

866: /*@
867:   PCBDDCGetDirichletBoundaries - Get parallel `IS` for Dirichlet boundaries

869:   Collective

871:   Input Parameter:
872: . pc - the preconditioning context

874:   Output Parameter:
875: . DirichletBoundaries - index set defining the Dirichlet boundaries

877:   Level: intermediate

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

882: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDirichletBoundaries()`
883: @*/
884: PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc, IS *DirichletBoundaries)
885: {
886:   PetscFunctionBegin;
888:   PetscUseMethod(pc, "PCBDDCGetDirichletBoundaries_C", (PC, IS *), (pc, DirichletBoundaries));
889:   PetscFunctionReturn(PETSC_SUCCESS);
890: }

892: static PetscErrorCode PCBDDCGetDirichletBoundariesLocal_BDDC(PC pc, IS *DirichletBoundaries)
893: {
894:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

896:   PetscFunctionBegin;
897:   *DirichletBoundaries = pcbddc->DirichletBoundariesLocal;
898:   PetscFunctionReturn(PETSC_SUCCESS);
899: }

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

904:   Collective

906:   Input Parameter:
907: . pc - the preconditioning context

909:   Output Parameter:
910: . DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries

912:   Level: intermediate

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

919: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCGetDirichletBoundaries()`, `PCBDDCSetDirichletBoundaries()`
920: @*/
921: PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc, IS *DirichletBoundaries)
922: {
923:   PetscFunctionBegin;
925:   PetscUseMethod(pc, "PCBDDCGetDirichletBoundariesLocal_C", (PC, IS *), (pc, DirichletBoundaries));
926:   PetscFunctionReturn(PETSC_SUCCESS);
927: }

929: static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc, IS *NeumannBoundaries)
930: {
931:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

933:   PetscFunctionBegin;
934:   *NeumannBoundaries = pcbddc->NeumannBoundaries;
935:   PetscFunctionReturn(PETSC_SUCCESS);
936: }

938: /*@
939:   PCBDDCGetNeumannBoundaries - Get parallel `IS` for Neumann boundaries

941:   Not Collective

943:   Input Parameter:
944: . pc - the preconditioning context

946:   Output Parameter:
947: . NeumannBoundaries - index set defining the Neumann boundaries

949:   Level: intermediate

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

954: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCGetDirichletBoundaries()`, `PCBDDCSetDirichletBoundaries()`
955: @*/
956: PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc, IS *NeumannBoundaries)
957: {
958:   PetscFunctionBegin;
960:   PetscUseMethod(pc, "PCBDDCGetNeumannBoundaries_C", (PC, IS *), (pc, NeumannBoundaries));
961:   PetscFunctionReturn(PETSC_SUCCESS);
962: }

964: static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc, IS *NeumannBoundaries)
965: {
966:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

968:   PetscFunctionBegin;
969:   *NeumannBoundaries = pcbddc->NeumannBoundariesLocal;
970:   PetscFunctionReturn(PETSC_SUCCESS);
971: }

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

976:   Not Collective

978:   Input Parameter:
979: . pc - the preconditioning context

981:   Output Parameter:
982: . NeumannBoundaries - index set defining the subdomain part of Neumann boundaries

984:   Level: intermediate

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

991: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCSetNeumannBoundariesLocal)`, `PCBDDCGetNeumannBoundaries()`
992: @*/
993: PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc, IS *NeumannBoundaries)
994: {
995:   PetscFunctionBegin;
997:   PetscUseMethod(pc, "PCBDDCGetNeumannBoundariesLocal_C", (PC, IS *), (pc, NeumannBoundaries));
998:   PetscFunctionReturn(PETSC_SUCCESS);
999: }

1001: static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs, const PetscInt xadj[], const PetscInt adjncy[], PetscCopyMode copymode)
1002: {
1003:   PC_BDDC    *pcbddc    = (PC_BDDC *)pc->data;
1004:   PCBDDCGraph mat_graph = pcbddc->mat_graph;
1005:   PetscBool   same_data = PETSC_FALSE;

1007:   PetscFunctionBegin;
1008:   if (!nvtxs) {
1009:     if (copymode == PETSC_OWN_POINTER) {
1010:       PetscCall(PetscFree(xadj));
1011:       PetscCall(PetscFree(adjncy));
1012:     }
1013:     PetscCall(PCBDDCGraphResetCSR(mat_graph));
1014:     PetscFunctionReturn(PETSC_SUCCESS);
1015:   }
1016:   if (mat_graph->nvtxs == nvtxs && mat_graph->freecsr) { /* we own the data */
1017:     if (mat_graph->xadj == xadj && mat_graph->adjncy == adjncy) same_data = PETSC_TRUE;
1018:     if (!same_data && mat_graph->xadj[nvtxs] == xadj[nvtxs]) {
1019:       PetscCall(PetscArraycmp(xadj, mat_graph->xadj, nvtxs + 1, &same_data));
1020:       if (same_data) PetscCall(PetscArraycmp(adjncy, mat_graph->adjncy, xadj[nvtxs], &same_data));
1021:     }
1022:   }
1023:   if (!same_data) {
1024:     /* free old CSR */
1025:     PetscCall(PCBDDCGraphResetCSR(mat_graph));
1026:     /* get CSR into graph structure */
1027:     if (copymode == PETSC_COPY_VALUES) {
1028:       PetscCall(PetscMalloc1(nvtxs + 1, &mat_graph->xadj));
1029:       PetscCall(PetscMalloc1(xadj[nvtxs], &mat_graph->adjncy));
1030:       PetscCall(PetscArraycpy(mat_graph->xadj, xadj, nvtxs + 1));
1031:       PetscCall(PetscArraycpy(mat_graph->adjncy, adjncy, xadj[nvtxs]));
1032:       mat_graph->freecsr = PETSC_TRUE;
1033:     } else if (copymode == PETSC_OWN_POINTER) {
1034:       mat_graph->xadj    = (PetscInt *)xadj;
1035:       mat_graph->adjncy  = (PetscInt *)adjncy;
1036:       mat_graph->freecsr = PETSC_TRUE;
1037:     } else if (copymode == PETSC_USE_POINTER) {
1038:       mat_graph->xadj    = (PetscInt *)xadj;
1039:       mat_graph->adjncy  = (PetscInt *)adjncy;
1040:       mat_graph->freecsr = PETSC_FALSE;
1041:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported copy mode %d", copymode);
1042:     mat_graph->nvtxs_csr         = nvtxs;
1043:     pcbddc->recompute_topography = PETSC_TRUE;
1044:   }
1045:   PetscFunctionReturn(PETSC_SUCCESS);
1046: }

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

1051:   Not collective

1053:   Input Parameters:
1054: + pc       - the preconditioning context.
1055: . nvtxs    - number of local vertices of the graph (i.e., the number of local dofs).
1056: . xadj     - CSR format row pointers for the connectivity of the dofs
1057: . adjncy   - CSR format column pointers for the connectivity of the dofs
1058: - copymode - supported modes are `PETSC_COPY_VALUES`, `PETSC_USE_POINTER` or `PETSC_OWN_POINTER`.

1060:   Level: intermediate

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

1065: .seealso: [](ch_ksp), `PCBDDC`, `PetscCopyMode`
1066: @*/
1067: PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc, PetscInt nvtxs, const PetscInt xadj[], const PetscInt adjncy[], PetscCopyMode copymode)
1068: {
1069:   void (*f)(void) = NULL;

1071:   PetscFunctionBegin;
1073:   if (nvtxs) {
1074:     PetscAssertPointer(xadj, 3);
1075:     if (xadj[nvtxs]) PetscAssertPointer(adjncy, 4);
1076:   }
1077:   PetscTryMethod(pc, "PCBDDCSetLocalAdjacencyGraph_C", (PC, PetscInt, const PetscInt[], const PetscInt[], PetscCopyMode), (pc, nvtxs, xadj, adjncy, copymode));
1078:   /* free arrays if PCBDDC is not the PC type */
1079:   PetscCall(PetscObjectQueryFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", &f));
1080:   if (!f && copymode == PETSC_OWN_POINTER) {
1081:     PetscCall(PetscFree(xadj));
1082:     PetscCall(PetscFree(adjncy));
1083:   }
1084:   PetscFunctionReturn(PETSC_SUCCESS);
1085: }

1087: static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc, PetscInt n_is, IS ISForDofs[])
1088: {
1089:   PC_BDDC  *pcbddc = (PC_BDDC *)pc->data;
1090:   PetscInt  i;
1091:   PetscBool isequal = PETSC_FALSE;

1093:   PetscFunctionBegin;
1094:   if (pcbddc->n_ISForDofsLocal == n_is) {
1095:     for (i = 0; i < n_is; i++) {
1096:       PetscBool isequalt;
1097:       PetscCall(ISEqual(ISForDofs[i], pcbddc->ISForDofsLocal[i], &isequalt));
1098:       if (!isequalt) break;
1099:     }
1100:     if (i == n_is) isequal = PETSC_TRUE;
1101:   }
1102:   for (i = 0; i < n_is; i++) PetscCall(PetscObjectReference((PetscObject)ISForDofs[i]));
1103:   /* Destroy ISes if they were already set */
1104:   for (i = 0; i < pcbddc->n_ISForDofsLocal; i++) PetscCall(ISDestroy(&pcbddc->ISForDofsLocal[i]));
1105:   PetscCall(PetscFree(pcbddc->ISForDofsLocal));
1106:   /* last user setting takes precedence -> destroy any other customization */
1107:   for (i = 0; i < pcbddc->n_ISForDofs; i++) PetscCall(ISDestroy(&pcbddc->ISForDofs[i]));
1108:   PetscCall(PetscFree(pcbddc->ISForDofs));
1109:   pcbddc->n_ISForDofs = 0;
1110:   /* allocate space then set */
1111:   if (n_is) PetscCall(PetscMalloc1(n_is, &pcbddc->ISForDofsLocal));
1112:   for (i = 0; i < n_is; i++) pcbddc->ISForDofsLocal[i] = ISForDofs[i];
1113:   pcbddc->n_ISForDofsLocal = n_is;
1114:   if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
1115:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1116:   PetscFunctionReturn(PETSC_SUCCESS);
1117: }

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

1122:   Collective

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

1129:   Level: intermediate

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

1134: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDofsSplitting()`
1135: @*/
1136: PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc, PetscInt n_is, IS ISForDofs[])
1137: {
1138:   PetscInt i;

1140:   PetscFunctionBegin;
1143:   for (i = 0; i < n_is; i++) {
1144:     PetscCheckSameComm(pc, 1, ISForDofs[i], 3);
1146:   }
1147:   PetscTryMethod(pc, "PCBDDCSetDofsSplittingLocal_C", (PC, PetscInt, IS[]), (pc, n_is, ISForDofs));
1148:   PetscFunctionReturn(PETSC_SUCCESS);
1149: }

1151: static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc, PetscInt n_is, IS ISForDofs[])
1152: {
1153:   PC_BDDC  *pcbddc = (PC_BDDC *)pc->data;
1154:   PetscInt  i;
1155:   PetscBool isequal = PETSC_FALSE;

1157:   PetscFunctionBegin;
1158:   if (pcbddc->n_ISForDofs == n_is) {
1159:     for (i = 0; i < n_is; i++) {
1160:       PetscBool isequalt;
1161:       PetscCall(ISEqual(ISForDofs[i], pcbddc->ISForDofs[i], &isequalt));
1162:       if (!isequalt) break;
1163:     }
1164:     if (i == n_is) isequal = PETSC_TRUE;
1165:   }
1166:   for (i = 0; i < n_is; i++) PetscCall(PetscObjectReference((PetscObject)ISForDofs[i]));
1167:   /* Destroy ISes if they were already set */
1168:   for (i = 0; i < pcbddc->n_ISForDofs; i++) PetscCall(ISDestroy(&pcbddc->ISForDofs[i]));
1169:   PetscCall(PetscFree(pcbddc->ISForDofs));
1170:   /* last user setting takes precedence -> destroy any other customization */
1171:   for (i = 0; i < pcbddc->n_ISForDofsLocal; i++) PetscCall(ISDestroy(&pcbddc->ISForDofsLocal[i]));
1172:   PetscCall(PetscFree(pcbddc->ISForDofsLocal));
1173:   pcbddc->n_ISForDofsLocal = 0;
1174:   /* allocate space then set */
1175:   if (n_is) PetscCall(PetscMalloc1(n_is, &pcbddc->ISForDofs));
1176:   for (i = 0; i < n_is; i++) pcbddc->ISForDofs[i] = ISForDofs[i];
1177:   pcbddc->n_ISForDofs = n_is;
1178:   if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
1179:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1180:   PetscFunctionReturn(PETSC_SUCCESS);
1181: }

1183: /*@
1184:   PCBDDCSetDofsSplitting - Set the `IS` defining fields of the global matrix

1186:   Collective

1188:   Input Parameters:
1189: + pc        - the preconditioning context
1190: . n_is      - number of index sets defining the fields
1191: - ISForDofs - array of `IS` describing the fields in global ordering

1193:   Level: intermediate

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

1198: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCSetDofsSplittingLocal()`
1199: @*/
1200: PetscErrorCode PCBDDCSetDofsSplitting(PC pc, PetscInt n_is, IS ISForDofs[])
1201: {
1202:   PetscInt i;

1204:   PetscFunctionBegin;
1207:   for (i = 0; i < n_is; i++) {
1209:     PetscCheckSameComm(pc, 1, ISForDofs[i], 3);
1210:   }
1211:   PetscTryMethod(pc, "PCBDDCSetDofsSplitting_C", (PC, PetscInt, IS[]), (pc, n_is, ISForDofs));
1212:   PetscFunctionReturn(PETSC_SUCCESS);
1213: }

1215: /*
1216:    PCPreSolve_BDDC - Changes the right hand side and (if necessary) the initial
1217:                      guess if a transformation of basis approach has been selected.

1219:    Input Parameter:
1220: +  pc - the preconditioner context

1222:    Note:
1223:      The interface routine PCPreSolve() is not usually called directly by
1224:    the user, but instead is called by KSPSolve().
1225: */
1226: static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
1227: {
1228:   PC_BDDC  *pcbddc = (PC_BDDC *)pc->data;
1229:   PC_IS    *pcis   = (PC_IS *)(pc->data);
1230:   Vec       used_vec;
1231:   PetscBool iscg, save_rhs = PETSC_TRUE, benign_correction_computed;

1233:   PetscFunctionBegin;
1234:   /* if we are working with CG, one dirichlet solve can be avoided during Krylov iterations */
1235:   if (ksp) {
1236:     PetscCall(PetscObjectTypeCompareAny((PetscObject)ksp, &iscg, KSPCG, KSPGROPPCG, KSPPIPECG, KSPPIPELCG, KSPPIPECGRR, ""));
1237:     if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || !iscg || pc->mat != pc->pmat) PetscCall(PCBDDCSetUseExactDirichlet(pc, PETSC_FALSE));
1238:   }
1239:   if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || pc->mat != pc->pmat) PetscCall(PCBDDCSetUseExactDirichlet(pc, PETSC_FALSE));

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

1245:   pcbddc->temp_solution_used = PETSC_FALSE;
1246:   if (x) {
1247:     PetscCall(PetscObjectReference((PetscObject)x));
1248:     used_vec = x;
1249:   } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */
1250:     PetscCall(PetscObjectReference((PetscObject)pcbddc->temp_solution));
1251:     used_vec = pcbddc->temp_solution;
1252:     PetscCall(VecSet(used_vec, 0.0));
1253:     pcbddc->temp_solution_used = PETSC_TRUE;
1254:     PetscCall(VecCopy(rhs, pcbddc->original_rhs));
1255:     save_rhs                  = PETSC_FALSE;
1256:     pcbddc->eliminate_dirdofs = PETSC_TRUE;
1257:   }

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

1266:   pcbddc->rhs_change = PETSC_FALSE;
1267:   /* Take into account zeroed rows -> change rhs and store solution removed */
1268:   if (rhs && pcbddc->eliminate_dirdofs) {
1269:     IS dirIS = NULL;

1271:     /* DirichletBoundariesLocal may not be consistent among neighbours; gets a dirichlet dofs IS from graph (may be cached) */
1272:     PetscCall(PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph, &dirIS));
1273:     if (dirIS) {
1274:       Mat_IS            *matis = (Mat_IS *)pc->pmat->data;
1275:       PetscInt           dirsize, i, *is_indices;
1276:       PetscScalar       *array_x;
1277:       const PetscScalar *array_diagonal;

1279:       PetscCall(MatGetDiagonal(pc->pmat, pcis->vec1_global));
1280:       PetscCall(VecPointwiseDivide(pcis->vec1_global, rhs, pcis->vec1_global));
1281:       PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_global, pcis->vec2_N, INSERT_VALUES, SCATTER_FORWARD));
1282:       PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_global, pcis->vec2_N, INSERT_VALUES, SCATTER_FORWARD));
1283:       PetscCall(VecScatterBegin(matis->rctx, used_vec, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD));
1284:       PetscCall(VecScatterEnd(matis->rctx, used_vec, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD));
1285:       PetscCall(ISGetLocalSize(dirIS, &dirsize));
1286:       PetscCall(VecGetArray(pcis->vec1_N, &array_x));
1287:       PetscCall(VecGetArrayRead(pcis->vec2_N, &array_diagonal));
1288:       PetscCall(ISGetIndices(dirIS, (const PetscInt **)&is_indices));
1289:       for (i = 0; i < dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]];
1290:       PetscCall(ISRestoreIndices(dirIS, (const PetscInt **)&is_indices));
1291:       PetscCall(VecRestoreArrayRead(pcis->vec2_N, &array_diagonal));
1292:       PetscCall(VecRestoreArray(pcis->vec1_N, &array_x));
1293:       PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_N, used_vec, INSERT_VALUES, SCATTER_REVERSE));
1294:       PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_N, used_vec, INSERT_VALUES, SCATTER_REVERSE));
1295:       pcbddc->rhs_change = PETSC_TRUE;
1296:       PetscCall(ISDestroy(&dirIS));
1297:     }
1298:   }

1300:   /* remove the computed solution or the initial guess from the rhs */
1301:   if (pcbddc->rhs_change || (ksp && pcbddc->ksp_guess_nonzero)) {
1302:     /* save the original rhs */
1303:     if (save_rhs) {
1304:       PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1305:       save_rhs = PETSC_FALSE;
1306:     }
1307:     pcbddc->rhs_change = PETSC_TRUE;
1308:     PetscCall(VecScale(used_vec, -1.0));
1309:     PetscCall(MatMultAdd(pc->mat, used_vec, pcbddc->original_rhs, rhs));
1310:     PetscCall(VecScale(used_vec, -1.0));
1311:     PetscCall(VecCopy(used_vec, pcbddc->temp_solution));
1312:     pcbddc->temp_solution_used = PETSC_TRUE;
1313:     if (ksp) PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_FALSE));
1314:   }
1315:   PetscCall(VecDestroy(&used_vec));

1317:   /* compute initial vector in benign space if needed
1318:      and remove non-benign solution from the rhs */
1319:   benign_correction_computed = PETSC_FALSE;
1320:   if (rhs && pcbddc->benign_compute_correction && (pcbddc->benign_have_null || pcbddc->benign_apply_coarse_only)) {
1321:     /* compute u^*_h using ideas similar to those in Xuemin Tu's PhD thesis (see Section 4.8.1)
1322:        Recursively apply BDDC in the multilevel case */
1323:     if (!pcbddc->benign_vec) PetscCall(VecDuplicate(rhs, &pcbddc->benign_vec));
1324:     /* keep applying coarse solver unless we no longer have benign subdomains */
1325:     pcbddc->benign_apply_coarse_only = pcbddc->benign_have_null ? PETSC_TRUE : PETSC_FALSE;
1326:     if (!pcbddc->benign_skip_correction) {
1327:       PetscCall(PCApply_BDDC(pc, rhs, pcbddc->benign_vec));
1328:       benign_correction_computed = PETSC_TRUE;
1329:       if (pcbddc->temp_solution_used) PetscCall(VecAXPY(pcbddc->temp_solution, 1.0, pcbddc->benign_vec));
1330:       PetscCall(VecScale(pcbddc->benign_vec, -1.0));
1331:       /* store the original rhs if not done earlier */
1332:       if (save_rhs) PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1333:       if (pcbddc->rhs_change) {
1334:         PetscCall(MatMultAdd(pc->mat, pcbddc->benign_vec, rhs, rhs));
1335:       } else {
1336:         PetscCall(MatMultAdd(pc->mat, pcbddc->benign_vec, pcbddc->original_rhs, rhs));
1337:       }
1338:       pcbddc->rhs_change = PETSC_TRUE;
1339:     }
1340:     pcbddc->benign_apply_coarse_only = PETSC_FALSE;
1341:   } else {
1342:     PetscCall(VecDestroy(&pcbddc->benign_vec));
1343:   }

1345:   /* dbg output */
1346:   if (pcbddc->dbg_flag && benign_correction_computed) {
1347:     Vec v;

1349:     PetscCall(VecDuplicate(pcis->vec1_global, &v));
1350:     if (pcbddc->ChangeOfBasisMatrix) {
1351:       PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, rhs, v));
1352:     } else {
1353:       PetscCall(VecCopy(rhs, v));
1354:     }
1355:     PetscCall(PCBDDCBenignGetOrSetP0(pc, v, PETSC_TRUE));
1356:     PetscCall(PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "LEVEL %" PetscInt_FMT ": is the correction benign?\n", pcbddc->current_level));
1357:     PetscCall(PetscScalarView(pcbddc->benign_n, pcbddc->benign_p0, pcbddc->dbg_viewer));
1358:     PetscCall(PetscViewerFlush(pcbddc->dbg_viewer));
1359:     PetscCall(VecDestroy(&v));
1360:   }

1362:   /* set initial guess if using PCG */
1363:   pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1364:   if (x && pcbddc->use_exact_dirichlet_trick) {
1365:     PetscCall(VecSet(x, 0.0));
1366:     if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) {
1367:       if (benign_correction_computed) { /* we have already saved the changed rhs */
1368:         PetscCall(VecLockReadPop(pcis->vec1_global));
1369:       } else {
1370:         PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, rhs, pcis->vec1_global));
1371:       }
1372:       PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_global, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1373:       PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_global, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1374:     } else {
1375:       PetscCall(VecScatterBegin(pcis->global_to_D, rhs, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1376:       PetscCall(VecScatterEnd(pcis->global_to_D, rhs, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1377:     }
1378:     PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1379:     PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
1380:     PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1381:     PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
1382:     if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) {
1383:       PetscCall(VecSet(pcis->vec1_global, 0.));
1384:       PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, pcis->vec1_global, INSERT_VALUES, SCATTER_REVERSE));
1385:       PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, pcis->vec1_global, INSERT_VALUES, SCATTER_REVERSE));
1386:       PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcis->vec1_global, x));
1387:     } else {
1388:       PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, x, INSERT_VALUES, SCATTER_REVERSE));
1389:       PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, x, INSERT_VALUES, SCATTER_REVERSE));
1390:     }
1391:     if (ksp) PetscCall(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
1392:     pcbddc->exact_dirichlet_trick_app = PETSC_TRUE;
1393:   } else if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior && benign_correction_computed && pcbddc->use_exact_dirichlet_trick) {
1394:     PetscCall(VecLockReadPop(pcis->vec1_global));
1395:   }
1396:   PetscFunctionReturn(PETSC_SUCCESS);
1397: }

1399: /*
1400:    PCPostSolve_BDDC - Changes the computed solution if a transformation of basis
1401:                      approach has been selected. Also, restores rhs to its original state.

1403:    Input Parameter:
1404: +  pc - the preconditioner context

1406:    Application Interface Routine: PCPostSolve()

1408:    Note:
1409:      The interface routine PCPostSolve() is not usually called directly by
1410:      the user, but instead is called by KSPSolve().
1411: */
1412: static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
1413: {
1414:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

1416:   PetscFunctionBegin;
1417:   /* add solution removed in presolve */
1418:   if (x && pcbddc->rhs_change) {
1419:     if (pcbddc->temp_solution_used) {
1420:       PetscCall(VecAXPY(x, 1.0, pcbddc->temp_solution));
1421:     } else if (pcbddc->benign_compute_correction && pcbddc->benign_vec) {
1422:       PetscCall(VecAXPY(x, -1.0, pcbddc->benign_vec));
1423:     }
1424:     /* restore to original state (not for FETI-DP) */
1425:     if (ksp) pcbddc->temp_solution_used = PETSC_FALSE;
1426:   }

1428:   /* restore rhs to its original state (not needed for FETI-DP) */
1429:   if (rhs && pcbddc->rhs_change) {
1430:     PetscCall(VecSwap(rhs, pcbddc->original_rhs));
1431:     pcbddc->rhs_change = PETSC_FALSE;
1432:   }
1433:   /* restore ksp guess state */
1434:   if (ksp) {
1435:     PetscCall(KSPSetInitialGuessNonzero(ksp, pcbddc->ksp_guess_nonzero));
1436:     /* reset flag for exact dirichlet trick */
1437:     pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1438:   }
1439:   PetscFunctionReturn(PETSC_SUCCESS);
1440: }

1442: // PetscClangLinter pragma disable: -fdoc-sowing-chars
1443: /*
1444:   PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner
1445:   by setting data structures and options.

1447:   Input Parameter:
1448: . pc - the preconditioner context

1450:    Application Interface Routine: PCSetUp()

1452: */
1453: static PetscErrorCode PCSetUp_BDDC(PC pc)
1454: {
1455:   PC_BDDC        *pcbddc = (PC_BDDC *)pc->data;
1456:   PCBDDCSubSchurs sub_schurs;
1457:   Mat_IS         *matis;
1458:   MatNullSpace    nearnullspace;
1459:   Mat             lA;
1460:   IS              lP, zerodiag = NULL;
1461:   PetscInt        nrows, ncols;
1462:   PetscMPIInt     size;
1463:   PetscBool       computesubschurs;
1464:   PetscBool       computeconstraintsmatrix;
1465:   PetscBool       new_nearnullspace_provided, ismatis, rl;
1466:   PetscBool       isset, issym, isspd;

1468:   PetscFunctionBegin;
1469:   PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATIS, &ismatis));
1470:   PetscCheck(ismatis, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PCBDDC preconditioner requires matrix of type MATIS");
1471:   PetscCall(MatGetSize(pc->pmat, &nrows, &ncols));
1472:   PetscCheck(nrows == ncols, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCBDDC preconditioner requires a square preconditioning matrix");
1473:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));

1475:   matis = (Mat_IS *)pc->pmat->data;
1476:   /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */
1477:   /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup
1478:      Also, BDDC builds its own KSP for the Dirichlet problem */
1479:   rl = pcbddc->recompute_topography;
1480:   if (!pc->setupcalled || pc->flag == DIFFERENT_NONZERO_PATTERN) rl = PETSC_TRUE;
1481:   PetscCall(MPIU_Allreduce(&rl, &pcbddc->recompute_topography, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject)pc)));
1482:   if (pcbddc->recompute_topography) {
1483:     pcbddc->graphanalyzed    = PETSC_FALSE;
1484:     computeconstraintsmatrix = PETSC_TRUE;
1485:   } else {
1486:     computeconstraintsmatrix = PETSC_FALSE;
1487:   }

1489:   /* check parameters' compatibility */
1490:   if (!pcbddc->use_deluxe_scaling) pcbddc->deluxe_zerorows = PETSC_FALSE;
1491:   pcbddc->adaptive_selection   = (PetscBool)(pcbddc->adaptive_threshold[0] != 0.0 || pcbddc->adaptive_threshold[1] != 0.0);
1492:   pcbddc->use_deluxe_scaling   = (PetscBool)(pcbddc->use_deluxe_scaling && size > 1);
1493:   pcbddc->adaptive_selection   = (PetscBool)(pcbddc->adaptive_selection && size > 1);
1494:   pcbddc->adaptive_userdefined = (PetscBool)(pcbddc->adaptive_selection && pcbddc->adaptive_userdefined);
1495:   if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE;

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

1499:   /* activate all connected components if the netflux has been requested */
1500:   if (pcbddc->compute_nonetflux) {
1501:     pcbddc->use_vertices = PETSC_TRUE;
1502:     pcbddc->use_edges    = PETSC_TRUE;
1503:     pcbddc->use_faces    = PETSC_TRUE;
1504:   }

1506:   /* Get stdout for dbg */
1507:   if (pcbddc->dbg_flag) {
1508:     if (!pcbddc->dbg_viewer) pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc));
1509:     PetscCall(PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer));
1510:     PetscCall(PetscViewerASCIIAddTab(pcbddc->dbg_viewer, 2 * pcbddc->current_level));
1511:   }

1513:   /* process topology information */
1514:   PetscCall(PetscLogEventBegin(PC_BDDC_Topology[pcbddc->current_level], pc, 0, 0, 0));
1515:   if (pcbddc->recompute_topography) {
1516:     PetscCall(PCBDDCComputeLocalTopologyInfo(pc));
1517:     if (pcbddc->discretegradient) PetscCall(PCBDDCNedelecSupport(pc));
1518:   }
1519:   if (pcbddc->corner_selected) pcbddc->use_vertices = PETSC_TRUE;

1521:   /* change basis if requested by the user */
1522:   if (pcbddc->user_ChangeOfBasisMatrix) {
1523:     /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */
1524:     pcbddc->use_change_of_basis = PETSC_FALSE;
1525:     PetscCall(PCBDDCComputeLocalMatrix(pc, pcbddc->user_ChangeOfBasisMatrix));
1526:   } else {
1527:     PetscCall(MatDestroy(&pcbddc->local_mat));
1528:     PetscCall(PetscObjectReference((PetscObject)matis->A));
1529:     pcbddc->local_mat = matis->A;
1530:   }

1532:   /*
1533:      Compute change of basis on local pressures (aka zerodiag dofs) with the benign trick
1534:      This should come earlier then PCISSetUp for extracting the correct subdomain matrices
1535:   */
1536:   PetscCall(PCBDDCBenignShellMat(pc, PETSC_TRUE));
1537:   if (pcbddc->benign_saddle_point) {
1538:     PC_IS *pcis = (PC_IS *)pc->data;

1540:     if (pcbddc->user_ChangeOfBasisMatrix || pcbddc->use_change_of_basis || !computesubschurs) pcbddc->benign_change_explicit = PETSC_TRUE;
1541:     /* detect local saddle point and change the basis in pcbddc->local_mat */
1542:     PetscCall(PCBDDCBenignDetectSaddlePoint(pc, (PetscBool)(!pcbddc->recompute_topography), &zerodiag));
1543:     /* pop B0 mat from local mat */
1544:     PetscCall(PCBDDCBenignPopOrPushB0(pc, PETSC_TRUE));
1545:     /* give pcis a hint to not reuse submatrices during PCISCreate */
1546:     if (pc->flag == SAME_NONZERO_PATTERN && pcis->reusesubmatrices == PETSC_TRUE) {
1547:       if (pcbddc->benign_n && (pcbddc->benign_change_explicit || pcbddc->dbg_flag)) {
1548:         pcis->reusesubmatrices = PETSC_FALSE;
1549:       } else {
1550:         pcis->reusesubmatrices = PETSC_TRUE;
1551:       }
1552:     } else {
1553:       pcis->reusesubmatrices = PETSC_FALSE;
1554:     }
1555:   }

1557:   /* propagate relevant information */
1558:   PetscCall(MatIsSymmetricKnown(matis->A, &isset, &issym));
1559:   if (isset) PetscCall(MatSetOption(pcbddc->local_mat, MAT_SYMMETRIC, issym));
1560:   PetscCall(MatIsSPDKnown(matis->A, &isset, &isspd));
1561:   if (isset) PetscCall(MatSetOption(pcbddc->local_mat, MAT_SPD, isspd));

1563:   /* Set up all the "iterative substructuring" common block without computing solvers */
1564:   {
1565:     Mat temp_mat;

1567:     temp_mat = matis->A;
1568:     matis->A = pcbddc->local_mat;
1569:     PetscCall(PCISSetUp(pc, PETSC_TRUE, PETSC_FALSE));
1570:     pcbddc->local_mat = matis->A;
1571:     matis->A          = temp_mat;
1572:   }

1574:   /* Analyze interface */
1575:   if (!pcbddc->graphanalyzed) {
1576:     PetscCall(PCBDDCAnalyzeInterface(pc));
1577:     computeconstraintsmatrix = PETSC_TRUE;
1578:     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");
1579:     if (pcbddc->compute_nonetflux) {
1580:       MatNullSpace nnfnnsp;

1582:       PetscCheck(pcbddc->divudotp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Missing divudotp operator");
1583:       PetscCall(PCBDDCComputeNoNetFlux(pc->pmat, pcbddc->divudotp, pcbddc->divudotp_trans, pcbddc->divudotp_vl2l, pcbddc->mat_graph, &nnfnnsp));
1584:       /* TODO what if a nearnullspace is already attached? */
1585:       if (nnfnnsp) {
1586:         PetscCall(MatSetNearNullSpace(pc->pmat, nnfnnsp));
1587:         PetscCall(MatNullSpaceDestroy(&nnfnnsp));
1588:       }
1589:     }
1590:   }
1591:   PetscCall(PetscLogEventEnd(PC_BDDC_Topology[pcbddc->current_level], pc, 0, 0, 0));

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

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

1604:   /* finish setup solvers and do adaptive selection of constraints */
1605:   sub_schurs = pcbddc->sub_schurs;
1606:   if (sub_schurs && sub_schurs->schur_explicit) {
1607:     if (computesubschurs) PetscCall(PCBDDCSetUpSubSchurs(pc));
1608:     PetscCall(PCBDDCSetUpLocalSolvers(pc, PETSC_TRUE, PETSC_FALSE));
1609:   } else {
1610:     PetscCall(PCBDDCSetUpLocalSolvers(pc, PETSC_TRUE, PETSC_FALSE));
1611:     if (computesubschurs) PetscCall(PCBDDCSetUpSubSchurs(pc));
1612:   }
1613:   if (pcbddc->adaptive_selection) {
1614:     PetscCall(PCBDDCAdaptiveSelection(pc));
1615:     computeconstraintsmatrix = PETSC_TRUE;
1616:   }

1618:   /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */
1619:   new_nearnullspace_provided = PETSC_FALSE;
1620:   PetscCall(MatGetNearNullSpace(pc->pmat, &nearnullspace));
1621:   if (pcbddc->onearnullspace) { /* already used nearnullspace */
1622:     if (!nearnullspace) {       /* near null space attached to mat has been destroyed */
1623:       new_nearnullspace_provided = PETSC_TRUE;
1624:     } else {
1625:       /* determine if the two nullspaces are different (should be lightweight) */
1626:       if (nearnullspace != pcbddc->onearnullspace) {
1627:         new_nearnullspace_provided = PETSC_TRUE;
1628:       } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */
1629:         PetscInt         i;
1630:         const Vec       *nearnullvecs;
1631:         PetscObjectState state;
1632:         PetscInt         nnsp_size;
1633:         PetscCall(MatNullSpaceGetVecs(nearnullspace, NULL, &nnsp_size, &nearnullvecs));
1634:         for (i = 0; i < nnsp_size; i++) {
1635:           PetscCall(PetscObjectStateGet((PetscObject)nearnullvecs[i], &state));
1636:           if (pcbddc->onearnullvecs_state[i] != state) {
1637:             new_nearnullspace_provided = PETSC_TRUE;
1638:             break;
1639:           }
1640:         }
1641:       }
1642:     }
1643:   } else {
1644:     if (!nearnullspace) { /* both nearnullspaces are null */
1645:       new_nearnullspace_provided = PETSC_FALSE;
1646:     } else { /* nearnullspace attached later */
1647:       new_nearnullspace_provided = PETSC_TRUE;
1648:     }
1649:   }

1651:   /* Setup constraints and related work vectors */
1652:   /* reset primal space flags */
1653:   PetscCall(PetscLogEventBegin(PC_BDDC_LocalWork[pcbddc->current_level], pc, 0, 0, 0));
1654:   pcbddc->new_primal_space       = PETSC_FALSE;
1655:   pcbddc->new_primal_space_local = PETSC_FALSE;
1656:   if (computeconstraintsmatrix || new_nearnullspace_provided) {
1657:     /* It also sets the primal space flags */
1658:     PetscCall(PCBDDCConstraintsSetUp(pc));
1659:   }
1660:   /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */
1661:   PetscCall(PCBDDCSetUpLocalWorkVectors(pc));

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

1666:     PetscCall(PCBDDCComputeLocalMatrix(pc, pcbddc->ChangeOfBasisMatrix));
1667:     if (pcbddc->benign_change) {
1668:       PetscCall(MatDestroy(&pcbddc->benign_B0));
1669:       /* pop B0 from pcbddc->local_mat */
1670:       PetscCall(PCBDDCBenignPopOrPushB0(pc, PETSC_TRUE));
1671:     }
1672:     /* get submatrices */
1673:     PetscCall(MatDestroy(&pcis->A_IB));
1674:     PetscCall(MatDestroy(&pcis->A_BI));
1675:     PetscCall(MatDestroy(&pcis->A_BB));
1676:     PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_B_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &pcis->A_BB));
1677:     PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_I_local, pcis->is_B_local, MAT_INITIAL_MATRIX, &pcis->A_IB));
1678:     PetscCall(MatCreateSubMatrix(pcbddc->local_mat, pcis->is_B_local, pcis->is_I_local, MAT_INITIAL_MATRIX, &pcis->A_BI));
1679:     /* set flag in pcis to not reuse submatrices during PCISCreate */
1680:     pcis->reusesubmatrices = PETSC_FALSE;
1681:   } else if (!pcbddc->user_ChangeOfBasisMatrix && !pcbddc->benign_change) {
1682:     PetscCall(MatDestroy(&pcbddc->local_mat));
1683:     PetscCall(PetscObjectReference((PetscObject)matis->A));
1684:     pcbddc->local_mat = matis->A;
1685:   }

1687:   /* interface pressure block row for B_C */
1688:   PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_lP", (PetscObject *)&lP));
1689:   PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_lA", (PetscObject *)&lA));
1690:   if (lA && lP) {
1691:     PC_IS    *pcis = (PC_IS *)pc->data;
1692:     Mat       B_BI, B_BB, Bt_BI, Bt_BB;
1693:     PetscBool issym;

1695:     PetscCall(MatIsSymmetric(lA, PETSC_SMALL, &issym));
1696:     if (issym) {
1697:       PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_I_local, MAT_INITIAL_MATRIX, &B_BI));
1698:       PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_B_local, MAT_INITIAL_MATRIX, &B_BB));
1699:       PetscCall(MatCreateTranspose(B_BI, &Bt_BI));
1700:       PetscCall(MatCreateTranspose(B_BB, &Bt_BB));
1701:     } else {
1702:       PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_I_local, MAT_INITIAL_MATRIX, &B_BI));
1703:       PetscCall(MatCreateSubMatrix(lA, lP, pcis->is_B_local, MAT_INITIAL_MATRIX, &B_BB));
1704:       PetscCall(MatCreateSubMatrix(lA, pcis->is_I_local, lP, MAT_INITIAL_MATRIX, &Bt_BI));
1705:       PetscCall(MatCreateSubMatrix(lA, pcis->is_B_local, lP, MAT_INITIAL_MATRIX, &Bt_BB));
1706:     }
1707:     PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_B_BI", (PetscObject)B_BI));
1708:     PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_B_BB", (PetscObject)B_BB));
1709:     PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_Bt_BI", (PetscObject)Bt_BI));
1710:     PetscCall(PetscObjectCompose((PetscObject)pc, "__KSPFETIDP_Bt_BB", (PetscObject)Bt_BB));
1711:     PetscCall(MatDestroy(&B_BI));
1712:     PetscCall(MatDestroy(&B_BB));
1713:     PetscCall(MatDestroy(&Bt_BI));
1714:     PetscCall(MatDestroy(&Bt_BB));
1715:   }
1716:   PetscCall(PetscLogEventEnd(PC_BDDC_LocalWork[pcbddc->current_level], pc, 0, 0, 0));

1718:   /* SetUp coarse and local Neumann solvers */
1719:   PetscCall(PCBDDCSetUpSolvers(pc));
1720:   /* SetUp Scaling operator */
1721:   if (pcbddc->use_deluxe_scaling) PetscCall(PCBDDCScalingSetUp(pc));

1723:   /* mark topography as done */
1724:   pcbddc->recompute_topography = PETSC_FALSE;

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

1729:   if (pcbddc->dbg_flag) {
1730:     PetscCall(PetscViewerASCIISubtractTab(pcbddc->dbg_viewer, 2 * pcbddc->current_level));
1731:     PetscCall(PetscViewerASCIIPopSynchronized(pcbddc->dbg_viewer));
1732:   }
1733:   PetscFunctionReturn(PETSC_SUCCESS);
1734: }

1736: // PetscClangLinter pragma disable: -fdoc-sowing-chars
1737: /*
1738:    PCApply_BDDC - Applies the BDDC operator to a vector.

1740:    Input Parameters:
1741: +  pc - the preconditioner context
1742: -  r - input vector (global)

1744:    Output Parameter:
1745: .  z - output vector (global)

1747:    Application Interface Routine: PCApply()
1748:  */
1749: static PetscErrorCode PCApply_BDDC(PC pc, Vec r, Vec z)
1750: {
1751:   PC_IS            *pcis   = (PC_IS *)(pc->data);
1752:   PC_BDDC          *pcbddc = (PC_BDDC *)(pc->data);
1753:   Mat               lA     = NULL;
1754:   PetscInt          n_B = pcis->n_B, n_D = pcis->n - n_B;
1755:   const PetscScalar one   = 1.0;
1756:   const PetscScalar m_one = -1.0;
1757:   const PetscScalar zero  = 0.0;
1758:   /* This code is similar to that provided in nn.c for PCNN
1759:    NN interface preconditioner changed to BDDC
1760:    Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */

1762:   PetscFunctionBegin;
1763:   PetscCall(PetscCitationsRegister(citation, &cited));
1764:   if (pcbddc->switch_static) PetscCall(MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat, &lA));

1766:   if (pcbddc->ChangeOfBasisMatrix) {
1767:     Vec swap;

1769:     PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, r, pcbddc->work_change));
1770:     swap                = pcbddc->work_change;
1771:     pcbddc->work_change = r;
1772:     r                   = swap;
1773:     /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */
1774:     if (pcbddc->benign_apply_coarse_only && pcbddc->use_exact_dirichlet_trick && pcbddc->change_interior) {
1775:       PetscCall(VecCopy(r, pcis->vec1_global));
1776:       PetscCall(VecLockReadPush(pcis->vec1_global));
1777:     }
1778:   }
1779:   if (pcbddc->benign_have_null) { /* get p0 from r */
1780:     PetscCall(PCBDDCBenignGetOrSetP0(pc, r, PETSC_TRUE));
1781:   }
1782:   if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_DIRICHLET && !pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
1783:     PetscCall(VecCopy(r, z));
1784:     /* First Dirichlet solve */
1785:     PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1786:     PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1787:     /*
1788:       Assembling right hand side for BDDC operator
1789:       - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
1790:       - pcis->vec1_B the interface part of the global vector z
1791:     */
1792:     if (n_D) {
1793:       PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1794:       PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
1795:       PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1796:       PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
1797:       PetscCall(VecScale(pcis->vec2_D, m_one));
1798:       if (pcbddc->switch_static) {
1799:         PetscCall(VecSet(pcis->vec1_N, 0.));
1800:         PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1801:         PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1802:         if (!pcbddc->switch_static_change) {
1803:           PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
1804:         } else {
1805:           PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1806:           PetscCall(MatMult(lA, pcis->vec2_N, pcis->vec1_N));
1807:           PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1808:         }
1809:         PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
1810:         PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
1811:         PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1812:         PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1813:       } else {
1814:         PetscCall(MatMult(pcis->A_BI, pcis->vec2_D, pcis->vec1_B));
1815:       }
1816:     } else {
1817:       PetscCall(VecSet(pcis->vec1_B, zero));
1818:     }
1819:     PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
1820:     PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
1821:     PetscCall(PCBDDCScalingRestriction(pc, z, pcis->vec1_B));
1822:   } else {
1823:     if (!pcbddc->benign_apply_coarse_only) PetscCall(PCBDDCScalingRestriction(pc, r, pcis->vec1_B));
1824:   }
1825:   if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_LUMP) {
1826:     PetscCheck(pcbddc->switch_static, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "You forgot to pass -pc_bddc_switch_static");
1827:     PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1828:     PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1829:   }

1831:   /* Apply interface preconditioner
1832:      input/output vecs: pcis->vec1_B and pcis->vec1_D */
1833:   PetscCall(PCBDDCApplyInterfacePreconditioner(pc, PETSC_FALSE));

1835:   /* Apply transpose of partition of unity operator */
1836:   PetscCall(PCBDDCScalingExtension(pc, pcis->vec1_B, z));
1837:   if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_LUMP) {
1838:     PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_D, z, INSERT_VALUES, SCATTER_REVERSE));
1839:     PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_D, z, INSERT_VALUES, SCATTER_REVERSE));
1840:     PetscFunctionReturn(PETSC_SUCCESS);
1841:   }
1842:   /* Second Dirichlet solve and assembling of output */
1843:   PetscCall(VecScatterBegin(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1844:   PetscCall(VecScatterEnd(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1845:   if (n_B) {
1846:     if (pcbddc->switch_static) {
1847:       PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1848:       PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1849:       PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1850:       PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1851:       if (!pcbddc->switch_static_change) {
1852:         PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
1853:       } else {
1854:         PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1855:         PetscCall(MatMult(lA, pcis->vec2_N, pcis->vec1_N));
1856:         PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1857:       }
1858:       PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
1859:       PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
1860:     } else {
1861:       PetscCall(MatMult(pcis->A_IB, pcis->vec1_B, pcis->vec3_D));
1862:     }
1863:   } else if (pcbddc->switch_static) { /* n_B is zero */
1864:     if (!pcbddc->switch_static_change) {
1865:       PetscCall(MatMult(lA, pcis->vec1_D, pcis->vec3_D));
1866:     } else {
1867:       PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_D, pcis->vec1_N));
1868:       PetscCall(MatMult(lA, pcis->vec1_N, pcis->vec2_N));
1869:       PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec2_N, pcis->vec3_D));
1870:     }
1871:   }
1872:   PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1873:   PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec3_D, pcis->vec4_D));
1874:   PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1875:   PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec4_D));

1877:   if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
1878:     if (pcbddc->switch_static) {
1879:       PetscCall(VecAXPBYPCZ(pcis->vec2_D, m_one, one, m_one, pcis->vec4_D, pcis->vec1_D));
1880:     } else {
1881:       PetscCall(VecAXPBY(pcis->vec2_D, m_one, m_one, pcis->vec4_D));
1882:     }
1883:     PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
1884:     PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
1885:   } else {
1886:     if (pcbddc->switch_static) {
1887:       PetscCall(VecAXPBY(pcis->vec4_D, one, m_one, pcis->vec1_D));
1888:     } else {
1889:       PetscCall(VecScale(pcis->vec4_D, m_one));
1890:     }
1891:     PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
1892:     PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
1893:   }
1894:   if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */
1895:     if (pcbddc->benign_apply_coarse_only) PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
1896:     PetscCall(PCBDDCBenignGetOrSetP0(pc, z, PETSC_FALSE));
1897:   }

1899:   if (pcbddc->ChangeOfBasisMatrix) {
1900:     pcbddc->work_change = r;
1901:     PetscCall(VecCopy(z, pcbddc->work_change));
1902:     PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcbddc->work_change, z));
1903:   }
1904:   PetscFunctionReturn(PETSC_SUCCESS);
1905: }

1907: /*
1908:    PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector.

1910:    Input Parameters:
1911: +  pc - the preconditioner context
1912: -  r - input vector (global)

1914:    Output Parameter:
1915: .  z - output vector (global)

1917:    Application Interface Routine: PCApplyTranspose()
1918:  */
1919: static PetscErrorCode PCApplyTranspose_BDDC(PC pc, Vec r, Vec z)
1920: {
1921:   PC_IS            *pcis   = (PC_IS *)(pc->data);
1922:   PC_BDDC          *pcbddc = (PC_BDDC *)(pc->data);
1923:   Mat               lA     = NULL;
1924:   PetscInt          n_B = pcis->n_B, n_D = pcis->n - n_B;
1925:   const PetscScalar one   = 1.0;
1926:   const PetscScalar m_one = -1.0;
1927:   const PetscScalar zero  = 0.0;

1929:   PetscFunctionBegin;
1930:   PetscCall(PetscCitationsRegister(citation, &cited));
1931:   if (pcbddc->switch_static) PetscCall(MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat, &lA));
1932:   if (pcbddc->ChangeOfBasisMatrix) {
1933:     Vec swap;

1935:     PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, r, pcbddc->work_change));
1936:     swap                = pcbddc->work_change;
1937:     pcbddc->work_change = r;
1938:     r                   = swap;
1939:     /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */
1940:     if (pcbddc->benign_apply_coarse_only && pcbddc->exact_dirichlet_trick_app && pcbddc->change_interior) {
1941:       PetscCall(VecCopy(r, pcis->vec1_global));
1942:       PetscCall(VecLockReadPush(pcis->vec1_global));
1943:     }
1944:   }
1945:   if (pcbddc->benign_have_null) { /* get p0 from r */
1946:     PetscCall(PCBDDCBenignGetOrSetP0(pc, r, PETSC_TRUE));
1947:   }
1948:   if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
1949:     PetscCall(VecCopy(r, z));
1950:     /* First Dirichlet solve */
1951:     PetscCall(VecScatterBegin(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1952:     PetscCall(VecScatterEnd(pcis->global_to_D, r, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD));
1953:     /*
1954:       Assembling right hand side for BDDC operator
1955:       - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
1956:       - pcis->vec1_B the interface part of the global vector z
1957:     */
1958:     if (n_D) {
1959:       PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1960:       PetscCall(KSPSolveTranspose(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D));
1961:       PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
1962:       PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D));
1963:       PetscCall(VecScale(pcis->vec2_D, m_one));
1964:       if (pcbddc->switch_static) {
1965:         PetscCall(VecSet(pcis->vec1_N, 0.));
1966:         PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1967:         PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
1968:         if (!pcbddc->switch_static_change) {
1969:           PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
1970:         } else {
1971:           PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1972:           PetscCall(MatMultTranspose(lA, pcis->vec2_N, pcis->vec1_N));
1973:           PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
1974:         }
1975:         PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
1976:         PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec1_D, ADD_VALUES, SCATTER_FORWARD));
1977:         PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1978:         PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec2_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
1979:       } else {
1980:         PetscCall(MatMultTranspose(pcis->A_IB, pcis->vec2_D, pcis->vec1_B));
1981:       }
1982:     } else {
1983:       PetscCall(VecSet(pcis->vec1_B, zero));
1984:     }
1985:     PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
1986:     PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, z, ADD_VALUES, SCATTER_REVERSE));
1987:     PetscCall(PCBDDCScalingRestriction(pc, z, pcis->vec1_B));
1988:   } else {
1989:     PetscCall(PCBDDCScalingRestriction(pc, r, pcis->vec1_B));
1990:   }

1992:   /* Apply interface preconditioner
1993:      input/output vecs: pcis->vec1_B and pcis->vec1_D */
1994:   PetscCall(PCBDDCApplyInterfacePreconditioner(pc, PETSC_TRUE));

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

1999:   /* Second Dirichlet solve and assembling of output */
2000:   PetscCall(VecScatterBegin(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2001:   PetscCall(VecScatterEnd(pcis->global_to_B, z, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
2002:   if (n_B) {
2003:     if (pcbddc->switch_static) {
2004:       PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2005:       PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec1_D, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2006:       PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2007:       PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_B, pcis->vec1_N, INSERT_VALUES, SCATTER_REVERSE));
2008:       if (!pcbddc->switch_static_change) {
2009:         PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
2010:       } else {
2011:         PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2012:         PetscCall(MatMultTranspose(lA, pcis->vec2_N, pcis->vec1_N));
2013:         PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec1_N, pcis->vec2_N));
2014:       }
2015:       PetscCall(VecScatterBegin(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
2016:       PetscCall(VecScatterEnd(pcis->N_to_D, pcis->vec2_N, pcis->vec3_D, INSERT_VALUES, SCATTER_FORWARD));
2017:     } else {
2018:       PetscCall(MatMultTranspose(pcis->A_BI, pcis->vec1_B, pcis->vec3_D));
2019:     }
2020:   } else if (pcbddc->switch_static) { /* n_B is zero */
2021:     if (!pcbddc->switch_static_change) {
2022:       PetscCall(MatMultTranspose(lA, pcis->vec1_D, pcis->vec3_D));
2023:     } else {
2024:       PetscCall(MatMult(pcbddc->switch_static_change, pcis->vec1_D, pcis->vec1_N));
2025:       PetscCall(MatMultTranspose(lA, pcis->vec1_N, pcis->vec2_N));
2026:       PetscCall(MatMultTranspose(pcbddc->switch_static_change, pcis->vec2_N, pcis->vec3_D));
2027:     }
2028:   }
2029:   PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2030:   PetscCall(KSPSolveTranspose(pcbddc->ksp_D, pcis->vec3_D, pcis->vec4_D));
2031:   PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0));
2032:   PetscCall(KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec4_D));
2033:   if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
2034:     if (pcbddc->switch_static) {
2035:       PetscCall(VecAXPBYPCZ(pcis->vec2_D, m_one, one, m_one, pcis->vec4_D, pcis->vec1_D));
2036:     } else {
2037:       PetscCall(VecAXPBY(pcis->vec2_D, m_one, m_one, pcis->vec4_D));
2038:     }
2039:     PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
2040:     PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec2_D, z, INSERT_VALUES, SCATTER_REVERSE));
2041:   } else {
2042:     if (pcbddc->switch_static) {
2043:       PetscCall(VecAXPBY(pcis->vec4_D, one, m_one, pcis->vec1_D));
2044:     } else {
2045:       PetscCall(VecScale(pcis->vec4_D, m_one));
2046:     }
2047:     PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
2048:     PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec4_D, z, INSERT_VALUES, SCATTER_REVERSE));
2049:   }
2050:   if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */
2051:     PetscCall(PCBDDCBenignGetOrSetP0(pc, z, PETSC_FALSE));
2052:   }
2053:   if (pcbddc->ChangeOfBasisMatrix) {
2054:     pcbddc->work_change = r;
2055:     PetscCall(VecCopy(z, pcbddc->work_change));
2056:     PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, pcbddc->work_change, z));
2057:   }
2058:   PetscFunctionReturn(PETSC_SUCCESS);
2059: }

2061: static PetscErrorCode PCReset_BDDC(PC pc)
2062: {
2063:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
2064:   PC_IS   *pcis   = (PC_IS *)pc->data;
2065:   KSP      kspD, kspR, kspC;

2067:   PetscFunctionBegin;
2068:   /* free BDDC custom data  */
2069:   PetscCall(PCBDDCResetCustomization(pc));
2070:   /* destroy objects related to topography */
2071:   PetscCall(PCBDDCResetTopography(pc));
2072:   /* destroy objects for scaling operator */
2073:   PetscCall(PCBDDCScalingDestroy(pc));
2074:   /* free solvers stuff */
2075:   PetscCall(PCBDDCResetSolvers(pc));
2076:   /* free global vectors needed in presolve */
2077:   PetscCall(VecDestroy(&pcbddc->temp_solution));
2078:   PetscCall(VecDestroy(&pcbddc->original_rhs));
2079:   /* free data created by PCIS */
2080:   PetscCall(PCISReset(pc));

2082:   /* restore defaults */
2083:   kspD = pcbddc->ksp_D;
2084:   kspR = pcbddc->ksp_R;
2085:   kspC = pcbddc->coarse_ksp;
2086:   PetscCall(PetscMemzero(pc->data, sizeof(*pcbddc)));
2087:   pcis->n_neigh                     = -1;
2088:   pcis->scaling_factor              = 1.0;
2089:   pcis->reusesubmatrices            = PETSC_TRUE;
2090:   pcbddc->use_local_adj             = PETSC_TRUE;
2091:   pcbddc->use_vertices              = PETSC_TRUE;
2092:   pcbddc->use_edges                 = PETSC_TRUE;
2093:   pcbddc->symmetric_primal          = PETSC_TRUE;
2094:   pcbddc->vertex_size               = 1;
2095:   pcbddc->recompute_topography      = PETSC_TRUE;
2096:   pcbddc->coarse_size               = -1;
2097:   pcbddc->use_exact_dirichlet_trick = PETSC_TRUE;
2098:   pcbddc->coarsening_ratio          = 8;
2099:   pcbddc->coarse_eqs_per_proc       = 1;
2100:   pcbddc->benign_compute_correction = PETSC_TRUE;
2101:   pcbddc->nedfield                  = -1;
2102:   pcbddc->nedglobal                 = PETSC_TRUE;
2103:   pcbddc->graphmaxcount             = PETSC_MAX_INT;
2104:   pcbddc->sub_schurs_layers         = -1;
2105:   pcbddc->ksp_D                     = kspD;
2106:   pcbddc->ksp_R                     = kspR;
2107:   pcbddc->coarse_ksp                = kspC;
2108:   PetscFunctionReturn(PETSC_SUCCESS);
2109: }

2111: static PetscErrorCode PCDestroy_BDDC(PC pc)
2112: {
2113:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

2115:   PetscFunctionBegin;
2116:   PetscCall(PCReset_BDDC(pc));
2117:   PetscCall(KSPDestroy(&pcbddc->ksp_D));
2118:   PetscCall(KSPDestroy(&pcbddc->ksp_R));
2119:   PetscCall(KSPDestroy(&pcbddc->coarse_ksp));
2120:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDiscreteGradient_C", NULL));
2121:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDivergenceMat_C", NULL));
2122:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetChangeOfBasisMat_C", NULL));
2123:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesLocalIS_C", NULL));
2124:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesIS_C", NULL));
2125:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesLocalIS_C", NULL));
2126:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesIS_C", NULL));
2127:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetCoarseningRatio_C", NULL));
2128:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevel_C", NULL));
2129:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetUseExactDirichlet_C", NULL));
2130:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevels_C", NULL));
2131:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundaries_C", NULL));
2132:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundariesLocal_C", NULL));
2133:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundaries_C", NULL));
2134:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundariesLocal_C", NULL));
2135:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundaries_C", NULL));
2136:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundariesLocal_C", NULL));
2137:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundaries_C", NULL));
2138:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundariesLocal_C", NULL));
2139:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplitting_C", NULL));
2140:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplittingLocal_C", NULL));
2141:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", NULL));
2142:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCCreateFETIDPOperators_C", NULL));
2143:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetRHS_C", NULL));
2144:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetSolution_C", NULL));
2145:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", NULL));
2146:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
2147:   PetscCall(PetscFree(pc->data));
2148:   PetscFunctionReturn(PETSC_SUCCESS);
2149: }

2151: static PetscErrorCode PCSetCoordinates_BDDC(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
2152: {
2153:   PC_BDDC    *pcbddc    = (PC_BDDC *)pc->data;
2154:   PCBDDCGraph mat_graph = pcbddc->mat_graph;

2156:   PetscFunctionBegin;
2157:   PetscCall(PetscFree(mat_graph->coords));
2158:   PetscCall(PetscMalloc1(nloc * dim, &mat_graph->coords));
2159:   PetscCall(PetscArraycpy(mat_graph->coords, coords, nloc * dim));
2160:   mat_graph->cnloc = nloc;
2161:   mat_graph->cdim  = dim;
2162:   mat_graph->cloc  = PETSC_FALSE;
2163:   /* flg setup */
2164:   pcbddc->recompute_topography = PETSC_TRUE;
2165:   pcbddc->corner_selected      = PETSC_FALSE;
2166:   PetscFunctionReturn(PETSC_SUCCESS);
2167: }

2169: static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool *change)
2170: {
2171:   PetscFunctionBegin;
2172:   *change = PETSC_TRUE;
2173:   PetscFunctionReturn(PETSC_SUCCESS);
2174: }

2176: static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
2177: {
2178:   FETIDPMat_ctx mat_ctx;
2179:   Vec           work;
2180:   PC_IS        *pcis;
2181:   PC_BDDC      *pcbddc;

2183:   PetscFunctionBegin;
2184:   PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2185:   pcis   = (PC_IS *)mat_ctx->pc->data;
2186:   pcbddc = (PC_BDDC *)mat_ctx->pc->data;

2188:   PetscCall(VecSet(fetidp_flux_rhs, 0.0));
2189:   /* copy rhs since we may change it during PCPreSolve_BDDC */
2190:   if (!pcbddc->original_rhs) PetscCall(VecDuplicate(pcis->vec1_global, &pcbddc->original_rhs));
2191:   if (mat_ctx->rhs_flip) {
2192:     PetscCall(VecPointwiseMult(pcbddc->original_rhs, standard_rhs, mat_ctx->rhs_flip));
2193:   } else {
2194:     PetscCall(VecCopy(standard_rhs, pcbddc->original_rhs));
2195:   }
2196:   if (mat_ctx->g2g_p) {
2197:     /* interface pressure rhs */
2198:     PetscCall(VecScatterBegin(mat_ctx->g2g_p, fetidp_flux_rhs, pcbddc->original_rhs, INSERT_VALUES, SCATTER_REVERSE));
2199:     PetscCall(VecScatterEnd(mat_ctx->g2g_p, fetidp_flux_rhs, pcbddc->original_rhs, INSERT_VALUES, SCATTER_REVERSE));
2200:     PetscCall(VecScatterBegin(mat_ctx->g2g_p, standard_rhs, fetidp_flux_rhs, INSERT_VALUES, SCATTER_FORWARD));
2201:     PetscCall(VecScatterEnd(mat_ctx->g2g_p, standard_rhs, fetidp_flux_rhs, INSERT_VALUES, SCATTER_FORWARD));
2202:     if (!mat_ctx->rhs_flip) PetscCall(VecScale(fetidp_flux_rhs, -1.));
2203:   }
2204:   /*
2205:      change of basis for physical rhs if needed
2206:      It also changes the rhs in case of dirichlet boundaries
2207:   */
2208:   PetscCall(PCPreSolve_BDDC(mat_ctx->pc, NULL, pcbddc->original_rhs, NULL));
2209:   if (pcbddc->ChangeOfBasisMatrix) {
2210:     PetscCall(MatMultTranspose(pcbddc->ChangeOfBasisMatrix, pcbddc->original_rhs, pcbddc->work_change));
2211:     work = pcbddc->work_change;
2212:   } else {
2213:     work = pcbddc->original_rhs;
2214:   }
2215:   /* store vectors for computation of fetidp final solution */
2216:   PetscCall(VecScatterBegin(pcis->global_to_D, work, mat_ctx->temp_solution_D, INSERT_VALUES, SCATTER_FORWARD));
2217:   PetscCall(VecScatterEnd(pcis->global_to_D, work, mat_ctx->temp_solution_D, INSERT_VALUES, SCATTER_FORWARD));
2218:   /* scale rhs since it should be unassembled */
2219:   /* TODO use counter scaling? (also below) */
2220:   PetscCall(VecScatterBegin(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2221:   PetscCall(VecScatterEnd(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2222:   /* Apply partition of unity */
2223:   PetscCall(VecPointwiseMult(mat_ctx->temp_solution_B, pcis->D, mat_ctx->temp_solution_B));
2224:   /* PetscCall(PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B)); */
2225:   if (!pcbddc->switch_static) {
2226:     /* compute partially subassembled Schur complement right-hand side */
2227:     PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2228:     PetscCall(KSPSolve(pcbddc->ksp_D, mat_ctx->temp_solution_D, pcis->vec1_D));
2229:     PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2230:     /* Cannot propagate up error in KSPSolve() because there is no access to the PC */
2231:     PetscCall(MatMult(pcis->A_BI, pcis->vec1_D, pcis->vec1_B));
2232:     PetscCall(VecAXPY(mat_ctx->temp_solution_B, -1.0, pcis->vec1_B));
2233:     PetscCall(VecSet(work, 0.0));
2234:     PetscCall(VecScatterBegin(pcis->global_to_B, mat_ctx->temp_solution_B, work, ADD_VALUES, SCATTER_REVERSE));
2235:     PetscCall(VecScatterEnd(pcis->global_to_B, mat_ctx->temp_solution_B, work, ADD_VALUES, SCATTER_REVERSE));
2236:     /* PetscCall(PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B)); */
2237:     PetscCall(VecScatterBegin(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2238:     PetscCall(VecScatterEnd(pcis->global_to_B, work, mat_ctx->temp_solution_B, INSERT_VALUES, SCATTER_FORWARD));
2239:     PetscCall(VecPointwiseMult(mat_ctx->temp_solution_B, pcis->D, mat_ctx->temp_solution_B));
2240:   }
2241:   /* BDDC rhs */
2242:   PetscCall(VecCopy(mat_ctx->temp_solution_B, pcis->vec1_B));
2243:   if (pcbddc->switch_static) PetscCall(VecCopy(mat_ctx->temp_solution_D, pcis->vec1_D));
2244:   /* apply BDDC */
2245:   PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2246:   PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, PETSC_FALSE));
2247:   PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));

2249:   /* Application of B_delta and assembling of rhs for fetidp fluxes */
2250:   PetscCall(MatMult(mat_ctx->B_delta, pcis->vec1_B, mat_ctx->lambda_local));
2251:   PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, mat_ctx->lambda_local, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2252:   PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, mat_ctx->lambda_local, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2253:   /* Add contribution to interface pressures */
2254:   if (mat_ctx->l2g_p) {
2255:     PetscCall(MatMult(mat_ctx->B_BB, pcis->vec1_B, mat_ctx->vP));
2256:     if (pcbddc->switch_static) PetscCall(MatMultAdd(mat_ctx->B_BI, pcis->vec1_D, mat_ctx->vP, mat_ctx->vP));
2257:     PetscCall(VecScatterBegin(mat_ctx->l2g_p, mat_ctx->vP, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2258:     PetscCall(VecScatterEnd(mat_ctx->l2g_p, mat_ctx->vP, fetidp_flux_rhs, ADD_VALUES, SCATTER_FORWARD));
2259:   }
2260:   PetscFunctionReturn(PETSC_SUCCESS);
2261: }

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

2266:   Collective

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

2272:   Output Parameter:
2273: . fetidp_flux_rhs - the right-hand side for the FETI-DP linear system

2275:   Level: developer

2277: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCCreateFETIDPOperators()`, `PCBDDCMatFETIDPGetSolution()`
2278: @*/
2279: PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
2280: {
2281:   FETIDPMat_ctx mat_ctx;

2283:   PetscFunctionBegin;
2287:   PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2288:   PetscUseMethod(mat_ctx->pc, "PCBDDCMatFETIDPGetRHS_C", (Mat, Vec, Vec), (fetidp_mat, standard_rhs, fetidp_flux_rhs));
2289:   PetscFunctionReturn(PETSC_SUCCESS);
2290: }

2292: static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
2293: {
2294:   FETIDPMat_ctx mat_ctx;
2295:   PC_IS        *pcis;
2296:   PC_BDDC      *pcbddc;
2297:   Vec           work;

2299:   PetscFunctionBegin;
2300:   PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2301:   pcis   = (PC_IS *)mat_ctx->pc->data;
2302:   pcbddc = (PC_BDDC *)mat_ctx->pc->data;

2304:   /* apply B_delta^T */
2305:   PetscCall(VecSet(pcis->vec1_B, 0.));
2306:   PetscCall(VecScatterBegin(mat_ctx->l2g_lambda, fetidp_flux_sol, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
2307:   PetscCall(VecScatterEnd(mat_ctx->l2g_lambda, fetidp_flux_sol, mat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
2308:   PetscCall(MatMultTranspose(mat_ctx->B_delta, mat_ctx->lambda_local, pcis->vec1_B));
2309:   if (mat_ctx->l2g_p) {
2310:     PetscCall(VecScatterBegin(mat_ctx->l2g_p, fetidp_flux_sol, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
2311:     PetscCall(VecScatterEnd(mat_ctx->l2g_p, fetidp_flux_sol, mat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
2312:     PetscCall(MatMultAdd(mat_ctx->Bt_BB, mat_ctx->vP, pcis->vec1_B, pcis->vec1_B));
2313:   }

2315:   /* compute rhs for BDDC application */
2316:   PetscCall(VecAYPX(pcis->vec1_B, -1.0, mat_ctx->temp_solution_B));
2317:   if (pcbddc->switch_static) {
2318:     PetscCall(VecCopy(mat_ctx->temp_solution_D, pcis->vec1_D));
2319:     if (mat_ctx->l2g_p) {
2320:       PetscCall(VecScale(mat_ctx->vP, -1.));
2321:       PetscCall(MatMultAdd(mat_ctx->Bt_BI, mat_ctx->vP, pcis->vec1_D, pcis->vec1_D));
2322:     }
2323:   }

2325:   /* apply BDDC */
2326:   PetscCall(PetscArrayzero(pcbddc->benign_p0, pcbddc->benign_n));
2327:   PetscCall(PCBDDCApplyInterfacePreconditioner(mat_ctx->pc, PETSC_FALSE));

2329:   /* put values into global vector */
2330:   if (pcbddc->ChangeOfBasisMatrix) work = pcbddc->work_change;
2331:   else work = standard_sol;
2332:   PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_B, work, INSERT_VALUES, SCATTER_REVERSE));
2333:   PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_B, work, INSERT_VALUES, SCATTER_REVERSE));
2334:   if (!pcbddc->switch_static) {
2335:     /* compute values into the interior if solved for the partially subassembled Schur complement */
2336:     PetscCall(MatMult(pcis->A_IB, pcis->vec1_B, pcis->vec1_D));
2337:     PetscCall(VecAYPX(pcis->vec1_D, -1.0, mat_ctx->temp_solution_D));
2338:     PetscCall(PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2339:     PetscCall(KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec1_D));
2340:     PetscCall(PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], mat_ctx->pc, 0, 0, 0));
2341:     /* Cannot propagate up error in KSPSolve() because there is no access to the PC */
2342:   }

2344:   PetscCall(VecScatterBegin(pcis->global_to_D, pcis->vec1_D, work, INSERT_VALUES, SCATTER_REVERSE));
2345:   PetscCall(VecScatterEnd(pcis->global_to_D, pcis->vec1_D, work, INSERT_VALUES, SCATTER_REVERSE));
2346:   /* add p0 solution to final solution */
2347:   PetscCall(PCBDDCBenignGetOrSetP0(mat_ctx->pc, work, PETSC_FALSE));
2348:   if (pcbddc->ChangeOfBasisMatrix) PetscCall(MatMult(pcbddc->ChangeOfBasisMatrix, work, standard_sol));
2349:   PetscCall(PCPostSolve_BDDC(mat_ctx->pc, NULL, NULL, standard_sol));
2350:   if (mat_ctx->g2g_p) {
2351:     PetscCall(VecScatterBegin(mat_ctx->g2g_p, fetidp_flux_sol, standard_sol, INSERT_VALUES, SCATTER_REVERSE));
2352:     PetscCall(VecScatterEnd(mat_ctx->g2g_p, fetidp_flux_sol, standard_sol, INSERT_VALUES, SCATTER_REVERSE));
2353:   }
2354:   PetscFunctionReturn(PETSC_SUCCESS);
2355: }

2357: static PetscErrorCode PCView_BDDCIPC(PC pc, PetscViewer viewer)
2358: {
2359:   BDDCIPC_ctx bddcipc_ctx;
2360:   PetscBool   isascii;

2362:   PetscFunctionBegin;
2363:   PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2364:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2365:   if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "BDDC interface preconditioner\n"));
2366:   PetscCall(PetscViewerASCIIPushTab(viewer));
2367:   PetscCall(PCView(bddcipc_ctx->bddc, viewer));
2368:   PetscCall(PetscViewerASCIIPopTab(viewer));
2369:   PetscFunctionReturn(PETSC_SUCCESS);
2370: }

2372: static PetscErrorCode PCSetUp_BDDCIPC(PC pc)
2373: {
2374:   BDDCIPC_ctx bddcipc_ctx;
2375:   PetscBool   isbddc;
2376:   Vec         vv;
2377:   IS          is;
2378:   PC_IS      *pcis;

2380:   PetscFunctionBegin;
2381:   PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2382:   PetscCall(PetscObjectTypeCompare((PetscObject)bddcipc_ctx->bddc, PCBDDC, &isbddc));
2383:   PetscCheck(isbddc, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Invalid type %s. Must be of type bddc", ((PetscObject)bddcipc_ctx->bddc)->type_name);
2384:   PetscCall(PCSetUp(bddcipc_ctx->bddc));

2386:   /* create interface scatter */
2387:   pcis = (PC_IS *)(bddcipc_ctx->bddc->data);
2388:   PetscCall(VecScatterDestroy(&bddcipc_ctx->g2l));
2389:   PetscCall(MatCreateVecs(pc->pmat, &vv, NULL));
2390:   PetscCall(ISRenumber(pcis->is_B_global, NULL, NULL, &is));
2391:   PetscCall(VecScatterCreate(vv, is, pcis->vec1_B, NULL, &bddcipc_ctx->g2l));
2392:   PetscCall(ISDestroy(&is));
2393:   PetscCall(VecDestroy(&vv));
2394:   PetscFunctionReturn(PETSC_SUCCESS);
2395: }

2397: static PetscErrorCode PCApply_BDDCIPC(PC pc, Vec r, Vec x)
2398: {
2399:   BDDCIPC_ctx bddcipc_ctx;
2400:   PC_IS      *pcis;
2401:   VecScatter  tmps;

2403:   PetscFunctionBegin;
2404:   PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2405:   pcis              = (PC_IS *)(bddcipc_ctx->bddc->data);
2406:   tmps              = pcis->global_to_B;
2407:   pcis->global_to_B = bddcipc_ctx->g2l;
2408:   PetscCall(PCBDDCScalingRestriction(bddcipc_ctx->bddc, r, pcis->vec1_B));
2409:   PetscCall(PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc, PETSC_FALSE));
2410:   PetscCall(PCBDDCScalingExtension(bddcipc_ctx->bddc, pcis->vec1_B, x));
2411:   pcis->global_to_B = tmps;
2412:   PetscFunctionReturn(PETSC_SUCCESS);
2413: }

2415: static PetscErrorCode PCApplyTranspose_BDDCIPC(PC pc, Vec r, Vec x)
2416: {
2417:   BDDCIPC_ctx bddcipc_ctx;
2418:   PC_IS      *pcis;
2419:   VecScatter  tmps;

2421:   PetscFunctionBegin;
2422:   PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2423:   pcis              = (PC_IS *)(bddcipc_ctx->bddc->data);
2424:   tmps              = pcis->global_to_B;
2425:   pcis->global_to_B = bddcipc_ctx->g2l;
2426:   PetscCall(PCBDDCScalingRestriction(bddcipc_ctx->bddc, r, pcis->vec1_B));
2427:   PetscCall(PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc, PETSC_TRUE));
2428:   PetscCall(PCBDDCScalingExtension(bddcipc_ctx->bddc, pcis->vec1_B, x));
2429:   pcis->global_to_B = tmps;
2430:   PetscFunctionReturn(PETSC_SUCCESS);
2431: }

2433: static PetscErrorCode PCDestroy_BDDCIPC(PC pc)
2434: {
2435:   BDDCIPC_ctx bddcipc_ctx;

2437:   PetscFunctionBegin;
2438:   PetscCall(PCShellGetContext(pc, &bddcipc_ctx));
2439:   PetscCall(PCDestroy(&bddcipc_ctx->bddc));
2440:   PetscCall(VecScatterDestroy(&bddcipc_ctx->g2l));
2441:   PetscCall(PetscFree(bddcipc_ctx));
2442:   PetscFunctionReturn(PETSC_SUCCESS);
2443: }

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

2448:   Collective

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

2454:   Output Parameter:
2455: . standard_sol - the solution defined on the physical domain

2457:   Level: developer

2459: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCCreateFETIDPOperators()`, `PCBDDCMatFETIDPGetRHS()`
2460: @*/
2461: PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
2462: {
2463:   FETIDPMat_ctx mat_ctx;

2465:   PetscFunctionBegin;
2469:   PetscCall(MatShellGetContext(fetidp_mat, &mat_ctx));
2470:   PetscUseMethod(mat_ctx->pc, "PCBDDCMatFETIDPGetSolution_C", (Mat, Vec, Vec), (fetidp_mat, fetidp_flux_sol, standard_sol));
2471:   PetscFunctionReturn(PETSC_SUCCESS);
2472: }

2474: static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc)
2475: {
2476:   FETIDPMat_ctx fetidpmat_ctx;
2477:   Mat           newmat;
2478:   FETIDPPC_ctx  fetidppc_ctx;
2479:   PC            newpc;
2480:   MPI_Comm      comm;

2482:   PetscFunctionBegin;
2483:   PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
2484:   /* FETI-DP matrix */
2485:   PetscCall(PCBDDCCreateFETIDPMatContext(pc, &fetidpmat_ctx));
2486:   fetidpmat_ctx->fully_redundant = fully_redundant;
2487:   PetscCall(PCBDDCSetupFETIDPMatContext(fetidpmat_ctx));
2488:   PetscCall(MatCreateShell(comm, fetidpmat_ctx->n, fetidpmat_ctx->n, fetidpmat_ctx->N, fetidpmat_ctx->N, fetidpmat_ctx, &newmat));
2489:   PetscCall(PetscObjectSetName((PetscObject)newmat, !fetidpmat_ctx->l2g_lambda_only ? "F" : "G"));
2490:   PetscCall(MatShellSetOperation(newmat, MATOP_MULT, (void (*)(void))FETIDPMatMult));
2491:   PetscCall(MatShellSetOperation(newmat, MATOP_MULT_TRANSPOSE, (void (*)(void))FETIDPMatMultTranspose));
2492:   PetscCall(MatShellSetOperation(newmat, MATOP_DESTROY, (void (*)(void))PCBDDCDestroyFETIDPMat));
2493:   /* propagate MatOptions */
2494:   {
2495:     PC_BDDC  *pcbddc = (PC_BDDC *)fetidpmat_ctx->pc->data;
2496:     PetscBool isset, issym;

2498:     PetscCall(MatIsSymmetricKnown(pc->mat, &isset, &issym));
2499:     if ((isset && issym) || pcbddc->symmetric_primal) PetscCall(MatSetOption(newmat, MAT_SYMMETRIC, PETSC_TRUE));
2500:   }
2501:   PetscCall(MatSetOptionsPrefix(newmat, prefix));
2502:   PetscCall(MatAppendOptionsPrefix(newmat, "fetidp_"));
2503:   PetscCall(MatSetUp(newmat));
2504:   /* FETI-DP preconditioner */
2505:   PetscCall(PCBDDCCreateFETIDPPCContext(pc, &fetidppc_ctx));
2506:   PetscCall(PCBDDCSetupFETIDPPCContext(newmat, fetidppc_ctx));
2507:   PetscCall(PCCreate(comm, &newpc));
2508:   PetscCall(PCSetOperators(newpc, newmat, newmat));
2509:   PetscCall(PCSetOptionsPrefix(newpc, prefix));
2510:   PetscCall(PCAppendOptionsPrefix(newpc, "fetidp_"));
2511:   PetscCall(PCSetErrorIfFailure(newpc, pc->erroriffailure));
2512:   if (!fetidpmat_ctx->l2g_lambda_only) { /* standard FETI-DP */
2513:     PetscCall(PCSetType(newpc, PCSHELL));
2514:     PetscCall(PCShellSetName(newpc, "FETI-DP multipliers"));
2515:     PetscCall(PCShellSetContext(newpc, fetidppc_ctx));
2516:     PetscCall(PCShellSetApply(newpc, FETIDPPCApply));
2517:     PetscCall(PCShellSetApplyTranspose(newpc, FETIDPPCApplyTranspose));
2518:     PetscCall(PCShellSetView(newpc, FETIDPPCView));
2519:     PetscCall(PCShellSetDestroy(newpc, PCBDDCDestroyFETIDPPC));
2520:   } else { /* saddle-point FETI-DP */
2521:     Mat       M;
2522:     PetscInt  psize;
2523:     PetscBool fake = PETSC_FALSE, isfieldsplit;

2525:     PetscCall(ISViewFromOptions(fetidpmat_ctx->lagrange, NULL, "-lag_view"));
2526:     PetscCall(ISViewFromOptions(fetidpmat_ctx->pressure, NULL, "-press_view"));
2527:     PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_PPmat", (PetscObject *)&M));
2528:     PetscCall(PCSetType(newpc, PCFIELDSPLIT));
2529:     PetscCall(PCFieldSplitSetIS(newpc, "lag", fetidpmat_ctx->lagrange));
2530:     PetscCall(PCFieldSplitSetIS(newpc, "p", fetidpmat_ctx->pressure));
2531:     PetscCall(PCFieldSplitSetType(newpc, PC_COMPOSITE_SCHUR));
2532:     PetscCall(PCFieldSplitSetSchurFactType(newpc, PC_FIELDSPLIT_SCHUR_FACT_DIAG));
2533:     PetscCall(ISGetSize(fetidpmat_ctx->pressure, &psize));
2534:     if (psize != M->rmap->N) {
2535:       Mat      M2;
2536:       PetscInt lpsize;

2538:       fake = PETSC_TRUE;
2539:       PetscCall(ISGetLocalSize(fetidpmat_ctx->pressure, &lpsize));
2540:       PetscCall(MatCreate(comm, &M2));
2541:       PetscCall(MatSetType(M2, MATAIJ));
2542:       PetscCall(MatSetSizes(M2, lpsize, lpsize, psize, psize));
2543:       PetscCall(MatSetUp(M2));
2544:       PetscCall(MatAssemblyBegin(M2, MAT_FINAL_ASSEMBLY));
2545:       PetscCall(MatAssemblyEnd(M2, MAT_FINAL_ASSEMBLY));
2546:       PetscCall(PCFieldSplitSetSchurPre(newpc, PC_FIELDSPLIT_SCHUR_PRE_USER, M2));
2547:       PetscCall(MatDestroy(&M2));
2548:     } else {
2549:       PetscCall(PCFieldSplitSetSchurPre(newpc, PC_FIELDSPLIT_SCHUR_PRE_USER, M));
2550:     }
2551:     PetscCall(PCFieldSplitSetSchurScale(newpc, 1.0));

2553:     /* we need to setfromoptions and setup here to access the blocks */
2554:     PetscCall(PCSetFromOptions(newpc));
2555:     PetscCall(PCSetUp(newpc));

2557:     /* user may have changed the type (e.g. -fetidp_pc_type none) */
2558:     PetscCall(PetscObjectTypeCompare((PetscObject)newpc, PCFIELDSPLIT, &isfieldsplit));
2559:     if (isfieldsplit) {
2560:       KSP      *ksps;
2561:       PC        ppc, lagpc;
2562:       PetscInt  nn;
2563:       PetscBool ismatis, matisok = PETSC_FALSE, check = PETSC_FALSE;

2565:       /* set the solver for the (0,0) block */
2566:       PetscCall(PCFieldSplitSchurGetSubKSP(newpc, &nn, &ksps));
2567:       if (!nn) { /* not of type PC_COMPOSITE_SCHUR */
2568:         PetscCall(PCFieldSplitGetSubKSP(newpc, &nn, &ksps));
2569:         if (!fake) { /* pass pmat to the pressure solver */
2570:           Mat F;

2572:           PetscCall(KSPGetOperators(ksps[1], &F, NULL));
2573:           PetscCall(KSPSetOperators(ksps[1], F, M));
2574:         }
2575:       } else {
2576:         PetscBool issym, isset;
2577:         Mat       S;

2579:         PetscCall(PCFieldSplitSchurGetS(newpc, &S));
2580:         PetscCall(MatIsSymmetricKnown(newmat, &isset, &issym));
2581:         if (isset) PetscCall(MatSetOption(S, MAT_SYMMETRIC, issym));
2582:       }
2583:       PetscCall(KSPGetPC(ksps[0], &lagpc));
2584:       PetscCall(PCSetType(lagpc, PCSHELL));
2585:       PetscCall(PCShellSetName(lagpc, "FETI-DP multipliers"));
2586:       PetscCall(PCShellSetContext(lagpc, fetidppc_ctx));
2587:       PetscCall(PCShellSetApply(lagpc, FETIDPPCApply));
2588:       PetscCall(PCShellSetApplyTranspose(lagpc, FETIDPPCApplyTranspose));
2589:       PetscCall(PCShellSetView(lagpc, FETIDPPCView));
2590:       PetscCall(PCShellSetDestroy(lagpc, PCBDDCDestroyFETIDPPC));

2592:       /* Olof's idea: interface Schur complement preconditioner for the mass matrix */
2593:       PetscCall(KSPGetPC(ksps[1], &ppc));
2594:       if (fake) {
2595:         BDDCIPC_ctx    bddcipc_ctx;
2596:         PetscContainer c;

2598:         matisok = PETSC_TRUE;

2600:         /* create inner BDDC solver */
2601:         PetscCall(PetscNew(&bddcipc_ctx));
2602:         PetscCall(PCCreate(comm, &bddcipc_ctx->bddc));
2603:         PetscCall(PCSetType(bddcipc_ctx->bddc, PCBDDC));
2604:         PetscCall(PCSetOperators(bddcipc_ctx->bddc, M, M));
2605:         PetscCall(PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_pCSR", (PetscObject *)&c));
2606:         PetscCall(PetscObjectTypeCompare((PetscObject)M, MATIS, &ismatis));
2607:         if (c && ismatis) {
2608:           Mat       lM;
2609:           PetscInt *csr, n;

2611:           PetscCall(MatISGetLocalMat(M, &lM));
2612:           PetscCall(MatGetSize(lM, &n, NULL));
2613:           PetscCall(PetscContainerGetPointer(c, (void **)&csr));
2614:           PetscCall(PCBDDCSetLocalAdjacencyGraph(bddcipc_ctx->bddc, n, csr, csr + (n + 1), PETSC_COPY_VALUES));
2615:           PetscCall(MatISRestoreLocalMat(M, &lM));
2616:         }
2617:         PetscCall(PCSetOptionsPrefix(bddcipc_ctx->bddc, ((PetscObject)ksps[1])->prefix));
2618:         PetscCall(PCSetErrorIfFailure(bddcipc_ctx->bddc, pc->erroriffailure));
2619:         PetscCall(PCSetFromOptions(bddcipc_ctx->bddc));

2621:         /* wrap the interface application */
2622:         PetscCall(PCSetType(ppc, PCSHELL));
2623:         PetscCall(PCShellSetName(ppc, "FETI-DP pressure"));
2624:         PetscCall(PCShellSetContext(ppc, bddcipc_ctx));
2625:         PetscCall(PCShellSetSetUp(ppc, PCSetUp_BDDCIPC));
2626:         PetscCall(PCShellSetApply(ppc, PCApply_BDDCIPC));
2627:         PetscCall(PCShellSetApplyTranspose(ppc, PCApplyTranspose_BDDCIPC));
2628:         PetscCall(PCShellSetView(ppc, PCView_BDDCIPC));
2629:         PetscCall(PCShellSetDestroy(ppc, PCDestroy_BDDCIPC));
2630:       }

2632:       /* determine if we need to assemble M to construct a preconditioner */
2633:       if (!matisok) {
2634:         PetscCall(PetscObjectTypeCompare((PetscObject)M, MATIS, &ismatis));
2635:         PetscCall(PetscObjectTypeCompareAny((PetscObject)ppc, &matisok, PCBDDC, PCJACOBI, PCNONE, PCMG, ""));
2636:         if (ismatis && !matisok) PetscCall(MatConvert(M, MATAIJ, MAT_INPLACE_MATRIX, &M));
2637:       }

2639:       /* run the subproblems to check convergence */
2640:       PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)newmat)->prefix, "-check_saddlepoint", &check, NULL));
2641:       if (check) {
2642:         PetscInt i;

2644:         for (i = 0; i < nn; i++) {
2645:           KSP       kspC;
2646:           PC        npc;
2647:           Mat       F, pF;
2648:           Vec       x, y;
2649:           PetscBool isschur, prec = PETSC_TRUE;

2651:           PetscCall(KSPCreate(PetscObjectComm((PetscObject)ksps[i]), &kspC));
2652:           PetscCall(KSPSetNestLevel(kspC, pc->kspnestlevel));
2653:           PetscCall(KSPSetOptionsPrefix(kspC, ((PetscObject)ksps[i])->prefix));
2654:           PetscCall(KSPAppendOptionsPrefix(kspC, "check_"));
2655:           PetscCall(KSPGetOperators(ksps[i], &F, &pF));
2656:           PetscCall(PetscObjectTypeCompare((PetscObject)F, MATSCHURCOMPLEMENT, &isschur));
2657:           if (isschur) {
2658:             KSP  kspS, kspS2;
2659:             Mat  A00, pA00, A10, A01, A11;
2660:             char prefix[256];

2662:             PetscCall(MatSchurComplementGetKSP(F, &kspS));
2663:             PetscCall(MatSchurComplementGetSubMatrices(F, &A00, &pA00, &A01, &A10, &A11));
2664:             PetscCall(MatCreateSchurComplement(A00, pA00, A01, A10, A11, &F));
2665:             PetscCall(MatSchurComplementGetKSP(F, &kspS2));
2666:             PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sschur_", ((PetscObject)kspC)->prefix));
2667:             PetscCall(KSPSetOptionsPrefix(kspS2, prefix));
2668:             PetscCall(KSPGetPC(kspS2, &npc));
2669:             PetscCall(PCSetType(npc, PCKSP));
2670:             PetscCall(PCKSPSetKSP(npc, kspS));
2671:             PetscCall(KSPSetFromOptions(kspS2));
2672:             PetscCall(KSPGetPC(kspS2, &npc));
2673:             PetscCall(PCSetUseAmat(npc, PETSC_TRUE));
2674:           } else {
2675:             PetscCall(PetscObjectReference((PetscObject)F));
2676:           }
2677:           PetscCall(KSPSetFromOptions(kspC));
2678:           PetscCall(PetscOptionsGetBool(NULL, ((PetscObject)kspC)->prefix, "-preconditioned", &prec, NULL));
2679:           if (prec) {
2680:             PetscCall(KSPGetPC(ksps[i], &npc));
2681:             PetscCall(KSPSetPC(kspC, npc));
2682:           }
2683:           PetscCall(KSPSetOperators(kspC, F, pF));
2684:           PetscCall(MatCreateVecs(F, &x, &y));
2685:           PetscCall(VecSetRandom(x, NULL));
2686:           PetscCall(MatMult(F, x, y));
2687:           PetscCall(KSPSolve(kspC, y, x));
2688:           PetscCall(KSPCheckSolve(kspC, npc, x));
2689:           PetscCall(KSPDestroy(&kspC));
2690:           PetscCall(MatDestroy(&F));
2691:           PetscCall(VecDestroy(&x));
2692:           PetscCall(VecDestroy(&y));
2693:         }
2694:       }
2695:       PetscCall(PetscFree(ksps));
2696:     }
2697:   }
2698:   /* return pointers for objects created */
2699:   *fetidp_mat = newmat;
2700:   *fetidp_pc  = newpc;
2701:   PetscFunctionReturn(PETSC_SUCCESS);
2702: }

2704: /*@C
2705:   PCBDDCCreateFETIDPOperators - Create FETI-DP operators

2707:   Collective

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

2714:   Output Parameters:
2715: + fetidp_mat - shell FETI-DP matrix object
2716: - fetidp_pc  - shell Dirichlet preconditioner for FETI-DP matrix

2718:   Level: developer

2720:   Note:
2721:   Currently the only operations provided for FETI-DP matrix are `MatMult()` and `MatMultTranspose()`

2723: .seealso: [](ch_ksp), `PCBDDC`, `PCBDDCMatFETIDPGetRHS()`, `PCBDDCMatFETIDPGetSolution()`
2724: @*/
2725: PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc)
2726: {
2727:   PetscFunctionBegin;
2729:   PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "You must call PCSetup_BDDC() first");
2730:   PetscUseMethod(pc, "PCBDDCCreateFETIDPOperators_C", (PC, PetscBool, const char *, Mat *, PC *), (pc, fully_redundant, prefix, fetidp_mat, fetidp_pc));
2731:   PetscFunctionReturn(PETSC_SUCCESS);
2732: }

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

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

2739:    It also works with unsymmetric and indefinite problems.

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

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

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

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

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

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

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

2761:    Adaptive selection of primal constraints is supported for SPD systems with high-contrast in the coefficients if MUMPS or MKL_PARDISO are present.
2762:    Future versions of the code will also consider using PASTIX.

2764:    An experimental interface to the FETI-DP method is available. FETI-DP operators could be created using `PCBDDCCreateFETIDPOperators()`.
2765:     A stand-alone class for the FETI-DP method will be provided in the next releases.

2767:    Options Database Keys:
2768: +    -pc_bddc_use_vertices <true>         - use or not vertices in primal space
2769: .    -pc_bddc_use_edges <true>            - use or not edges in primal space
2770: .    -pc_bddc_use_faces <false>           - use or not faces in primal space
2771: .    -pc_bddc_symmetric <true>            - symmetric computation of primal basis functions. Specify false for unsymmetric problems
2772: .    -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only)
2773: .    -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested
2774: .    -pc_bddc_switch_static <false>       - switches from M_2 (default) to M_3 operator (see reference article [1])
2775: .    -pc_bddc_levels <0>                  - maximum number of levels for multilevel
2776: .    -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)
2777: .    -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)
2778: .    -pc_bddc_use_deluxe_scaling <false>  - use deluxe scaling
2779: .    -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)
2780: .    -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)
2781: -    -pc_bddc_check_level <0>             - set verbosity level of debugging output

2783:    Options for Dirichlet, Neumann or coarse solver can be set using the appropriate options prefix
2784: .vb
2785:       -pc_bddc_dirichlet_
2786:       -pc_bddc_neumann_
2787:       -pc_bddc_coarse_
2788: .ve
2789:    e.g. -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. `PCBDDC` uses by default `KSPPREONLY` and `PCLU`.

2791:    When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified using the options prefix
2792: .vb
2793:       -pc_bddc_dirichlet_lN_
2794:       -pc_bddc_neumann_lN_
2795:       -pc_bddc_coarse_lN_
2796: .ve
2797:    Note that level number ranges from the finest (0) to the coarsest (N).
2798:    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
2799:    to the option, e.g.
2800: .vb
2801:      -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3
2802: .ve
2803:    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

2805:    Level: intermediate

2807:    Contributed by:
2808:    Stefano Zampini

2810: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `MATIS`, `PCLU`, `PCGAMG`, `PC`, `PCBDDCSetLocalAdjacencyGraph()`, `PCBDDCSetDofsSplitting()`,
2811:           `PCBDDCSetDirichletBoundaries()`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCSetPrimalVerticesIS()`, `MatNullSpace`, `MatSetNearNullSpace()`,
2812:           `PCBDDCSetChangeOfBasisMat()`, `PCBDDCCreateFETIDPOperators()`, `PCNN`
2813: M*/

2815: PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc)
2816: {
2817:   PC_BDDC *pcbddc;

2819:   PetscFunctionBegin;
2820:   PetscCall(PetscNew(&pcbddc));
2821:   pc->data = pcbddc;

2823:   PetscCall(PCISInitialize(pc));

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

2828:   /* BDDC nonzero defaults */
2829:   pcbddc->use_nnsp                  = PETSC_TRUE;
2830:   pcbddc->use_local_adj             = PETSC_TRUE;
2831:   pcbddc->use_vertices              = PETSC_TRUE;
2832:   pcbddc->use_edges                 = PETSC_TRUE;
2833:   pcbddc->symmetric_primal          = PETSC_TRUE;
2834:   pcbddc->vertex_size               = 1;
2835:   pcbddc->recompute_topography      = PETSC_TRUE;
2836:   pcbddc->coarse_size               = -1;
2837:   pcbddc->use_exact_dirichlet_trick = PETSC_TRUE;
2838:   pcbddc->coarsening_ratio          = 8;
2839:   pcbddc->coarse_eqs_per_proc       = 1;
2840:   pcbddc->benign_compute_correction = PETSC_TRUE;
2841:   pcbddc->nedfield                  = -1;
2842:   pcbddc->nedglobal                 = PETSC_TRUE;
2843:   pcbddc->graphmaxcount             = PETSC_MAX_INT;
2844:   pcbddc->sub_schurs_layers         = -1;
2845:   pcbddc->adaptive_threshold[0]     = 0.0;
2846:   pcbddc->adaptive_threshold[1]     = 0.0;

2848:   /* function pointers */
2849:   pc->ops->apply               = PCApply_BDDC;
2850:   pc->ops->applytranspose      = PCApplyTranspose_BDDC;
2851:   pc->ops->setup               = PCSetUp_BDDC;
2852:   pc->ops->destroy             = PCDestroy_BDDC;
2853:   pc->ops->setfromoptions      = PCSetFromOptions_BDDC;
2854:   pc->ops->view                = PCView_BDDC;
2855:   pc->ops->applyrichardson     = NULL;
2856:   pc->ops->applysymmetricleft  = NULL;
2857:   pc->ops->applysymmetricright = NULL;
2858:   pc->ops->presolve            = PCPreSolve_BDDC;
2859:   pc->ops->postsolve           = PCPostSolve_BDDC;
2860:   pc->ops->reset               = PCReset_BDDC;

2862:   /* composing function */
2863:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDiscreteGradient_C", PCBDDCSetDiscreteGradient_BDDC));
2864:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDivergenceMat_C", PCBDDCSetDivergenceMat_BDDC));
2865:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetChangeOfBasisMat_C", PCBDDCSetChangeOfBasisMat_BDDC));
2866:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesLocalIS_C", PCBDDCSetPrimalVerticesLocalIS_BDDC));
2867:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesIS_C", PCBDDCSetPrimalVerticesIS_BDDC));
2868:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesLocalIS_C", PCBDDCGetPrimalVerticesLocalIS_BDDC));
2869:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesIS_C", PCBDDCGetPrimalVerticesIS_BDDC));
2870:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetCoarseningRatio_C", PCBDDCSetCoarseningRatio_BDDC));
2871:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevel_C", PCBDDCSetLevel_BDDC));
2872:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetUseExactDirichlet_C", PCBDDCSetUseExactDirichlet_BDDC));
2873:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevels_C", PCBDDCSetLevels_BDDC));
2874:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundaries_C", PCBDDCSetDirichletBoundaries_BDDC));
2875:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundariesLocal_C", PCBDDCSetDirichletBoundariesLocal_BDDC));
2876:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundaries_C", PCBDDCSetNeumannBoundaries_BDDC));
2877:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundariesLocal_C", PCBDDCSetNeumannBoundariesLocal_BDDC));
2878:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundaries_C", PCBDDCGetDirichletBoundaries_BDDC));
2879:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundariesLocal_C", PCBDDCGetDirichletBoundariesLocal_BDDC));
2880:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundaries_C", PCBDDCGetNeumannBoundaries_BDDC));
2881:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundariesLocal_C", PCBDDCGetNeumannBoundariesLocal_BDDC));
2882:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplitting_C", PCBDDCSetDofsSplitting_BDDC));
2883:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplittingLocal_C", PCBDDCSetDofsSplittingLocal_BDDC));
2884:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", PCBDDCSetLocalAdjacencyGraph_BDDC));
2885:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCCreateFETIDPOperators_C", PCBDDCCreateFETIDPOperators_BDDC));
2886:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetRHS_C", PCBDDCMatFETIDPGetRHS_BDDC));
2887:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetSolution_C", PCBDDCMatFETIDPGetSolution_BDDC));
2888:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", PCPreSolveChangeRHS_BDDC));
2889:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_BDDC));
2890:   PetscFunctionReturn(PETSC_SUCCESS);
2891: }

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

2897:   Level: developer

2899: .seealso: [](ch_ksp), `PetscInitialize()`, `PCBDDCFinalizePackage()`
2900: @*/
2901: PetscErrorCode PCBDDCInitializePackage(void)
2902: {
2903:   int i;

2905:   PetscFunctionBegin;
2906:   if (PCBDDCPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
2907:   PCBDDCPackageInitialized = PETSC_TRUE;
2908:   PetscCall(PetscRegisterFinalize(PCBDDCFinalizePackage));

2910:   /* general events */
2911:   PetscCall(PetscLogEventRegister("PCBDDCTopo", PC_CLASSID, &PC_BDDC_Topology[0]));
2912:   PetscCall(PetscLogEventRegister("PCBDDCLKSP", PC_CLASSID, &PC_BDDC_LocalSolvers[0]));
2913:   PetscCall(PetscLogEventRegister("PCBDDCLWor", PC_CLASSID, &PC_BDDC_LocalWork[0]));
2914:   PetscCall(PetscLogEventRegister("PCBDDCCorr", PC_CLASSID, &PC_BDDC_CorrectionSetUp[0]));
2915:   PetscCall(PetscLogEventRegister("PCBDDCASet", PC_CLASSID, &PC_BDDC_ApproxSetUp[0]));
2916:   PetscCall(PetscLogEventRegister("PCBDDCAApp", PC_CLASSID, &PC_BDDC_ApproxApply[0]));
2917:   PetscCall(PetscLogEventRegister("PCBDDCCSet", PC_CLASSID, &PC_BDDC_CoarseSetUp[0]));
2918:   PetscCall(PetscLogEventRegister("PCBDDCCKSP", PC_CLASSID, &PC_BDDC_CoarseSolver[0]));
2919:   PetscCall(PetscLogEventRegister("PCBDDCAdap", PC_CLASSID, &PC_BDDC_AdaptiveSetUp[0]));
2920:   PetscCall(PetscLogEventRegister("PCBDDCScal", PC_CLASSID, &PC_BDDC_Scaling[0]));
2921:   PetscCall(PetscLogEventRegister("PCBDDCSchr", PC_CLASSID, &PC_BDDC_Schurs[0]));
2922:   PetscCall(PetscLogEventRegister("PCBDDCDirS", PC_CLASSID, &PC_BDDC_Solves[0][0]));
2923:   PetscCall(PetscLogEventRegister("PCBDDCNeuS", PC_CLASSID, &PC_BDDC_Solves[0][1]));
2924:   PetscCall(PetscLogEventRegister("PCBDDCCoaS", PC_CLASSID, &PC_BDDC_Solves[0][2]));
2925:   for (i = 1; i < PETSC_PCBDDC_MAXLEVELS; i++) {
2926:     char ename[32];

2928:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCTopo l%02d", i));
2929:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Topology[i]));
2930:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCLKSP l%02d", i));
2931:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_LocalSolvers[i]));
2932:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCLWor l%02d", i));
2933:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_LocalWork[i]));
2934:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCorr l%02d", i));
2935:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CorrectionSetUp[i]));
2936:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCASet l%02d", i));
2937:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_ApproxSetUp[i]));
2938:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCAApp l%02d", i));
2939:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_ApproxApply[i]));
2940:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCSet l%02d", i));
2941:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CoarseSetUp[i]));
2942:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCKSP l%02d", i));
2943:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_CoarseSolver[i]));
2944:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCAdap l%02d", i));
2945:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_AdaptiveSetUp[i]));
2946:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCScal l%02d", i));
2947:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Scaling[i]));
2948:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCSchr l%02d", i));
2949:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Schurs[i]));
2950:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCDirS l%02d", i));
2951:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][0]));
2952:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCNeuS l%02d", i));
2953:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][1]));
2954:     PetscCall(PetscSNPrintf(ename, sizeof(ename), "PCBDDCCoaS l%02d", i));
2955:     PetscCall(PetscLogEventRegister(ename, PC_CLASSID, &PC_BDDC_Solves[i][2]));
2956:   }
2957:   PetscFunctionReturn(PETSC_SUCCESS);
2958: }

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

2964:   Level: developer

2966: .seealso: [](ch_ksp), `PetscFinalize()`, `PCBDDCInitializePackage()`
2967: @*/
2968: PetscErrorCode PCBDDCFinalizePackage(void)
2969: {
2970:   PetscFunctionBegin;
2971:   PCBDDCPackageInitialized = PETSC_FALSE;
2972:   PetscFunctionReturn(PETSC_SUCCESS);
2973: }