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: PetscErrorCode PCApply_BDDC(PC, Vec, Vec);

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

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

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

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

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

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

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

240:     MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);

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

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

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

284:   PetscObjectReference((PetscObject)G);
285:   MatDestroy(&pcbddc->discretegradient);
286:   pcbddc->discretegradient = G;
287:   pcbddc->nedorder         = order > 0 ? order : -order;
288:   pcbddc->nedfield         = field;
289:   pcbddc->nedglobal        = global;
290:   pcbddc->conforming       = conforming;
291:   return 0;
292: }

294: /*@
295:   PCBDDCSetDiscreteGradient - Sets the discrete gradient

297:    Collective

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

307:    Level: advanced

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

317: .seealso: `PCBDDC`, `PCBDDCSetDofsSplitting()`, `PCBDDCSetDofsSplittingLocal()`, `MATAIJ`, `PCBDDCSetDivergenceMat()`
318: @*/
319: PetscErrorCode PCBDDCSetDiscreteGradient(PC pc, Mat G, PetscInt order, PetscInt field, PetscBool global, PetscBool conforming)
320: {
328:   PetscTryMethod(pc, "PCBDDCSetDiscreteGradient_C", (PC, Mat, PetscInt, PetscInt, PetscBool, PetscBool), (pc, G, order, field, global, conforming));
329:   return 0;
330: }

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

336:   PetscObjectReference((PetscObject)divudotp);
337:   MatDestroy(&pcbddc->divudotp);
338:   pcbddc->divudotp          = divudotp;
339:   pcbddc->divudotp_trans    = trans;
340:   pcbddc->compute_nonetflux = PETSC_TRUE;
341:   if (vl2l) {
342:     PetscObjectReference((PetscObject)vl2l);
343:     ISDestroy(&pcbddc->divudotp_vl2l);
344:     pcbddc->divudotp_vl2l = vl2l;
345:   }
346:   return 0;
347: }

349: /*@
350:   PCBDDCSetDivergenceMat - Sets the linear operator representing \int_\Omega \div {\bf u} \cdot p dx

352:    Collective

354:    Input Parameters:
355: +  pc - the preconditioning context
356: .  divudotp - the matrix (must be of type `MATIS`)
357: .  trans - if trans if false (resp. true), then pressures are in the test (trial) space and velocities are in the trial (test) space.
358: -  vl2l - optional index set describing the local (wrt the local matrix in divudotp) to local (wrt the local matrix
359:    in the preconditioning matrix) map for the velocities

361:    Level: advanced

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

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

368: .seealso: `PCBDDC`, `PCBDDCSetDiscreteGradient()`
369: @*/
370: PetscErrorCode PCBDDCSetDivergenceMat(PC pc, Mat divudotp, PetscBool trans, IS vl2l)
371: {
372:   PetscBool ismatis;

379:   PetscObjectTypeCompare((PetscObject)divudotp, MATIS, &ismatis);
381:   PetscTryMethod(pc, "PCBDDCSetDivergenceMat_C", (PC, Mat, PetscBool, IS), (pc, divudotp, trans, vl2l));
382:   return 0;
383: }

385: static PetscErrorCode PCBDDCSetChangeOfBasisMat_BDDC(PC pc, Mat change, PetscBool interior)
386: {
387:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

389:   PetscObjectReference((PetscObject)change);
390:   MatDestroy(&pcbddc->user_ChangeOfBasisMatrix);
391:   pcbddc->user_ChangeOfBasisMatrix = change;
392:   pcbddc->change_interior          = interior;
393:   return 0;
394: }

396: /*@
397:  PCBDDCSetChangeOfBasisMat - Set user defined change of basis for dofs

399:    Collective

401:    Input Parameters:
402: +  pc - the preconditioning context
403: .  change - the change of basis matrix
404: -  interior - whether or not the change of basis modifies interior dofs

406:    Level: intermediate

408: .seealso: `PCBDDC`
409: @*/
410: PetscErrorCode PCBDDCSetChangeOfBasisMat(PC pc, Mat change, PetscBool interior)
411: {
415:   if (pc->mat) {
416:     PetscInt rows_c, cols_c, rows, cols;
417:     MatGetSize(pc->mat, &rows, &cols);
418:     MatGetSize(change, &rows_c, &cols_c);
421:     MatGetLocalSize(pc->mat, &rows, &cols);
422:     MatGetLocalSize(change, &rows_c, &cols_c);
425:   }
426:   PetscTryMethod(pc, "PCBDDCSetChangeOfBasisMat_C", (PC, Mat, PetscBool), (pc, change, interior));
427:   return 0;
428: }

430: static PetscErrorCode PCBDDCSetPrimalVerticesIS_BDDC(PC pc, IS PrimalVertices)
431: {
432:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
433:   PetscBool isequal = PETSC_FALSE;

435:   PetscObjectReference((PetscObject)PrimalVertices);
436:   if (pcbddc->user_primal_vertices) ISEqual(PrimalVertices, pcbddc->user_primal_vertices, &isequal);
437:   ISDestroy(&pcbddc->user_primal_vertices);
438:   ISDestroy(&pcbddc->user_primal_vertices_local);
439:   pcbddc->user_primal_vertices = PrimalVertices;
440:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
441:   return 0;
442: }

444: /*@
445:  PCBDDCSetPrimalVerticesIS - Set additional user defined primal vertices in `PCBDDC`

447:    Collective

449:    Input Parameters:
450: +  pc - the preconditioning context
451: -  PrimalVertices - index set of primal vertices in global numbering (can be empty)

453:    Level: intermediate

455:    Note:
456:    Any process can list any global node

458: .seealso: `PCBDDC`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
459: @*/
460: PetscErrorCode PCBDDCSetPrimalVerticesIS(PC pc, IS PrimalVertices)
461: {
465:   PetscTryMethod(pc, "PCBDDCSetPrimalVerticesIS_C", (PC, IS), (pc, PrimalVertices));
466:   return 0;
467: }

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

473:   *is = pcbddc->user_primal_vertices;
474:   return 0;
475: }

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

480:    Collective

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

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

488:    Level: intermediate

