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 <../src/ksp/pc/impls/bddc/bddc.h>
     23: #include <../src/ksp/pc/impls/bddc/bddcprivate.h>
     24: #include <petscblaslapack.h>

     26: static PetscBool PCBDDCPackageInitialized = PETSC_FALSE;

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

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

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

     58: PetscErrorCode PCApply_BDDC(PC,Vec,Vec);

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

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

    135: static PetscErrorCode PCView_BDDC(PC pc,PetscViewer viewer)
    136: {
    137:   PC_BDDC              *pcbddc = (PC_BDDC*)pc->data;
    138:   PC_IS                *pcis = (PC_IS*)pc->data;
    139:   PetscErrorCode       ierr;
    140:   PetscBool            isascii;
    141:   PetscSubcomm         subcomm;
    142:   PetscViewer          subviewer;

    145:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
    146:   /* ASCII viewer */
    147:   if (isascii) {
    148:     PetscMPIInt   color,rank,size;
    149:     PetscInt64    loc[7],gsum[6],gmax[6],gmin[6],totbenign;
    150:     PetscScalar   interface_size;
    151:     PetscReal     ratio1=0.,ratio2=0.;
    152:     Vec           counter;

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

    207:     /* compute interface size */
    208:     VecSet(pcis->vec1_B,1.0);
    209:     MatCreateVecs(pc->pmat,&counter,NULL);
    210:     VecSet(counter,0.0);
    211:     VecScatterBegin(pcis->global_to_B,pcis->vec1_B,counter,INSERT_VALUES,SCATTER_REVERSE);
    212:     VecScatterEnd(pcis->global_to_B,pcis->vec1_B,counter,INSERT_VALUES,SCATTER_REVERSE);
    213:     VecSum(counter,&interface_size);
    214:     VecDestroy(&counter);

    216:     /* compute some statistics on the domain decomposition */
    217:     gsum[0] = 1;
    218:     gsum[1] = gsum[2] = gsum[3] = gsum[4] = gsum[5] = 0;
    219:     loc[0]  = !!pcis->n;
    220:     loc[1]  = pcis->n - pcis->n_B;
    221:     loc[2]  = pcis->n_B;
    222:     loc[3]  = pcbddc->local_primal_size;
    223:     loc[4]  = pcis->n;
    224:     loc[5]  = pcbddc->n_local_subs > 0 ? pcbddc->n_local_subs : (pcis->n ? 1 : 0);
    225:     loc[6]  = pcbddc->benign_n;
    226:     MPI_Reduce(loc,gsum,6,MPIU_INT64,MPI_SUM,0,PetscObjectComm((PetscObject)pc));
    227:     if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = loc[5] = -1;
    228:     MPI_Reduce(loc,gmax,6,MPIU_INT64,MPI_MAX,0,PetscObjectComm((PetscObject)pc));
    229:     if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = loc[5] = PETSC_MAX_INT;
    230:     MPI_Reduce(loc,gmin,6,MPIU_INT64,MPI_MIN,0,PetscObjectComm((PetscObject)pc));
    231:     MPI_Reduce(&loc[6],&totbenign,1,MPIU_INT64,MPI_SUM,0,PetscObjectComm((PetscObject)pc));
    232:     if (pcbddc->coarse_size) {
    233:       ratio1 = pc->pmat->rmap->N/(1.*pcbddc->coarse_size);
    234:       ratio2 = PetscRealPart(interface_size)/pcbddc->coarse_size;
    235:     }
    236:     PetscViewerASCIIPrintf(viewer,"********************************** STATISTICS AT LEVEL %d **********************************\n",pcbddc->current_level);
    237:     PetscViewerASCIIPrintf(viewer,"  Global dofs sizes: all %D interface %D coarse %D\n",pc->pmat->rmap->N,(PetscInt)PetscRealPart(interface_size),pcbddc->coarse_size);
    238:     PetscViewerASCIIPrintf(viewer,"  Coarsening ratios: all/coarse %D interface/coarse %D\n",(PetscInt)ratio1,(PetscInt)ratio2);
    239:     PetscViewerASCIIPrintf(viewer,"  Active processes : %D\n",(PetscInt)gsum[0]);
    240:     PetscViewerASCIIPrintf(viewer,"  Total subdomains : %D\n",(PetscInt)gsum[5]);
    241:     if (pcbddc->benign_have_null) {
    242:       PetscViewerASCIIPrintf(viewer,"  Benign subs      : %D\n",(PetscInt)totbenign);
    243:     }
    244:     PetscViewerASCIIPrintf(viewer,"  Dofs type        :\tMIN\tMAX\tMEAN\n");
    245:     PetscViewerASCIIPrintf(viewer,"  Interior  dofs   :\t%D\t%D\t%D\n",(PetscInt)gmin[1],(PetscInt)gmax[1],(PetscInt)(gsum[1]/gsum[0]));
    246:     PetscViewerASCIIPrintf(viewer,"  Interface dofs   :\t%D\t%D\t%D\n",(PetscInt)gmin[2],(PetscInt)gmax[2],(PetscInt)(gsum[2]/gsum[0]));
    247:     PetscViewerASCIIPrintf(viewer,"  Primal    dofs   :\t%D\t%D\t%D\n",(PetscInt)gmin[3],(PetscInt)gmax[3],(PetscInt)(gsum[3]/gsum[0]));
    248:     PetscViewerASCIIPrintf(viewer,"  Local     dofs   :\t%D\t%D\t%D\n",(PetscInt)gmin[4],(PetscInt)gmax[4],(PetscInt)(gsum[4]/gsum[0]));
    249:     PetscViewerASCIIPrintf(viewer,"  Local     subs   :\t%D\t%D\n"    ,(PetscInt)gmin[5],(PetscInt)gmax[5]);
    250:     PetscViewerFlush(viewer);

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

    254:     /* local solvers */
    255:     PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)pcbddc->ksp_D),&subviewer);
    256:     if (!rank) {
    257:       PetscViewerASCIIPrintf(subviewer,"--- Interior solver (rank 0)\n");
    258:       PetscViewerASCIIPushTab(subviewer);
    259:       KSPView(pcbddc->ksp_D,subviewer);
    260:       PetscViewerASCIIPopTab(subviewer);
    261:       PetscViewerASCIIPrintf(subviewer,"--- Correction solver (rank 0)\n");
    262:       PetscViewerASCIIPushTab(subviewer);
    263:       KSPView(pcbddc->ksp_R,subviewer);
    264:       PetscViewerASCIIPopTab(subviewer);
    265:       PetscViewerFlush(subviewer);
    266:     }
    267:     PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)pcbddc->ksp_D),&subviewer);
    268:     PetscViewerFlush(viewer);

    270:     /* the coarse problem can be handled by a different communicator */
    271:     if (pcbddc->coarse_ksp) color = 1;
    272:     else color = 0;
    273:     MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
    274:     PetscSubcommCreate(PetscObjectComm((PetscObject)pc),&subcomm);
    275:     PetscSubcommSetNumber(subcomm,PetscMin(size,2));
    276:     PetscSubcommSetTypeGeneral(subcomm,color,rank);
    277:     PetscViewerGetSubViewer(viewer,PetscSubcommChild(subcomm),&subviewer);
    278:     if (color == 1) {
    279:       PetscViewerASCIIPrintf(subviewer,"--- Coarse solver\n");
    280:       PetscViewerASCIIPushTab(subviewer);
    281:       KSPView(pcbddc->coarse_ksp,subviewer);
    282:       PetscViewerASCIIPopTab(subviewer);
    283:       PetscViewerFlush(subviewer);
    284:     }
    285:     PetscViewerRestoreSubViewer(viewer,PetscSubcommChild(subcomm),&subviewer);
    286:     PetscSubcommDestroy(&subcomm);
    287:     PetscViewerFlush(viewer);
    288:   }
    289:   return(0);
    290: }

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

    298:   PetscObjectReference((PetscObject)G);
    299:   MatDestroy(&pcbddc->discretegradient);
    300:   pcbddc->discretegradient = G;
    301:   pcbddc->nedorder         = order > 0 ? order : -order;
    302:   pcbddc->nedfield         = field;
    303:   pcbddc->nedglobal        = global;
    304:   pcbddc->conforming       = conforming;
    305:   return(0);
    306: }

    308: /*@
    309:  PCBDDCSetDiscreteGradient - Sets the discrete gradient

    311:    Collective on PC

    313:    Input Parameters:
    314: +  pc         - the preconditioning context
    315: .  G          - the discrete gradient matrix (should be in AIJ format)
    316: .  order      - the order of the Nedelec space (1 for the lowest order)
    317: .  field      - the field id of the Nedelec dofs (not used if the fields have not been specified)
    318: .  global     - the type of global ordering for the rows of G
    319: -  conforming - whether the mesh is conforming or not

    321:    Level: advanced

    323:    Notes:
    324:     The discrete gradient matrix G is used to analyze the subdomain edges, and it should not contain any zero entry.
    325:           For variable order spaces, the order should be set to zero.
    326:           If global is true, the rows of G should be given in global ordering for the whole dofs;
    327:           if false, the ordering should be global for the Nedelec field.
    328:           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
    329:           and geid the one for the Nedelec field.

    331: .seealso: PCBDDC,PCBDDCSetDofsSplitting(),PCBDDCSetDofsSplittingLocal()
    332: @*/
    333: PetscErrorCode PCBDDCSetDiscreteGradient(PC pc, Mat G, PetscInt order, PetscInt field, PetscBool global, PetscBool conforming)
    334: {

    345:   PetscTryMethod(pc,"PCBDDCSetDiscreteGradient_C",(PC,Mat,PetscInt,PetscInt,PetscBool,PetscBool),(pc,G,order,field,global,conforming));
    346:   return(0);
    347: }

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

    355:   PetscObjectReference((PetscObject)divudotp);
    356:   MatDestroy(&pcbddc->divudotp);
    357:   pcbddc->divudotp = divudotp;
    358:   pcbddc->divudotp_trans = trans;
    359:   pcbddc->compute_nonetflux = PETSC_TRUE;
    360:   if (vl2l) {
    361:     PetscObjectReference((PetscObject)vl2l);
    362:     ISDestroy(&pcbddc->divudotp_vl2l);
    363:     pcbddc->divudotp_vl2l = vl2l;
    364:   }
    365:   return(0);
    366: }

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

    371:    Collective on PC

    373:    Input Parameters:
    374: +  pc - the preconditioning context
    375: .  divudotp - the matrix (must be of type MATIS)
    376: .  trans - if trans if false (resp. true), then pressures are in the test (trial) space and velocities are in the trial (test) space.
    377: -  vl2l - optional index set describing the local (wrt the local matrix in divudotp) to local (wrt the local matrix in the preconditioning matrix) map for the velocities

    379:    Level: advanced

    381:    Notes:
    382:     This auxiliary matrix is used to compute quadrature weights representing the net-flux across subdomain boundaries
    383:           If vl2l is NULL, the local ordering for velocities in divudotp should match that of the preconditioning matrix

    385: .seealso: PCBDDC
    386: @*/
    387: PetscErrorCode PCBDDCSetDivergenceMat(PC pc, Mat divudotp, PetscBool trans, IS vl2l)
    388: {
    389:   PetscBool      ismatis;

    398:   PetscObjectTypeCompare((PetscObject)divudotp,MATIS,&ismatis);
    399:   if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Divergence matrix needs to be of type MATIS");
    400:   PetscTryMethod(pc,"PCBDDCSetDivergenceMat_C",(PC,Mat,PetscBool,IS),(pc,divudotp,trans,vl2l));
    401:   return(0);
    402: }

    404: static PetscErrorCode PCBDDCSetChangeOfBasisMat_BDDC(PC pc, Mat change, PetscBool interior)
    405: {
    406:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;

    410:   PetscObjectReference((PetscObject)change);
    411:   MatDestroy(&pcbddc->user_ChangeOfBasisMatrix);
    412:   pcbddc->user_ChangeOfBasisMatrix = change;
    413:   pcbddc->change_interior = interior;
    414:   return(0);
    415: }
    416: /*@
    417:  PCBDDCSetChangeOfBasisMat - Set user defined change of basis for dofs

    419:    Collective on PC

    421:    Input Parameters:
    422: +  pc - the preconditioning context
    423: .  change - the change of basis matrix
    424: -  interior - whether or not the change of basis modifies interior dofs

    426:    Level: intermediate

    428:    Notes:

    430: .seealso: PCBDDC
    431: @*/
    432: PetscErrorCode PCBDDCSetChangeOfBasisMat(PC pc, Mat change, PetscBool interior)
    433: {

    440:   if (pc->mat) {
    441:     PetscInt rows_c,cols_c,rows,cols;
    442:     MatGetSize(pc->mat,&rows,&cols);
    443:     MatGetSize(change,&rows_c,&cols_c);
    444:     if (rows_c != rows) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of rows for change of basis matrix! %D != %D",rows_c,rows);
    445:     if (cols_c != cols) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of columns for change of basis matrix! %D != %D",cols_c,cols);
    446:     MatGetLocalSize(pc->mat,&rows,&cols);
    447:     MatGetLocalSize(change,&rows_c,&cols_c);
    448:     if (rows_c != rows) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of local rows for change of basis matrix! %D != %D",rows_c,rows);
    449:     if (cols_c != cols) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of local columns for change of basis matrix! %D != %D",cols_c,cols);
    450:   }
    451:   PetscTryMethod(pc,"PCBDDCSetChangeOfBasisMat_C",(PC,Mat,PetscBool),(pc,change,interior));
    452:   return(0);
    453: }

    455: static PetscErrorCode PCBDDCSetPrimalVerticesIS_BDDC(PC pc, IS PrimalVertices)
    456: {
    457:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
    458:   PetscBool      isequal = PETSC_FALSE;

    462:   PetscObjectReference((PetscObject)PrimalVertices);
    463:   if (pcbddc->user_primal_vertices) {
    464:     ISEqual(PrimalVertices,pcbddc->user_primal_vertices,&isequal);
    465:   }
    466:   ISDestroy(&pcbddc->user_primal_vertices);
    467:   ISDestroy(&pcbddc->user_primal_vertices_local);
    468:   pcbddc->user_primal_vertices = PrimalVertices;
    469:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
    470:   return(0);
    471: }

    473: /*@
    474:  PCBDDCSetPrimalVerticesIS - Set additional user defined primal vertices in PCBDDC

    476:    Collective

    478:    Input Parameters:
    479: +  pc - the preconditioning context
    480: -  PrimalVertices - index set of primal vertices in global numbering (can be empty)

    482:    Level: intermediate

    484:    Notes:
    485:      Any process can list any global node

    487: .seealso: PCBDDC, PCBDDCGetPrimalVerticesIS(), PCBDDCSetPrimalVerticesLocalIS(), PCBDDCGetPrimalVerticesLocalIS()
    488: @*/
    489: PetscErrorCode PCBDDCSetPrimalVerticesIS(PC pc, IS PrimalVertices)
    490: {

    497:   PetscTryMethod(pc,"PCBDDCSetPrimalVerticesIS_C",(PC,IS),(pc,PrimalVertices));
    498:   return(0);
    499: }

    501: static PetscErrorCode PCBDDCGetPrimalVerticesIS_BDDC(PC pc, IS *is)
    502: {
    503:   PC_BDDC *pcbddc = (PC_BDDC*)pc->data;

    506:   *is = pcbddc->user_primal_vertices;
    507:   return(0);
    508: }

    510: /*@
    511:  PCBDDCGetPrimalVerticesIS - Get user defined primal vertices set with PCBDDCSetPrimalVerticesIS()

    513:    Collective

    515:    Input Parameters:
    516: .  pc - the preconditioning context

    518:    Output Parameters:
    519: .  is - index set of primal vertices in global numbering (NULL if not set)

    521:    Level: intermediate

    523:    Notes:

    525: .seealso: PCBDDC, PCBDDCSetPrimalVerticesIS(), PCBDDCSetPrimalVerticesLocalIS(), PCBDDCGetPrimalVerticesLocalIS()
    526: @*/
    527: PetscErrorCode PCBDDCGetPrimalVerticesIS(PC pc, IS *is)
    528: {

    534:   PetscUseMethod(pc,"PCBDDCGetPrimalVerticesIS_C",(PC,IS*),(pc,is));
    535:   return(0);
    536: }

    538: static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices)
    539: {
    540:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
    541:   PetscBool      isequal = PETSC_FALSE;

    545:   PetscObjectReference((PetscObject)PrimalVertices);
    546:   if (pcbddc->user_primal_vertices_local) {
    547:     ISEqual(PrimalVertices,pcbddc->user_primal_vertices_local,&isequal);
    548:   }
    549:   ISDestroy(&pcbddc->user_primal_vertices);
    550:   ISDestroy(&pcbddc->user_primal_vertices_local);
    551:   pcbddc->user_primal_vertices_local = PrimalVertices;
    552:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
    553:   return(0);
    554: }

    556: /*@
    557:  PCBDDCSetPrimalVerticesLocalIS - Set additional user defined primal vertices in PCBDDC

    559:    Collective

    561:    Input Parameters:
    562: +  pc - the preconditioning context
    563: -  PrimalVertices - index set of primal vertices in local numbering (can be empty)

    565:    Level: intermediate

    567:    Notes:

    569: .seealso: PCBDDC, PCBDDCSetPrimalVerticesIS(), PCBDDCGetPrimalVerticesIS(), PCBDDCGetPrimalVerticesLocalIS()
    570: @*/
    571: PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices)
    572: {

    579:   PetscTryMethod(pc,"PCBDDCSetPrimalVerticesLocalIS_C",(PC,IS),(pc,PrimalVertices));
    580:   return(0);
    581: }

    583: static PetscErrorCode PCBDDCGetPrimalVerticesLocalIS_BDDC(PC pc, IS *is)
    584: {
    585:   PC_BDDC *pcbddc = (PC_BDDC*)pc->data;

    588:   *is = pcbddc->user_primal_vertices_local;
    589:   return(0);
    590: }

    592: /*@
    593:  PCBDDCGetPrimalVerticesLocalIS - Get user defined primal vertices set with PCBDDCSetPrimalVerticesLocalIS()

    595:    Collective

    597:    Input Parameters:
    598: .  pc - the preconditioning context

    600:    Output Parameters:
    601: .  is - index set of primal vertices in local numbering (NULL if not set)

    603:    Level: intermediate

    605:    Notes:

    607: .seealso: PCBDDC, PCBDDCSetPrimalVerticesIS(), PCBDDCGetPrimalVerticesIS(), PCBDDCSetPrimalVerticesLocalIS()
    608: @*/
    609: PetscErrorCode PCBDDCGetPrimalVerticesLocalIS(PC pc, IS *is)
    610: {

    616:   PetscUseMethod(pc,"PCBDDCGetPrimalVerticesLocalIS_C",(PC,IS*),(pc,is));
    617:   return(0);
    618: }

    620: static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc,PetscInt k)
    621: {
    622:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

    625:   pcbddc->coarsening_ratio = k;
    626:   return(0);
    627: }

    629: /*@
    630:  PCBDDCSetCoarseningRatio - Set coarsening ratio used in multilevel

    632:    Logically collective on PC

    634:    Input Parameters:
    635: +  pc - the preconditioning context
    636: -  k - coarsening ratio (H/h at the coarser level)

    638:    Options Database Keys:
    639: .    -pc_bddc_coarsening_ratio

    641:    Level: intermediate

    643:    Notes:
    644:      Approximatively k subdomains at the finer level will be aggregated into a single subdomain at the coarser level

    646: .seealso: PCBDDC, PCBDDCSetLevels()
    647: @*/
    648: PetscErrorCode PCBDDCSetCoarseningRatio(PC pc,PetscInt k)
    649: {

    655:   PetscTryMethod(pc,"PCBDDCSetCoarseningRatio_C",(PC,PetscInt),(pc,k));
    656:   return(0);
    657: }

    659: /* The following functions (PCBDDCSetUseExactDirichlet PCBDDCSetLevel) are not public */
    660: static PetscErrorCode PCBDDCSetUseExactDirichlet_BDDC(PC pc,PetscBool flg)
    661: {
    662:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

    665:   pcbddc->use_exact_dirichlet_trick = flg;
    666:   return(0);
    667: }

    669: PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool flg)
    670: {

    676:   PetscTryMethod(pc,"PCBDDCSetUseExactDirichlet_C",(PC,PetscBool),(pc,flg));
    677:   return(0);
    678: }

    680: static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc,PetscInt level)
    681: {
    682:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

    685:   pcbddc->current_level = level;
    686:   return(0);
    687: }

    689: PetscErrorCode PCBDDCSetLevel(PC pc,PetscInt level)
    690: {

    696:   PetscTryMethod(pc,"PCBDDCSetLevel_C",(PC,PetscInt),(pc,level));
    697:   return(0);
    698: }

    700: static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc,PetscInt levels)
    701: {
    702:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

    705:   if (levels > PETSC_PCBDDC_MAXLEVELS-1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Maximum number of additional levels for BDDC is %d",PETSC_PCBDDC_MAXLEVELS-1);
    706:   pcbddc->max_levels = levels;
    707:   return(0);
    708: }

    710: /*@
    711:  PCBDDCSetLevels - Sets the maximum number of additional levels allowed for multilevel BDDC

    713:    Logically collective on PC

    715:    Input Parameters:
    716: +  pc - the preconditioning context
    717: -  levels - the maximum number of levels

    719:    Options Database Keys:
    720: .    -pc_bddc_levels

    722:    Level: intermediate

    724:    Notes:
    725:      The default value is 0, that gives the classical two-levels BDDC

    727: .seealso: PCBDDC, PCBDDCSetCoarseningRatio()
    728: @*/
    729: PetscErrorCode PCBDDCSetLevels(PC pc,PetscInt levels)
    730: {

    736:   PetscTryMethod(pc,"PCBDDCSetLevels_C",(PC,PetscInt),(pc,levels));
    737:   return(0);
    738: }

    740: static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc,IS DirichletBoundaries)
    741: {
    742:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
    743:   PetscBool      isequal = PETSC_FALSE;

    747:   PetscObjectReference((PetscObject)DirichletBoundaries);
    748:   if (pcbddc->DirichletBoundaries) {
    749:     ISEqual(DirichletBoundaries,pcbddc->DirichletBoundaries,&isequal);
    750:   }
    751:   /* last user setting takes precendence -> destroy any other customization */
    752:   ISDestroy(&pcbddc->DirichletBoundariesLocal);
    753:   ISDestroy(&pcbddc->DirichletBoundaries);
    754:   pcbddc->DirichletBoundaries = DirichletBoundaries;
    755:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
    756:   return(0);
    757: }

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

    762:    Collective

    764:    Input Parameters:
    765: +  pc - the preconditioning context
    766: -  DirichletBoundaries - parallel IS defining the Dirichlet boundaries

    768:    Level: intermediate

    770:    Notes:
    771:      Provide the information if you used MatZeroRows/Columns routines. Any process can list any global node

    773: .seealso: PCBDDC, PCBDDCSetDirichletBoundariesLocal(), MatZeroRows(), MatZeroRowsColumns()
    774: @*/
    775: PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc,IS DirichletBoundaries)
    776: {

    783:   PetscTryMethod(pc,"PCBDDCSetDirichletBoundaries_C",(PC,IS),(pc,DirichletBoundaries));
    784:   return(0);
    785: }

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

    794:   PetscObjectReference((PetscObject)DirichletBoundaries);
    795:   if (pcbddc->DirichletBoundariesLocal) {
    796:     ISEqual(DirichletBoundaries,pcbddc->DirichletBoundariesLocal,&isequal);
    797:   }
    798:   /* last user setting takes precendence -> destroy any other customization */
    799:   ISDestroy(&pcbddc->DirichletBoundariesLocal);
    800:   ISDestroy(&pcbddc->DirichletBoundaries);
    801:   pcbddc->DirichletBoundariesLocal = DirichletBoundaries;
    802:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
    803:   return(0);
    804: }

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

    809:    Collective

    811:    Input Parameters:
    812: +  pc - the preconditioning context
    813: -  DirichletBoundaries - parallel IS defining the Dirichlet boundaries (in local ordering)

    815:    Level: intermediate

    817:    Notes:

    819: .seealso: PCBDDC, PCBDDCSetDirichletBoundaries(), MatZeroRows(), MatZeroRowsColumns()
    820: @*/
    821: PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc,IS DirichletBoundaries)
    822: {

    829:   PetscTryMethod(pc,"PCBDDCSetDirichletBoundariesLocal_C",(PC,IS),(pc,DirichletBoundaries));
    830:   return(0);
    831: }

    833: static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc,IS NeumannBoundaries)
    834: {
    835:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
    836:   PetscBool      isequal = PETSC_FALSE;

    840:   PetscObjectReference((PetscObject)NeumannBoundaries);
    841:   if (pcbddc->NeumannBoundaries) {
    842:     ISEqual(NeumannBoundaries,pcbddc->NeumannBoundaries,&isequal);
    843:   }
    844:   /* last user setting takes precendence -> destroy any other customization */
    845:   ISDestroy(&pcbddc->NeumannBoundariesLocal);
    846:   ISDestroy(&pcbddc->NeumannBoundaries);
    847:   pcbddc->NeumannBoundaries = NeumannBoundaries;
    848:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
    849:   return(0);
    850: }

    852: /*@
    853:  PCBDDCSetNeumannBoundaries - Set IS defining Neumann boundaries for the global problem.

    855:    Collective

    857:    Input Parameters:
    858: +  pc - the preconditioning context
    859: -  NeumannBoundaries - parallel IS defining the Neumann boundaries

    861:    Level: intermediate

    863:    Notes:
    864:      Any process can list any global node

    866: .seealso: PCBDDC, PCBDDCSetNeumannBoundariesLocal()
    867: @*/
    868: PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc,IS NeumannBoundaries)
    869: {

    876:   PetscTryMethod(pc,"PCBDDCSetNeumannBoundaries_C",(PC,IS),(pc,NeumannBoundaries));
    877:   return(0);
    878: }

    880: static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc,IS NeumannBoundaries)
    881: {
    882:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
    883:   PetscBool      isequal = PETSC_FALSE;

    887:   PetscObjectReference((PetscObject)NeumannBoundaries);
    888:   if (pcbddc->NeumannBoundariesLocal) {
    889:     ISEqual(NeumannBoundaries,pcbddc->NeumannBoundariesLocal,&isequal);
    890:   }
    891:   /* last user setting takes precendence -> destroy any other customization */
    892:   ISDestroy(&pcbddc->NeumannBoundariesLocal);
    893:   ISDestroy(&pcbddc->NeumannBoundaries);
    894:   pcbddc->NeumannBoundariesLocal = NeumannBoundaries;
    895:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
    896:   return(0);
    897: }

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

    902:    Collective

    904:    Input Parameters:
    905: +  pc - the preconditioning context
    906: -  NeumannBoundaries - parallel IS defining the subdomain part of Neumann boundaries (in local ordering)

    908:    Level: intermediate

    910:    Notes:

    912: .seealso: PCBDDC, PCBDDCSetNeumannBoundaries()
    913: @*/
    914: PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc,IS NeumannBoundaries)
    915: {

    922:   PetscTryMethod(pc,"PCBDDCSetNeumannBoundariesLocal_C",(PC,IS),(pc,NeumannBoundaries));
    923:   return(0);
    924: }

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

    931:   *DirichletBoundaries = pcbddc->DirichletBoundaries;
    932:   return(0);
    933: }

    935: /*@
    936:  PCBDDCGetDirichletBoundaries - Get parallel IS for Dirichlet boundaries

    938:    Collective

    940:    Input Parameters:
    941: .  pc - the preconditioning context

    943:    Output Parameters:
    944: .  DirichletBoundaries - index set defining the Dirichlet boundaries

    946:    Level: intermediate

    948:    Notes:
    949:      The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetDirichletBoundaries

    951: .seealso: PCBDDC
    952: @*/
    953: PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc,IS *DirichletBoundaries)
    954: {

    959:   PetscUseMethod(pc,"PCBDDCGetDirichletBoundaries_C",(PC,IS*),(pc,DirichletBoundaries));
    960:   return(0);
    961: }

    963: static PetscErrorCode PCBDDCGetDirichletBoundariesLocal_BDDC(PC pc,IS *DirichletBoundaries)
    964: {
    965:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

    968:   *DirichletBoundaries = pcbddc->DirichletBoundariesLocal;
    969:   return(0);
    970: }

    972: /*@
    973:  PCBDDCGetDirichletBoundariesLocal - Get parallel IS for Dirichlet boundaries (in local ordering)

    975:    Collective

    977:    Input Parameters:
    978: .  pc - the preconditioning context

    980:    Output Parameters:
    981: .  DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries

    983:    Level: intermediate

    985:    Notes:
    986:      The IS returned could be the same passed in earlier by the user (if provided with PCBDDCSetDirichletBoundariesLocal) or a global-to-local map of the global IS (if provided with PCBDDCSetDirichletBoundaries).
    987:           In the latter case, the IS will be available after PCSetUp.

    989: .seealso: PCBDDC
    990: @*/
    991: PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc,IS *DirichletBoundaries)
    992: {

    997:   PetscUseMethod(pc,"PCBDDCGetDirichletBoundariesLocal_C",(PC,IS*),(pc,DirichletBoundaries));
    998:   return(0);
    999: }

    1001: static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc,IS *NeumannBoundaries)
    1002: {
    1003:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

    1006:   *NeumannBoundaries = pcbddc->NeumannBoundaries;
    1007:   return(0);
    1008: }

    1010: /*@
    1011:  PCBDDCGetNeumannBoundaries - Get parallel IS for Neumann boundaries

    1013:    Collective

    1015:    Input Parameters:
    1016: .  pc - the preconditioning context

    1018:    Output Parameters:
    1019: .  NeumannBoundaries - index set defining the Neumann boundaries

    1021:    Level: intermediate

    1023:    Notes:
    1024:      The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetNeumannBoundaries

    1026: .seealso: PCBDDC
    1027: @*/
    1028: PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc,IS *NeumannBoundaries)
    1029: {

    1034:   PetscUseMethod(pc,"PCBDDCGetNeumannBoundaries_C",(PC,IS*),(pc,NeumannBoundaries));
    1035:   return(0);
    1036: }

    1038: static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc,IS *NeumannBoundaries)
    1039: {
    1040:   PC_BDDC  *pcbddc = (PC_BDDC*)pc->data;

    1043:   *NeumannBoundaries = pcbddc->NeumannBoundariesLocal;
    1044:   return(0);
    1045: }

    1047: /*@
    1048:  PCBDDCGetNeumannBoundariesLocal - Get parallel IS for Neumann boundaries (in local ordering)

    1050:    Collective

    1052:    Input Parameters:
    1053: .  pc - the preconditioning context

    1055:    Output Parameters:
    1056: .  NeumannBoundaries - index set defining the subdomain part of Neumann boundaries

    1058:    Level: intermediate

    1060:    Notes:
    1061:      The IS returned could be the same passed in earlier by the user (if provided with PCBDDCSetNeumannBoundariesLocal) or a global-to-local map of the global IS (if provided with PCBDDCSetNeumannBoundaries).
    1062:           In the latter case, the IS will be available after PCSetUp.

    1064: .seealso: PCBDDC
    1065: @*/
    1066: PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc,IS *NeumannBoundaries)
    1067: {

    1072:   PetscUseMethod(pc,"PCBDDCGetNeumannBoundariesLocal_C",(PC,IS*),(pc,NeumannBoundaries));
    1073:   return(0);
    1074: }

    1076: static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode)
    1077: {
    1078:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
    1079:   PCBDDCGraph    mat_graph = pcbddc->mat_graph;
    1080:   PetscBool      same_data = PETSC_FALSE;

    1084:   if (!nvtxs) {
    1085:     if (copymode == PETSC_OWN_POINTER) {
    1086:       PetscFree(xadj);
    1087:       PetscFree(adjncy);
    1088:     }
    1089:     PCBDDCGraphResetCSR(mat_graph);
    1090:     return(0);
    1091:   }
    1092:   if (mat_graph->nvtxs == nvtxs && mat_graph->freecsr) { /* we own the data */
    1093:     if (mat_graph->xadj == xadj && mat_graph->adjncy == adjncy) same_data = PETSC_TRUE;
    1094:     if (!same_data && mat_graph->xadj[nvtxs] == xadj[nvtxs]) {
    1095:       PetscArraycmp(xadj,mat_graph->xadj,nvtxs+1,&same_data);
    1096:       if (same_data) {
    1097:         PetscArraycmp(adjncy,mat_graph->adjncy,xadj[nvtxs],&same_data);
    1098:       }
    1099:     }
    1100:   }
    1101:   if (!same_data) {
    1102:     /* free old CSR */
    1103:     PCBDDCGraphResetCSR(mat_graph);
    1104:     /* get CSR into graph structure */
    1105:     if (copymode == PETSC_COPY_VALUES) {
    1106:       PetscMalloc1(nvtxs+1,&mat_graph->xadj);
    1107:       PetscMalloc1(xadj[nvtxs],&mat_graph->adjncy);
    1108:       PetscArraycpy(mat_graph->xadj,xadj,nvtxs+1);
    1109:       PetscArraycpy(mat_graph->adjncy,adjncy,xadj[nvtxs]);
    1110:       mat_graph->freecsr = PETSC_TRUE;
    1111:     } else if (copymode == PETSC_OWN_POINTER) {
    1112:       mat_graph->xadj    = (PetscInt*)xadj;
    1113:       mat_graph->adjncy  = (PetscInt*)adjncy;
    1114:       mat_graph->freecsr = PETSC_TRUE;
    1115:     } else if (copymode == PETSC_USE_POINTER) {
    1116:       mat_graph->xadj    = (PetscInt*)xadj;
    1117:       mat_graph->adjncy  = (PetscInt*)adjncy;
    1118:       mat_graph->freecsr = PETSC_FALSE;
    1119:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported copy mode %D",copymode);
    1120:     mat_graph->nvtxs_csr = nvtxs;
    1121:     pcbddc->recompute_topography = PETSC_TRUE;
    1122:   }
    1123:   return(0);
    1124: }

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

    1129:    Not collective

    1131:    Input Parameters:
    1132: +  pc - the preconditioning context.
    1133: .  nvtxs - number of local vertices of the graph (i.e., the number of local dofs).
    1134: .  xadj, adjncy - the connectivity of the dofs in CSR format.
    1135: -  copymode - supported modes are PETSC_COPY_VALUES, PETSC_USE_POINTER or PETSC_OWN_POINTER.

    1137:    Level: intermediate

    1139:    Notes:
    1140:     A dof is considered connected with all local dofs if xadj[dof+1]-xadj[dof] == 1 and adjncy[xadj[dof]] is negative.

    1142: .seealso: PCBDDC,PetscCopyMode
    1143: @*/
    1144: PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc,PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode)
    1145: {
    1146:   void (*f)(void) = NULL;

    1151:   if (nvtxs) {
    1154:   }
    1155:   PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));
    1156:   /* free arrays if PCBDDC is not the PC type */
    1157:   PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);
    1158:   if (!f && copymode == PETSC_OWN_POINTER) {
    1159:     PetscFree(xadj);
    1160:     PetscFree(adjncy);
    1161:   }
    1162:   return(0);
    1163: }

    1165: static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc,PetscInt n_is, IS ISForDofs[])
    1166: {
    1167:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
    1168:   PetscInt       i;
    1169:   PetscBool      isequal = PETSC_FALSE;

    1173:   if (pcbddc->n_ISForDofsLocal == n_is) {
    1174:     for (i=0;i<n_is;i++) {
    1175:       PetscBool isequalt;
    1176:       ISEqual(ISForDofs[i],pcbddc->ISForDofsLocal[i],&isequalt);
    1177:       if (!isequalt) break;
    1178:     }
    1179:     if (i == n_is) isequal = PETSC_TRUE;
    1180:   }
    1181:   for (i=0;i<n_is;i++) {
    1182:     PetscObjectReference((PetscObject)ISForDofs[i]);
    1183:   }
    1184:   /* Destroy ISes if they were already set */
    1185:   for (i=0;i<pcbddc->n_ISForDofsLocal;i++) {
    1186:     ISDestroy(&pcbddc->ISForDofsLocal[i]);
    1187:   }
    1188:   PetscFree(pcbddc->ISForDofsLocal);
    1189:   /* last user setting takes precendence -> destroy any other customization */
    1190:   for (i=0;i<pcbddc->n_ISForDofs;i++) {
    1191:     ISDestroy(&pcbddc->ISForDofs[i]);
    1192:   }
    1193:   PetscFree(pcbddc->ISForDofs);
    1194:   pcbddc->n_ISForDofs = 0;
    1195:   /* allocate space then set */
    1196:   if (n_is) {
    1197:     PetscMalloc1(n_is,&pcbddc->ISForDofsLocal);
    1198:   }
    1199:   for (i=0;i<n_is;i++) {
    1200:     pcbddc->ISForDofsLocal[i] = ISForDofs[i];
    1201:   }
    1202:   pcbddc->n_ISForDofsLocal = n_is;
    1203:   if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
    1204:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
    1205:   return(0);
    1206: }

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

    1211:    Collective

    1213:    Input Parameters:
    1214: +  pc - the preconditioning context
    1215: .  n_is - number of index sets defining the fields
    1216: -  ISForDofs - array of IS describing the fields in local ordering

    1218:    Level: intermediate

    1220:    Notes:
    1221:      n_is should be the same among processes. Not all nodes need to be listed: unlisted nodes will belong to the complement field.

    1223: .seealso: PCBDDC
    1224: @*/
    1225: PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc,PetscInt n_is, IS ISForDofs[])
    1226: {
    1227:   PetscInt       i;

    1233:   for (i=0;i<n_is;i++) {
    1236:   }
    1237:   PetscTryMethod(pc,"PCBDDCSetDofsSplittingLocal_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));
    1238:   return(0);
    1239: }

    1241: static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc,PetscInt n_is, IS ISForDofs[])
    1242: {
    1243:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
    1244:   PetscInt       i;
    1245:   PetscBool      isequal = PETSC_FALSE;

    1249:   if (pcbddc->n_ISForDofs == n_is) {
    1250:     for (i=0;i<n_is;i++) {
    1251:       PetscBool isequalt;
    1252:       ISEqual(ISForDofs[i],pcbddc->ISForDofs[i],&isequalt);
    1253:       if (!isequalt) break;
    1254:     }
    1255:     if (i == n_is) isequal = PETSC_TRUE;
    1256:   }
    1257:   for (i=0;i<n_is;i++) {
    1258:     PetscObjectReference((PetscObject)ISForDofs[i]);
    1259:   }
    1260:   /* Destroy ISes if they were already set */
    1261:   for (i=0;i<pcbddc->n_ISForDofs;i++) {
    1262:     ISDestroy(&pcbddc->ISForDofs[i]);
    1263:   }
    1264:   PetscFree(pcbddc->ISForDofs);
    1265:   /* last user setting takes precendence -> destroy any other customization */
    1266:   for (i=0;i<pcbddc->n_ISForDofsLocal;i++) {
    1267:     ISDestroy(&pcbddc->ISForDofsLocal[i]);
    1268:   }
    1269:   PetscFree(pcbddc->ISForDofsLocal);
    1270:   pcbddc->n_ISForDofsLocal = 0;
    1271:   /* allocate space then set */
    1272:   if (n_is) {
    1273:     PetscMalloc1(n_is,&pcbddc->ISForDofs);
    1274:   }
    1275:   for (i=0;i<n_is;i++) {
    1276:     pcbddc->ISForDofs[i] = ISForDofs[i];
    1277:   }
    1278:   pcbddc->n_ISForDofs = n_is;
    1279:   if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE;
    1280:   if (!isequal) pcbddc->recompute_topography = PETSC_TRUE;
    1281:   return(0);
    1282: }

    1284: /*@
    1285:  PCBDDCSetDofsSplitting - Set index sets defining fields of the global matrix

    1287:    Collective

    1289:    Input Parameters:
    1290: +  pc - the preconditioning context
    1291: .  n_is - number of index sets defining the fields
    1292: -  ISForDofs - array of IS describing the fields in global ordering

    1294:    Level: intermediate

    1296:    Notes:
    1297:      Any process can list any global node. Not all nodes need to be listed: unlisted nodes will belong to the complement field.

    1299: .seealso: PCBDDC
    1300: @*/
    1301: PetscErrorCode PCBDDCSetDofsSplitting(PC pc,PetscInt n_is, IS ISForDofs[])
    1302: {
    1303:   PetscInt       i;

    1309:   for (i=0;i<n_is;i++) {
    1312:   }
    1313:   PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));
    1314:   return(0);
    1315: }

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

    1321:    Input Parameter:
    1322: +  pc - the preconditioner context

    1324:    Application Interface Routine: PCPreSolve()

    1326:    Notes:
    1327:      The interface routine PCPreSolve() is not usually called directly by
    1328:    the user, but instead is called by KSPSolve().
    1329: */
    1330: static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
    1331: {
    1333:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
    1334:   PC_IS          *pcis = (PC_IS*)(pc->data);
    1335:   Vec            used_vec;
    1336:   PetscBool      iscg = PETSC_FALSE, save_rhs = PETSC_TRUE, benign_correction_computed;

    1339:   /* if we are working with CG, one dirichlet solve can be avoided during Krylov iterations */
    1340:   if (ksp) {
    1341:     PetscBool isgroppcg, ispipecg, ispipelcg, ispipecgrr;

    1343:     PetscObjectTypeCompare((PetscObject)ksp,KSPCG,&iscg);
    1344:     PetscObjectTypeCompare((PetscObject)ksp,KSPGROPPCG,&isgroppcg);
    1345:     PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECG,&ispipecg);
    1346:     PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECG,&ispipelcg);
    1347:     PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECGRR,&ispipecgrr);
    1348:     iscg = (PetscBool)(iscg || isgroppcg || ispipecg || ispipelcg || ispipecgrr);
    1349:     if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || !iscg || pc->mat != pc->pmat) {
    1350:       PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);
    1351:     }
    1352:   }
    1353:   if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || pc->mat != pc->pmat) {
    1354:     PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);
    1355:   }

    1357:   /* Creates parallel work vectors used in presolve */
    1358:   if (!pcbddc->original_rhs) {
    1359:     VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);
    1360:   }
    1361:   if (!pcbddc->temp_solution) {
    1362:     VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);
    1363:   }

    1365:   pcbddc->temp_solution_used = PETSC_FALSE;
    1366:   if (x) {
    1367:     PetscObjectReference((PetscObject)x);
    1368:     used_vec = x;
    1369:   } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */
    1370:     PetscObjectReference((PetscObject)pcbddc->temp_solution);
    1371:     used_vec = pcbddc->temp_solution;
    1372:     VecSet(used_vec,0.0);
    1373:     pcbddc->temp_solution_used = PETSC_TRUE;
    1374:     VecCopy(rhs,pcbddc->original_rhs);
    1375:     save_rhs = PETSC_FALSE;
    1376:     pcbddc->eliminate_dirdofs = PETSC_TRUE;
    1377:   }

    1379:   /* hack into ksp data structure since PCPreSolve comes earlier than setting to zero the guess in src/ksp/ksp/interface/itfunc.c */
    1380:   if (ksp) {
    1381:     /* store the flag for the initial guess since it will be restored back during PCPostSolve_BDDC */
    1382:     KSPGetInitialGuessNonzero(ksp,&pcbddc->ksp_guess_nonzero);
    1383:     if (!pcbddc->ksp_guess_nonzero) {
    1384:       VecSet(used_vec,0.0);
    1385:     }
    1386:   }

    1388:   pcbddc->rhs_change = PETSC_FALSE;
    1389:   /* Take into account zeroed rows -> change rhs and store solution removed */
    1390:   if (rhs && pcbddc->eliminate_dirdofs) {
    1391:     IS dirIS = NULL;

    1393:     /* DirichletBoundariesLocal may not be consistent among neighbours; gets a dirichlet dofs IS from graph (may be cached) */
    1394:     PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph,&dirIS);
    1395:     if (dirIS) {
    1396:       Mat_IS            *matis = (Mat_IS*)pc->pmat->data;
    1397:       PetscInt          dirsize,i,*is_indices;
    1398:       PetscScalar       *array_x;
    1399:       const PetscScalar *array_diagonal;

    1401:       MatGetDiagonal(pc->pmat,pcis->vec1_global);
    1402:       VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);
    1403:       VecScatterBegin(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);
    1404:       VecScatterEnd(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);
    1405:       VecScatterBegin(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
    1406:       VecScatterEnd(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);
    1407:       ISGetLocalSize(dirIS,&dirsize);
    1408:       VecGetArray(pcis->vec1_N,&array_x);
    1409:       VecGetArrayRead(pcis->vec2_N,&array_diagonal);
    1410:       ISGetIndices(dirIS,(const PetscInt**)&is_indices);
    1411:       for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]];
    1412:       ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);
    1413:       VecRestoreArrayRead(pcis->vec2_N,&array_diagonal);
    1414:       VecRestoreArray(pcis->vec1_N,&array_x);
    1415:       VecScatterBegin(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);
    1416:       VecScatterEnd(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);
    1417:       pcbddc->rhs_change = PETSC_TRUE;
    1418:       ISDestroy(&dirIS);
    1419:     }
    1420:   }

    1422:   /* remove the computed solution or the initial guess from the rhs */
    1423:   if (pcbddc->rhs_change || (ksp && pcbddc->ksp_guess_nonzero)) {
    1424:     /* save the original rhs */
    1425:     if (save_rhs) {
    1426:       VecSwap(rhs,pcbddc->original_rhs);
    1427:       save_rhs = PETSC_FALSE;
    1428:     }
    1429:     pcbddc->rhs_change = PETSC_TRUE;
    1430:     VecScale(used_vec,-1.0);
    1431:     MatMultAdd(pc->mat,used_vec,pcbddc->original_rhs,rhs);
    1432:     VecScale(used_vec,-1.0);
    1433:     VecCopy(used_vec,pcbddc->temp_solution);
    1434:     pcbddc->temp_solution_used = PETSC_TRUE;
    1435:     if (ksp) {
    1436:       KSPSetInitialGuessNonzero(ksp,PETSC_FALSE);
    1437:     }
    1438:   }
    1439:   VecDestroy(&used_vec);

    1441:   /* compute initial vector in benign space if needed
    1442:      and remove non-benign solution from the rhs */
    1443:   benign_correction_computed = PETSC_FALSE;
    1444:   if (rhs && pcbddc->benign_compute_correction && (pcbddc->benign_have_null || pcbddc->benign_apply_coarse_only)) {
    1445:     /* compute u^*_h using ideas similar to those in Xuemin Tu's PhD thesis (see Section 4.8.1)
    1446:        Recursively apply BDDC in the multilevel case */
    1447:     if (!pcbddc->benign_vec) {
    1448:       VecDuplicate(rhs,&pcbddc->benign_vec);
    1449:     }
    1450:     /* keep applying coarse solver unless we no longer have benign subdomains */
    1451:     pcbddc->benign_apply_coarse_only = pcbddc->benign_have_null ? PETSC_TRUE : PETSC_FALSE;
    1452:     if (!pcbddc->benign_skip_correction) {
    1453:       PCApply_BDDC(pc,rhs,pcbddc->benign_vec);
    1454:       benign_correction_computed = PETSC_TRUE;
    1455:       if (pcbddc->temp_solution_used) {
    1456:         VecAXPY(pcbddc->temp_solution,1.0,pcbddc->benign_vec);
    1457:       }
    1458:       VecScale(pcbddc->benign_vec,-1.0);
    1459:       /* store the original rhs if not done earlier */
    1460:       if (save_rhs) {
    1461:         VecSwap(rhs,pcbddc->original_rhs);
    1462:       }
    1463:       if (pcbddc->rhs_change) {
    1464:         MatMultAdd(pc->mat,pcbddc->benign_vec,rhs,rhs);
    1465:       } else {
    1466:         MatMultAdd(pc->mat,pcbddc->benign_vec,pcbddc->original_rhs,rhs);
    1467:       }
    1468:       pcbddc->rhs_change = PETSC_TRUE;
    1469:     }
    1470:     pcbddc->benign_apply_coarse_only = PETSC_FALSE;
    1471:   } else {
    1472:     VecDestroy(&pcbddc->benign_vec);
    1473:   }

    1475:   /* dbg output */
    1476:   if (pcbddc->dbg_flag && benign_correction_computed) {
    1477:     Vec v;

    1479:     VecDuplicate(pcis->vec1_global,&v);
    1480:     if (pcbddc->ChangeOfBasisMatrix) {
    1481:       MatMultTranspose(pcbddc->ChangeOfBasisMatrix,rhs,v);
    1482:     } else {
    1483:       VecCopy(rhs,v);
    1484:     }
    1485:     PCBDDCBenignGetOrSetP0(pc,v,PETSC_TRUE);
    1486:     PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"LEVEL %D: is the correction benign?\n",pcbddc->current_level);
    1487:     PetscScalarView(pcbddc->benign_n,pcbddc->benign_p0,pcbddc->dbg_viewer);
    1488:     PetscViewerFlush(pcbddc->dbg_viewer);
    1489:     VecDestroy(&v);
    1490:   }

    1492:   /* set initial guess if using PCG */
    1493:   pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
    1494:   if (x && pcbddc->use_exact_dirichlet_trick) {
    1495:     VecSet(x,0.0);
    1496:     if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) {
    1497:       if (benign_correction_computed) { /* we have already saved the changed rhs */
    1498:         VecLockReadPop(pcis->vec1_global);
    1499:       } else {
    1500:         MatMultTranspose(pcbddc->ChangeOfBasisMatrix,rhs,pcis->vec1_global);
    1501:       }
    1502:       VecScatterBegin(pcis->global_to_D,pcis->vec1_global,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
    1503:       VecScatterEnd(pcis->global_to_D,pcis->vec1_global,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
    1504:     } else {
    1505:       VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
    1506:       VecScatterEnd(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
    1507:     }
    1508:     PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);
    1509:     KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);
    1510:     PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);
    1511:     KSPCheckSolve(pcbddc->ksp_D,pc,pcis->vec2_D);
    1512:     if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) {
    1513:       VecSet(pcis->vec1_global,0.);
    1514:       VecScatterBegin(pcis->global_to_D,pcis->vec2_D,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);
    1515:       VecScatterEnd(pcis->global_to_D,pcis->vec2_D,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);
    1516:       MatMult(pcbddc->ChangeOfBasisMatrix,pcis->vec1_global,x);
    1517:     } else {
    1518:       VecScatterBegin(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);
    1519:       VecScatterEnd(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);
    1520:     }
    1521:     if (ksp) {
    1522:       KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
    1523:     }
    1524:     pcbddc->exact_dirichlet_trick_app = PETSC_TRUE;
    1525:   } else if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior && benign_correction_computed && pcbddc->use_exact_dirichlet_trick) {
    1526:     VecLockReadPop(pcis->vec1_global);
    1527:   }
    1528:   return(0);
    1529: }

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

    1535:    Input Parameter:
    1536: +  pc - the preconditioner context

    1538:    Application Interface Routine: PCPostSolve()

    1540:    Notes:
    1541:      The interface routine PCPostSolve() is not usually called directly by
    1542:      the user, but instead is called by KSPSolve().
    1543: */
    1544: static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x)
    1545: {
    1547:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;

    1550:   /* add solution removed in presolve */
    1551:   if (x && pcbddc->rhs_change) {
    1552:     if (pcbddc->temp_solution_used) {
    1553:       VecAXPY(x,1.0,pcbddc->temp_solution);
    1554:     } else if (pcbddc->benign_compute_correction && pcbddc->benign_vec) {
    1555:       VecAXPY(x,-1.0,pcbddc->benign_vec);
    1556:     }
    1557:     /* restore to original state (not for FETI-DP) */
    1558:     if (ksp) pcbddc->temp_solution_used = PETSC_FALSE;
    1559:   }

    1561:   /* restore rhs to its original state (not needed for FETI-DP) */
    1562:   if (rhs && pcbddc->rhs_change) {
    1563:     VecSwap(rhs,pcbddc->original_rhs);
    1564:     pcbddc->rhs_change = PETSC_FALSE;
    1565:   }
    1566:   /* restore ksp guess state */
    1567:   if (ksp) {
    1568:     KSPSetInitialGuessNonzero(ksp,pcbddc->ksp_guess_nonzero);
    1569:     /* reset flag for exact dirichlet trick */
    1570:     pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
    1571:   }
    1572:   return(0);
    1573: }

    1575: /*
    1576:    PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner
    1577:                   by setting data structures and options.

    1579:    Input Parameter:
    1580: +  pc - the preconditioner context

    1582:    Application Interface Routine: PCSetUp()

    1584:    Notes:
    1585:      The interface routine PCSetUp() is not usually called directly by
    1586:      the user, but instead is called by PCApply() if necessary.
    1587: */
    1588: PetscErrorCode PCSetUp_BDDC(PC pc)
    1589: {
    1590:   PC_BDDC*        pcbddc = (PC_BDDC*)pc->data;
    1591:   PCBDDCSubSchurs sub_schurs;
    1592:   Mat_IS*         matis;
    1593:   MatNullSpace    nearnullspace;
    1594:   Mat             lA;
    1595:   IS              lP,zerodiag = NULL;
    1596:   PetscInt        nrows,ncols;
    1597:   PetscMPIInt     size;
    1598:   PetscBool       computesubschurs;
    1599:   PetscBool       computeconstraintsmatrix;
    1600:   PetscBool       new_nearnullspace_provided,ismatis,rl;
    1601:   PetscErrorCode  ierr;

    1604:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&ismatis);
    1605:   if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PCBDDC preconditioner requires matrix of type MATIS");
    1606:   MatGetSize(pc->pmat,&nrows,&ncols);
    1607:   if (nrows != ncols) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCBDDC preconditioner requires a square preconditioning matrix");
    1608:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);

    1610:   matis = (Mat_IS*)pc->pmat->data;
    1611:   /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */
    1612:   /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup
    1613:      Also, BDDC builds its own KSP for the Dirichlet problem */
    1614:   rl = pcbddc->recompute_topography;
    1615:   if (!pc->setupcalled || pc->flag == DIFFERENT_NONZERO_PATTERN) rl = PETSC_TRUE;
    1616:   MPIU_Allreduce(&rl,&pcbddc->recompute_topography,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)pc));
    1617:   if (pcbddc->recompute_topography) {
    1618:     pcbddc->graphanalyzed    = PETSC_FALSE;
    1619:     computeconstraintsmatrix = PETSC_TRUE;
    1620:   } else {
    1621:     computeconstraintsmatrix = PETSC_FALSE;
    1622:   }

    1624:   /* check parameters' compatibility */
    1625:   if (!pcbddc->use_deluxe_scaling) pcbddc->deluxe_zerorows = PETSC_FALSE;
    1626:   pcbddc->adaptive_selection   = (PetscBool)(pcbddc->adaptive_threshold[0] != 0.0 || pcbddc->adaptive_threshold[1] != 0.0);
    1627:   pcbddc->use_deluxe_scaling   = (PetscBool)(pcbddc->use_deluxe_scaling && size > 1);
    1628:   pcbddc->adaptive_selection   = (PetscBool)(pcbddc->adaptive_selection && size > 1);
    1629:   pcbddc->adaptive_userdefined = (PetscBool)(pcbddc->adaptive_selection && pcbddc->adaptive_userdefined);
    1630:   if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE;

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

    1634:   /* activate all connected components if the netflux has been requested */
    1635:   if (pcbddc->compute_nonetflux) {
    1636:     pcbddc->use_vertices = PETSC_TRUE;
    1637:     pcbddc->use_edges    = PETSC_TRUE;
    1638:     pcbddc->use_faces    = PETSC_TRUE;
    1639:   }

    1641:   /* Get stdout for dbg */
    1642:   if (pcbddc->dbg_flag) {
    1643:     if (!pcbddc->dbg_viewer) {
    1644:       pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc));
    1645:     }
    1646:     PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer);
    1647:     PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);
    1648:   }

    1650:   /* process topology information */
    1651:   PetscLogEventBegin(PC_BDDC_Topology[pcbddc->current_level],pc,0,0,0);
    1652:   if (pcbddc->recompute_topography) {
    1653:     PCBDDCComputeLocalTopologyInfo(pc);
    1654:     if (pcbddc->discretegradient) {
    1655:       PCBDDCNedelecSupport(pc);
    1656:     }
    1657:   }
    1658:   if (pcbddc->corner_selected) pcbddc->use_vertices = PETSC_TRUE;

    1660:   /* change basis if requested by the user */
    1661:   if (pcbddc->user_ChangeOfBasisMatrix) {
    1662:     /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */
    1663:     pcbddc->use_change_of_basis = PETSC_FALSE;
    1664:     PCBDDCComputeLocalMatrix(pc,pcbddc->user_ChangeOfBasisMatrix);
    1665:   } else {
    1666:     MatDestroy(&pcbddc->local_mat);
    1667:     PetscObjectReference((PetscObject)matis->A);
    1668:     pcbddc->local_mat = matis->A;
    1669:   }

    1671:   /*
    1672:      Compute change of basis on local pressures (aka zerodiag dofs) with the benign trick
    1673:      This should come earlier then PCISSetUp for extracting the correct subdomain matrices
    1674:   */
    1675:   PCBDDCBenignShellMat(pc,PETSC_TRUE);
    1676:   if (pcbddc->benign_saddle_point) {
    1677:     PC_IS* pcis = (PC_IS*)pc->data;

    1679:     if (pcbddc->user_ChangeOfBasisMatrix || pcbddc->use_change_of_basis || !computesubschurs) pcbddc->benign_change_explicit = PETSC_TRUE;
    1680:     /* detect local saddle point and change the basis in pcbddc->local_mat */
    1681:     PCBDDCBenignDetectSaddlePoint(pc,(PetscBool)(!pcbddc->recompute_topography),&zerodiag);
    1682:     /* pop B0 mat from local mat */
    1683:     PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);
    1684:     /* give pcis a hint to not reuse submatrices during PCISCreate */
    1685:     if (pc->flag == SAME_NONZERO_PATTERN && pcis->reusesubmatrices == PETSC_TRUE) {
    1686:       if (pcbddc->benign_n && (pcbddc->benign_change_explicit || pcbddc->dbg_flag)) {
    1687:         pcis->reusesubmatrices = PETSC_FALSE;
    1688:       } else {
    1689:         pcis->reusesubmatrices = PETSC_TRUE;
    1690:       }
    1691:     } else {
    1692:       pcis->reusesubmatrices = PETSC_FALSE;
    1693:     }
    1694:   }

    1696:   /* propagate relevant information */
    1697:   if (matis->A->symmetric_set) {
    1698:     MatSetOption(pcbddc->local_mat,MAT_SYMMETRIC,matis->A->symmetric);
    1699:   }
    1700:   if (matis->A->spd_set) {
    1701:     MatSetOption(pcbddc->local_mat,MAT_SPD,matis->A->spd);
    1702:   }

    1704:   /* Set up all the "iterative substructuring" common block without computing solvers */
    1705:   {
    1706:     Mat temp_mat;

    1708:     temp_mat = matis->A;
    1709:     matis->A = pcbddc->local_mat;
    1710:     PCISSetUp(pc,PETSC_TRUE,PETSC_FALSE);
    1711:     pcbddc->local_mat = matis->A;
    1712:     matis->A = temp_mat;
    1713:   }

    1715:   /* Analyze interface */
    1716:   if (!pcbddc->graphanalyzed) {
    1717:     PCBDDCAnalyzeInterface(pc);
    1718:     computeconstraintsmatrix = PETSC_TRUE;
    1719:     if (pcbddc->adaptive_selection && !pcbddc->use_deluxe_scaling && !pcbddc->mat_graph->twodim) {
    1720:       SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot compute the adaptive primal space for a problem with 3D edges without deluxe scaling");
    1721:     }
    1722:     if (pcbddc->compute_nonetflux) {
    1723:       MatNullSpace nnfnnsp;

    1725:       if (!pcbddc->divudotp) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Missing divudotp operator");
    1726:       PCBDDCComputeNoNetFlux(pc->pmat,pcbddc->divudotp,pcbddc->divudotp_trans,pcbddc->divudotp_vl2l,pcbddc->mat_graph,&nnfnnsp);
    1727:       /* TODO what if a nearnullspace is already attached? */
    1728:       if (nnfnnsp) {
    1729:         MatSetNearNullSpace(pc->pmat,nnfnnsp);
    1730:         MatNullSpaceDestroy(&nnfnnsp);
    1731:       }
    1732:     }
    1733:   }
    1734:   PetscLogEventEnd(PC_BDDC_Topology[pcbddc->current_level],pc,0,0,0);

    1736:   /* check existence of a divergence free extension, i.e.
    1737:      b(v_I,p_0) = 0 for all v_I (raise error if not).
    1738:      Also, check that PCBDDCBenignGetOrSetP0 works */
    1739:   if (pcbddc->benign_saddle_point && pcbddc->dbg_flag > 1) {
    1740:     PCBDDCBenignCheck(pc,zerodiag);
    1741:   }
    1742:   ISDestroy(&zerodiag);

    1744:   /* Setup local dirichlet solver ksp_D and sub_schurs solvers */
    1745:   if (computesubschurs && pcbddc->recompute_topography) {
    1746:     PCBDDCInitSubSchurs(pc);
    1747:   }
    1748:   /* SetUp Scaling operator (scaling matrices could be needed in SubSchursSetUp)*/
    1749:   if (!pcbddc->use_deluxe_scaling) {
    1750:     PCBDDCScalingSetUp(pc);
    1751:   }

    1753:   /* finish setup solvers and do adaptive selection of constraints */
    1754:   sub_schurs = pcbddc->sub_schurs;
    1755:   if (sub_schurs && sub_schurs->schur_explicit) {
    1756:     if (computesubschurs) {
    1757:       PCBDDCSetUpSubSchurs(pc);
    1758:     }
    1759:     PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);
    1760:   } else {
    1761:     PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);
    1762:     if (computesubschurs) {
    1763:       PCBDDCSetUpSubSchurs(pc);
    1764:     }
    1765:   }
    1766:   if (pcbddc->adaptive_selection) {
    1767:     PCBDDCAdaptiveSelection(pc);
    1768:     computeconstraintsmatrix = PETSC_TRUE;
    1769:   }

    1771:   /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */
    1772:   new_nearnullspace_provided = PETSC_FALSE;
    1773:   MatGetNearNullSpace(pc->pmat,&nearnullspace);
    1774:   if (pcbddc->onearnullspace) { /* already used nearnullspace */
    1775:     if (!nearnullspace) { /* near null space attached to mat has been destroyed */
    1776:       new_nearnullspace_provided = PETSC_TRUE;
    1777:     } else {
    1778:       /* determine if the two nullspaces are different (should be lightweight) */
    1779:       if (nearnullspace != pcbddc->onearnullspace) {
    1780:         new_nearnullspace_provided = PETSC_TRUE;
    1781:       } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */
    1782:         PetscInt         i;
    1783:         const Vec        *nearnullvecs;
    1784:         PetscObjectState state;
    1785:         PetscInt         nnsp_size;
    1786:         MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);
    1787:         for (i=0;i<nnsp_size;i++) {
    1788:           PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);
    1789:           if (pcbddc->onearnullvecs_state[i] != state) {
    1790:             new_nearnullspace_provided = PETSC_TRUE;
    1791:             break;
    1792:           }
    1793:         }
    1794:       }
    1795:     }
    1796:   } else {
    1797:     if (!nearnullspace) { /* both nearnullspaces are null */
    1798:       new_nearnullspace_provided = PETSC_FALSE;
    1799:     } else { /* nearnullspace attached later */
    1800:       new_nearnullspace_provided = PETSC_TRUE;
    1801:     }
    1802:   }

    1804:   /* Setup constraints and related work vectors */
    1805:   /* reset primal space flags */
    1806:   PetscLogEventBegin(PC_BDDC_LocalWork[pcbddc->current_level],pc,0,0,0);
    1807:   pcbddc->new_primal_space = PETSC_FALSE;
    1808:   pcbddc->new_primal_space_local = PETSC_FALSE;
    1809:   if (computeconstraintsmatrix || new_nearnullspace_provided) {
    1810:     /* It also sets the primal space flags */
    1811:     PCBDDCConstraintsSetUp(pc);
    1812:   }
    1813:   /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */
    1814:   PCBDDCSetUpLocalWorkVectors(pc);

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

    1819:     PCBDDCComputeLocalMatrix(pc,pcbddc->ChangeOfBasisMatrix);
    1820:     if (pcbddc->benign_change) {
    1821:       MatDestroy(&pcbddc->benign_B0);
    1822:       /* pop B0 from pcbddc->local_mat */
    1823:       PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);
    1824:     }
    1825:     /* get submatrices */
    1826:     MatDestroy(&pcis->A_IB);
    1827:     MatDestroy(&pcis->A_BI);
    1828:     MatDestroy(&pcis->A_BB);
    1829:     MatCreateSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);
    1830:     MatCreateSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);
    1831:     MatCreateSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);
    1832:     /* set flag in pcis to not reuse submatrices during PCISCreate */
    1833:     pcis->reusesubmatrices = PETSC_FALSE;
    1834:   } else if (!pcbddc->user_ChangeOfBasisMatrix && !pcbddc->benign_change) {
    1835:     MatDestroy(&pcbddc->local_mat);
    1836:     PetscObjectReference((PetscObject)matis->A);
    1837:     pcbddc->local_mat = matis->A;
    1838:   }

    1840:   /* interface pressure block row for B_C */
    1841:   PetscObjectQuery((PetscObject)pc,"__KSPFETIDP_lP" ,(PetscObject*)&lP);
    1842:   PetscObjectQuery((PetscObject)pc,"__KSPFETIDP_lA" ,(PetscObject*)&lA);
    1843:   if (lA && lP) {
    1844:     PC_IS*    pcis = (PC_IS*)pc->data;
    1845:     Mat       B_BI,B_BB,Bt_BI,Bt_BB;
    1846:     PetscBool issym;
    1847:     MatIsSymmetric(lA,PETSC_SMALL,&issym);
    1848:     if (issym) {
    1849:       MatCreateSubMatrix(lA,lP,pcis->is_I_local,MAT_INITIAL_MATRIX,&B_BI);
    1850:       MatCreateSubMatrix(lA,lP,pcis->is_B_local,MAT_INITIAL_MATRIX,&B_BB);
    1851:       MatCreateTranspose(B_BI,&Bt_BI);
    1852:       MatCreateTranspose(B_BB,&Bt_BB);
    1853:     } else {
    1854:       MatCreateSubMatrix(lA,lP,pcis->is_I_local,MAT_INITIAL_MATRIX,&B_BI);
    1855:       MatCreateSubMatrix(lA,lP,pcis->is_B_local,MAT_INITIAL_MATRIX,&B_BB);
    1856:       MatCreateSubMatrix(lA,pcis->is_I_local,lP,MAT_INITIAL_MATRIX,&Bt_BI);
    1857:       MatCreateSubMatrix(lA,pcis->is_B_local,lP,MAT_INITIAL_MATRIX,&Bt_BB);
    1858:     }
    1859:     PetscObjectCompose((PetscObject)pc,"__KSPFETIDP_B_BI",(PetscObject)B_BI);
    1860:     PetscObjectCompose((PetscObject)pc,"__KSPFETIDP_B_BB",(PetscObject)B_BB);
    1861:     PetscObjectCompose((PetscObject)pc,"__KSPFETIDP_Bt_BI",(PetscObject)Bt_BI);
    1862:     PetscObjectCompose((PetscObject)pc,"__KSPFETIDP_Bt_BB",(PetscObject)Bt_BB);
    1863:     MatDestroy(&B_BI);
    1864:     MatDestroy(&B_BB);
    1865:     MatDestroy(&Bt_BI);
    1866:     MatDestroy(&Bt_BB);
    1867:   }
    1868:   PetscLogEventEnd(PC_BDDC_LocalWork[pcbddc->current_level],pc,0,0,0);

    1870:   /* SetUp coarse and local Neumann solvers */
    1871:   PCBDDCSetUpSolvers(pc);
    1872:   /* SetUp Scaling operator */
    1873:   if (pcbddc->use_deluxe_scaling) {
    1874:     PCBDDCScalingSetUp(pc);
    1875:   }

    1877:   /* mark topography as done */
    1878:   pcbddc->recompute_topography = PETSC_FALSE;

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

    1883:   if (pcbddc->dbg_flag) {
    1884:     PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);
    1885:     PetscViewerASCIIPopSynchronized(pcbddc->dbg_viewer);
    1886:   }
    1887:   return(0);
    1888: }

    1890: /*
    1891:    PCApply_BDDC - Applies the BDDC operator to a vector.

    1893:    Input Parameters:
    1894: +  pc - the preconditioner context
    1895: -  r - input vector (global)

    1897:    Output Parameter:
    1898: .  z - output vector (global)

    1900:    Application Interface Routine: PCApply()
    1901:  */
    1902: PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z)
    1903: {
    1904:   PC_IS             *pcis = (PC_IS*)(pc->data);
    1905:   PC_BDDC           *pcbddc = (PC_BDDC*)(pc->data);
    1906:   Mat               lA = NULL;
    1907:   PetscInt          n_B = pcis->n_B, n_D = pcis->n - n_B;
    1908:   PetscErrorCode    ierr;
    1909:   const PetscScalar one = 1.0;
    1910:   const PetscScalar m_one = -1.0;
    1911:   const PetscScalar zero = 0.0;
    1912: /* This code is similar to that provided in nn.c for PCNN
    1913:    NN interface preconditioner changed to BDDC
    1914:    Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */

    1917:   PetscCitationsRegister(citation,&cited);
    1918:   if (pcbddc->switch_static) {
    1919:     MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat,&lA);
    1920:   }

    1922:   if (pcbddc->ChangeOfBasisMatrix) {
    1923:     Vec swap;

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

    1989:   /* Apply interface preconditioner
    1990:      input/output vecs: pcis->vec1_B and pcis->vec1_D */
    1991:   PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);

    1993:   /* Apply transpose of partition of unity operator */
    1994:   PCBDDCScalingExtension(pc,pcis->vec1_B,z);
    1995:   if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_LUMP) {
    1996:     VecScatterBegin(pcis->global_to_D,pcis->vec1_D,z,INSERT_VALUES,SCATTER_REVERSE);
    1997:     VecScatterEnd(pcis->global_to_D,pcis->vec1_D,z,INSERT_VALUES,SCATTER_REVERSE);
    1998:     return(0);
    1999:   }
    2000:   /* Second Dirichlet solve and assembling of output */
    2001:   VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
    2002:   VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
    2003:   if (n_B) {
    2004:     if (pcbddc->switch_static) {
    2005:       VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);
    2006:       VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);
    2007:       VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);
    2008:       VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);
    2009:       if (!pcbddc->switch_static_change) {
    2010:         MatMult(lA,pcis->vec1_N,pcis->vec2_N);
    2011:       } else {
    2012:         MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);
    2013:         MatMult(lA,pcis->vec2_N,pcis->vec1_N);
    2014:         MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);
    2015:       }
    2016:       VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);
    2017:       VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);
    2018:     } else {
    2019:       MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);
    2020:     }
    2021:   } else if (pcbddc->switch_static) { /* n_B is zero */
    2022:     if (!pcbddc->switch_static_change) {
    2023:       MatMult(lA,pcis->vec1_D,pcis->vec3_D);
    2024:     } else {
    2025:       MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);
    2026:       MatMult(lA,pcis->vec1_N,pcis->vec2_N);
    2027:       MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);
    2028:     }
    2029:   }
    2030:   PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);
    2031:   KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);
    2032:   PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);
    2033:   KSPCheckSolve(pcbddc->ksp_D,pc,pcis->vec4_D);

    2035:   if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
    2036:     if (pcbddc->switch_static) {
    2037:       VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);
    2038:     } else {
    2039:       VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);
    2040:     }
    2041:     VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);
    2042:     VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);
    2043:   } else {
    2044:     if (pcbddc->switch_static) {
    2045:       VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);
    2046:     } else {
    2047:       VecScale(pcis->vec4_D,m_one);
    2048:     }
    2049:     VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);
    2050:     VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);
    2051:   }
    2052:   if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */
    2053:     if (pcbddc->benign_apply_coarse_only) {
    2054:       PetscArrayzero(pcbddc->benign_p0,pcbddc->benign_n);
    2055:     }
    2056:     PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);
    2057:   }

    2059:   if (pcbddc->ChangeOfBasisMatrix) {
    2060:     pcbddc->work_change = r;
    2061:     VecCopy(z,pcbddc->work_change);
    2062:     MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);
    2063:   }
    2064:   return(0);
    2065: }

    2067: /*
    2068:    PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector.

    2070:    Input Parameters:
    2071: +  pc - the preconditioner context
    2072: -  r - input vector (global)

    2074:    Output Parameter:
    2075: .  z - output vector (global)

    2077:    Application Interface Routine: PCApplyTranspose()
    2078:  */
    2079: PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z)
    2080: {
    2081:   PC_IS             *pcis = (PC_IS*)(pc->data);
    2082:   PC_BDDC           *pcbddc = (PC_BDDC*)(pc->data);
    2083:   Mat               lA = NULL;
    2084:   PetscInt          n_B = pcis->n_B, n_D = pcis->n - n_B;
    2085:   PetscErrorCode    ierr;
    2086:   const PetscScalar one = 1.0;
    2087:   const PetscScalar m_one = -1.0;
    2088:   const PetscScalar zero = 0.0;

    2091:   PetscCitationsRegister(citation,&cited);
    2092:   if (pcbddc->switch_static) {
    2093:     MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat,&lA);
    2094:   }
    2095:   if (pcbddc->ChangeOfBasisMatrix) {
    2096:     Vec swap;

    2098:     MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);
    2099:     swap = pcbddc->work_change;
    2100:     pcbddc->work_change = r;
    2101:     r = swap;
    2102:     /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */
    2103:     if (pcbddc->benign_apply_coarse_only && pcbddc->exact_dirichlet_trick_app && pcbddc->change_interior) {
    2104:       VecCopy(r,pcis->vec1_global);
    2105:       VecLockReadPush(pcis->vec1_global);
    2106:     }
    2107:   }
    2108:   if (pcbddc->benign_have_null) { /* get p0 from r */
    2109:     PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);
    2110:   }
    2111:   if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
    2112:     VecCopy(r,z);
    2113:     /* First Dirichlet solve */
    2114:     VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
    2115:     VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);
    2116:     /*
    2117:       Assembling right hand side for BDDC operator
    2118:       - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE)
    2119:       - pcis->vec1_B the interface part of the global vector z
    2120:     */
    2121:     if (n_D) {
    2122:       PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);
    2123:       KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);
    2124:       PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);
    2125:       KSPCheckSolve(pcbddc->ksp_D,pc,pcis->vec2_D);
    2126:       VecScale(pcis->vec2_D,m_one);
    2127:       if (pcbddc->switch_static) {
    2128:         VecSet(pcis->vec1_N,0.);
    2129:         VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);
    2130:         VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);
    2131:         if (!pcbddc->switch_static_change) {
    2132:           MatMultTranspose(lA,pcis->vec1_N,pcis->vec2_N);
    2133:         } else {
    2134:           MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);
    2135:           MatMultTranspose(lA,pcis->vec2_N,pcis->vec1_N);
    2136:           MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);
    2137:         }
    2138:         VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);
    2139:         VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);
    2140:         VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
    2141:         VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
    2142:       } else {
    2143:         MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);
    2144:       }
    2145:     } else {
    2146:       VecSet(pcis->vec1_B,zero);
    2147:     }
    2148:     VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);
    2149:     VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);
    2150:     PCBDDCScalingRestriction(pc,z,pcis->vec1_B);
    2151:   } else {
    2152:     PCBDDCScalingRestriction(pc,r,pcis->vec1_B);
    2153:   }

    2155:   /* Apply interface preconditioner
    2156:      input/output vecs: pcis->vec1_B and pcis->vec1_D */
    2157:   PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);

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

    2162:   /* Second Dirichlet solve and assembling of output */
    2163:   VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
    2164:   VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);
    2165:   if (n_B) {
    2166:     if (pcbddc->switch_static) {
    2167:       VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);
    2168:       VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);
    2169:       VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);
    2170:       VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);
    2171:       if (!pcbddc->switch_static_change) {
    2172:         MatMultTranspose(lA,pcis->vec1_N,pcis->vec2_N);
    2173:       } else {
    2174:         MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);
    2175:         MatMultTranspose(lA,pcis->vec2_N,pcis->vec1_N);
    2176:         MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);
    2177:       }
    2178:       VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);
    2179:       VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);
    2180:     } else {
    2181:       MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);
    2182:     }
    2183:   } else if (pcbddc->switch_static) { /* n_B is zero */
    2184:     if (!pcbddc->switch_static_change) {
    2185:       MatMultTranspose(lA,pcis->vec1_D,pcis->vec3_D);
    2186:     } else {
    2187:       MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);
    2188:       MatMultTranspose(lA,pcis->vec1_N,pcis->vec2_N);
    2189:       MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);
    2190:     }
    2191:   }
    2192:   PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);
    2193:   KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);
    2194:   PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);
    2195:   KSPCheckSolve(pcbddc->ksp_D,pc,pcis->vec4_D);
    2196:   if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) {
    2197:     if (pcbddc->switch_static) {
    2198:       VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);
    2199:     } else {
    2200:       VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);
    2201:     }
    2202:     VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);
    2203:     VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);
    2204:   } else {
    2205:     if (pcbddc->switch_static) {
    2206:       VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);
    2207:     } else {
    2208:       VecScale(pcis->vec4_D,m_one);
    2209:     }
    2210:     VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);
    2211:     VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);
    2212:   }
    2213:   if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */
    2214:     PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);
    2215:   }
    2216:   if (pcbddc->ChangeOfBasisMatrix) {
    2217:     pcbddc->work_change = r;
    2218:     VecCopy(z,pcbddc->work_change);
    2219:     MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);
    2220:   }
    2221:   return(0);
    2222: }

    2224: PetscErrorCode PCReset_BDDC(PC pc)
    2225: {
    2226:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
    2227:   PC_IS          *pcis = (PC_IS*)pc->data;
    2228:   KSP            kspD,kspR,kspC;

    2232:   /* free BDDC custom data  */
    2233:   PCBDDCResetCustomization(pc);
    2234:   /* destroy objects related to topography */
    2235:   PCBDDCResetTopography(pc);
    2236:   /* destroy objects for scaling operator */
    2237:   PCBDDCScalingDestroy(pc);
    2238:   /* free solvers stuff */
    2239:   PCBDDCResetSolvers(pc);
    2240:   /* free global vectors needed in presolve */
    2241:   VecDestroy(&pcbddc->temp_solution);
    2242:   VecDestroy(&pcbddc->original_rhs);
    2243:   /* free data created by PCIS */
    2244:   PCISDestroy(pc);

    2246:   /* restore defaults */
    2247:   kspD = pcbddc->ksp_D;
    2248:   kspR = pcbddc->ksp_R;
    2249:   kspC = pcbddc->coarse_ksp;
    2250:   PetscMemzero(pc->data,sizeof(*pcbddc));
    2251:   pcis->n_neigh                     = -1;
    2252:   pcis->scaling_factor              = 1.0;
    2253:   pcis->reusesubmatrices            = PETSC_TRUE;
    2254:   pcbddc->use_local_adj             = PETSC_TRUE;
    2255:   pcbddc->use_vertices              = PETSC_TRUE;
    2256:   pcbddc->use_edges                 = PETSC_TRUE;
    2257:   pcbddc->symmetric_primal          = PETSC_TRUE;
    2258:   pcbddc->vertex_size               = 1;
    2259:   pcbddc->recompute_topography      = PETSC_TRUE;
    2260:   pcbddc->coarse_size               = -1;
    2261:   pcbddc->use_exact_dirichlet_trick = PETSC_TRUE;
    2262:   pcbddc->coarsening_ratio          = 8;
    2263:   pcbddc->coarse_eqs_per_proc       = 1;
    2264:   pcbddc->benign_compute_correction = PETSC_TRUE;
    2265:   pcbddc->nedfield                  = -1;
    2266:   pcbddc->nedglobal                 = PETSC_TRUE;
    2267:   pcbddc->graphmaxcount             = PETSC_MAX_INT;
    2268:   pcbddc->sub_schurs_layers         = -1;
    2269:   pcbddc->ksp_D                     = kspD;
    2270:   pcbddc->ksp_R                     = kspR;
    2271:   pcbddc->coarse_ksp                = kspC;
    2272:   return(0);
    2273: }

    2275: PetscErrorCode PCDestroy_BDDC(PC pc)
    2276: {
    2277:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;

    2281:   PCReset_BDDC(pc);
    2282:   KSPDestroy(&pcbddc->ksp_D);
    2283:   KSPDestroy(&pcbddc->ksp_R);
    2284:   KSPDestroy(&pcbddc->coarse_ksp);
    2285:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDiscreteGradient_C",NULL);
    2286:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",NULL);
    2287:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);
    2288:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);
    2289:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",NULL);
    2290:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);
    2291:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);
    2292:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);
    2293:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);
    2294:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);
    2295:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);
    2296:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);
    2297:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);
    2298:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);
    2299:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);
    2300:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);
    2301:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);
    2302:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);
    2303:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);
    2304:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);
    2305:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);
    2306:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);
    2307:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);
    2308:   PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",NULL);
    2309:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL);
    2310:   PetscFree(pc->data);
    2311:   return(0);
    2312: }

    2314: static PetscErrorCode PCSetCoordinates_BDDC(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
    2315: {
    2316:   PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
    2317:   PCBDDCGraph    mat_graph = pcbddc->mat_graph;

    2321:   PetscFree(mat_graph->coords);
    2322:   PetscMalloc1(nloc*dim,&mat_graph->coords);
    2323:   PetscArraycpy(mat_graph->coords,coords,nloc*dim);
    2324:   mat_graph->cnloc = nloc;
    2325:   mat_graph->cdim  = dim;
    2326:   mat_graph->cloc  = PETSC_FALSE;
    2327:   /* flg setup */
    2328:   pcbddc->recompute_topography = PETSC_TRUE;
    2329:   pcbddc->corner_selected = PETSC_FALSE;
    2330:   return(0);
    2331: }

    2333: static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool* change)
    2334: {
    2336:   *change = PETSC_TRUE;
    2337:   return(0);
    2338: }

    2340: static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
    2341: {
    2342:   FETIDPMat_ctx  mat_ctx;
    2343:   Vec            work;
    2344:   PC_IS*         pcis;
    2345:   PC_BDDC*       pcbddc;

    2349:   MatShellGetContext(fetidp_mat,&mat_ctx);
    2350:   pcis = (PC_IS*)mat_ctx->pc->data;
    2351:   pcbddc = (PC_BDDC*)mat_ctx->pc->data;

    2353:   VecSet(fetidp_flux_rhs,0.0);
    2354:   /* copy rhs since we may change it during PCPreSolve_BDDC */
    2355:   if (!pcbddc->original_rhs) {
    2356:     VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);
    2357:   }
    2358:   if (mat_ctx->rhs_flip) {
    2359:     VecPointwiseMult(pcbddc->original_rhs,standard_rhs,mat_ctx->rhs_flip);
    2360:   } else {
    2361:     VecCopy(standard_rhs,pcbddc->original_rhs);
    2362:   }
    2363:   if (mat_ctx->g2g_p) {
    2364:     /* interface pressure rhs */
    2365:     VecScatterBegin(mat_ctx->g2g_p,fetidp_flux_rhs,pcbddc->original_rhs,INSERT_VALUES,SCATTER_REVERSE);
    2366:     VecScatterEnd(mat_ctx->g2g_p,fetidp_flux_rhs,pcbddc->original_rhs,INSERT_VALUES,SCATTER_REVERSE);
    2367:     VecScatterBegin(mat_ctx->g2g_p,standard_rhs,fetidp_flux_rhs,INSERT_VALUES,SCATTER_FORWARD);
    2368:     VecScatterEnd(mat_ctx->g2g_p,standard_rhs,fetidp_flux_rhs,INSERT_VALUES,SCATTER_FORWARD);
    2369:     if (!mat_ctx->rhs_flip) {
    2370:       VecScale(fetidp_flux_rhs,-1.);
    2371:     }
    2372:   }
    2373:   /*
    2374:      change of basis for physical rhs if needed
    2375:      It also changes the rhs in case of dirichlet boundaries
    2376:   */
    2377:   PCPreSolve_BDDC(mat_ctx->pc,NULL,pcbddc->original_rhs,NULL);
    2378:   if (pcbddc->ChangeOfBasisMatrix) {
    2379:     MatMultTranspose(pcbddc->ChangeOfBasisMatrix,pcbddc->original_rhs,pcbddc->work_change);
    2380:     work = pcbddc->work_change;
    2381:    } else {
    2382:     work = pcbddc->original_rhs;
    2383:   }
    2384:   /* store vectors for computation of fetidp final solution */
    2385:   VecScatterBegin(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);
    2386:   VecScatterEnd(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);
    2387:   /* scale rhs since it should be unassembled */
    2388:   /* TODO use counter scaling? (also below) */
    2389:   VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);
    2390:   VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);
    2391:   /* Apply partition of unity */
    2392:   VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);
    2393:   /* PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B); */
    2394:   if (!pcbddc->switch_static) {
    2395:     /* compute partially subassembled Schur complement right-hand side */
    2396:     PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],mat_ctx->pc,0,0,0);
    2397:     KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);
    2398:     PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],mat_ctx->pc,0,0,0);
    2399:     /* Cannot propagate up error in KSPSolve() because there is no access to the PC */
    2400:     MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);
    2401:     VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);
    2402:     VecSet(work,0.0);
    2403:     VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);
    2404:     VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);
    2405:     /* PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B); */
    2406:     VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);
    2407:     VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);
    2408:     VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);
    2409:   }
    2410:   /* BDDC rhs */
    2411:   VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);
    2412:   if (pcbddc->switch_static) {
    2413:     VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);
    2414:   }
    2415:   /* apply BDDC */
    2416:   PetscArrayzero(pcbddc->benign_p0,pcbddc->benign_n);
    2417:   PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);
    2418:   PetscArrayzero(pcbddc->benign_p0,pcbddc->benign_n);

    2420:   /* Application of B_delta and assembling of rhs for fetidp fluxes */
    2421:   MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);
    2422:   VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);
    2423:   VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);
    2424:   /* Add contribution to interface pressures */
    2425:   if (mat_ctx->l2g_p) {
    2426:     MatMult(mat_ctx->B_BB,pcis->vec1_B,mat_ctx->vP);
    2427:     if (pcbddc->switch_static) {
    2428:       MatMultAdd(mat_ctx->B_BI,pcis->vec1_D,mat_ctx->vP,mat_ctx->vP);
    2429:     }
    2430:     VecScatterBegin(mat_ctx->l2g_p,mat_ctx->vP,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);
    2431:     VecScatterEnd(mat_ctx->l2g_p,mat_ctx->vP,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);
    2432:   }
    2433:   return(0);
    2434: }

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

    2439:    Collective

    2441:    Input Parameters:
    2442: +  fetidp_mat      - the FETI-DP matrix object obtained by a call to PCBDDCCreateFETIDPOperators
    2443: -  standard_rhs    - the right-hand side of the original linear system

    2445:    Output Parameters:
    2446: .  fetidp_flux_rhs - the right-hand side for the FETI-DP linear system

    2448:    Level: developer

    2450:    Notes:

    2452: .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetSolution
    2453: @*/
    2454: PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs)
    2455: {
    2456:   FETIDPMat_ctx  mat_ctx;

    2463:   MatShellGetContext(fetidp_mat,&mat_ctx);
    2464:   PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));
    2465:   return(0);
    2466: }

    2468: static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
    2469: {
    2470:   FETIDPMat_ctx  mat_ctx;
    2471:   PC_IS*         pcis;
    2472:   PC_BDDC*       pcbddc;
    2474:   Vec            work;

    2477:   MatShellGetContext(fetidp_mat,&mat_ctx);
    2478:   pcis = (PC_IS*)mat_ctx->pc->data;
    2479:   pcbddc = (PC_BDDC*)mat_ctx->pc->data;

    2481:   /* apply B_delta^T */
    2482:   VecSet(pcis->vec1_B,0.);
    2483:   VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);
    2484:   VecScatterEnd(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);
    2485:   MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);
    2486:   if (mat_ctx->l2g_p) {
    2487:     VecScatterBegin(mat_ctx->l2g_p,fetidp_flux_sol,mat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);
    2488:     VecScatterEnd(mat_ctx->l2g_p,fetidp_flux_sol,mat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);
    2489:     MatMultAdd(mat_ctx->Bt_BB,mat_ctx->vP,pcis->vec1_B,pcis->vec1_B);
    2490:   }

    2492:   /* compute rhs for BDDC application */
    2493:   VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);
    2494:   if (pcbddc->switch_static) {
    2495:     VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);
    2496:     if (mat_ctx->l2g_p) {
    2497:       VecScale(mat_ctx->vP,-1.);
    2498:       MatMultAdd(mat_ctx->Bt_BI,mat_ctx->vP,pcis->vec1_D,pcis->vec1_D);
    2499:     }
    2500:   }

    2502:   /* apply BDDC */
    2503:   PetscArrayzero(pcbddc->benign_p0,pcbddc->benign_n);
    2504:   PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);

    2506:   /* put values into global vector */
    2507:   if (pcbddc->ChangeOfBasisMatrix) work = pcbddc->work_change;
    2508:   else work = standard_sol;
    2509:   VecScatterBegin(pcis->global_to_B,pcis->vec1_B,work,INSERT_VALUES,SCATTER_REVERSE);
    2510:   VecScatterEnd(pcis->global_to_B,pcis->vec1_B,work,INSERT_VALUES,SCATTER_REVERSE);
    2511:   if (!pcbddc->switch_static) {
    2512:     /* compute values into the interior if solved for the partially subassembled Schur complement */
    2513:     MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);
    2514:     VecAYPX(pcis->vec1_D,-1.0,mat_ctx->temp_solution_D);
    2515:     PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],mat_ctx->pc,0,0,0);
    2516:     KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec1_D);
    2517:     PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],mat_ctx->pc,0,0,0);
    2518:     /* Cannot propagate up error in KSPSolve() because there is no access to the PC */
    2519:   }

    2521:   VecScatterBegin(pcis->global_to_D,pcis->vec1_D,work,INSERT_VALUES,SCATTER_REVERSE);
    2522:   VecScatterEnd(pcis->global_to_D,pcis->vec1_D,work,INSERT_VALUES,SCATTER_REVERSE);
    2523:   /* add p0 solution to final solution */
    2524:   PCBDDCBenignGetOrSetP0(mat_ctx->pc,work,PETSC_FALSE);
    2525:   if (pcbddc->ChangeOfBasisMatrix) {
    2526:     MatMult(pcbddc->ChangeOfBasisMatrix,work,standard_sol);
    2527:   }
    2528:   PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);
    2529:   if (mat_ctx->g2g_p) {
    2530:     VecScatterBegin(mat_ctx->g2g_p,fetidp_flux_sol,standard_sol,INSERT_VALUES,SCATTER_REVERSE);
    2531:     VecScatterEnd(mat_ctx->g2g_p,fetidp_flux_sol,standard_sol,INSERT_VALUES,SCATTER_REVERSE);
    2532:   }
    2533:   return(0);
    2534: }

    2536: static PetscErrorCode PCView_BDDCIPC(PC pc, PetscViewer viewer)
    2537: {
    2539:   BDDCIPC_ctx    bddcipc_ctx;
    2540:   PetscBool      isascii;

    2543:   PCShellGetContext(pc,(void **)&bddcipc_ctx);
    2544:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
    2545:   if (isascii) {
    2546:     PetscViewerASCIIPrintf(viewer,"BDDC interface preconditioner\n");
    2547:   }
    2548:   PetscViewerASCIIPushTab(viewer);
    2549:   PCView(bddcipc_ctx->bddc,viewer);
    2550:   PetscViewerASCIIPopTab(viewer);
    2551:   return(0);
    2552: }

    2554: static PetscErrorCode PCSetUp_BDDCIPC(PC pc)
    2555: {
    2557:   BDDCIPC_ctx    bddcipc_ctx;
    2558:   PetscBool      isbddc;
    2559:   Vec            vv;
    2560:   IS             is;
    2561:   PC_IS          *pcis;

    2564:   PCShellGetContext(pc,(void **)&bddcipc_ctx);
    2565:   PetscObjectTypeCompare((PetscObject)bddcipc_ctx->bddc,PCBDDC,&isbddc);
    2566:   if (!isbddc) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid type %s. Must be of type bddc",((PetscObject)bddcipc_ctx->bddc)->type_name);
    2567:   PCSetUp(bddcipc_ctx->bddc);

    2569:   /* create interface scatter */
    2570:   pcis = (PC_IS*)(bddcipc_ctx->bddc->data);
    2571:   VecScatterDestroy(&bddcipc_ctx->g2l);
    2572:   MatCreateVecs(pc->pmat,&vv,NULL);
    2573:   ISRenumber(pcis->is_B_global,NULL,NULL,&is);
    2574:   VecScatterCreate(vv,is,pcis->vec1_B,NULL,&bddcipc_ctx->g2l);
    2575:   ISDestroy(&is);
    2576:   VecDestroy(&vv);
    2577:   return(0);
    2578: }

    2580: static PetscErrorCode PCApply_BDDCIPC(PC pc, Vec r, Vec x)
    2581: {
    2583:   BDDCIPC_ctx    bddcipc_ctx;
    2584:   PC_IS          *pcis;
    2585:   VecScatter     tmps;

    2588:   PCShellGetContext(pc,(void **)&bddcipc_ctx);
    2589:   pcis = (PC_IS*)(bddcipc_ctx->bddc->data);
    2590:   tmps = pcis->global_to_B;
    2591:   pcis->global_to_B = bddcipc_ctx->g2l;
    2592:   PCBDDCScalingRestriction(bddcipc_ctx->bddc,r,pcis->vec1_B);
    2593:   PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc,PETSC_FALSE);
    2594:   PCBDDCScalingExtension(bddcipc_ctx->bddc,pcis->vec1_B,x);
    2595:   pcis->global_to_B = tmps;
    2596:   return(0);
    2597: }

    2599: static PetscErrorCode PCApplyTranspose_BDDCIPC(PC pc, Vec r, Vec x)
    2600: {
    2602:   BDDCIPC_ctx    bddcipc_ctx;
    2603:   PC_IS          *pcis;
    2604:   VecScatter     tmps;

    2607:   PCShellGetContext(pc,(void **)&bddcipc_ctx);
    2608:   pcis = (PC_IS*)(bddcipc_ctx->bddc->data);
    2609:   tmps = pcis->global_to_B;
    2610:   pcis->global_to_B = bddcipc_ctx->g2l;
    2611:   PCBDDCScalingRestriction(bddcipc_ctx->bddc,r,pcis->vec1_B);
    2612:   PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc,PETSC_TRUE);
    2613:   PCBDDCScalingExtension(bddcipc_ctx->bddc,pcis->vec1_B,x);
    2614:   pcis->global_to_B = tmps;
    2615:   return(0);
    2616: }

    2618: static PetscErrorCode PCDestroy_BDDCIPC(PC pc)
    2619: {
    2621:   BDDCIPC_ctx    bddcipc_ctx;

    2624:   PCShellGetContext(pc,(void **)&bddcipc_ctx);
    2625:   PCDestroy(&bddcipc_ctx->bddc);
    2626:   VecScatterDestroy(&bddcipc_ctx->g2l);
    2627:   PetscFree(bddcipc_ctx);
    2628:   return(0);
    2629: }

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

    2634:    Collective

    2636:    Input Parameters:
    2637: +  fetidp_mat      - the FETI-DP matrix obtained by a call to PCBDDCCreateFETIDPOperators
    2638: -  fetidp_flux_sol - the solution of the FETI-DP linear system

    2640:    Output Parameters:
    2641: .  standard_sol    - the solution defined on the physical domain

    2643:    Level: developer

    2645:    Notes:

    2647: .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetRHS
    2648: @*/
    2649: PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol)
    2650: {
    2651:   FETIDPMat_ctx  mat_ctx;

    2658:   MatShellGetContext(fetidp_mat,&mat_ctx);
    2659:   PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));
    2660:   return(0);
    2661: }

    2663: static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, PetscBool fully_redundant, const char* prefix, Mat *fetidp_mat, PC *fetidp_pc)
    2664: {

    2666:   FETIDPMat_ctx  fetidpmat_ctx;
    2667:   Mat            newmat;
    2668:   FETIDPPC_ctx   fetidppc_ctx;
    2669:   PC             newpc;
    2670:   MPI_Comm       comm;

    2674:   PetscObjectGetComm((PetscObject)pc,&comm);
    2675:   /* FETI-DP matrix */
    2676:   PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);
    2677:   fetidpmat_ctx->fully_redundant = fully_redundant;
    2678:   PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);
    2679:   MatCreateShell(comm,fetidpmat_ctx->n,fetidpmat_ctx->n,fetidpmat_ctx->N,fetidpmat_ctx->N,fetidpmat_ctx,&newmat);
    2680:   PetscObjectSetName((PetscObject)newmat,!fetidpmat_ctx->l2g_lambda_only ? "F" : "G");
    2681:   MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);
    2682:   MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);
    2683:   MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);
    2684:   /* propagate MatOptions */
    2685:   {
    2686:     PC_BDDC   *pcbddc = (PC_BDDC*)fetidpmat_ctx->pc->data;
    2687:     PetscBool issym;

    2689:     MatGetOption(pc->mat,MAT_SYMMETRIC,&issym);
    2690:     if (issym || pcbddc->symmetric_primal) {
    2691:       MatSetOption(newmat,MAT_SYMMETRIC,PETSC_TRUE);
    2692:     }
    2693:   }
    2694:   MatSetOptionsPrefix(newmat,prefix);
    2695:   MatAppendOptionsPrefix(newmat,"fetidp_");
    2696:   MatSetUp(newmat);
    2697:   /* FETI-DP preconditioner */
    2698:   PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);
    2699:   PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);
    2700:   PCCreate(comm,&newpc);
    2701:   PCSetOperators(newpc,newmat,newmat);
    2702:   PCSetOptionsPrefix(newpc,prefix);
    2703:   PCAppendOptionsPrefix(newpc,"fetidp_");
    2704:   PCSetErrorIfFailure(newpc,pc->erroriffailure);
    2705:   if (!fetidpmat_ctx->l2g_lambda_only) { /* standard FETI-DP */
    2706:     PCSetType(newpc,PCSHELL);
    2707:     PCShellSetName(newpc,"FETI-DP multipliers");
    2708:     PCShellSetContext(newpc,fetidppc_ctx);
    2709:     PCShellSetApply(newpc,FETIDPPCApply);
    2710:     PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);
    2711:     PCShellSetView(newpc,FETIDPPCView);
    2712:     PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);
    2713:   } else { /* saddle-point FETI-DP */
    2714:     Mat       M;
    2715:     PetscInt  psize;
    2716:     PetscBool fake = PETSC_FALSE, isfieldsplit;

    2718:     ISViewFromOptions(fetidpmat_ctx->lagrange,NULL,"-lag_view");
    2719:     ISViewFromOptions(fetidpmat_ctx->pressure,NULL,"-press_view");
    2720:     PetscObjectQuery((PetscObject)pc,"__KSPFETIDP_PPmat",(PetscObject*)&M);
    2721:     PCSetType(newpc,PCFIELDSPLIT);
    2722:     PCFieldSplitSetIS(newpc,"lag",fetidpmat_ctx->lagrange);
    2723:     PCFieldSplitSetIS(newpc,"p",fetidpmat_ctx->pressure);
    2724:     PCFieldSplitSetType(newpc,PC_COMPOSITE_SCHUR);
    2725:     PCFieldSplitSetSchurFactType(newpc,PC_FIELDSPLIT_SCHUR_FACT_DIAG);
    2726:     ISGetSize(fetidpmat_ctx->pressure,&psize);
    2727:     if (psize != M->rmap->N) {
    2728:       Mat      M2;
    2729:       PetscInt lpsize;

    2731:       fake = PETSC_TRUE;
    2732:       ISGetLocalSize(fetidpmat_ctx->pressure,&lpsize);
    2733:       MatCreate(comm,&M2);
    2734:       MatSetType(M2,MATAIJ);
    2735:       MatSetSizes(M2,lpsize,lpsize,psize,psize);
    2736:       MatSetUp(M2);
    2737:       MatAssemblyBegin(M2,MAT_FINAL_ASSEMBLY);
    2738:       MatAssemblyEnd(M2,MAT_FINAL_ASSEMBLY);
    2739:       PCFieldSplitSetSchurPre(newpc,PC_FIELDSPLIT_SCHUR_PRE_USER,M2);
    2740:       MatDestroy(&M2);
    2741:     } else {
    2742:       PCFieldSplitSetSchurPre(newpc,PC_FIELDSPLIT_SCHUR_PRE_USER,M);
    2743:     }
    2744:     PCFieldSplitSetSchurScale(newpc,1.0);

    2746:     /* we need to setfromoptions and setup here to access the blocks */
    2747:     PCSetFromOptions(newpc);
    2748:     PCSetUp(newpc);

    2750:     /* user may have changed the type (e.g. -fetidp_pc_type none) */
    2751:     PetscObjectTypeCompare((PetscObject)newpc,PCFIELDSPLIT,&isfieldsplit);
    2752:     if (isfieldsplit) {
    2753:       KSP       *ksps;
    2754:       PC        ppc,lagpc;
    2755:       PetscInt  nn;
    2756:       PetscBool ismatis,matisok = PETSC_FALSE,check = PETSC_FALSE;

    2758:       /* set the solver for the (0,0) block */
    2759:       PCFieldSplitSchurGetSubKSP(newpc,&nn,&ksps);
    2760:       if (!nn) { /* not of type PC_COMPOSITE_SCHUR */
    2761:         PCFieldSplitGetSubKSP(newpc,&nn,&ksps);
    2762:         if (!fake) { /* pass pmat to the pressure solver */
    2763:           Mat F;

    2765:           KSPGetOperators(ksps[1],&F,NULL);
    2766:           KSPSetOperators(ksps[1],F,M);
    2767:         }
    2768:       } else {
    2769:         PetscBool issym;
    2770:         Mat       S;

    2772:         PCFieldSplitSchurGetS(newpc,&S);

    2774:         MatGetOption(newmat,MAT_SYMMETRIC,&issym);
    2775:         if (issym) {
    2776:           MatSetOption(S,MAT_SYMMETRIC,PETSC_TRUE);
    2777:         }
    2778:       }
    2779:       KSPGetPC(ksps[0],&lagpc);
    2780:       PCSetType(lagpc,PCSHELL);
    2781:       PCShellSetName(lagpc,"FETI-DP multipliers");
    2782:       PCShellSetContext(lagpc,fetidppc_ctx);
    2783:       PCShellSetApply(lagpc,FETIDPPCApply);
    2784:       PCShellSetApplyTranspose(lagpc,FETIDPPCApplyTranspose);
    2785:       PCShellSetView(lagpc,FETIDPPCView);
    2786:       PCShellSetDestroy(lagpc,PCBDDCDestroyFETIDPPC);

    2788:       /* Olof's idea: interface Schur complement preconditioner for the mass matrix */
    2789:       KSPGetPC(ksps[1],&ppc);
    2790:       if (fake) {
    2791:         BDDCIPC_ctx    bddcipc_ctx;
    2792:         PetscContainer c;

    2794:         matisok = PETSC_TRUE;

    2796:         /* create inner BDDC solver */
    2797:         PetscNew(&bddcipc_ctx);
    2798:         PCCreate(comm,&bddcipc_ctx->bddc);
    2799:         PCSetType(bddcipc_ctx->bddc,PCBDDC);
    2800:         PCSetOperators(bddcipc_ctx->bddc,M,M);
    2801:         PetscObjectQuery((PetscObject)pc,"__KSPFETIDP_pCSR",(PetscObject*)&c);
    2802:         PetscObjectTypeCompare((PetscObject)M,MATIS,&ismatis);
    2803:         if (c && ismatis) {
    2804:           Mat      lM;
    2805:           PetscInt *csr,n;

    2807:           MatISGetLocalMat(M,&lM);
    2808:           MatGetSize(lM,&n,NULL);
    2809:           PetscContainerGetPointer(c,(void**)&csr);
    2810:           PCBDDCSetLocalAdjacencyGraph(bddcipc_ctx->bddc,n,csr,csr + (n + 1),PETSC_COPY_VALUES);
    2811:           MatISRestoreLocalMat(M,&lM);
    2812:         }
    2813:         PCSetOptionsPrefix(bddcipc_ctx->bddc,((PetscObject)ksps[1])->prefix);
    2814:         PCSetErrorIfFailure(bddcipc_ctx->bddc,pc->erroriffailure);
    2815:         PCSetFromOptions(bddcipc_ctx->bddc);

    2817:         /* wrap the interface application */
    2818:         PCSetType(ppc,PCSHELL);
    2819:         PCShellSetName(ppc,"FETI-DP pressure");
    2820:         PCShellSetContext(ppc,bddcipc_ctx);
    2821:         PCShellSetSetUp(ppc,PCSetUp_BDDCIPC);
    2822:         PCShellSetApply(ppc,PCApply_BDDCIPC);
    2823:         PCShellSetApplyTranspose(ppc,PCApplyTranspose_BDDCIPC);
    2824:         PCShellSetView(ppc,PCView_BDDCIPC);
    2825:         PCShellSetDestroy(ppc,PCDestroy_BDDCIPC);
    2826:       }

    2828:       /* determine if we need to assemble M to construct a preconditioner */
    2829:       if (!matisok) {
    2830:         PetscObjectTypeCompare((PetscObject)M,MATIS,&ismatis);
    2831:         PetscObjectTypeCompareAny((PetscObject)ppc,&matisok,PCBDDC,PCJACOBI,PCNONE,PCMG,"");
    2832:         if (ismatis && !matisok) {
    2833:           MatConvert(M,MATAIJ,MAT_INPLACE_MATRIX,&M);
    2834:         }
    2835:       }

    2837:       /* run the subproblems to check convergence */
    2838:       PetscOptionsGetBool(NULL,((PetscObject)newmat)->prefix,"-check_saddlepoint",&check,NULL);
    2839:       if (check) {
    2840:         PetscInt i;

    2842:         for (i=0;i<nn;i++) {
    2843:           KSP       kspC;
    2844:           PC        pc;
    2845:           Mat       F,pF;
    2846:           Vec       x,y;
    2847:           PetscBool isschur,prec = PETSC_TRUE;

    2849:           KSPCreate(PetscObjectComm((PetscObject)ksps[i]),&kspC);
    2850:           KSPSetOptionsPrefix(kspC,((PetscObject)ksps[i])->prefix);
    2851:           KSPAppendOptionsPrefix(kspC,"check_");
    2852:           KSPGetOperators(ksps[i],&F,&pF);
    2853:           PetscObjectTypeCompare((PetscObject)F,MATSCHURCOMPLEMENT,&isschur);
    2854:           if (isschur) {
    2855:             KSP  kspS,kspS2;
    2856:             Mat  A00,pA00,A10,A01,A11;
    2857:             char prefix[256];

    2859:             MatSchurComplementGetKSP(F,&kspS);
    2860:             MatSchurComplementGetSubMatrices(F,&A00,&pA00,&A01,&A10,&A11);
    2861:             MatCreateSchurComplement(A00,pA00,A01,A10,A11,&F);
    2862:             MatSchurComplementGetKSP(F,&kspS2);
    2863:             PetscSNPrintf(prefix,sizeof(prefix),"%sschur_",((PetscObject)kspC)->prefix);
    2864:             KSPSetOptionsPrefix(kspS2,prefix);
    2865:             KSPGetPC(kspS2,&pc);
    2866:             PCSetType(pc,PCKSP);
    2867:             PCKSPSetKSP(pc,kspS);
    2868:             KSPSetFromOptions(kspS2);
    2869:             KSPGetPC(kspS2,&pc);
    2870:             PCSetUseAmat(pc,PETSC_TRUE);
    2871:           } else {
    2872:             PetscObjectReference((PetscObject)F);
    2873:           }
    2874:           KSPSetFromOptions(kspC);
    2875:           PetscOptionsGetBool(NULL,((PetscObject)kspC)->prefix,"-preconditioned",&prec,NULL);
    2876:           if (prec)  {
    2877:             KSPGetPC(ksps[i],&pc);
    2878:             KSPSetPC(kspC,pc);
    2879:           }
    2880:           KSPSetOperators(kspC,F,pF);
    2881:           MatCreateVecs(F,&x,&y);
    2882:           VecSetRandom(x,NULL);
    2883:           MatMult(F,x,y);
    2884:           KSPSolve(kspC,y,x);
    2885:           KSPCheckSolve(kspC,pc,x);
    2886:           KSPDestroy(&kspC);
    2887:           MatDestroy(&F);
    2888:           VecDestroy(&x);
    2889:           VecDestroy(&y);
    2890:         }
    2891:       }
    2892:       PetscFree(ksps);
    2893:     }
    2894:   }
    2895:   /* return pointers for objects created */
    2896:   *fetidp_mat = newmat;
    2897:   *fetidp_pc  = newpc;
    2898:   return(0);
    2899: }

    2901: /*@C
    2902:  PCBDDCCreateFETIDPOperators - Create FETI-DP operators

    2904:    Collective

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

    2911:    Output Parameters:
    2912: +  fetidp_mat - shell FETI-DP matrix object
    2913: -  fetidp_pc  - shell Dirichlet preconditioner for FETI-DP matrix

    2915:    Level: developer

    2917:    Notes:
    2918:      Currently the only operations provided for FETI-DP matrix are MatMult and MatMultTranspose

    2920: .seealso: PCBDDC, PCBDDCMatFETIDPGetRHS, PCBDDCMatFETIDPGetSolution
    2921: @*/
    2922: PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc)
    2923: {

    2928:   if (pc->setupcalled) {
    2929:     PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,PetscBool,const char*,Mat*,PC*),(pc,fully_redundant,prefix,fetidp_mat,fetidp_pc));
    2930:   } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"You must call PCSetup_BDDC() first");
    2931:   return(0);
    2932: }
    2933: /* -------------------------------------------------------------------------- */
    2934: /*MC
    2935:    PCBDDC - Balancing Domain Decomposition by Constraints.

    2937:    An implementation of the BDDC preconditioner based on the bibliography found below.

    2939:    The matrix to be preconditioned (Pmat) must be of type MATIS.

    2941:    Currently works with MATIS matrices with local matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers.

    2943:    It also works with unsymmetric and indefinite problems.

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

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

    2949:    Boundary nodes are split in vertices, edges and faces classes using information from the local to global mapping of dofs and the local connectivity graph of nodes. The latter can be customized by using PCBDDCSetLocalAdjacencyGraph()
    2950:    Additional information on dofs can be provided by using PCBDDCSetDofsSplitting(), PCBDDCSetDirichletBoundaries(), PCBDDCSetNeumannBoundaries(), and PCBDDCSetPrimalVerticesIS() and their local counterparts.

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

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

    2957:    The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using a MatPartitioning object.

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

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

    2963:    Options Database Keys (some of them, run with -help for a complete list):

    2965: +    -pc_bddc_use_vertices <true> - use or not vertices in primal space
    2966: .    -pc_bddc_use_edges <true> - use or not edges in primal space
    2967: .    -pc_bddc_use_faces <false> - use or not faces in primal space
    2968: .    -pc_bddc_symmetric <true> - symmetric computation of primal basis functions. Specify false for unsymmetric problems
    2969: .    -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only)
    2970: .    -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested
    2971: .    -pc_bddc_switch_static <false> - switches from M_2 (default) to M_3 operator (see reference article [1])
    2972: .    -pc_bddc_levels <0> - maximum number of levels for multilevel
    2973: .    -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)
    2974: .    -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)
    2975: .    -pc_bddc_use_deluxe_scaling <false> - use deluxe scaling
    2976: .    -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)
    2977: .    -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)
    2978: -    -pc_bddc_check_level <0> - set verbosity level of debugging output

    2980:    Options for Dirichlet, Neumann or coarse solver can be set with
    2981: .vb
    2982:       -pc_bddc_dirichlet_
    2983:       -pc_bddc_neumann_
    2984:       -pc_bddc_coarse_
    2985: .ve
    2986:    e.g. -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. PCBDDC uses by default KSPPREONLY and PCLU.

    2988:    When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified as
    2989: .vb
    2990:       -pc_bddc_dirichlet_lN_
    2991:       -pc_bddc_neumann_lN_
    2992:       -pc_bddc_coarse_lN_
    2993: .ve
    2994:    Note that level number ranges from the finest (0) to the coarsest (N).
    2995:    In order to specify options for the BDDC operators at the coarser levels (and not for the solvers), prepend -pc_bddc_coarse_ or -pc_bddc_coarse_l to the option, e.g.
    2996: .vb
    2997:      -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3
    2998: .ve
    2999:    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

    3001:    References:
    3002: +   [1] - C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149--168, March 2007
    3003: .   [2] - 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
    3004: .   [3] - J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", Computing Volume 83, Issue 2--3, pages 55--85, November 2008
    3005: -   [4] - C. Pechstein and C. R. Dohrmann. "Modern domain decomposition methods BDDC, deluxe scaling, and an algebraic approach", Seminar talk, Linz, December 2013, http://people.ricam.oeaw.ac.at/c.pechstein/pechstein-bddc2013.pdf

    3007:    Level: intermediate

    3009:    Developer Notes:

    3011:    Contributed by Stefano Zampini

    3013: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,  MATIS
    3014: M*/

    3016: PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc)
    3017: {
    3018:   PetscErrorCode      ierr;
    3019:   PC_BDDC             *pcbddc;

    3022:   PetscNewLog(pc,&pcbddc);
    3023:   pc->data = (void*)pcbddc;

    3025:   /* create PCIS data structure */
    3026:   PCISCreate(pc);

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

    3031:   /* BDDC nonzero defaults */
    3032:   pcbddc->use_nnsp                  = PETSC_TRUE;
    3033:   pcbddc->use_local_adj             = PETSC_TRUE;
    3034:   pcbddc->use_vertices              = PETSC_TRUE;
    3035:   pcbddc->use_edges                 = PETSC_TRUE;
    3036:   pcbddc->symmetric_primal          = PETSC_TRUE;
    3037:   pcbddc->vertex_size               = 1;
    3038:   pcbddc->recompute_topography      = PETSC_TRUE;
    3039:   pcbddc->coarse_size               = -1;
    3040:   pcbddc->use_exact_dirichlet_trick = PETSC_TRUE;
    3041:   pcbddc->coarsening_ratio          = 8;
    3042:   pcbddc->coarse_eqs_per_proc       = 1;
    3043:   pcbddc->benign_compute_correction = PETSC_TRUE;
    3044:   pcbddc->nedfield                  = -1;
    3045:   pcbddc->nedglobal                 = PETSC_TRUE;
    3046:   pcbddc->graphmaxcount             = PETSC_MAX_INT;
    3047:   pcbddc->sub_schurs_layers         = -1;
    3048:   pcbddc->adaptive_threshold[0]     = 0.0;
    3049:   pcbddc->adaptive_threshold[1]     = 0.0;

    3051:   /* function pointers */
    3052:   pc->ops->apply               = PCApply_BDDC;
    3053:   pc->ops->applytranspose      = PCApplyTranspose_BDDC;
    3054:   pc->ops->setup               = PCSetUp_BDDC;
    3055:   pc->ops->destroy             = PCDestroy_BDDC;
    3056:   pc->ops->setfromoptions      = PCSetFromOptions_BDDC;
    3057:   pc->ops->view                = PCView_BDDC;
    3058:   pc->ops->applyrichardson     = NULL;
    3059:   pc->ops->applysymmetricleft  = NULL;
    3060:   pc->ops->applysymmetricright = NULL;
    3061:   pc->ops->presolve            = PCPreSolve_BDDC;
    3062:   pc->ops->postsolve           = PCPostSolve_BDDC;
    3063:   pc->ops->reset               = PCReset_BDDC;

    3065:   /* composing function */
    3066:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDiscreteGradient_C",PCBDDCSetDiscreteGradient_BDDC);
    3067:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",PCBDDCSetDivergenceMat_BDDC);
    3068:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);
    3069:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);
    3070:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",PCBDDCSetPrimalVerticesIS_BDDC);
    3071:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetPrimalVerticesLocalIS_C",PCBDDCGetPrimalVerticesLocalIS_BDDC);
    3072:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetPrimalVerticesIS_C",PCBDDCGetPrimalVerticesIS_BDDC);
    3073:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);
    3074:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);
    3075:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);
    3076:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);
    3077:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);
    3078:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);
    3079:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);
    3080:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);
    3081:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);
    3082:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);
    3083:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);
    3084:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);
    3085:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);
    3086:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);
    3087:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);
    3088:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);
    3089:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);
    3090:   PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);
    3091:   PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_BDDC);
    3092:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_BDDC);
    3093:   return(0);
    3094: }

    3096: /*@C
    3097:  PCBDDCInitializePackage - This function initializes everything in the PCBDDC package. It is called
    3098:     from PCInitializePackage().

    3100:  Level: developer

    3102:  .seealso: PetscInitialize()
    3103: @*/
    3104: PetscErrorCode PCBDDCInitializePackage(void)
    3105: {
    3107:   int            i;

    3110:   if (PCBDDCPackageInitialized) return(0);
    3111:   PCBDDCPackageInitialized = PETSC_TRUE;
    3112:   PetscRegisterFinalize(PCBDDCFinalizePackage);

    3114:   /* general events */
    3115:   PetscLogEventRegister("PCBDDCTopo",PC_CLASSID,&PC_BDDC_Topology[0]);
    3116:   PetscLogEventRegister("PCBDDCLKSP",PC_CLASSID,&PC_BDDC_LocalSolvers[0]);
    3117:   PetscLogEventRegister("PCBDDCLWor",PC_CLASSID,&PC_BDDC_LocalWork[0]);
    3118:   PetscLogEventRegister("PCBDDCCorr",PC_CLASSID,&PC_BDDC_CorrectionSetUp[0]);
    3119:   PetscLogEventRegister("PCBDDCASet",PC_CLASSID,&PC_BDDC_ApproxSetUp[0]);
    3120:   PetscLogEventRegister("PCBDDCAApp",PC_CLASSID,&PC_BDDC_ApproxApply[0]);
    3121:   PetscLogEventRegister("PCBDDCCSet",PC_CLASSID,&PC_BDDC_CoarseSetUp[0]);
    3122:   PetscLogEventRegister("PCBDDCCKSP",PC_CLASSID,&PC_BDDC_CoarseSolver[0]);
    3123:   PetscLogEventRegister("PCBDDCAdap",PC_CLASSID,&PC_BDDC_AdaptiveSetUp[0]);
    3124:   PetscLogEventRegister("PCBDDCScal",PC_CLASSID,&PC_BDDC_Scaling[0]);
    3125:   PetscLogEventRegister("PCBDDCSchr",PC_CLASSID,&PC_BDDC_Schurs[0]);
    3126:   PetscLogEventRegister("PCBDDCDirS",PC_CLASSID,&PC_BDDC_Solves[0][0]);
    3127:   PetscLogEventRegister("PCBDDCNeuS",PC_CLASSID,&PC_BDDC_Solves[0][1]);
    3128:   PetscLogEventRegister("PCBDDCCoaS",PC_CLASSID,&PC_BDDC_Solves[0][2]);
    3129:   for (i=1;i<PETSC_PCBDDC_MAXLEVELS;i++) {
    3130:     char ename[32];

    3132:     PetscSNPrintf(ename,sizeof(ename),"PCBDDCTopo l%02d",i);
    3133:     PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Topology[i]);
    3134:     PetscSNPrintf(ename,sizeof(ename),"PCBDDCLKSP l%02d",i);
    3135:     PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_LocalSolvers[i]);
    3136:     PetscSNPrintf(ename,sizeof(ename),"PCBDDCLWor l%02d",i);
    3137:     PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_LocalWork[i]);
    3138:     PetscSNPrintf(ename,sizeof(ename),"PCBDDCCorr l%02d",i);
    3139:     PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_CorrectionSetUp[i]);
    3140:     PetscSNPrintf(ename,sizeof(ename),"PCBDDCASet l%02d",i);
    3141:     PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_ApproxSetUp[i]);
    3142:     PetscSNPrintf(ename,sizeof(ename),"PCBDDCAApp l%02d",i);
    3143:     PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_ApproxApply[i]);
    3144:     PetscSNPrintf(ename,sizeof(ename),"PCBDDCCSet l%02d",i);
    3145:     PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_CoarseSetUp[i]);
    3146:     PetscSNPrintf(ename,sizeof(ename),"PCBDDCCKSP l%02d",i);
    3147:     PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_CoarseSolver[i]);
    3148:     PetscSNPrintf(ename,sizeof(ename),"PCBDDCAdap l%02d",i);
    3149:     PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_AdaptiveSetUp[i]);
    3150:     PetscSNPrintf(ename,sizeof(ename),"PCBDDCScal l%02d",i);
    3151:     PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Scaling[i]);
    3152:     PetscSNPrintf(ename,sizeof(ename),"PCBDDCSchr l%02d",i);
    3153:     PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Schurs[i]);
    3154:     PetscSNPrintf(ename,sizeof(ename),"PCBDDCDirS l%02d",i);
    3155:     PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Solves[i][0]);
    3156:     PetscSNPrintf(ename,sizeof(ename),"PCBDDCNeuS l%02d",i);
    3157:     PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Solves[i][1]);
    3158:     PetscSNPrintf(ename,sizeof(ename),"PCBDDCCoaS l%02d",i);
    3159:     PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Solves[i][2]);
    3160:   }
    3161:   return(0);
    3162: }

    3164: /*@C
    3165:  PCBDDCFinalizePackage - This function frees everything from the PCBDDC package. It is
    3166:     called from PetscFinalize() automatically.

    3168:  Level: developer

    3170:  .seealso: PetscFinalize()
    3171: @*/
    3172: PetscErrorCode PCBDDCFinalizePackage(void)
    3173: {
    3175:   PCBDDCPackageInitialized = PETSC_FALSE;
    3176:   return(0);
    3177: }