490: .seealso: `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
491: @*/
492: PetscErrorCode PCBDDCGetPrimalVerticesIS(PC pc, IS *is)
493: {
496:   PetscUseMethod(pc, "PCBDDCGetPrimalVerticesIS_C", (PC, IS *), (pc, is));
497:   return 0;
498: }

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

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

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

517:    Collective

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

523:    Level: intermediate

525: .seealso: `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesLocalIS()`
526: @*/
527: PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices)
528: {
532:   PetscTryMethod(pc, "PCBDDCSetPrimalVerticesLocalIS_C", (PC, IS), (pc, PrimalVertices));
533:   return 0;
534: }

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

540:   *is = pcbddc->user_primal_vertices_local;
541:   return 0;
542: }

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

547:    Collective

549:    Input Parameter:
550: .  pc - the preconditioning context

552:    Output Parameter:
553: .  is - index set of primal vertices in local numbering (NULL if not set)

555:    Level: intermediate

557: .seealso: `PCBDDC`, `PCBDDCSetPrimalVerticesIS()`, `PCBDDCGetPrimalVerticesIS()`, `PCBDDCSetPrimalVerticesLocalIS()`
558: @*/
559: PetscErrorCode PCBDDCGetPrimalVerticesLocalIS(PC pc, IS *is)
560: {
563:   PetscUseMethod(pc, "PCBDDCGetPrimalVerticesLocalIS_C", (PC, IS *), (pc, is));
564:   return 0;
565: }

567: static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc, PetscInt k)
568: {
569:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

571:   pcbddc->coarsening_ratio = k;
572:   return 0;
573: }

575: /*@
576:   PCBDDCSetCoarseningRatio - Set coarsening ratio used in multilevel version

578:    Logically collective

580:    Input Parameters:
581: +  pc - the preconditioning context
582: -  k - coarsening ratio (H/h at the coarser level)

584:    Options Database Key:
585: .    -pc_bddc_coarsening_ratio <int> - Set coarsening ratio used in multilevel coarsening

587:    Level: intermediate

589:    Note:
590:    Approximately k subdomains at the finer level will be aggregated into a single subdomain at the coarser level

592: .seealso: `PCBDDC`, `PCBDDCSetLevels()`
593: @*/
594: PetscErrorCode PCBDDCSetCoarseningRatio(PC pc, PetscInt k)
595: {
598:   PetscTryMethod(pc, "PCBDDCSetCoarseningRatio_C", (PC, PetscInt), (pc, k));
599:   return 0;
600: }

602: /* The following functions (PCBDDCSetUseExactDirichlet PCBDDCSetLevel) are not public */
603: static PetscErrorCode PCBDDCSetUseExactDirichlet_BDDC(PC pc, PetscBool flg)
604: {
605:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

607:   pcbddc->use_exact_dirichlet_trick = flg;
608:   return 0;
609: }

611: PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc, PetscBool flg)
612: {
615:   PetscTryMethod(pc, "PCBDDCSetUseExactDirichlet_C", (PC, PetscBool), (pc, flg));
616:   return 0;
617: }

619: static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc, PetscInt level)
620: {
621:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

623:   pcbddc->current_level = level;
624:   return 0;
625: }

627: PetscErrorCode PCBDDCSetLevel(PC pc, PetscInt level)
628: {
631:   PetscTryMethod(pc, "PCBDDCSetLevel_C", (PC, PetscInt), (pc, level));
632:   return 0;
633: }

635: static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc, PetscInt levels)
636: {
637:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

640:   pcbddc->max_levels = levels;
641:   return 0;
642: }

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

647:    Logically collective

649:    Input Parameters:
650: +  pc - the preconditioning context
651: -  levels - the maximum number of levels

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

656:    Level: intermediate

658:    Note:
659:    The default value is 0, that gives the classical two-levels BDDC

661: .seealso: `PCBDDC`, `PCBDDCSetCoarseningRatio()`
662: @*/
663: PetscErrorCode PCBDDCSetLevels(PC pc, PetscInt levels)
664: {
667:   PetscTryMethod(pc, "PCBDDCSetLevels_C", (PC, PetscInt), (pc, levels));
668:   return 0;
669: }

671: static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc, IS DirichletBoundaries)
672: {
673:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
674:   PetscBool isequal = PETSC_FALSE;

676:   PetscObjectReference((PetscObject)DirichletBoundaries);
677:   if (pcbddc->DirichletBoundaries) ISEqual(DirichletBoundaries, pcbddc->DirichletBoundaries, &isequal);
678:   /* last user setting takes precedence -> destroy any other customization */
679:   ISDestroy(&pcbddc->DirichletBoundariesLocal);
680:   ISDestroy(&pcbddc->DirichletBoundaries);
681:   pcbddc->DirichletBoundaries = DirichletBoundaries;
682:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
683:   return 0;
684: }

686: /*@
687:  PCBDDCSetDirichletBoundaries - Set IS defining Dirichlet boundaries for the global problem.

689:    Collective

691:    Input Parameters:
692: +  pc - the preconditioning context
693: -  DirichletBoundaries - parallel IS defining the Dirichlet boundaries

695:    Level: intermediate

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

700: .seealso: `PCBDDC`, `PCBDDCSetDirichletBoundariesLocal()`, `MatZeroRows()`, `MatZeroRowsColumns()`
701: @*/
702: PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc, IS DirichletBoundaries)
703: {
707:   PetscTryMethod(pc, "PCBDDCSetDirichletBoundaries_C", (PC, IS), (pc, DirichletBoundaries));
708:   return 0;
709: }

711: static PetscErrorCode PCBDDCSetDirichletBoundariesLocal_BDDC(PC pc, IS DirichletBoundaries)
712: {
713:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
714:   PetscBool isequal = PETSC_FALSE;

716:   PetscObjectReference((PetscObject)DirichletBoundaries);
717:   if (pcbddc->DirichletBoundariesLocal) ISEqual(DirichletBoundaries, pcbddc->DirichletBoundariesLocal, &isequal);
718:   /* last user setting takes precedence -> destroy any other customization */
719:   ISDestroy(&pcbddc->DirichletBoundariesLocal);
720:   ISDestroy(&pcbddc->DirichletBoundaries);
721:   pcbddc->DirichletBoundariesLocal = DirichletBoundaries;
722:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
723:   return 0;
724: }

726: /*@
727:  PCBDDCSetDirichletBoundariesLocal - Set IS defining Dirichlet boundaries for the global problem in local ordering.

729:    Collective

731:    Input Parameters:
732: +  pc - the preconditioning context
733: -  DirichletBoundaries - parallel IS defining the Dirichlet boundaries (in local ordering)

735:    Level: intermediate

737: .seealso: `PCBDDC`, `PCBDDCSetDirichletBoundaries()`, `MatZeroRows()`, `MatZeroRowsColumns()`
738: @*/
739: PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc, IS DirichletBoundaries)
740: {
744:   PetscTryMethod(pc, "PCBDDCSetDirichletBoundariesLocal_C", (PC, IS), (pc, DirichletBoundaries));
745:   return 0;
746: }

748: static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc, IS NeumannBoundaries)
749: {
750:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
751:   PetscBool isequal = PETSC_FALSE;

753:   PetscObjectReference((PetscObject)NeumannBoundaries);
754:   if (pcbddc->NeumannBoundaries) ISEqual(NeumannBoundaries, pcbddc->NeumannBoundaries, &isequal);
755:   /* last user setting takes precedence -> destroy any other customization */
756:   ISDestroy(&pcbddc->NeumannBoundariesLocal);
757:   ISDestroy(&pcbddc->NeumannBoundaries);
758:   pcbddc->NeumannBoundaries = NeumannBoundaries;
759:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
760:   return 0;
761: }

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

766:    Collective

768:    Input Parameters:
769: +  pc - the preconditioning context
770: -  NeumannBoundaries - parallel IS defining the Neumann boundaries

772:    Level: intermediate

774:    Note:
775:    Any process can list any global node

777: .seealso: `PCBDDC`, `PCBDDCSetNeumannBoundariesLocal()`
778: @*/
779: PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc, IS NeumannBoundaries)
780: {
784:   PetscTryMethod(pc, "PCBDDCSetNeumannBoundaries_C", (PC, IS), (pc, NeumannBoundaries));
785:   return 0;
786: }

788: static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc, IS NeumannBoundaries)
789: {
790:   PC_BDDC  *pcbddc  = (PC_BDDC *)pc->data;
791:   PetscBool isequal = PETSC_FALSE;

793:   PetscObjectReference((PetscObject)NeumannBoundaries);
794:   if (pcbddc->NeumannBoundariesLocal) ISEqual(NeumannBoundaries, pcbddc->NeumannBoundariesLocal, &isequal);
795:   /* last user setting takes precedence -> destroy any other customization */
796:   ISDestroy(&pcbddc->NeumannBoundariesLocal);
797:   ISDestroy(&pcbddc->NeumannBoundaries);
798:   pcbddc->NeumannBoundariesLocal = NeumannBoundaries;
799:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
800:   return 0;
801: }

803: /*@
804:  PCBDDCSetNeumannBoundariesLocal - Set IS defining Neumann boundaries for the global problem in local ordering.

806:    Collective

808:    Input Parameters:
809: +  pc - the preconditioning context
810: -  NeumannBoundaries - parallel IS defining the subdomain part of Neumann boundaries (in local ordering)

812:    Level: intermediate

814: .seealso: `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCGetDirichletBoundaries()`
815: @*/
816: PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc, IS NeumannBoundaries)
817: {
821:   PetscTryMethod(pc, "PCBDDCSetNeumannBoundariesLocal_C", (PC, IS), (pc, NeumannBoundaries));
822:   return 0;
823: }

825: static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc, IS *DirichletBoundaries)
826: {
827:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

829:   *DirichletBoundaries = pcbddc->DirichletBoundaries;
830:   return 0;
831: }

833: /*@
834:    PCBDDCGetDirichletBoundaries - Get parallel `IS` for Dirichlet boundaries

836:    Collective

838:    Input Parameter:
839: .  pc - the preconditioning context

841:    Output Parameter:
842: .  DirichletBoundaries - index set defining the Dirichlet boundaries

844:    Level: intermediate

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

849: .seealso: `PCBDDC`, `PCBDDCSetDirichletBoundaries()`
850: @*/
851: PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc, IS *DirichletBoundaries)
852: {
854:   PetscUseMethod(pc, "PCBDDCGetDirichletBoundaries_C", (PC, IS *), (pc, DirichletBoundaries));
855:   return 0;
856: }

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

862:   *DirichletBoundaries = pcbddc->DirichletBoundariesLocal;
863:   return 0;
864: }

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

869:    Collective

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

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

877:    Level: intermediate

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

884: .seealso: `PCBDDC`, `PCBDDCGetDirichletBoundariesLocal()`, `PCBDDCGetDirichletBoundaries()`, `PCBDDCSetDirichletBoundaries()`
885: @*/
886: PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc, IS *DirichletBoundaries)
887: {
889:   PetscUseMethod(pc, "PCBDDCGetDirichletBoundariesLocal_C", (PC, IS *), (pc, DirichletBoundaries));
890:   return 0;
891: }

893: static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc, IS *NeumannBoundaries)
894: {
895:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

897:   *NeumannBoundaries = pcbddc->NeumannBoundaries;
898:   return 0;
899: }

901: /*@
902:    PCBDDCGetNeumannBoundaries - Get parallel `IS` for Neumann boundaries

904:    Not Collective

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

909:    Output Parameter:
910: .  NeumannBoundaries - index set defining the Neumann boundaries

912:    Level: intermediate

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

917: .seealso: `PCBDDC`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCGetDirichletBoundaries()`, `PCBDDCSetDirichletBoundaries()`
918: @*/
919: PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc, IS *NeumannBoundaries)
920: {
922:   PetscUseMethod(pc, "PCBDDCGetNeumannBoundaries_C", (PC, IS *), (pc, NeumannBoundaries));
923:   return 0;
924: }

926: static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc, IS *NeumannBoundaries)
927: {
928:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

930:   *NeumannBoundaries = pcbddc->NeumannBoundariesLocal;
931:   return 0;
932: }

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

937:    Not Collective

939:    Input Parameter:
940: .  pc - the preconditioning context

942:    Output Parameter:
943: .  NeumannBoundaries - index set defining the subdomain part of Neumann boundaries

945:    Level: intermediate

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

952: .seealso: `PCBDDC``PCBDDCSetNeumannBoundaries()`, `PCBDDCSetNeumannBoundariesLocal)`, `PCBDDCGetNeumannBoundaries()`
953: @*/
954: PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc, IS *NeumannBoundaries)
955: {
957:   PetscUseMethod(pc, "PCBDDCGetNeumannBoundariesLocal_C", (PC, IS *), (pc, NeumannBoundaries));
958:   return 0;
959: }

961: static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs, const PetscInt xadj[], const PetscInt adjncy[], PetscCopyMode copymode)
962: {
963:   PC_BDDC    *pcbddc    = (PC_BDDC *)pc->data;
964:   PCBDDCGraph mat_graph = pcbddc->mat_graph;
965:   PetscBool   same_data = PETSC_FALSE;

967:   if (!nvtxs) {
968:     if (copymode == PETSC_OWN_POINTER) {
969:       PetscFree(xadj);
970:       PetscFree(adjncy);
971:     }
972:     PCBDDCGraphResetCSR(mat_graph);
973:     return 0;
974:   }
975:   if (mat_graph->nvtxs == nvtxs && mat_graph->freecsr) { /* we own the data */
976:     if (mat_graph->xadj == xadj && mat_graph->adjncy == adjncy) same_data = PETSC_TRUE;
977:     if (!same_data && mat_graph->xadj[nvtxs] == xadj[nvtxs]) {
978:       PetscArraycmp(xadj, mat_graph->xadj, nvtxs + 1, &same_data);
979:       if (same_data) PetscArraycmp(adjncy, mat_graph->adjncy, xadj[nvtxs], &same_data);
980:     }
981:   }
982:   if (!same_data) {
983:     /* free old CSR */
984:     PCBDDCGraphResetCSR(mat_graph);
985:     /* get CSR into graph structure */
986:     if (copymode == PETSC_COPY_VALUES) {
987:       PetscMalloc1(nvtxs + 1, &mat_graph->xadj);
988:       PetscMalloc1(xadj[nvtxs], &mat_graph->adjncy);
989:       PetscArraycpy(mat_graph->xadj, xadj, nvtxs + 1);
990:       PetscArraycpy(mat_graph->adjncy, adjncy, xadj[nvtxs]);
991:       mat_graph->freecsr = PETSC_TRUE;
992:     } else if (copymode == PETSC_OWN_POINTER) {
993:       mat_graph->xadj    = (PetscInt *)xadj;
994:       mat_graph->adjncy  = (PetscInt *)adjncy;
995:       mat_graph->freecsr = PETSC_TRUE;
996:     } else if (copymode == PETSC_USE_POINTER) {
997:       mat_graph->xadj    = (PetscInt *)xadj;
998:       mat_graph->adjncy  = (PetscInt *)adjncy;
999:       mat_graph->freecsr = PETSC_FALSE;
1000:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported copy mode %d", copymode);
1001:     mat_graph->nvtxs_csr         = nvtxs;
1002:     pcbddc->recompute_topography = PETSC_TRUE;
1003:   }
1004:   return 0;
1005: }

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

1010:    Not collective

1012:    Input Parameters:
1013: +  pc - the preconditioning context.
1014: .  nvtxs - number of local vertices of the graph (i.e., the number of local dofs).
1015: .  xadj, adjncy - the connectivity of the dofs in CSR format.
1016: -  copymode - supported modes are `PETSC_COPY_VALUES`, `PETSC_USE_POINTER` or `PETSC_OWN_POINTER`.

1018:    Level: intermediate

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

1023: .seealso: `PCBDDC`, `PetscCopyMode`
1024: @*/
1025: PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc, PetscInt nvtxs, const PetscInt xadj[], const PetscInt adjncy[], PetscCopyMode copymode)
1026: {
1027:   void (*f)(void) = NULL;

1030:   if (nvtxs) {
1033:   }
1034:   PetscTryMethod(pc, "PCBDDCSetLocalAdjacencyGraph_C", (PC, PetscInt, const PetscInt[], const PetscInt[], PetscCopyMode), (pc, nvtxs, xadj, adjncy, copymode));
1035:   /* free arrays if PCBDDC is not the PC type */
1036:   PetscObjectQueryFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", &f);
1037:   if (!f && copymode == PETSC_OWN_POINTER) {
1038:     PetscFree(xadj);
1039:     PetscFree(adjncy);
1040:   }
1041:   return 0;
1042: }

1044: static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc, PetscInt n_is, IS ISForDofs[])
1045: {
1046:   PC_BDDC  *pcbddc = (PC_BDDC *)pc->data;
1047:   PetscInt  i;
1048:   PetscBool isequal = PETSC_FALSE;

1050:   if (pcbddc->n_ISForDofsLocal == n_is) {
1051:     for (i = 0; i < n_is; i++) {
1052:       PetscBool isequalt;
1053:       ISEqual(ISForDofs[i], pcbddc->ISForDofsLocal[i], &isequalt);
1054:       if (!isequalt) break;
1055:     }
1056:     if (i == n_is) isequal = PETSC_TRUE;
1057:   }
1058:   for (i = 0; i < n_is; i++) PetscObjectReference((PetscObject)ISForDofs[i]);
1059:   /* Destroy ISes if they were already set */
1060:   for (i = 0; i < pcbddc->n_ISForDofsLocal; i++) ISDestroy(&pcbddc->ISForDofsLocal[i]);
1061:   PetscFree(pcbddc->ISForDofsLocal);
1062:   /* last user setting takes precedence -> destroy any other customization */
1063:   for (i = 0; i < pcbddc->n_ISForDofs; i++) ISDestroy(&pcbddc->ISForDofs[i]);
1064:   PetscFree(pcbddc->ISForDofs);
1065:   pcbddc->n_ISForDofs = 0;
1066:   /* allocate space then set */
1067:   if (n_is) PetscMalloc1(n_is, &pcbddc->ISForDofsLocal);
1068:   for (i = 0; i < n_is; i++) pcbddc->ISForDofsLocal[i] = ISForDofs[i];
1069:   pcbddc->n_ISForDofsLocal = n_is;
1070:   if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
1071:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1072:   return 0;
1073: }

1075: /*@
1076:    PCBDDCSetDofsSplittingLocal - Set index sets defining fields of the local subdomain matrix

1078:    Collective

1080:    Input Parameters:
1081: +  pc - the preconditioning context
1082: .  n_is - number of index sets defining the fields, must be the same on all MPI ranks
1083: -  ISForDofs - array of `IS` describing the fields in local ordering

1085:    Level: intermediate

1087:    Note:
1088:    Not all nodes need to be listed: unlisted nodes will belong to the complement field.

1090: .seealso: `PCBDDC`, `PCBDDCSetDofsSplitting()`
1091: @*/
1092: PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc, PetscInt n_is, IS ISForDofs[])
1093: {
1094:   PetscInt i;

1098:   for (i = 0; i < n_is; i++) {
1101:   }
1102:   PetscTryMethod(pc, "PCBDDCSetDofsSplittingLocal_C", (PC, PetscInt, IS[]), (pc, n_is, ISForDofs));
1103:   return 0;
1104: }

1106: static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc, PetscInt n_is, IS ISForDofs[])
1107: {
1108:   PC_BDDC  *pcbddc = (PC_BDDC *)pc->data;
1109:   PetscInt  i;
1110:   PetscBool isequal = PETSC_FALSE;

1112:   if (pcbddc->n_ISForDofs == n_is) {
1113:     for (i = 0; i < n_is; i++) {
1114:       PetscBool isequalt;
1115:       ISEqual(ISForDofs[i], pcbddc->ISForDofs[i], &isequalt);
1116:       if (!isequalt) break;
1117:     }
1118:     if (i == n_is) isequal = PETSC_TRUE;
1119:   }
1120:   for (i = 0; i < n_is; i++) PetscObjectReference((PetscObject)ISForDofs[i]);
1121:   /* Destroy ISes if they were already set */
1122:   for (i = 0; i < pcbddc->n_ISForDofs; i++) ISDestroy(&pcbddc->ISForDofs[i]);
1123:   PetscFree(pcbddc->ISForDofs);
1124:   /* last user setting takes precedence -> destroy any other customization */
1125:   for (i = 0; i < pcbddc->n_ISForDofsLocal; i++) ISDestroy(&pcbddc->ISForDofsLocal[i]);
1126:   PetscFree(pcbddc->ISForDofsLocal);
1127:   pcbddc->n_ISForDofsLocal = 0;
1128:   /* allocate space then set */
1129:   if (n_is) PetscMalloc1(n_is, &pcbddc->ISForDofs);
1130:   for (i = 0; i < n_is; i++) pcbddc->ISForDofs[i] = ISForDofs[i];
1131:   pcbddc->n_ISForDofs = n_is;
1132:   if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
1133:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
1134:   return 0;
1135: }

1137: /*@
1138:    PCBDDCSetDofsSplitting - Set index sets defining fields of the global matrix

1140:    Collective

1142:    Input Parameters:
1143: +  pc - the preconditioning context
1144: .  n_is - number of index sets defining the fields
1145: -  ISForDofs - array of IS describing the fields in global ordering

1147:    Level: intermediate

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

1152: .seealso: `PCBDDC`, `PCBDDCSetDofsSplittingLocal()`
1153: @*/
1154: PetscErrorCode PCBDDCSetDofsSplitting(PC pc, PetscInt n_is, IS ISForDofs[])
1155: {
1156:   PetscInt i;

1160:   for (i = 0; i < n_is; i++) {
1163:   }
1164:   PetscTryMethod(pc, "PCBDDCSetDofsSplitting_C", (PC, PetscInt, IS[]), (pc, n_is, ISForDofs));
1165:   return 0;
1166: }

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

1172:    Input Parameter:
1173: +  pc - the preconditioner context

1175:    Note:
1176:      The interface routine PCPreSolve() is not usually called directly by
1177:    the user, but instead is called by KSPSolve().
1178: */
1179: static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
1180: {
1181:   PC_BDDC  *pcbddc = (PC_BDDC *)pc->data;
1182:   PC_IS    *pcis   = (PC_IS *)(pc->data);
1183:   Vec       used_vec;
1184:   PetscBool iscg, save_rhs = PETSC_TRUE, benign_correction_computed;

1186:   /* if we are working with CG, one dirichlet solve can be avoided during Krylov iterations */
1187:   if (ksp) {
1188:     PetscObjectTypeCompareAny((PetscObject)ksp, &iscg, KSPCG, KSPGROPPCG, KSPPIPECG, KSPPIPELCG, KSPPIPECGRR, "");
1189:     if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || !iscg || pc->mat != pc->pmat) PCBDDCSetUseExactDirichlet(pc, PETSC_FALSE);
1190:   }
1191:   if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || pc->mat != pc->pmat) PCBDDCSetUseExactDirichlet(pc, PETSC_FALSE);

1193:   /* Creates parallel work vectors used in presolve */
1194:   if (!pcbddc->original_rhs) VecDuplicate(pcis->vec1_global, &pcbddc->original_rhs);
1195:   if (!pcbddc->temp_solution) VecDuplicate(pcis->vec1_global, &pcbddc->temp_solution);

1197:   pcbddc->temp_solution_used = PETSC_FALSE;
1198:   if (x) {
1199:     PetscObjectReference((PetscObject)x);
1200:     used_vec = x;
1201:   } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */
1202:     PetscObjectReference((PetscObject)pcbddc->temp_solution);
1203:     used_vec = pcbddc->temp_solution;
1204:     VecSet(used_vec, 0.0);
1205:     pcbddc->temp_solution_used = PETSC_TRUE;
1206:     VecCopy(rhs, pcbddc->original_rhs);
1207:     save_rhs                  = PETSC_FALSE;
1208:     pcbddc->eliminate_dirdofs = PETSC_TRUE;
1209:   }

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

1218:   pcbddc->rhs_change = PETSC_FALSE;
1219:   /* Take into account zeroed rows -> change rhs and store solution removed */
1220:   if (rhs && pcbddc->eliminate_dirdofs) {
1221:     IS dirIS = NULL;

1223:     /* DirichletBoundariesLocal may not be consistent among neighbours; gets a dirichlet dofs IS from graph (may be cached) */
1224:     PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph, &dirIS);
1225:     if (dirIS) {
1226:       Mat_IS            *matis = (Mat_IS *)pc->pmat->data;
1227:       PetscInt           dirsize, i, *is_indices;
1228:       PetscScalar       *array_x;
1229:       const PetscScalar *array_diagonal;

1231:       MatGetDiagonal(pc->pmat, pcis->vec1_global);
1232:       VecPointwiseDivide(pcis->vec1_global, rhs, pcis->vec1_global);
1233:       VecScatterBegin(matis->rctx, pcis->vec1_global, pcis->vec2_N, INSERT_VALUES, SCATTER_FORWARD);
1234:       VecScatterEnd(matis->rctx, pcis->vec1_global, pcis->vec2_N, INSERT_VALUES, SCATTER_FORWARD);
1235:       VecScatterBegin(matis->rctx, used_vec, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD);
1236:       VecScatterEnd(matis->rctx, used_vec, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD);
1237:       ISGetLocalSize(dirIS, &dirsize);
1238:       VecGetArray(pcis->vec1_N, &array_x);
1239:       VecGetArrayRead(pcis->vec2_N, &array_diagonal);
1240:       ISGetIndices(dirIS, (const PetscInt **)&is_indices);
1241:       for (i = 0; i < dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]];
1242:       ISRestoreIndices(dirIS, (const PetscInt **)&is_indices);
1243:       VecRestoreArrayRead(pcis->vec2_N, &array_diagonal);
1244:       VecRestoreArray(pcis->vec1_N, &array_x);
1245:       VecScatterBegin(matis->rctx, pcis->vec1_N, used_vec, INSERT_VALUES, SCATTER_REVERSE);
1246:       VecScatterEnd(matis->rctx, pcis->vec1_N, used_vec, INSERT_VALUES, SCATTER_REVERSE);
1247:       pcbddc->rhs_change = PETSC_TRUE;
1248:       ISDestroy(&dirIS);
1249:     }
1250:   }

1252:   /* remove the computed solution or the initial guess from the rhs */
1253:   if (pcbddc->rhs_change || (ksp && pcbddc->ksp_guess_nonzero)) {
1254:     /* save the original rhs */
1255:     if (save_rhs) {
1256:       VecSwap(rhs, pcbddc->original_rhs);
1257:       save_rhs = PETSC_FALSE;
1258:     }
1259:     pcbddc->rhs_change = PETSC_TRUE;
1260:     VecScale(used_vec, -1.0);
1261:     MatMultAdd(pc->mat, used_vec, pcbddc->original_rhs, rhs);
1262:     VecScale(used_vec, -1.0);
1263:     VecCopy(used_vec, pcbddc->temp_solution);
1264:     pcbddc->temp_solution_used = PETSC_TRUE;
1265:     if (ksp) KSPSetInitialGuessNonzero(ksp, PETSC_FALSE);
1266:   }
1267:   VecDestroy(&used_vec);

1269:   /* compute initial vector in benign space if needed
1270:      and remove non-benign solution from the rhs */
1271:   benign_correction_computed = PETSC_FALSE;
1272:   if (rhs && pcbddc->benign_compute_correction && (pcbddc->benign_have_null || pcbddc->benign_apply_coarse_only)) {
1273:     /* compute u^*_h using ideas similar to those in Xuemin Tu's PhD thesis (see Section 4.8.1)
1274:        Recursively apply BDDC in the multilevel case */
1275:     if (!pcbddc->benign_vec) VecDuplicate(rhs, &pcbddc->benign_vec);
1276:     /* keep applying coarse solver unless we no longer have benign subdomains */
1277:     pcbddc->benign_apply_coarse_only = pcbddc->benign_have_null ? PETSC_TRUE : PETSC_FALSE;
1278:     if (!pcbddc->benign_skip_correction) {
1279:       PCApply_BDDC(pc, rhs, pcbddc->benign_vec);
1280:       benign_correction_computed = PETSC_TRUE;
1281:       if (pcbddc->temp_solution_used) VecAXPY(pcbddc->temp_solution, 1.0, pcbddc->benign_vec);
1282:       VecScale(pcbddc->benign_vec, -1.0);
1283:       /* store the original rhs if not done earlier */
1284:       if (save_rhs) VecSwap(rhs, pcbddc->original_rhs);
1285:       if (pcbddc->rhs_change) {
1286:         MatMultAdd(pc->mat, pcbddc->benign_vec, rhs, rhs);
1287:       } else {
1288:         MatMultAdd(pc->mat, pcbddc->benign_vec, pcbddc->original_rhs, rhs);
1289:       }
1290:       pcbddc->rhs_change = PETSC_TRUE;
1291:     }
1292:     pcbddc->benign_apply_coarse_only = PETSC_FALSE;
1293:   } else {
1294:     VecDestroy(&pcbddc->benign_vec);
1295:   }

1297:   /* dbg output */
1298:   if (pcbddc->dbg_flag && benign_correction_computed) {
1299:     Vec v;

1301:     VecDuplicate(pcis->vec1_global, &v);
1302:     if (pcbddc->ChangeOfBasisMatrix) {
1303:       MatMultTranspose(pcbddc->ChangeOfBasisMatrix, rhs, v);
1304:     } else {
1305:       VecCopy(rhs, v);
1306:     }
1307:     PCBDDCBenignGetOrSetP0(pc, v, PETSC_TRUE);
1308:     PetscViewerASCIIPrintf(pcbddc->dbg_viewer, "LEVEL %" PetscInt_FMT ": is the correction benign?\n", pcbddc->current_level);
1309:     PetscScalarView(pcbddc->benign_n, pcbddc->benign_p0, pcbddc->dbg_viewer);
1310:     PetscViewerFlush(pcbddc->dbg_viewer);
1311:     VecDestroy(&v);
1312:   }

1314:   /* set initial guess if using PCG */
1315:   pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1316:   if (x && pcbddc->use_exact_dirichlet_trick) {
1317:     VecSet(x, 0.0);
1318:     if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) {
1319:       if (benign_correction_computed) { /* we have already saved the changed rhs */
1320:         VecLockReadPop(pcis->vec1_global);
1321:       } else {
1322:         MatMultTranspose(pcbddc->ChangeOfBasisMatrix, rhs, pcis->vec1_global);
1323:       }
1324:       VecScatterBegin(pcis->global_to_D, pcis->vec1_global, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD);
1325:       VecScatterEnd(pcis->global_to_D, pcis->vec1_global, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD);
1326:     } else {
1327:       VecScatterBegin(pcis->global_to_D, rhs, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD);
1328:       VecScatterEnd(pcis->global_to_D, rhs, pcis->vec1_D, INSERT_VALUES, SCATTER_FORWARD);
1329:     }
1330:     PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0);
1331:     KSPSolve(pcbddc->ksp_D, pcis->vec1_D, pcis->vec2_D);
1332:     PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0], pc, 0, 0, 0);
1333:     KSPCheckSolve(pcbddc->ksp_D, pc, pcis->vec2_D);
1334:     if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) {
1335:       VecSet(pcis->vec1_global, 0.);
1336:       VecScatterBegin(pcis->global_to_D, pcis->vec2_D, pcis->vec1_global, INSERT_VALUES, SCATTER_REVERSE);
1337:       VecScatterEnd(pcis->global_to_D, pcis->vec2_D, pcis->vec1_global, INSERT_VALUES, SCATTER_REVERSE);
1338:       MatMult(pcbddc->ChangeOfBasisMatrix, pcis->vec1_global, x);
1339:     } else {
1340:       VecScatterBegin(pcis->global_to_D, pcis->vec2_D, x, INSERT_VALUES, SCATTER_REVERSE);
1341:       VecScatterEnd(pcis->global_to_D, pcis->vec2_D, x, INSERT_VALUES, SCATTER_REVERSE);
1342:     }
1343:     if (ksp) KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
1344:     pcbddc->exact_dirichlet_trick_app = PETSC_TRUE;
1345:   } else if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior && benign_correction_computed && pcbddc->use_exact_dirichlet_trick) {
1346:     VecLockReadPop(pcis->vec1_global);
1347:   }
1348:   return 0;
1349: }

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

1355:    Input Parameter:
1356: +  pc - the preconditioner context

1358:    Application Interface Routine: PCPostSolve()

1360:    Note:
1361:      The interface routine PCPostSolve() is not usually called directly by
1362:      the user, but instead is called by KSPSolve().
1363: */
1364: static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
1365: {
1366:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

1368:   /* add solution removed in presolve */
1369:   if (x && pcbddc->rhs_change) {
1370:     if (pcbddc->temp_solution_used) {
1371:       VecAXPY(x, 1.0, pcbddc->temp_solution);
1372:     } else if (pcbddc->benign_compute_correction && pcbddc->benign_vec) {
1373:       VecAXPY(x, -1.0, pcbddc->benign_vec);
1374:     }
1375:     /* restore to original state (not for FETI-DP) */
1376:     if (ksp) pcbddc->temp_solution_used = PETSC_FALSE;
1377:   }

1379:   /* restore rhs to its original state (not needed for FETI-DP) */
1380:   if (rhs && pcbddc->rhs_change) {
1381:     VecSwap(rhs, pcbddc->original_rhs);
1382:     pcbddc->rhs_change = PETSC_FALSE;
1383:   }
1384:   /* restore ksp guess state */
1385:   if (ksp) {
1386:     KSPSetInitialGuessNonzero(ksp, pcbddc->ksp_guess_nonzero);
1387:     /* reset flag for exact dirichlet trick */
1388:     pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
1389:   }
1390:   return 0;
1391: }

1393: /*
1394:    PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner
1395:                   by setting data structures and options.

1397:    Input Parameter:
1398: +  pc - the preconditioner context

1400:    Application Interface Routine: PCSetUp()

1402:    Note:
1403:      The interface routine PCSetUp() is not usually called directly by
1404:      the user, but instead is called by PCApply() if necessary.
1405: */
1406: PetscErrorCode PCSetUp_BDDC(PC pc)
1407: {
1408:   PC_BDDC        *pcbddc = (PC_BDDC *)pc->data;
1409:   PCBDDCSubSchurs sub_schurs;
1410:   Mat_IS         *matis;
1411:   MatNullSpace    nearnullspace;
1412:   Mat             lA;
1413:   IS              lP, zerodiag = NULL;
1414:   PetscInt        nrows, ncols;
1415:   PetscMPIInt     size;
1416:   PetscBool       computesubschurs;
1417:   PetscBool       computeconstraintsmatrix;
1418:   PetscBool       new_nearnullspace_provided, ismatis, rl;
1419:   PetscBool       isset, issym, isspd;

1421:   PetscObjectTypeCompare((PetscObject)pc->pmat, MATIS, &ismatis);
1423:   MatGetSize(pc->pmat, &nrows, &ncols);
1425:   MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1683:   if (pcbddc->dbg_flag) {
1684:     PetscViewerASCIISubtractTab(pcbddc->dbg_viewer, 2 * pcbddc->current_level);
1685:     PetscViewerASCIIPopSynchronized(pcbddc->dbg_viewer);
1686:   }
1687:   return 0;
1688: }

1690: /*
1691:    PCApply_BDDC - Applies the BDDC operator to a vector.

1693:    Input Parameters:
1694: +  pc - the preconditioner context
1695: -  r - input vector (global)

1697:    Output Parameter:
1698: .  z - output vector (global)

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

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

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

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

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

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

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

1851:   if (pcbddc->ChangeOfBasisMatrix) {
1852:     pcbddc->work_change = r;
1853:     VecCopy(z, pcbddc->work_change);
1854:     MatMult(pcbddc->ChangeOfBasisMatrix, pcbddc->work_change, z);
1855:   }
1856:   return 0;
1857: }

1859: /*
1860:    PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector.

1862:    Input Parameters:
1863: +  pc - the preconditioner context
1864: -  r - input vector (global)

1866:    Output Parameter:
1867: .  z - output vector (global)

1869:    Application Interface Routine: PCApplyTranspose()
1870:  */
1871: PetscErrorCode PCApplyTranspose_BDDC(PC pc, Vec r, Vec z)
1872: {
1873:   PC_IS            *pcis   = (PC_IS *)(pc->data);
1874:   PC_BDDC          *pcbddc = (PC_BDDC *)(pc->data);
1875:   Mat               lA     = NULL;
1876:   PetscInt          n_B = pcis->n_B, n_D = pcis->n - n_B;
1877:   const PetscScalar one   = 1.0;
1878:   const PetscScalar m_one = -1.0;
1879:   const PetscScalar zero  = 0.0;

1881:   PetscCitationsRegister(citation, &cited);
1882:   if (pcbddc->switch_static) MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat, &lA);
1883:   if (pcbddc->ChangeOfBasisMatrix) {
1884:     Vec swap;

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

1943:   /* Apply interface preconditioner
1944:      input/output vecs: pcis->vec1_B and pcis->vec1_D */
1945:   PCBDDCApplyInterfacePreconditioner(pc, PETSC_TRUE);

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

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

2012: PetscErrorCode PCReset_BDDC(PC pc)
2013: {
2014:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;
2015:   PC_IS   *pcis   = (PC_IS *)pc->data;
2016:   KSP      kspD, kspR, kspC;

2018:   /* free BDDC custom data  */
2019:   PCBDDCResetCustomization(pc);
2020:   /* destroy objects related to topography */
2021:   PCBDDCResetTopography(pc);
2022:   /* destroy objects for scaling operator */
2023:   PCBDDCScalingDestroy(pc);
2024:   /* free solvers stuff */
2025:   PCBDDCResetSolvers(pc);
2026:   /* free global vectors needed in presolve */
2027:   VecDestroy(&pcbddc->temp_solution);
2028:   VecDestroy(&pcbddc->original_rhs);
2029:   /* free data created by PCIS */
2030:   PCISDestroy(pc);

2032:   /* restore defaults */
2033:   kspD = pcbddc->ksp_D;
2034:   kspR = pcbddc->ksp_R;
2035:   kspC = pcbddc->coarse_ksp;
2036:   PetscMemzero(pc->data, sizeof(*pcbddc));
2037:   pcis->n_neigh                     = -1;
2038:   pcis->scaling_factor              = 1.0;
2039:   pcis->reusesubmatrices            = PETSC_TRUE;
2040:   pcbddc->use_local_adj             = PETSC_TRUE;
2041:   pcbddc->use_vertices              = PETSC_TRUE;
2042:   pcbddc->use_edges                 = PETSC_TRUE;
2043:   pcbddc->symmetric_primal          = PETSC_TRUE;
2044:   pcbddc->vertex_size               = 1;
2045:   pcbddc->recompute_topography      = PETSC_TRUE;
2046:   pcbddc->coarse_size               = -1;
2047:   pcbddc->use_exact_dirichlet_trick = PETSC_TRUE;
2048:   pcbddc->coarsening_ratio          = 8;
2049:   pcbddc->coarse_eqs_per_proc       = 1;
2050:   pcbddc->benign_compute_correction = PETSC_TRUE;
2051:   pcbddc->nedfield                  = -1;
2052:   pcbddc->nedglobal                 = PETSC_TRUE;
2053:   pcbddc->graphmaxcount             = PETSC_MAX_INT;
2054:   pcbddc->sub_schurs_layers         = -1;
2055:   pcbddc->ksp_D                     = kspD;
2056:   pcbddc->ksp_R                     = kspR;
2057:   pcbddc->coarse_ksp                = kspC;
2058:   return 0;
2059: }

2061: PetscErrorCode PCDestroy_BDDC(PC pc)
2062: {
2063:   PC_BDDC *pcbddc = (PC_BDDC *)pc->data;

2065:   PCReset_BDDC(pc);
2066:   KSPDestroy(&pcbddc->ksp_D);
2067:   KSPDestroy(&pcbddc->ksp_R);
2068:   KSPDestroy(&pcbddc->coarse_ksp);
2069:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDiscreteGradient_C", NULL);
2070:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDivergenceMat_C", NULL);
2071:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetChangeOfBasisMat_C", NULL);
2072:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesLocalIS_C", NULL);
2073:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetPrimalVerticesIS_C", NULL);
2074:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesLocalIS_C", NULL);
2075:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetPrimalVerticesIS_C", NULL);
2076:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetCoarseningRatio_C", NULL);
2077:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevel_C", NULL);
2078:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetUseExactDirichlet_C", NULL);
2079:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLevels_C", NULL);
2080:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundaries_C", NULL);
2081:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDirichletBoundariesLocal_C", NULL);
2082:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundaries_C", NULL);
2083:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetNeumannBoundariesLocal_C", NULL);
2084:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundaries_C", NULL);
2085:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetDirichletBoundariesLocal_C", NULL);
2086:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundaries_C", NULL);
2087:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCGetNeumannBoundariesLocal_C", NULL);
2088:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplitting_C", NULL);
2089:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetDofsSplittingLocal_C", NULL);
2090:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCSetLocalAdjacencyGraph_C", NULL);
2091:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCCreateFETIDPOperators_C", NULL);
2092:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetRHS_C", NULL);
2093:   PetscObjectComposeFunction((PetscObject)pc, "PCBDDCMatFETIDPGetSolution_C", NULL);
2094:   PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", NULL);
2095:   PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL);
2096:   PetscFree(pc->data);
2097:   return 0;
2098: }

2100: static PetscErrorCode PCSetCoordinates_BDDC(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
2101: {
2102:   PC_BDDC    *pcbddc    = (PC_BDDC *)pc->data;
2103:   PCBDDCGraph mat_graph = pcbddc->mat_graph;

2105:   PetscFree(mat_graph->coords);
2106:   PetscMalloc1(nloc * dim, &mat_graph->coords);
2107:   PetscArraycpy(mat_graph->coords, coords, nloc * dim);
2108:   mat_graph->cnloc = nloc;
2109:   mat_graph->cdim  = dim;
2110:   mat_graph->cloc  = PETSC_FALSE;
2111:   /* flg setup */
2112:   pcbddc->recompute_topography = PETSC_TRUE;
2113:   pcbddc->corner_selected      = PETSC_FALSE;
2114:   return 0;
2115: }

2117: static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool *change)
2118: {
2119:   *change = PETSC_TRUE;
2120:   return 0;
2121: }

2123: static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
2124: {
2125:   FETIDPMat_ctx mat_ctx;
2126:   Vec           work;
2127:   PC_IS        *pcis;
2128:   PC_BDDC      *pcbddc;

2130:   MatShellGetContext(fetidp_mat, &mat_ctx);
2131:   pcis   = (PC_IS *)mat_ctx->pc->data;
2132:   pcbddc = (PC_BDDC *)mat_ctx->pc->data;

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

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

2209: /*@
2210:    PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETI-DP linear system using the physical right-hand side

2212:    Collective

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

2218:    Output Parameter:
2219: .  fetidp_flux_rhs - the right-hand side for the FETI-DP linear system

2221:    Level: developer

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

2232:   MatShellGetContext(fetidp_mat, &mat_ctx);
2233:   PetscUseMethod(mat_ctx->pc, "PCBDDCMatFETIDPGetRHS_C", (Mat, Vec, Vec), (fetidp_mat, standard_rhs, fetidp_flux_rhs));
2234:   return 0;
2235: }

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

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

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

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

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

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

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

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

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

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

2323:   PCShellGetContext(pc, &bddcipc_ctx);
2324:   PetscObjectTypeCompare((PetscObject)bddcipc_ctx->bddc, PCBDDC, &isbddc);
2326:   PCSetUp(bddcipc_ctx->bddc);

2328:   /* create interface scatter */
2329:   pcis = (PC_IS *)(bddcipc_ctx->bddc->data);
2330:   VecScatterDestroy(&bddcipc_ctx->g2l);
2331:   MatCreateVecs(pc->pmat, &vv, NULL);
2332:   ISRenumber(pcis->is_B_global, NULL, NULL, &is);
2333:   VecScatterCreate(vv, is, pcis->vec1_B, NULL, &bddcipc_ctx->g2l);
2334:   ISDestroy(&is);
2335:   VecDestroy(&vv);
2336:   return 0;
2337: }

2339: static PetscErrorCode PCApply_BDDCIPC(PC pc, Vec r, Vec x)
2340: {
2341:   BDDCIPC_ctx bddcipc_ctx;
2342:   PC_IS      *pcis;
2343:   VecScatter  tmps;

2345:   PCShellGetContext(pc, &bddcipc_ctx);
2346:   pcis              = (PC_IS *)(bddcipc_ctx->bddc->data);
2347:   tmps              = pcis->global_to_B;
2348:   pcis->global_to_B = bddcipc_ctx->g2l;
2349:   PCBDDCScalingRestriction(bddcipc_ctx->bddc, r, pcis->vec1_B);
2350:   PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc, PETSC_FALSE);
2351:   PCBDDCScalingExtension(bddcipc_ctx->bddc, pcis->vec1_B, x);
2352:   pcis->global_to_B = tmps;
2353:   return 0;
2354: }

2356: static PetscErrorCode PCApplyTranspose_BDDCIPC(PC pc, Vec r, Vec x)
2357: {
2358:   BDDCIPC_ctx bddcipc_ctx;
2359:   PC_IS      *pcis;
2360:   VecScatter  tmps;

2362:   PCShellGetContext(pc, &bddcipc_ctx);
2363:   pcis              = (PC_IS *)(bddcipc_ctx->bddc->data);
2364:   tmps              = pcis->global_to_B;
2365:   pcis->global_to_B = bddcipc_ctx->g2l;
2366:   PCBDDCScalingRestriction(bddcipc_ctx->bddc, r, pcis->vec1_B);
2367:   PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc, PETSC_TRUE);
2368:   PCBDDCScalingExtension(bddcipc_ctx->bddc, pcis->vec1_B, x);
2369:   pcis->global_to_B = tmps;
2370:   return 0;
2371: }

2373: static PetscErrorCode PCDestroy_BDDCIPC(PC pc)
2374: {
2375:   BDDCIPC_ctx bddcipc_ctx;

2377:   PCShellGetContext(pc, &bddcipc_ctx);
2378:   PCDestroy(&bddcipc_ctx->bddc);
2379:   VecScatterDestroy(&bddcipc_ctx->g2l);
2380:   PetscFree(bddcipc_ctx);
2381:   return 0;
2382: }

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

2387:    Collective

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

2393:    Output Parameter:
2394: .  standard_sol    - the solution defined on the physical domain

2396:    Level: developer

2398: .seealso: `PCBDDC`, `PCBDDCCreateFETIDPOperators()`, `PCBDDCMatFETIDPGetRHS()`
2399: @*/
2400: PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
2401: {
2402:   FETIDPMat_ctx mat_ctx;

2407:   MatShellGetContext(fetidp_mat, &mat_ctx);
2408:   PetscUseMethod(mat_ctx->pc, "PCBDDCMatFETIDPGetSolution_C", (Mat, Vec, Vec), (fetidp_mat, fetidp_flux_sol, standard_sol));
2409:   return 0;
2410: }

2412: static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc)
2413: {
2414:   FETIDPMat_ctx fetidpmat_ctx;
2415:   Mat           newmat;
2416:   FETIDPPC_ctx  fetidppc_ctx;
2417:   PC            newpc;
2418:   MPI_Comm      comm;

2420:   PetscObjectGetComm((PetscObject)pc, &comm);
2421:   /* FETI-DP matrix */
2422:   PCBDDCCreateFETIDPMatContext(pc, &fetidpmat_ctx);
2423:   fetidpmat_ctx->fully_redundant = fully_redundant;
2424:   PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);
2425:   MatCreateShell(comm, fetidpmat_ctx->n, fetidpmat_ctx->n, fetidpmat_ctx->N, fetidpmat_ctx->N, fetidpmat_ctx, &newmat);
2426:   PetscObjectSetName((PetscObject)newmat, !fetidpmat_ctx->l2g_lambda_only ? "F" : "G");
2427:   MatShellSetOperation(newmat, MATOP_MULT, (void (*)(void))FETIDPMatMult);
2428:   MatShellSetOperation(newmat, MATOP_MULT_TRANSPOSE, (void (*)(void))FETIDPMatMultTranspose);
2429:   MatShellSetOperation(newmat, MATOP_DESTROY, (void (*)(void))PCBDDCDestroyFETIDPMat);
2430:   /* propagate MatOptions */
2431:   {
2432:     PC_BDDC  *pcbddc = (PC_BDDC *)fetidpmat_ctx->pc->data;
2433:     PetscBool isset, issym;

2435:     MatIsSymmetricKnown(pc->mat, &isset, &issym);
2436:     if ((isset && issym) || pcbddc->symmetric_primal) MatSetOption(newmat, MAT_SYMMETRIC, PETSC_TRUE);
2437:   }
2438:   MatSetOptionsPrefix(newmat, prefix);
2439:   MatAppendOptionsPrefix(newmat, "fetidp_");
2440:   MatSetUp(newmat);
2441:   /* FETI-DP preconditioner */
2442:   PCBDDCCreateFETIDPPCContext(pc, &fetidppc_ctx);
2443:   PCBDDCSetupFETIDPPCContext(newmat, fetidppc_ctx);
2444:   PCCreate(comm, &newpc);
2445:   PCSetOperators(newpc, newmat, newmat);
2446:   PCSetOptionsPrefix(newpc, prefix);
2447:   PCAppendOptionsPrefix(newpc, "fetidp_");
2448:   PCSetErrorIfFailure(newpc, pc->erroriffailure);
2449:   if (!fetidpmat_ctx->l2g_lambda_only) { /* standard FETI-DP */
2450:     PCSetType(newpc, PCSHELL);
2451:     PCShellSetName(newpc, "FETI-DP multipliers");
2452:     PCShellSetContext(newpc, fetidppc_ctx);
2453:     PCShellSetApply(newpc, FETIDPPCApply);
2454:     PCShellSetApplyTranspose(newpc, FETIDPPCApplyTranspose);
2455:     PCShellSetView(newpc, FETIDPPCView);
2456:     PCShellSetDestroy(newpc, PCBDDCDestroyFETIDPPC);
2457:   } else { /* saddle-point FETI-DP */
2458:     Mat       M;
2459:     PetscInt  psize;
2460:     PetscBool fake = PETSC_FALSE, isfieldsplit;

2462:     ISViewFromOptions(fetidpmat_ctx->lagrange, NULL, "-lag_view");
2463:     ISViewFromOptions(fetidpmat_ctx->pressure, NULL, "-press_view");
2464:     PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_PPmat", (PetscObject *)&M);
2465:     PCSetType(newpc, PCFIELDSPLIT);
2466:     PCFieldSplitSetIS(newpc, "lag", fetidpmat_ctx->lagrange);
2467:     PCFieldSplitSetIS(newpc, "p", fetidpmat_ctx->pressure);
2468:     PCFieldSplitSetType(newpc, PC_COMPOSITE_SCHUR);
2469:     PCFieldSplitSetSchurFactType(newpc, PC_FIELDSPLIT_SCHUR_FACT_DIAG);
2470:     ISGetSize(fetidpmat_ctx->pressure, &psize);
2471:     if (psize != M->rmap->N) {
2472:       Mat      M2;
2473:       PetscInt lpsize;

2475:       fake = PETSC_TRUE;
2476:       ISGetLocalSize(fetidpmat_ctx->pressure, &lpsize);
2477:       MatCreate(comm, &M2);
2478:       MatSetType(M2, MATAIJ);
2479:       MatSetSizes(M2, lpsize, lpsize, psize, psize);
2480:       MatSetUp(M2);
2481:       MatAssemblyBegin(M2, MAT_FINAL_ASSEMBLY);
2482:       MatAssemblyEnd(M2, MAT_FINAL_ASSEMBLY);
2483:       PCFieldSplitSetSchurPre(newpc, PC_FIELDSPLIT_SCHUR_PRE_USER, M2);
2484:       MatDestroy(&M2);
2485:     } else {
2486:       PCFieldSplitSetSchurPre(newpc, PC_FIELDSPLIT_SCHUR_PRE_USER, M);
2487:     }
2488:     PCFieldSplitSetSchurScale(newpc, 1.0);

2490:     /* we need to setfromoptions and setup here to access the blocks */
2491:     PCSetFromOptions(newpc);
2492:     PCSetUp(newpc);

2494:     /* user may have changed the type (e.g. -fetidp_pc_type none) */
2495:     PetscObjectTypeCompare((PetscObject)newpc, PCFIELDSPLIT, &isfieldsplit);
2496:     if (isfieldsplit) {
2497:       KSP      *ksps;
2498:       PC        ppc, lagpc;
2499:       PetscInt  nn;
2500:       PetscBool ismatis, matisok = PETSC_FALSE, check = PETSC_FALSE;

2502:       /* set the solver for the (0,0) block */
2503:       PCFieldSplitSchurGetSubKSP(newpc, &nn, &ksps);
2504:       if (!nn) { /* not of type PC_COMPOSITE_SCHUR */
2505:         PCFieldSplitGetSubKSP(newpc, &nn, &ksps);
2506:         if (!fake) { /* pass pmat to the pressure solver */
2507:           Mat F;

2509:           KSPGetOperators(ksps[1], &F, NULL);
2510:           KSPSetOperators(ksps[1], F, M);
2511:         }
2512:       } else {
2513:         PetscBool issym, isset;
2514:         Mat       S;

2516:         PCFieldSplitSchurGetS(newpc, &S);
2517:         MatIsSymmetricKnown(newmat, &isset, &issym);
2518:         if (isset) MatSetOption(S, MAT_SYMMETRIC, issym);
2519:       }
2520:       KSPGetPC(ksps[0], &lagpc);
2521:       PCSetType(lagpc, PCSHELL);
2522:       PCShellSetName(lagpc, "FETI-DP multipliers");
2523:       PCShellSetContext(lagpc, fetidppc_ctx);
2524:       PCShellSetApply(lagpc, FETIDPPCApply);
2525:       PCShellSetApplyTranspose(lagpc, FETIDPPCApplyTranspose);
2526:       PCShellSetView(lagpc, FETIDPPCView);
2527:       PCShellSetDestroy(lagpc, PCBDDCDestroyFETIDPPC);

2529:       /* Olof's idea: interface Schur complement preconditioner for the mass matrix */
2530:       KSPGetPC(ksps[1], &ppc);
2531:       if (fake) {
2532:         BDDCIPC_ctx    bddcipc_ctx;
2533:         PetscContainer c;

2535:         matisok = PETSC_TRUE;

2537:         /* create inner BDDC solver */
2538:         PetscNew(&bddcipc_ctx);
2539:         PCCreate(comm, &bddcipc_ctx->bddc);
2540:         PCSetType(bddcipc_ctx->bddc, PCBDDC);
2541:         PCSetOperators(bddcipc_ctx->bddc, M, M);
2542:         PetscObjectQuery((PetscObject)pc, "__KSPFETIDP_pCSR", (PetscObject *)&c);
2543:         PetscObjectTypeCompare((PetscObject)M, MATIS, &ismatis);
2544:         if (c && ismatis) {
2545:           Mat       lM;
2546:           PetscInt *csr, n;

2548:           MatISGetLocalMat(M, &lM);
2549:           MatGetSize(lM, &n, NULL);
2550:           PetscContainerGetPointer(c, (void **)&csr);
2551:           PCBDDCSetLocalAdjacencyGraph(bddcipc_ctx->bddc, n, csr, csr + (n + 1), PETSC_COPY_VALUES);
2552:           MatISRestoreLocalMat(M, &lM);
2553:         }
2554:         PCSetOptionsPrefix(bddcipc_ctx->bddc, ((PetscObject)ksps[1])->prefix);
2555:         PCSetErrorIfFailure(bddcipc_ctx->bddc, pc->erroriffailure);
2556:         PCSetFromOptions(bddcipc_ctx->bddc);

2558:         /* wrap the interface application */
2559:         PCSetType(ppc, PCSHELL);
2560:         PCShellSetName(ppc, "FETI-DP pressure");
2561:         PCShellSetContext(ppc, bddcipc_ctx);
2562:         PCShellSetSetUp(ppc, PCSetUp_BDDCIPC);
2563:         PCShellSetApply(ppc, PCApply_BDDCIPC);
2564:         PCShellSetApplyTranspose(ppc, PCApplyTranspose_BDDCIPC);
2565:         PCShellSetView(ppc, PCView_BDDCIPC);
2566:         PCShellSetDestroy(ppc, PCDestroy_BDDCIPC);
2567:       }

2569:       /* determine if we need to assemble M to construct a preconditioner */
2570:       if (!matisok) {
2571:         PetscObjectTypeCompare((PetscObject)M, MATIS, &ismatis);
2572:         PetscObjectTypeCompareAny((PetscObject)ppc, &matisok, PCBDDC, PCJACOBI, PCNONE, PCMG, "");
2573:         if (ismatis && !matisok) MatConvert(M, MATAIJ, MAT_INPLACE_MATRIX, &M);
2574:       }

2576:       /* run the subproblems to check convergence */
2577:       PetscOptionsGetBool(NULL, ((PetscObject)newmat)->prefix, "-check_saddlepoint", &check, NULL);
2578:       if (check) {
2579:         PetscInt i;

2581:         for (i = 0; i < nn; i++) {
2582:           KSP       kspC;
2583:           PC        pc;
2584:           Mat       F, pF;
2585:           Vec       x, y;
2586:           PetscBool isschur, prec = PETSC_TRUE;

2588:           KSPCreate(PetscObjectComm((PetscObject)ksps[i]), &kspC);
2589:           KSPSetOptionsPrefix(kspC, ((PetscObject)ksps[i])->prefix);
2590:           KSPAppendOptionsPrefix(kspC, "check_");
2591:           KSPGetOperators(ksps[i], &F, &pF);
2592:           PetscObjectTypeCompare((PetscObject)F, MATSCHURCOMPLEMENT, &isschur);
2593:           if (isschur) {
2594:             KSP  kspS, kspS2;
2595:             Mat  A00, pA00, A10, A01, A11;
2596:             char prefix[256];

2598:             MatSchurComplementGetKSP(F, &kspS);
2599:             MatSchurComplementGetSubMatrices(F, &A00, &pA00, &A01, &A10, &A11);
2600:             MatCreateSchurComplement(A00, pA00, A01, A10, A11, &F);
2601:             MatSchurComplementGetKSP(F, &kspS2);
2602:             PetscSNPrintf(prefix, sizeof(prefix), "%sschur_", ((PetscObject)kspC)->prefix);
2603:             KSPSetOptionsPrefix(kspS2, prefix);
2604:             KSPGetPC(kspS2, &pc);
2605:             PCSetType(pc, PCKSP);
2606:             PCKSPSetKSP(pc, kspS);
2607:             KSPSetFromOptions(kspS2);
2608:             KSPGetPC(kspS2, &pc);
2609:             PCSetUseAmat(pc, PETSC_TRUE);
2610:           } else {
2611:             PetscObjectReference((PetscObject)F);
2612:           }
2613:           KSPSetFromOptions(kspC);
2614:           PetscOptionsGetBool(NULL, ((PetscObject)kspC)->prefix, "-preconditioned", &prec, NULL);
2615:           if (prec) {
2616:             KSPGetPC(ksps[i], &pc);
2617:             KSPSetPC(kspC, pc);
2618:           }
2619:           KSPSetOperators(kspC, F, pF);
2620:           MatCreateVecs(F, &x, &y);
2621:           VecSetRandom(x, NULL);
2622:           MatMult(F, x, y);
2623:           KSPSolve(kspC, y, x);
2624:           KSPCheckSolve(kspC, pc, x);
2625:           KSPDestroy(&kspC);
2626:           MatDestroy(&F);
2627:           VecDestroy(&x);
2628:           VecDestroy(&y);
2629:         }
2630:       }
2631:       PetscFree(ksps);
2632:     }
2633:   }
2634:   /* return pointers for objects created */
2635:   *fetidp_mat = newmat;
2636:   *fetidp_pc  = newpc;
2637:   return 0;
2638: }

2640: /*@C
2641:  PCBDDCCreateFETIDPOperators - Create FETI-DP operators

2643:    Collective

2645:    Input Parameters:
2646: +  pc - the BDDC preconditioning context (setup should have been called before)
2647: .  fully_redundant - true for a fully redundant set of Lagrange multipliers
2648: -  prefix - optional options database prefix for the objects to be created (can be NULL)

2650:    Output Parameters:
2651: +  fetidp_mat - shell FETI-DP matrix object
2652: -  fetidp_pc  - shell Dirichlet preconditioner for FETI-DP matrix

2654:    Level: developer

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

2659: .seealso: `PCBDDC`, `PCBDDCMatFETIDPGetRHS()`, `PCBDDCMatFETIDPGetSolution()`
2660: @*/
2661: PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc)
2662: {
2664:   if (pc->setupcalled) {
2665:     PetscUseMethod(pc, "PCBDDCCreateFETIDPOperators_C", (PC, PetscBool, const char *, Mat *, PC *), (pc, fully_redundant, prefix, fetidp_mat, fetidp_pc));
2666:   } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "You must call PCSetup_BDDC() first");
2667:   return 0;
2668: }

2670: /*MC
2671:    PCBDDC - Balancing Domain Decomposition by Constraints preconditioners

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

2675:    It also works with unsymmetric and indefinite problems.

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

2680:    Approximate local solvers are automatically adapted (see [1]) if the user has attached a nullspace object to the subdomain matrices, and informed
2681:    `PCBDDC` of using approximate solvers (via the command line).

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

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

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

2691:    Change of basis is performed similarly to [2] when requested. When more than one constraint is present on a single connected component
2692:    (i.e. an edge or a face), a robust method based on local QR factorizations is used.
2693:    User defined change of basis can be passed to `PCBDDC` with `PCBDDCSetChangeOfBasisMat()`

2695:    The PETSc implementation also supports multilevel `PCBDDC` [3]. Coarse grids are partitioned using a `MatPartitioning` object.

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

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

2703:    Options Database Keys:
2704: +    -pc_bddc_use_vertices <true> - use or not vertices in primal space
2705: .    -pc_bddc_use_edges <true> - use or not edges in primal space
2706: .    -pc_bddc_use_faces <false> - use or not faces in primal space
2707: .    -pc_bddc_symmetric <true> - symmetric computation of primal basis functions. Specify false for unsymmetric problems
2708: .    -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only)
2709: .    -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested
2710: .    -pc_bddc_switch_static <false> - switches from M_2 (default) to M_3 operator (see reference article [1])
2711: .    -pc_bddc_levels <0> - maximum number of levels for multilevel
2712: .    -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)
2713: .    -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)
2714: .    -pc_bddc_use_deluxe_scaling <false> - use deluxe scaling
2715: .    -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)
2716: .    -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)
2717: -    -pc_bddc_check_level <0> - set verbosity level of debugging output

2719:    Options for Dirichlet, Neumann or coarse solver can be set using the appropriate options prefix
2720: .vb
2721:       -pc_bddc_dirichlet_
2722:       -pc_bddc_neumann_
2723:       -pc_bddc_coarse_
2724: .ve
2725:    e.g. -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. `PCBDDC` uses by default `KSPPREONLY` and `PCLU`.

2727:    When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified using the options prefix
2728: .vb
2729:       -pc_bddc_dirichlet_lN_
2730:       -pc_bddc_neumann_lN_
2731:       -pc_bddc_coarse_lN_
2732: .ve
2733:    Note that level number ranges from the finest (0) to the coarsest (N).
2734:    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
2735:    to the option, e.g.
2736: .vb
2737:      -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3
2738: .ve
2739:    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

2741:    References:
2742: +  * - C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149--168, March 2007
2743: .  * - A. Klawonn and O. B. Widlund. "Dual-Primal FETI Methods for Linear Elasticity", Communications on Pure and Applied Mathematics Volume 59, Issue 11, pages 1523--1572, November 2006
2744: .  * - J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", Computing Volume 83, Issue 2--3, pages 55--85, November 2008
2745: -  * - C. Pechstein and C. R. Dohrmann. "Modern domain decomposition methods BDDC, deluxe scaling, and an algebraic approach", Seminar talk, Linz, December 2013, http://people.ricam.oeaw.ac.at/c.pechstein/pechstein-bddc2013.pdf

2747:    Level: intermediate

2749:    Contributed by Stefano Zampini

2751:  .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `MATIS`, `PCLU`, `PGGAMG`, `PC`, `PCBDDCSetLocalAdjacencyGraph()`, `PCBDDCSetDofsSplitting()`,
2752:             `PCBDDCSetDirichletBoundaries()`, `PCBDDCSetNeumannBoundaries()`, `PCBDDCSetPrimalVerticesIS()`, `MatNullSpace`, `MatSetNearNullSpace()`,
2753:             `PCBDDCSetChangeOfBasisMat()`, `PCBDDCCreateFETIDPOperators()`, `PCNN`
2754: M*/

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

2760:   PetscNew(&pcbddc);
2761:   pc->data = pcbddc;

2763:   /* create PCIS data structure */
2764:   PCISCreate(pc);

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

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

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

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

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

2838:  Level: developer

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

2846:   if (PCBDDCPackageInitialized) return 0;
2847:   PCBDDCPackageInitialized = PETSC_TRUE;
2848:   PetscRegisterFinalize(PCBDDCFinalizePackage);

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

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

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

2904:     Level: developer

2906:  .seealso: `PetscFinalize()`, `PCBDDCInitializePackage()`
2907: @*/
2908: PetscErrorCode PCBDDCFinalizePackage(void)
2909: {
2910:   PCBDDCPackageInitialized = PETSC_FALSE;
2911:   return 0;
2912: }