1 /* TODOLIST 2 3 Solvers 4 - Add support for cholesky for coarse solver (similar to local solvers) 5 - Propagate ksp prefixes for solvers to mat objects? 6 7 User interface 8 - ** DM attached to pc? 9 10 Debugging output 11 - * Better management of verbosity levels of debugging output 12 13 Extra 14 - *** Is it possible to work with PCBDDCGraph on boundary indices only (less memory consumed)? 15 - BDDC with MG framework? 16 17 MATIS related operations contained in BDDC code 18 - Provide general case for subassembling 19 20 */ 21 22 #include <../src/ksp/pc/impls/bddc/bddc.h> /*I "petscpc.h" I*/ /* includes for fortran wrappers */ 23 #include <../src/ksp/pc/impls/bddc/bddcprivate.h> 24 #include <petscblaslapack.h> 25 26 static PetscBool PCBDDCPackageInitialized = PETSC_FALSE; 27 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"; 42 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]; 55 56 const char *const PCBDDCInterfaceExtTypes[] = {"DIRICHLET","LUMP","PCBDDCInterfaceExtType","PC_BDDC_INTERFACE_EXT_",NULL}; 57 58 PetscErrorCode PCApply_BDDC(PC,Vec,Vec); 59 60 PetscErrorCode PCSetFromOptions_BDDC(PetscOptionItems *PetscOptionsObject,PC pc) 61 { 62 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 63 PetscInt nt,i; 64 PetscErrorCode ierr; 65 66 PetscFunctionBegin; 67 ierr = PetscOptionsHead(PetscOptionsObject,"BDDC options");CHKERRQ(ierr); 68 /* Verbose debugging */ 69 ierr = PetscOptionsInt("-pc_bddc_check_level","Verbose output for PCBDDC (intended for debug)","none",pcbddc->dbg_flag,&pcbddc->dbg_flag,NULL);CHKERRQ(ierr); 70 /* Approximate solvers */ 71 ierr = 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);CHKERRQ(ierr); 72 if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_DIRICHLET) { 73 ierr = PetscOptionsBool("-pc_bddc_dirichlet_approximate","Inform PCBDDC that we are using approximate Dirichlet solvers","none",pcbddc->NullSpace_corr[0],&pcbddc->NullSpace_corr[0],NULL);CHKERRQ(ierr); 74 ierr = PetscOptionsBool("-pc_bddc_dirichlet_approximate_scale","Inform PCBDDC that we need to scale the Dirichlet solve","none",pcbddc->NullSpace_corr[1],&pcbddc->NullSpace_corr[1],NULL);CHKERRQ(ierr); 75 } else { 76 /* This flag is needed/implied by lumping */ 77 pcbddc->switch_static = PETSC_TRUE; 78 } 79 ierr = PetscOptionsBool("-pc_bddc_neumann_approximate","Inform PCBDDC that we are using approximate Neumann solvers","none",pcbddc->NullSpace_corr[2],&pcbddc->NullSpace_corr[2],NULL);CHKERRQ(ierr); 80 ierr = PetscOptionsBool("-pc_bddc_neumann_approximate_scale","Inform PCBDDC that we need to scale the Neumann solve","none",pcbddc->NullSpace_corr[3],&pcbddc->NullSpace_corr[3],NULL);CHKERRQ(ierr); 81 /* Primal space customization */ 82 ierr = PetscOptionsBool("-pc_bddc_use_local_mat_graph","Use or not adjacency graph of local mat for interface analysis","none",pcbddc->use_local_adj,&pcbddc->use_local_adj,NULL);CHKERRQ(ierr); 83 ierr = PetscOptionsInt("-pc_bddc_graph_maxcount","Maximum number of shared subdomains for a connected component","none",pcbddc->graphmaxcount,&pcbddc->graphmaxcount,NULL);CHKERRQ(ierr); 84 ierr = PetscOptionsBool("-pc_bddc_corner_selection","Activates face-based corner selection","none",pcbddc->corner_selection,&pcbddc->corner_selection,NULL);CHKERRQ(ierr); 85 ierr = PetscOptionsBool("-pc_bddc_use_vertices","Use or not corner dofs in coarse space","none",pcbddc->use_vertices,&pcbddc->use_vertices,NULL);CHKERRQ(ierr); 86 ierr = PetscOptionsBool("-pc_bddc_use_edges","Use or not edge constraints in coarse space","none",pcbddc->use_edges,&pcbddc->use_edges,NULL);CHKERRQ(ierr); 87 ierr = PetscOptionsBool("-pc_bddc_use_faces","Use or not face constraints in coarse space","none",pcbddc->use_faces,&pcbddc->use_faces,NULL);CHKERRQ(ierr); 88 ierr = 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);CHKERRQ(ierr); 89 ierr = PetscOptionsBool("-pc_bddc_use_nnsp","Use near null space attached to the matrix to compute constraints","none",pcbddc->use_nnsp,&pcbddc->use_nnsp,NULL);CHKERRQ(ierr); 90 ierr = 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);CHKERRQ(ierr); 91 ierr = 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);CHKERRQ(ierr); 92 /* Change of basis */ 93 ierr = PetscOptionsBool("-pc_bddc_use_change_of_basis","Use or not internal change of basis on local edge nodes","none",pcbddc->use_change_of_basis,&pcbddc->use_change_of_basis,NULL);CHKERRQ(ierr); 94 ierr = PetscOptionsBool("-pc_bddc_use_change_on_faces","Use or not internal change of basis on local face nodes","none",pcbddc->use_change_on_faces,&pcbddc->use_change_on_faces,NULL);CHKERRQ(ierr); 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 ierr = PetscOptionsBool("-pc_bddc_switch_static","Switch on static condensation ops around the interface preconditioner","none",pcbddc->switch_static,&pcbddc->switch_static,NULL);CHKERRQ(ierr); 100 ierr = 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);CHKERRQ(ierr); 101 i = pcbddc->coarsening_ratio; 102 ierr = PetscOptionsInt("-pc_bddc_coarsening_ratio","Set coarsening ratio used in multilevel coarsening","PCBDDCSetCoarseningRatio",i,&i,NULL);CHKERRQ(ierr); 103 ierr = PCBDDCSetCoarseningRatio(pc,i);CHKERRQ(ierr); 104 i = pcbddc->max_levels; 105 ierr = PetscOptionsInt("-pc_bddc_levels","Set maximum number of levels for multilevel","PCBDDCSetLevels",i,&i,NULL);CHKERRQ(ierr); 106 ierr = PCBDDCSetLevels(pc,i);CHKERRQ(ierr); 107 ierr = 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);CHKERRQ(ierr); 108 ierr = PetscOptionsBool("-pc_bddc_use_coarse_estimates","Use estimated eigenvalues for coarse problem","none",pcbddc->use_coarse_estimates,&pcbddc->use_coarse_estimates,NULL);CHKERRQ(ierr); 109 ierr = PetscOptionsBool("-pc_bddc_use_deluxe_scaling","Use deluxe scaling for BDDC","none",pcbddc->use_deluxe_scaling,&pcbddc->use_deluxe_scaling,NULL);CHKERRQ(ierr); 110 ierr = PetscOptionsBool("-pc_bddc_schur_rebuild","Whether or not the interface graph for Schur principal minors has to be rebuilt (i.e. define the interface without any adjacency)","none",pcbddc->sub_schurs_rebuild,&pcbddc->sub_schurs_rebuild,NULL);CHKERRQ(ierr); 111 ierr = PetscOptionsInt("-pc_bddc_schur_layers","Number of dofs' layers for the computation of principal minors (i.e. -1 uses all dofs)","none",pcbddc->sub_schurs_layers,&pcbddc->sub_schurs_layers,NULL);CHKERRQ(ierr); 112 ierr = PetscOptionsBool("-pc_bddc_schur_use_useradj","Whether or not the CSR graph specified by the user should be used for computing successive layers (default is to use adj of local mat)","none",pcbddc->sub_schurs_use_useradj,&pcbddc->sub_schurs_use_useradj,NULL);CHKERRQ(ierr); 113 ierr = 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);CHKERRQ(ierr); 114 ierr = PetscOptionsBool("-pc_bddc_deluxe_zerorows","Zero rows and columns of deluxe operators associated with primal dofs","none",pcbddc->deluxe_zerorows,&pcbddc->deluxe_zerorows,NULL);CHKERRQ(ierr); 115 ierr = PetscOptionsBool("-pc_bddc_deluxe_singlemat","Collapse deluxe operators","none",pcbddc->deluxe_singlemat,&pcbddc->deluxe_singlemat,NULL);CHKERRQ(ierr); 116 ierr = 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);CHKERRQ(ierr); 117 nt = 2; 118 ierr = PetscOptionsRealArray("-pc_bddc_adaptive_threshold","Thresholds to be used for adaptive selection of constraints","none",pcbddc->adaptive_threshold,&nt,NULL);CHKERRQ(ierr); 119 if (nt == 1) pcbddc->adaptive_threshold[1] = pcbddc->adaptive_threshold[0]; 120 ierr = PetscOptionsInt("-pc_bddc_adaptive_nmin","Minimum number of constraints per connected components","none",pcbddc->adaptive_nmin,&pcbddc->adaptive_nmin,NULL);CHKERRQ(ierr); 121 ierr = PetscOptionsInt("-pc_bddc_adaptive_nmax","Maximum number of constraints per connected components","none",pcbddc->adaptive_nmax,&pcbddc->adaptive_nmax,NULL);CHKERRQ(ierr); 122 ierr = PetscOptionsBool("-pc_bddc_symmetric","Symmetric computation of primal basis functions","none",pcbddc->symmetric_primal,&pcbddc->symmetric_primal,NULL);CHKERRQ(ierr); 123 ierr = PetscOptionsInt("-pc_bddc_coarse_adj","Number of processors where to map the coarse adjacency list","none",pcbddc->coarse_adj_red,&pcbddc->coarse_adj_red,NULL);CHKERRQ(ierr); 124 ierr = 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);CHKERRQ(ierr); 125 ierr = PetscOptionsBool("-pc_bddc_benign_change","Compute the pressure change of basis explicitly","none",pcbddc->benign_change_explicit,&pcbddc->benign_change_explicit,NULL);CHKERRQ(ierr); 126 ierr = PetscOptionsBool("-pc_bddc_benign_compute_correction","Compute the benign correction during PreSolve","none",pcbddc->benign_compute_correction,&pcbddc->benign_compute_correction,NULL);CHKERRQ(ierr); 127 ierr = PetscOptionsBool("-pc_bddc_nonetflux","Automatic computation of no-net-flux quadrature weights","none",pcbddc->compute_nonetflux,&pcbddc->compute_nonetflux,NULL);CHKERRQ(ierr); 128 ierr = PetscOptionsBool("-pc_bddc_detect_disconnected","Detects disconnected subdomains","none",pcbddc->detect_disconnected,&pcbddc->detect_disconnected,NULL);CHKERRQ(ierr); 129 ierr = 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);CHKERRQ(ierr); 130 ierr = PetscOptionsBool("-pc_bddc_eliminate_dirichlet","Whether or not we want to eliminate dirichlet dofs during presolve","none",pcbddc->eliminate_dirdofs,&pcbddc->eliminate_dirdofs,NULL);CHKERRQ(ierr); 131 ierr = PetscOptionsTail();CHKERRQ(ierr); 132 PetscFunctionReturn(0); 133 } 134 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; 143 144 PetscFunctionBegin; 145 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 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; 153 154 if (!pc->setupcalled) { 155 ierr = PetscViewerASCIIPrintf(viewer," Partial information available: preconditioner has not been setup yet\n");CHKERRQ(ierr); 156 } 157 ierr = PetscViewerASCIIPrintf(viewer," Use verbose output: %D\n",pcbddc->dbg_flag);CHKERRQ(ierr); 158 ierr = PetscViewerASCIIPrintf(viewer," Use user-defined CSR: %d\n",!!pcbddc->mat_graph->nvtxs_csr);CHKERRQ(ierr); 159 ierr = PetscViewerASCIIPrintf(viewer," Use local mat graph: %d\n",pcbddc->use_local_adj && !pcbddc->mat_graph->nvtxs_csr);CHKERRQ(ierr); 160 if (pcbddc->mat_graph->twodim) { 161 ierr = PetscViewerASCIIPrintf(viewer," Connectivity graph topological dimension: 2\n");CHKERRQ(ierr); 162 } else { 163 ierr = PetscViewerASCIIPrintf(viewer," Connectivity graph topological dimension: 3\n");CHKERRQ(ierr); 164 } 165 if (pcbddc->graphmaxcount != PETSC_MAX_INT) { 166 ierr = PetscViewerASCIIPrintf(viewer," Graph max count: %D\n",pcbddc->graphmaxcount);CHKERRQ(ierr); 167 } 168 ierr = PetscViewerASCIIPrintf(viewer," Use vertices: %d (vertex size %D)\n",pcbddc->use_vertices,pcbddc->vertex_size);CHKERRQ(ierr); 169 ierr = PetscViewerASCIIPrintf(viewer," Use edges: %d\n",pcbddc->use_edges);CHKERRQ(ierr); 170 ierr = PetscViewerASCIIPrintf(viewer," Use faces: %d\n",pcbddc->use_faces);CHKERRQ(ierr); 171 ierr = PetscViewerASCIIPrintf(viewer," Use true near null space: %d\n",pcbddc->use_nnsp_true);CHKERRQ(ierr); 172 ierr = PetscViewerASCIIPrintf(viewer," Use QR for single constraints on cc: %d\n",pcbddc->use_qr_single);CHKERRQ(ierr); 173 ierr = PetscViewerASCIIPrintf(viewer," Use change of basis on local edge nodes: %d\n",pcbddc->use_change_of_basis);CHKERRQ(ierr); 174 ierr = PetscViewerASCIIPrintf(viewer," Use change of basis on local face nodes: %d\n",pcbddc->use_change_on_faces);CHKERRQ(ierr); 175 ierr = PetscViewerASCIIPrintf(viewer," User defined change of basis matrix: %d\n",!!pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 176 ierr = PetscViewerASCIIPrintf(viewer," Has change of basis matrix: %d\n",!!pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 177 ierr = PetscViewerASCIIPrintf(viewer," Eliminate dirichlet boundary dofs: %d\n",pcbddc->eliminate_dirdofs);CHKERRQ(ierr); 178 ierr = PetscViewerASCIIPrintf(viewer," Switch on static condensation ops around the interface preconditioner: %d\n",pcbddc->switch_static);CHKERRQ(ierr); 179 ierr = PetscViewerASCIIPrintf(viewer," Use exact dirichlet trick: %d\n",pcbddc->use_exact_dirichlet_trick);CHKERRQ(ierr); 180 ierr = PetscViewerASCIIPrintf(viewer," Interface extension: %s\n",PCBDDCInterfaceExtTypes[pcbddc->interface_extension]);CHKERRQ(ierr); 181 ierr = PetscViewerASCIIPrintf(viewer," Multilevel max levels: %D\n",pcbddc->max_levels);CHKERRQ(ierr); 182 ierr = PetscViewerASCIIPrintf(viewer," Multilevel coarsening ratio: %D\n",pcbddc->coarsening_ratio);CHKERRQ(ierr); 183 ierr = PetscViewerASCIIPrintf(viewer," Use estimated eigs for coarse problem: %d\n",pcbddc->use_coarse_estimates);CHKERRQ(ierr); 184 ierr = PetscViewerASCIIPrintf(viewer," Use deluxe scaling: %d\n",pcbddc->use_deluxe_scaling);CHKERRQ(ierr); 185 ierr = PetscViewerASCIIPrintf(viewer," Use deluxe zerorows: %d\n",pcbddc->deluxe_zerorows);CHKERRQ(ierr); 186 ierr = PetscViewerASCIIPrintf(viewer," Use deluxe singlemat: %d\n",pcbddc->deluxe_singlemat);CHKERRQ(ierr); 187 ierr = PetscViewerASCIIPrintf(viewer," Rebuild interface graph for Schur principal minors: %d\n",pcbddc->sub_schurs_rebuild);CHKERRQ(ierr); 188 ierr = PetscViewerASCIIPrintf(viewer," Number of dofs' layers for the computation of principal minors: %D\n",pcbddc->sub_schurs_layers);CHKERRQ(ierr); 189 ierr = PetscViewerASCIIPrintf(viewer," Use user CSR graph to compute successive layers: %d\n",pcbddc->sub_schurs_use_useradj);CHKERRQ(ierr); 190 if (pcbddc->adaptive_threshold[1] != pcbddc->adaptive_threshold[0]) { 191 ierr = 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]);CHKERRQ(ierr); 192 } else { 193 ierr = PetscViewerASCIIPrintf(viewer," Adaptive constraint selection threshold (active %d, userdefined %d): %g\n",pcbddc->adaptive_selection,pcbddc->adaptive_userdefined,pcbddc->adaptive_threshold[0]);CHKERRQ(ierr); 194 } 195 ierr = PetscViewerASCIIPrintf(viewer," Min constraints / connected component: %D\n",pcbddc->adaptive_nmin);CHKERRQ(ierr); 196 ierr = PetscViewerASCIIPrintf(viewer," Max constraints / connected component: %D\n",pcbddc->adaptive_nmax);CHKERRQ(ierr); 197 ierr = PetscViewerASCIIPrintf(viewer," Invert exact Schur complement for adaptive selection: %d\n",pcbddc->sub_schurs_exact_schur);CHKERRQ(ierr); 198 ierr = PetscViewerASCIIPrintf(viewer," Symmetric computation of primal basis functions: %d\n",pcbddc->symmetric_primal);CHKERRQ(ierr); 199 ierr = PetscViewerASCIIPrintf(viewer," Num. Procs. to map coarse adjacency list: %D\n",pcbddc->coarse_adj_red);CHKERRQ(ierr); 200 ierr = PetscViewerASCIIPrintf(viewer," Coarse eqs per proc (significant at the coarsest level): %D\n",pcbddc->coarse_eqs_per_proc);CHKERRQ(ierr); 201 ierr = PetscViewerASCIIPrintf(viewer," Detect disconnected: %d (filter %d)\n",pcbddc->detect_disconnected,pcbddc->detect_disconnected_filter);CHKERRQ(ierr); 202 ierr = PetscViewerASCIIPrintf(viewer," Benign subspace trick: %d (change explicit %d)\n",pcbddc->benign_saddle_point,pcbddc->benign_change_explicit);CHKERRQ(ierr); 203 ierr = PetscViewerASCIIPrintf(viewer," Benign subspace trick is active: %d\n",pcbddc->benign_have_null);CHKERRQ(ierr); 204 ierr = PetscViewerASCIIPrintf(viewer," Algebraic computation of no-net-flux: %d\n",pcbddc->compute_nonetflux);CHKERRQ(ierr); 205 if (!pc->setupcalled) PetscFunctionReturn(0); 206 207 /* compute interface size */ 208 ierr = VecSet(pcis->vec1_B,1.0);CHKERRQ(ierr); 209 ierr = MatCreateVecs(pc->pmat,&counter,NULL);CHKERRQ(ierr); 210 ierr = VecSet(counter,0.0);CHKERRQ(ierr); 211 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,counter,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 212 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,counter,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 213 ierr = VecSum(counter,&interface_size);CHKERRQ(ierr); 214 ierr = VecDestroy(&counter);CHKERRQ(ierr); 215 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 ierr = MPI_Reduce(loc,gsum,6,MPIU_INT64,MPI_SUM,0,PetscObjectComm((PetscObject)pc));CHKERRMPI(ierr); 227 if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = loc[5] = -1; 228 ierr = MPI_Reduce(loc,gmax,6,MPIU_INT64,MPI_MAX,0,PetscObjectComm((PetscObject)pc));CHKERRMPI(ierr); 229 if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = loc[5] = PETSC_MAX_INT; 230 ierr = MPI_Reduce(loc,gmin,6,MPIU_INT64,MPI_MIN,0,PetscObjectComm((PetscObject)pc));CHKERRMPI(ierr); 231 ierr = MPI_Reduce(&loc[6],&totbenign,1,MPIU_INT64,MPI_SUM,0,PetscObjectComm((PetscObject)pc));CHKERRMPI(ierr); 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 ierr = PetscViewerASCIIPrintf(viewer,"********************************** STATISTICS AT LEVEL %d **********************************\n",pcbddc->current_level);CHKERRQ(ierr); 237 ierr = PetscViewerASCIIPrintf(viewer," Global dofs sizes: all %D interface %D coarse %D\n",pc->pmat->rmap->N,(PetscInt)PetscRealPart(interface_size),pcbddc->coarse_size);CHKERRQ(ierr); 238 ierr = PetscViewerASCIIPrintf(viewer," Coarsening ratios: all/coarse %D interface/coarse %D\n",(PetscInt)ratio1,(PetscInt)ratio2);CHKERRQ(ierr); 239 ierr = PetscViewerASCIIPrintf(viewer," Active processes : %D\n",(PetscInt)gsum[0]);CHKERRQ(ierr); 240 ierr = PetscViewerASCIIPrintf(viewer," Total subdomains : %D\n",(PetscInt)gsum[5]);CHKERRQ(ierr); 241 if (pcbddc->benign_have_null) { 242 ierr = PetscViewerASCIIPrintf(viewer," Benign subs : %D\n",(PetscInt)totbenign);CHKERRQ(ierr); 243 } 244 ierr = PetscViewerASCIIPrintf(viewer," Dofs type :\tMIN\tMAX\tMEAN\n");CHKERRQ(ierr); 245 ierr = PetscViewerASCIIPrintf(viewer," Interior dofs :\t%D\t%D\t%D\n",(PetscInt)gmin[1],(PetscInt)gmax[1],(PetscInt)(gsum[1]/gsum[0]));CHKERRQ(ierr); 246 ierr = PetscViewerASCIIPrintf(viewer," Interface dofs :\t%D\t%D\t%D\n",(PetscInt)gmin[2],(PetscInt)gmax[2],(PetscInt)(gsum[2]/gsum[0]));CHKERRQ(ierr); 247 ierr = PetscViewerASCIIPrintf(viewer," Primal dofs :\t%D\t%D\t%D\n",(PetscInt)gmin[3],(PetscInt)gmax[3],(PetscInt)(gsum[3]/gsum[0]));CHKERRQ(ierr); 248 ierr = PetscViewerASCIIPrintf(viewer," Local dofs :\t%D\t%D\t%D\n",(PetscInt)gmin[4],(PetscInt)gmax[4],(PetscInt)(gsum[4]/gsum[0]));CHKERRQ(ierr); 249 ierr = PetscViewerASCIIPrintf(viewer," Local subs :\t%D\t%D\n" ,(PetscInt)gmin[5],(PetscInt)gmax[5]);CHKERRQ(ierr); 250 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 251 252 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);CHKERRMPI(ierr); 253 254 /* local solvers */ 255 ierr = PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)pcbddc->ksp_D),&subviewer);CHKERRQ(ierr); 256 if (rank == 0) { 257 ierr = PetscViewerASCIIPrintf(subviewer,"--- Interior solver (rank 0)\n");CHKERRQ(ierr); 258 ierr = PetscViewerASCIIPushTab(subviewer);CHKERRQ(ierr); 259 ierr = KSPView(pcbddc->ksp_D,subviewer);CHKERRQ(ierr); 260 ierr = PetscViewerASCIIPopTab(subviewer);CHKERRQ(ierr); 261 ierr = PetscViewerASCIIPrintf(subviewer,"--- Correction solver (rank 0)\n");CHKERRQ(ierr); 262 ierr = PetscViewerASCIIPushTab(subviewer);CHKERRQ(ierr); 263 ierr = KSPView(pcbddc->ksp_R,subviewer);CHKERRQ(ierr); 264 ierr = PetscViewerASCIIPopTab(subviewer);CHKERRQ(ierr); 265 ierr = PetscViewerFlush(subviewer);CHKERRQ(ierr); 266 } 267 ierr = PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)pcbddc->ksp_D),&subviewer);CHKERRQ(ierr); 268 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 269 270 /* the coarse problem can be handled by a different communicator */ 271 if (pcbddc->coarse_ksp) color = 1; 272 else color = 0; 273 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRMPI(ierr); 274 ierr = PetscSubcommCreate(PetscObjectComm((PetscObject)pc),&subcomm);CHKERRQ(ierr); 275 ierr = PetscSubcommSetNumber(subcomm,PetscMin(size,2));CHKERRQ(ierr); 276 ierr = PetscSubcommSetTypeGeneral(subcomm,color,rank);CHKERRQ(ierr); 277 ierr = PetscViewerGetSubViewer(viewer,PetscSubcommChild(subcomm),&subviewer);CHKERRQ(ierr); 278 if (color == 1) { 279 ierr = PetscViewerASCIIPrintf(subviewer,"--- Coarse solver\n");CHKERRQ(ierr); 280 ierr = PetscViewerASCIIPushTab(subviewer);CHKERRQ(ierr); 281 ierr = KSPView(pcbddc->coarse_ksp,subviewer);CHKERRQ(ierr); 282 ierr = PetscViewerASCIIPopTab(subviewer);CHKERRQ(ierr); 283 ierr = PetscViewerFlush(subviewer);CHKERRQ(ierr); 284 } 285 ierr = PetscViewerRestoreSubViewer(viewer,PetscSubcommChild(subcomm),&subviewer);CHKERRQ(ierr); 286 ierr = PetscSubcommDestroy(&subcomm);CHKERRQ(ierr); 287 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 288 } 289 PetscFunctionReturn(0); 290 } 291 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; 295 PetscErrorCode ierr; 296 297 PetscFunctionBegin; 298 ierr = PetscObjectReference((PetscObject)G);CHKERRQ(ierr); 299 ierr = MatDestroy(&pcbddc->discretegradient);CHKERRQ(ierr); 300 pcbddc->discretegradient = G; 301 pcbddc->nedorder = order > 0 ? order : -order; 302 pcbddc->nedfield = field; 303 pcbddc->nedglobal = global; 304 pcbddc->conforming = conforming; 305 PetscFunctionReturn(0); 306 } 307 308 /*@ 309 PCBDDCSetDiscreteGradient - Sets the discrete gradient 310 311 Collective on PC 312 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 320 321 Level: advanced 322 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. 330 331 .seealso: PCBDDC,PCBDDCSetDofsSplitting(),PCBDDCSetDofsSplittingLocal() 332 @*/ 333 PetscErrorCode PCBDDCSetDiscreteGradient(PC pc, Mat G, PetscInt order, PetscInt field, PetscBool global, PetscBool conforming) 334 { 335 PetscErrorCode ierr; 336 337 PetscFunctionBegin; 338 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 339 PetscValidHeaderSpecific(G,MAT_CLASSID,2); 340 PetscValidLogicalCollectiveInt(pc,order,3); 341 PetscValidLogicalCollectiveInt(pc,field,4); 342 PetscValidLogicalCollectiveBool(pc,global,5); 343 PetscValidLogicalCollectiveBool(pc,conforming,6); 344 PetscCheckSameComm(pc,1,G,2); 345 ierr = PetscTryMethod(pc,"PCBDDCSetDiscreteGradient_C",(PC,Mat,PetscInt,PetscInt,PetscBool,PetscBool),(pc,G,order,field,global,conforming));CHKERRQ(ierr); 346 PetscFunctionReturn(0); 347 } 348 349 static PetscErrorCode PCBDDCSetDivergenceMat_BDDC(PC pc, Mat divudotp, PetscBool trans, IS vl2l) 350 { 351 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 352 PetscErrorCode ierr; 353 354 PetscFunctionBegin; 355 ierr = PetscObjectReference((PetscObject)divudotp);CHKERRQ(ierr); 356 ierr = MatDestroy(&pcbddc->divudotp);CHKERRQ(ierr); 357 pcbddc->divudotp = divudotp; 358 pcbddc->divudotp_trans = trans; 359 pcbddc->compute_nonetflux = PETSC_TRUE; 360 if (vl2l) { 361 ierr = PetscObjectReference((PetscObject)vl2l);CHKERRQ(ierr); 362 ierr = ISDestroy(&pcbddc->divudotp_vl2l);CHKERRQ(ierr); 363 pcbddc->divudotp_vl2l = vl2l; 364 } 365 PetscFunctionReturn(0); 366 } 367 368 /*@ 369 PCBDDCSetDivergenceMat - Sets the linear operator representing \int_\Omega \div {\bf u} \cdot p dx 370 371 Collective on PC 372 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 378 379 Level: advanced 380 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 384 385 .seealso: PCBDDC 386 @*/ 387 PetscErrorCode PCBDDCSetDivergenceMat(PC pc, Mat divudotp, PetscBool trans, IS vl2l) 388 { 389 PetscBool ismatis; 390 PetscErrorCode ierr; 391 392 PetscFunctionBegin; 393 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 394 PetscValidHeaderSpecific(divudotp,MAT_CLASSID,2); 395 PetscCheckSameComm(pc,1,divudotp,2); 396 PetscValidLogicalCollectiveBool(pc,trans,3); 397 if (vl2l) PetscValidHeaderSpecific(vl2l,IS_CLASSID,4); 398 ierr = PetscObjectTypeCompare((PetscObject)divudotp,MATIS,&ismatis);CHKERRQ(ierr); 399 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Divergence matrix needs to be of type MATIS"); 400 ierr = PetscTryMethod(pc,"PCBDDCSetDivergenceMat_C",(PC,Mat,PetscBool,IS),(pc,divudotp,trans,vl2l));CHKERRQ(ierr); 401 PetscFunctionReturn(0); 402 } 403 404 static PetscErrorCode PCBDDCSetChangeOfBasisMat_BDDC(PC pc, Mat change, PetscBool interior) 405 { 406 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 407 PetscErrorCode ierr; 408 409 PetscFunctionBegin; 410 ierr = PetscObjectReference((PetscObject)change);CHKERRQ(ierr); 411 ierr = MatDestroy(&pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 412 pcbddc->user_ChangeOfBasisMatrix = change; 413 pcbddc->change_interior = interior; 414 PetscFunctionReturn(0); 415 } 416 /*@ 417 PCBDDCSetChangeOfBasisMat - Set user defined change of basis for dofs 418 419 Collective on PC 420 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 425 426 Level: intermediate 427 428 Notes: 429 430 .seealso: PCBDDC 431 @*/ 432 PetscErrorCode PCBDDCSetChangeOfBasisMat(PC pc, Mat change, PetscBool interior) 433 { 434 PetscErrorCode ierr; 435 436 PetscFunctionBegin; 437 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 438 PetscValidHeaderSpecific(change,MAT_CLASSID,2); 439 PetscCheckSameComm(pc,1,change,2); 440 if (pc->mat) { 441 PetscInt rows_c,cols_c,rows,cols; 442 ierr = MatGetSize(pc->mat,&rows,&cols);CHKERRQ(ierr); 443 ierr = MatGetSize(change,&rows_c,&cols_c);CHKERRQ(ierr); 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 ierr = MatGetLocalSize(pc->mat,&rows,&cols);CHKERRQ(ierr); 447 ierr = MatGetLocalSize(change,&rows_c,&cols_c);CHKERRQ(ierr); 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 ierr = PetscTryMethod(pc,"PCBDDCSetChangeOfBasisMat_C",(PC,Mat,PetscBool),(pc,change,interior));CHKERRQ(ierr); 452 PetscFunctionReturn(0); 453 } 454 455 static PetscErrorCode PCBDDCSetPrimalVerticesIS_BDDC(PC pc, IS PrimalVertices) 456 { 457 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 458 PetscBool isequal = PETSC_FALSE; 459 PetscErrorCode ierr; 460 461 PetscFunctionBegin; 462 ierr = PetscObjectReference((PetscObject)PrimalVertices);CHKERRQ(ierr); 463 if (pcbddc->user_primal_vertices) { 464 ierr = ISEqual(PrimalVertices,pcbddc->user_primal_vertices,&isequal);CHKERRQ(ierr); 465 } 466 ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr); 467 ierr = ISDestroy(&pcbddc->user_primal_vertices_local);CHKERRQ(ierr); 468 pcbddc->user_primal_vertices = PrimalVertices; 469 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 470 PetscFunctionReturn(0); 471 } 472 473 /*@ 474 PCBDDCSetPrimalVerticesIS - Set additional user defined primal vertices in PCBDDC 475 476 Collective 477 478 Input Parameters: 479 + pc - the preconditioning context 480 - PrimalVertices - index set of primal vertices in global numbering (can be empty) 481 482 Level: intermediate 483 484 Notes: 485 Any process can list any global node 486 487 .seealso: PCBDDC, PCBDDCGetPrimalVerticesIS(), PCBDDCSetPrimalVerticesLocalIS(), PCBDDCGetPrimalVerticesLocalIS() 488 @*/ 489 PetscErrorCode PCBDDCSetPrimalVerticesIS(PC pc, IS PrimalVertices) 490 { 491 PetscErrorCode ierr; 492 493 PetscFunctionBegin; 494 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 495 PetscValidHeaderSpecific(PrimalVertices,IS_CLASSID,2); 496 PetscCheckSameComm(pc,1,PrimalVertices,2); 497 ierr = PetscTryMethod(pc,"PCBDDCSetPrimalVerticesIS_C",(PC,IS),(pc,PrimalVertices));CHKERRQ(ierr); 498 PetscFunctionReturn(0); 499 } 500 501 static PetscErrorCode PCBDDCGetPrimalVerticesIS_BDDC(PC pc, IS *is) 502 { 503 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 504 505 PetscFunctionBegin; 506 *is = pcbddc->user_primal_vertices; 507 PetscFunctionReturn(0); 508 } 509 510 /*@ 511 PCBDDCGetPrimalVerticesIS - Get user defined primal vertices set with PCBDDCSetPrimalVerticesIS() 512 513 Collective 514 515 Input Parameters: 516 . pc - the preconditioning context 517 518 Output Parameters: 519 . is - index set of primal vertices in global numbering (NULL if not set) 520 521 Level: intermediate 522 523 Notes: 524 525 .seealso: PCBDDC, PCBDDCSetPrimalVerticesIS(), PCBDDCSetPrimalVerticesLocalIS(), PCBDDCGetPrimalVerticesLocalIS() 526 @*/ 527 PetscErrorCode PCBDDCGetPrimalVerticesIS(PC pc, IS *is) 528 { 529 PetscErrorCode ierr; 530 531 PetscFunctionBegin; 532 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 533 PetscValidPointer(is,2); 534 ierr = PetscUseMethod(pc,"PCBDDCGetPrimalVerticesIS_C",(PC,IS*),(pc,is));CHKERRQ(ierr); 535 PetscFunctionReturn(0); 536 } 537 538 static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices) 539 { 540 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 541 PetscBool isequal = PETSC_FALSE; 542 PetscErrorCode ierr; 543 544 PetscFunctionBegin; 545 ierr = PetscObjectReference((PetscObject)PrimalVertices);CHKERRQ(ierr); 546 if (pcbddc->user_primal_vertices_local) { 547 ierr = ISEqual(PrimalVertices,pcbddc->user_primal_vertices_local,&isequal);CHKERRQ(ierr); 548 } 549 ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr); 550 ierr = ISDestroy(&pcbddc->user_primal_vertices_local);CHKERRQ(ierr); 551 pcbddc->user_primal_vertices_local = PrimalVertices; 552 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 553 PetscFunctionReturn(0); 554 } 555 556 /*@ 557 PCBDDCSetPrimalVerticesLocalIS - Set additional user defined primal vertices in PCBDDC 558 559 Collective 560 561 Input Parameters: 562 + pc - the preconditioning context 563 - PrimalVertices - index set of primal vertices in local numbering (can be empty) 564 565 Level: intermediate 566 567 Notes: 568 569 .seealso: PCBDDC, PCBDDCSetPrimalVerticesIS(), PCBDDCGetPrimalVerticesIS(), PCBDDCGetPrimalVerticesLocalIS() 570 @*/ 571 PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices) 572 { 573 PetscErrorCode ierr; 574 575 PetscFunctionBegin; 576 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 577 PetscValidHeaderSpecific(PrimalVertices,IS_CLASSID,2); 578 PetscCheckSameComm(pc,1,PrimalVertices,2); 579 ierr = PetscTryMethod(pc,"PCBDDCSetPrimalVerticesLocalIS_C",(PC,IS),(pc,PrimalVertices));CHKERRQ(ierr); 580 PetscFunctionReturn(0); 581 } 582 583 static PetscErrorCode PCBDDCGetPrimalVerticesLocalIS_BDDC(PC pc, IS *is) 584 { 585 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 586 587 PetscFunctionBegin; 588 *is = pcbddc->user_primal_vertices_local; 589 PetscFunctionReturn(0); 590 } 591 592 /*@ 593 PCBDDCGetPrimalVerticesLocalIS - Get user defined primal vertices set with PCBDDCSetPrimalVerticesLocalIS() 594 595 Collective 596 597 Input Parameters: 598 . pc - the preconditioning context 599 600 Output Parameters: 601 . is - index set of primal vertices in local numbering (NULL if not set) 602 603 Level: intermediate 604 605 Notes: 606 607 .seealso: PCBDDC, PCBDDCSetPrimalVerticesIS(), PCBDDCGetPrimalVerticesIS(), PCBDDCSetPrimalVerticesLocalIS() 608 @*/ 609 PetscErrorCode PCBDDCGetPrimalVerticesLocalIS(PC pc, IS *is) 610 { 611 PetscErrorCode ierr; 612 613 PetscFunctionBegin; 614 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 615 PetscValidPointer(is,2); 616 ierr = PetscUseMethod(pc,"PCBDDCGetPrimalVerticesLocalIS_C",(PC,IS*),(pc,is));CHKERRQ(ierr); 617 PetscFunctionReturn(0); 618 } 619 620 static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc,PetscInt k) 621 { 622 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 623 624 PetscFunctionBegin; 625 pcbddc->coarsening_ratio = k; 626 PetscFunctionReturn(0); 627 } 628 629 /*@ 630 PCBDDCSetCoarseningRatio - Set coarsening ratio used in multilevel 631 632 Logically collective on PC 633 634 Input Parameters: 635 + pc - the preconditioning context 636 - k - coarsening ratio (H/h at the coarser level) 637 638 Options Database Keys: 639 . -pc_bddc_coarsening_ratio 640 641 Level: intermediate 642 643 Notes: 644 Approximatively k subdomains at the finer level will be aggregated into a single subdomain at the coarser level 645 646 .seealso: PCBDDC, PCBDDCSetLevels() 647 @*/ 648 PetscErrorCode PCBDDCSetCoarseningRatio(PC pc,PetscInt k) 649 { 650 PetscErrorCode ierr; 651 652 PetscFunctionBegin; 653 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 654 PetscValidLogicalCollectiveInt(pc,k,2); 655 ierr = PetscTryMethod(pc,"PCBDDCSetCoarseningRatio_C",(PC,PetscInt),(pc,k));CHKERRQ(ierr); 656 PetscFunctionReturn(0); 657 } 658 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; 663 664 PetscFunctionBegin; 665 pcbddc->use_exact_dirichlet_trick = flg; 666 PetscFunctionReturn(0); 667 } 668 669 PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool flg) 670 { 671 PetscErrorCode ierr; 672 673 PetscFunctionBegin; 674 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 675 PetscValidLogicalCollectiveBool(pc,flg,2); 676 ierr = PetscTryMethod(pc,"PCBDDCSetUseExactDirichlet_C",(PC,PetscBool),(pc,flg));CHKERRQ(ierr); 677 PetscFunctionReturn(0); 678 } 679 680 static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc,PetscInt level) 681 { 682 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 683 684 PetscFunctionBegin; 685 pcbddc->current_level = level; 686 PetscFunctionReturn(0); 687 } 688 689 PetscErrorCode PCBDDCSetLevel(PC pc,PetscInt level) 690 { 691 PetscErrorCode ierr; 692 693 PetscFunctionBegin; 694 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 695 PetscValidLogicalCollectiveInt(pc,level,2); 696 ierr = PetscTryMethod(pc,"PCBDDCSetLevel_C",(PC,PetscInt),(pc,level));CHKERRQ(ierr); 697 PetscFunctionReturn(0); 698 } 699 700 static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc,PetscInt levels) 701 { 702 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 703 704 PetscFunctionBegin; 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 PetscFunctionReturn(0); 708 } 709 710 /*@ 711 PCBDDCSetLevels - Sets the maximum number of additional levels allowed for multilevel BDDC 712 713 Logically collective on PC 714 715 Input Parameters: 716 + pc - the preconditioning context 717 - levels - the maximum number of levels 718 719 Options Database Keys: 720 . -pc_bddc_levels 721 722 Level: intermediate 723 724 Notes: 725 The default value is 0, that gives the classical two-levels BDDC 726 727 .seealso: PCBDDC, PCBDDCSetCoarseningRatio() 728 @*/ 729 PetscErrorCode PCBDDCSetLevels(PC pc,PetscInt levels) 730 { 731 PetscErrorCode ierr; 732 733 PetscFunctionBegin; 734 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 735 PetscValidLogicalCollectiveInt(pc,levels,2); 736 ierr = PetscTryMethod(pc,"PCBDDCSetLevels_C",(PC,PetscInt),(pc,levels));CHKERRQ(ierr); 737 PetscFunctionReturn(0); 738 } 739 740 static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc,IS DirichletBoundaries) 741 { 742 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 743 PetscBool isequal = PETSC_FALSE; 744 PetscErrorCode ierr; 745 746 PetscFunctionBegin; 747 ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr); 748 if (pcbddc->DirichletBoundaries) { 749 ierr = ISEqual(DirichletBoundaries,pcbddc->DirichletBoundaries,&isequal);CHKERRQ(ierr); 750 } 751 /* last user setting takes precedence -> destroy any other customization */ 752 ierr = ISDestroy(&pcbddc->DirichletBoundariesLocal);CHKERRQ(ierr); 753 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 754 pcbddc->DirichletBoundaries = DirichletBoundaries; 755 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 756 PetscFunctionReturn(0); 757 } 758 759 /*@ 760 PCBDDCSetDirichletBoundaries - Set IS defining Dirichlet boundaries for the global problem. 761 762 Collective 763 764 Input Parameters: 765 + pc - the preconditioning context 766 - DirichletBoundaries - parallel IS defining the Dirichlet boundaries 767 768 Level: intermediate 769 770 Notes: 771 Provide the information if you used MatZeroRows/Columns routines. Any process can list any global node 772 773 .seealso: PCBDDC, PCBDDCSetDirichletBoundariesLocal(), MatZeroRows(), MatZeroRowsColumns() 774 @*/ 775 PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc,IS DirichletBoundaries) 776 { 777 PetscErrorCode ierr; 778 779 PetscFunctionBegin; 780 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 781 PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2); 782 PetscCheckSameComm(pc,1,DirichletBoundaries,2); 783 ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundaries_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr); 784 PetscFunctionReturn(0); 785 } 786 787 static PetscErrorCode PCBDDCSetDirichletBoundariesLocal_BDDC(PC pc,IS DirichletBoundaries) 788 { 789 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 790 PetscBool isequal = PETSC_FALSE; 791 PetscErrorCode ierr; 792 793 PetscFunctionBegin; 794 ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr); 795 if (pcbddc->DirichletBoundariesLocal) { 796 ierr = ISEqual(DirichletBoundaries,pcbddc->DirichletBoundariesLocal,&isequal);CHKERRQ(ierr); 797 } 798 /* last user setting takes precedence -> destroy any other customization */ 799 ierr = ISDestroy(&pcbddc->DirichletBoundariesLocal);CHKERRQ(ierr); 800 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 801 pcbddc->DirichletBoundariesLocal = DirichletBoundaries; 802 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 803 PetscFunctionReturn(0); 804 } 805 806 /*@ 807 PCBDDCSetDirichletBoundariesLocal - Set IS defining Dirichlet boundaries for the global problem in local ordering. 808 809 Collective 810 811 Input Parameters: 812 + pc - the preconditioning context 813 - DirichletBoundaries - parallel IS defining the Dirichlet boundaries (in local ordering) 814 815 Level: intermediate 816 817 Notes: 818 819 .seealso: PCBDDC, PCBDDCSetDirichletBoundaries(), MatZeroRows(), MatZeroRowsColumns() 820 @*/ 821 PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc,IS DirichletBoundaries) 822 { 823 PetscErrorCode ierr; 824 825 PetscFunctionBegin; 826 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 827 PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2); 828 PetscCheckSameComm(pc,1,DirichletBoundaries,2); 829 ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundariesLocal_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr); 830 PetscFunctionReturn(0); 831 } 832 833 static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc,IS NeumannBoundaries) 834 { 835 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 836 PetscBool isequal = PETSC_FALSE; 837 PetscErrorCode ierr; 838 839 PetscFunctionBegin; 840 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 841 if (pcbddc->NeumannBoundaries) { 842 ierr = ISEqual(NeumannBoundaries,pcbddc->NeumannBoundaries,&isequal);CHKERRQ(ierr); 843 } 844 /* last user setting takes precedence -> destroy any other customization */ 845 ierr = ISDestroy(&pcbddc->NeumannBoundariesLocal);CHKERRQ(ierr); 846 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 847 pcbddc->NeumannBoundaries = NeumannBoundaries; 848 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 849 PetscFunctionReturn(0); 850 } 851 852 /*@ 853 PCBDDCSetNeumannBoundaries - Set IS defining Neumann boundaries for the global problem. 854 855 Collective 856 857 Input Parameters: 858 + pc - the preconditioning context 859 - NeumannBoundaries - parallel IS defining the Neumann boundaries 860 861 Level: intermediate 862 863 Notes: 864 Any process can list any global node 865 866 .seealso: PCBDDC, PCBDDCSetNeumannBoundariesLocal() 867 @*/ 868 PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc,IS NeumannBoundaries) 869 { 870 PetscErrorCode ierr; 871 872 PetscFunctionBegin; 873 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 874 PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2); 875 PetscCheckSameComm(pc,1,NeumannBoundaries,2); 876 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundaries_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 877 PetscFunctionReturn(0); 878 } 879 880 static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc,IS NeumannBoundaries) 881 { 882 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 883 PetscBool isequal = PETSC_FALSE; 884 PetscErrorCode ierr; 885 886 PetscFunctionBegin; 887 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 888 if (pcbddc->NeumannBoundariesLocal) { 889 ierr = ISEqual(NeumannBoundaries,pcbddc->NeumannBoundariesLocal,&isequal);CHKERRQ(ierr); 890 } 891 /* last user setting takes precedence -> destroy any other customization */ 892 ierr = ISDestroy(&pcbddc->NeumannBoundariesLocal);CHKERRQ(ierr); 893 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 894 pcbddc->NeumannBoundariesLocal = NeumannBoundaries; 895 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 896 PetscFunctionReturn(0); 897 } 898 899 /*@ 900 PCBDDCSetNeumannBoundariesLocal - Set IS defining Neumann boundaries for the global problem in local ordering. 901 902 Collective 903 904 Input Parameters: 905 + pc - the preconditioning context 906 - NeumannBoundaries - parallel IS defining the subdomain part of Neumann boundaries (in local ordering) 907 908 Level: intermediate 909 910 Notes: 911 912 .seealso: PCBDDC, PCBDDCSetNeumannBoundaries() 913 @*/ 914 PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc,IS NeumannBoundaries) 915 { 916 PetscErrorCode ierr; 917 918 PetscFunctionBegin; 919 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 920 PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2); 921 PetscCheckSameComm(pc,1,NeumannBoundaries,2); 922 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundariesLocal_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 923 PetscFunctionReturn(0); 924 } 925 926 static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc,IS *DirichletBoundaries) 927 { 928 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 929 930 PetscFunctionBegin; 931 *DirichletBoundaries = pcbddc->DirichletBoundaries; 932 PetscFunctionReturn(0); 933 } 934 935 /*@ 936 PCBDDCGetDirichletBoundaries - Get parallel IS for Dirichlet boundaries 937 938 Collective 939 940 Input Parameters: 941 . pc - the preconditioning context 942 943 Output Parameters: 944 . DirichletBoundaries - index set defining the Dirichlet boundaries 945 946 Level: intermediate 947 948 Notes: 949 The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetDirichletBoundaries 950 951 .seealso: PCBDDC 952 @*/ 953 PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc,IS *DirichletBoundaries) 954 { 955 PetscErrorCode ierr; 956 957 PetscFunctionBegin; 958 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 959 ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundaries_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr); 960 PetscFunctionReturn(0); 961 } 962 963 static PetscErrorCode PCBDDCGetDirichletBoundariesLocal_BDDC(PC pc,IS *DirichletBoundaries) 964 { 965 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 966 967 PetscFunctionBegin; 968 *DirichletBoundaries = pcbddc->DirichletBoundariesLocal; 969 PetscFunctionReturn(0); 970 } 971 972 /*@ 973 PCBDDCGetDirichletBoundariesLocal - Get parallel IS for Dirichlet boundaries (in local ordering) 974 975 Collective 976 977 Input Parameters: 978 . pc - the preconditioning context 979 980 Output Parameters: 981 . DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries 982 983 Level: intermediate 984 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. 988 989 .seealso: PCBDDC 990 @*/ 991 PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc,IS *DirichletBoundaries) 992 { 993 PetscErrorCode ierr; 994 995 PetscFunctionBegin; 996 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 997 ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundariesLocal_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr); 998 PetscFunctionReturn(0); 999 } 1000 1001 static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc,IS *NeumannBoundaries) 1002 { 1003 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1004 1005 PetscFunctionBegin; 1006 *NeumannBoundaries = pcbddc->NeumannBoundaries; 1007 PetscFunctionReturn(0); 1008 } 1009 1010 /*@ 1011 PCBDDCGetNeumannBoundaries - Get parallel IS for Neumann boundaries 1012 1013 Collective 1014 1015 Input Parameters: 1016 . pc - the preconditioning context 1017 1018 Output Parameters: 1019 . NeumannBoundaries - index set defining the Neumann boundaries 1020 1021 Level: intermediate 1022 1023 Notes: 1024 The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetNeumannBoundaries 1025 1026 .seealso: PCBDDC 1027 @*/ 1028 PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc,IS *NeumannBoundaries) 1029 { 1030 PetscErrorCode ierr; 1031 1032 PetscFunctionBegin; 1033 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1034 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundaries_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 1035 PetscFunctionReturn(0); 1036 } 1037 1038 static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc,IS *NeumannBoundaries) 1039 { 1040 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1041 1042 PetscFunctionBegin; 1043 *NeumannBoundaries = pcbddc->NeumannBoundariesLocal; 1044 PetscFunctionReturn(0); 1045 } 1046 1047 /*@ 1048 PCBDDCGetNeumannBoundariesLocal - Get parallel IS for Neumann boundaries (in local ordering) 1049 1050 Collective 1051 1052 Input Parameters: 1053 . pc - the preconditioning context 1054 1055 Output Parameters: 1056 . NeumannBoundaries - index set defining the subdomain part of Neumann boundaries 1057 1058 Level: intermediate 1059 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. 1063 1064 .seealso: PCBDDC 1065 @*/ 1066 PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc,IS *NeumannBoundaries) 1067 { 1068 PetscErrorCode ierr; 1069 1070 PetscFunctionBegin; 1071 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1072 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundariesLocal_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 1073 PetscFunctionReturn(0); 1074 } 1075 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; 1081 PetscErrorCode ierr; 1082 1083 PetscFunctionBegin; 1084 if (!nvtxs) { 1085 if (copymode == PETSC_OWN_POINTER) { 1086 ierr = PetscFree(xadj);CHKERRQ(ierr); 1087 ierr = PetscFree(adjncy);CHKERRQ(ierr); 1088 } 1089 ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr); 1090 PetscFunctionReturn(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 ierr = PetscArraycmp(xadj,mat_graph->xadj,nvtxs+1,&same_data);CHKERRQ(ierr); 1096 if (same_data) { 1097 ierr = PetscArraycmp(adjncy,mat_graph->adjncy,xadj[nvtxs],&same_data);CHKERRQ(ierr); 1098 } 1099 } 1100 } 1101 if (!same_data) { 1102 /* free old CSR */ 1103 ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr); 1104 /* get CSR into graph structure */ 1105 if (copymode == PETSC_COPY_VALUES) { 1106 ierr = PetscMalloc1(nvtxs+1,&mat_graph->xadj);CHKERRQ(ierr); 1107 ierr = PetscMalloc1(xadj[nvtxs],&mat_graph->adjncy);CHKERRQ(ierr); 1108 ierr = PetscArraycpy(mat_graph->xadj,xadj,nvtxs+1);CHKERRQ(ierr); 1109 ierr = PetscArraycpy(mat_graph->adjncy,adjncy,xadj[nvtxs]);CHKERRQ(ierr); 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 PetscFunctionReturn(0); 1124 } 1125 1126 /*@ 1127 PCBDDCSetLocalAdjacencyGraph - Set adjacency structure (CSR graph) of the local degrees of freedom. 1128 1129 Not collective 1130 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. 1136 1137 Level: intermediate 1138 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. 1141 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; 1147 PetscErrorCode ierr; 1148 1149 PetscFunctionBegin; 1150 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1151 if (nvtxs) { 1152 PetscValidIntPointer(xadj,3); 1153 if (xadj[nvtxs]) PetscValidIntPointer(adjncy,4); 1154 } 1155 ierr = PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));CHKERRQ(ierr); 1156 /* free arrays if PCBDDC is not the PC type */ 1157 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);CHKERRQ(ierr); 1158 if (!f && copymode == PETSC_OWN_POINTER) { 1159 ierr = PetscFree(xadj);CHKERRQ(ierr); 1160 ierr = PetscFree(adjncy);CHKERRQ(ierr); 1161 } 1162 PetscFunctionReturn(0); 1163 } 1164 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; 1170 PetscErrorCode ierr; 1171 1172 PetscFunctionBegin; 1173 if (pcbddc->n_ISForDofsLocal == n_is) { 1174 for (i=0;i<n_is;i++) { 1175 PetscBool isequalt; 1176 ierr = ISEqual(ISForDofs[i],pcbddc->ISForDofsLocal[i],&isequalt);CHKERRQ(ierr); 1177 if (!isequalt) break; 1178 } 1179 if (i == n_is) isequal = PETSC_TRUE; 1180 } 1181 for (i=0;i<n_is;i++) { 1182 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 1183 } 1184 /* Destroy ISes if they were already set */ 1185 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 1186 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 1187 } 1188 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 1189 /* last user setting takes precedence -> destroy any other customization */ 1190 for (i=0;i<pcbddc->n_ISForDofs;i++) { 1191 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 1192 } 1193 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 1194 pcbddc->n_ISForDofs = 0; 1195 /* allocate space then set */ 1196 if (n_is) { 1197 ierr = PetscMalloc1(n_is,&pcbddc->ISForDofsLocal);CHKERRQ(ierr); 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 PetscFunctionReturn(0); 1206 } 1207 1208 /*@ 1209 PCBDDCSetDofsSplittingLocal - Set index sets defining fields of the local subdomain matrix 1210 1211 Collective 1212 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 1217 1218 Level: intermediate 1219 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. 1222 1223 .seealso: PCBDDC 1224 @*/ 1225 PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc,PetscInt n_is, IS ISForDofs[]) 1226 { 1227 PetscInt i; 1228 PetscErrorCode ierr; 1229 1230 PetscFunctionBegin; 1231 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1232 PetscValidLogicalCollectiveInt(pc,n_is,2); 1233 for (i=0;i<n_is;i++) { 1234 PetscCheckSameComm(pc,1,ISForDofs[i],3); 1235 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 1236 } 1237 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplittingLocal_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 1238 PetscFunctionReturn(0); 1239 } 1240 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; 1246 PetscErrorCode ierr; 1247 1248 PetscFunctionBegin; 1249 if (pcbddc->n_ISForDofs == n_is) { 1250 for (i=0;i<n_is;i++) { 1251 PetscBool isequalt; 1252 ierr = ISEqual(ISForDofs[i],pcbddc->ISForDofs[i],&isequalt);CHKERRQ(ierr); 1253 if (!isequalt) break; 1254 } 1255 if (i == n_is) isequal = PETSC_TRUE; 1256 } 1257 for (i=0;i<n_is;i++) { 1258 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 1259 } 1260 /* Destroy ISes if they were already set */ 1261 for (i=0;i<pcbddc->n_ISForDofs;i++) { 1262 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 1263 } 1264 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 1265 /* last user setting takes precedence -> destroy any other customization */ 1266 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 1267 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 1268 } 1269 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 1270 pcbddc->n_ISForDofsLocal = 0; 1271 /* allocate space then set */ 1272 if (n_is) { 1273 ierr = PetscMalloc1(n_is,&pcbddc->ISForDofs);CHKERRQ(ierr); 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 PetscFunctionReturn(0); 1282 } 1283 1284 /*@ 1285 PCBDDCSetDofsSplitting - Set index sets defining fields of the global matrix 1286 1287 Collective 1288 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 1293 1294 Level: intermediate 1295 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. 1298 1299 .seealso: PCBDDC 1300 @*/ 1301 PetscErrorCode PCBDDCSetDofsSplitting(PC pc,PetscInt n_is, IS ISForDofs[]) 1302 { 1303 PetscInt i; 1304 PetscErrorCode ierr; 1305 1306 PetscFunctionBegin; 1307 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1308 PetscValidLogicalCollectiveInt(pc,n_is,2); 1309 for (i=0;i<n_is;i++) { 1310 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 1311 PetscCheckSameComm(pc,1,ISForDofs[i],3); 1312 } 1313 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 1314 PetscFunctionReturn(0); 1315 } 1316 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. 1320 1321 Input Parameter: 1322 + pc - the preconditioner context 1323 1324 Application Interface Routine: PCPreSolve() 1325 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 { 1332 PetscErrorCode ierr; 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; 1337 1338 PetscFunctionBegin; 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; 1342 1343 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPCG,&iscg);CHKERRQ(ierr); 1344 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPGROPPCG,&isgroppcg);CHKERRQ(ierr); 1345 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECG,&ispipecg);CHKERRQ(ierr); 1346 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECG,&ispipelcg);CHKERRQ(ierr); 1347 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECGRR,&ispipecgrr);CHKERRQ(ierr); 1348 iscg = (PetscBool)(iscg || isgroppcg || ispipecg || ispipelcg || ispipecgrr); 1349 if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || !iscg || pc->mat != pc->pmat) { 1350 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 1351 } 1352 } 1353 if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || pc->mat != pc->pmat) { 1354 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 1355 } 1356 1357 /* Creates parallel work vectors used in presolve */ 1358 if (!pcbddc->original_rhs) { 1359 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 1360 } 1361 if (!pcbddc->temp_solution) { 1362 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);CHKERRQ(ierr); 1363 } 1364 1365 pcbddc->temp_solution_used = PETSC_FALSE; 1366 if (x) { 1367 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr); 1368 used_vec = x; 1369 } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */ 1370 ierr = PetscObjectReference((PetscObject)pcbddc->temp_solution);CHKERRQ(ierr); 1371 used_vec = pcbddc->temp_solution; 1372 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 1373 pcbddc->temp_solution_used = PETSC_TRUE; 1374 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1375 save_rhs = PETSC_FALSE; 1376 pcbddc->eliminate_dirdofs = PETSC_TRUE; 1377 } 1378 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 ierr = KSPGetInitialGuessNonzero(ksp,&pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 1383 if (!pcbddc->ksp_guess_nonzero) { 1384 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 1385 } 1386 } 1387 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; 1392 1393 /* DirichletBoundariesLocal may not be consistent among neighbours; gets a dirichlet dofs IS from graph (may be cached) */ 1394 ierr = PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph,&dirIS);CHKERRQ(ierr); 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; 1400 1401 ierr = MatGetDiagonal(pc->pmat,pcis->vec1_global);CHKERRQ(ierr); 1402 ierr = VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);CHKERRQ(ierr); 1403 ierr = VecScatterBegin(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1404 ierr = VecScatterEnd(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1405 ierr = VecScatterBegin(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1406 ierr = VecScatterEnd(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1407 ierr = ISGetLocalSize(dirIS,&dirsize);CHKERRQ(ierr); 1408 ierr = VecGetArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 1409 ierr = VecGetArrayRead(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 1410 ierr = ISGetIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1411 for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]]; 1412 ierr = ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1413 ierr = VecRestoreArrayRead(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 1414 ierr = VecRestoreArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 1415 ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1416 ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1417 pcbddc->rhs_change = PETSC_TRUE; 1418 ierr = ISDestroy(&dirIS);CHKERRQ(ierr); 1419 } 1420 } 1421 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 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1427 save_rhs = PETSC_FALSE; 1428 } 1429 pcbddc->rhs_change = PETSC_TRUE; 1430 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 1431 ierr = MatMultAdd(pc->mat,used_vec,pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1432 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 1433 ierr = VecCopy(used_vec,pcbddc->temp_solution);CHKERRQ(ierr); 1434 pcbddc->temp_solution_used = PETSC_TRUE; 1435 if (ksp) { 1436 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_FALSE);CHKERRQ(ierr); 1437 } 1438 } 1439 ierr = VecDestroy(&used_vec);CHKERRQ(ierr); 1440 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 ierr = VecDuplicate(rhs,&pcbddc->benign_vec);CHKERRQ(ierr); 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 ierr = PCApply_BDDC(pc,rhs,pcbddc->benign_vec);CHKERRQ(ierr); 1454 benign_correction_computed = PETSC_TRUE; 1455 if (pcbddc->temp_solution_used) { 1456 ierr = VecAXPY(pcbddc->temp_solution,1.0,pcbddc->benign_vec);CHKERRQ(ierr); 1457 } 1458 ierr = VecScale(pcbddc->benign_vec,-1.0);CHKERRQ(ierr); 1459 /* store the original rhs if not done earlier */ 1460 if (save_rhs) { 1461 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1462 } 1463 if (pcbddc->rhs_change) { 1464 ierr = MatMultAdd(pc->mat,pcbddc->benign_vec,rhs,rhs);CHKERRQ(ierr); 1465 } else { 1466 ierr = MatMultAdd(pc->mat,pcbddc->benign_vec,pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1467 } 1468 pcbddc->rhs_change = PETSC_TRUE; 1469 } 1470 pcbddc->benign_apply_coarse_only = PETSC_FALSE; 1471 } else { 1472 ierr = VecDestroy(&pcbddc->benign_vec);CHKERRQ(ierr); 1473 } 1474 1475 /* dbg output */ 1476 if (pcbddc->dbg_flag && benign_correction_computed) { 1477 Vec v; 1478 1479 ierr = VecDuplicate(pcis->vec1_global,&v);CHKERRQ(ierr); 1480 if (pcbddc->ChangeOfBasisMatrix) { 1481 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,rhs,v);CHKERRQ(ierr); 1482 } else { 1483 ierr = VecCopy(rhs,v);CHKERRQ(ierr); 1484 } 1485 ierr = PCBDDCBenignGetOrSetP0(pc,v,PETSC_TRUE);CHKERRQ(ierr); 1486 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"LEVEL %D: is the correction benign?\n",pcbddc->current_level);CHKERRQ(ierr); 1487 ierr = PetscScalarView(pcbddc->benign_n,pcbddc->benign_p0,pcbddc->dbg_viewer);CHKERRQ(ierr); 1488 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 1489 ierr = VecDestroy(&v);CHKERRQ(ierr); 1490 } 1491 1492 /* set initial guess if using PCG */ 1493 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 1494 if (x && pcbddc->use_exact_dirichlet_trick) { 1495 ierr = VecSet(x,0.0);CHKERRQ(ierr); 1496 if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) { 1497 if (benign_correction_computed) { /* we have already saved the changed rhs */ 1498 ierr = VecLockReadPop(pcis->vec1_global);CHKERRQ(ierr); 1499 } else { 1500 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,rhs,pcis->vec1_global);CHKERRQ(ierr); 1501 } 1502 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_global,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1503 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_global,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1504 } else { 1505 ierr = VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1506 ierr = VecScatterEnd(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1507 } 1508 ierr = PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 1509 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1510 ierr = PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 1511 ierr = KSPCheckSolve(pcbddc->ksp_D,pc,pcis->vec2_D);CHKERRQ(ierr); 1512 if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) { 1513 ierr = VecSet(pcis->vec1_global,0.);CHKERRQ(ierr); 1514 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1515 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1516 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcis->vec1_global,x);CHKERRQ(ierr); 1517 } else { 1518 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1519 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1520 } 1521 if (ksp) { 1522 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 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 ierr = VecLockReadPop(pcis->vec1_global);CHKERRQ(ierr); 1527 } 1528 PetscFunctionReturn(0); 1529 } 1530 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. 1534 1535 Input Parameter: 1536 + pc - the preconditioner context 1537 1538 Application Interface Routine: PCPostSolve() 1539 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 { 1546 PetscErrorCode ierr; 1547 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1548 1549 PetscFunctionBegin; 1550 /* add solution removed in presolve */ 1551 if (x && pcbddc->rhs_change) { 1552 if (pcbddc->temp_solution_used) { 1553 ierr = VecAXPY(x,1.0,pcbddc->temp_solution);CHKERRQ(ierr); 1554 } else if (pcbddc->benign_compute_correction && pcbddc->benign_vec) { 1555 ierr = VecAXPY(x,-1.0,pcbddc->benign_vec);CHKERRQ(ierr); 1556 } 1557 /* restore to original state (not for FETI-DP) */ 1558 if (ksp) pcbddc->temp_solution_used = PETSC_FALSE; 1559 } 1560 1561 /* restore rhs to its original state (not needed for FETI-DP) */ 1562 if (rhs && pcbddc->rhs_change) { 1563 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1564 pcbddc->rhs_change = PETSC_FALSE; 1565 } 1566 /* restore ksp guess state */ 1567 if (ksp) { 1568 ierr = KSPSetInitialGuessNonzero(ksp,pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 1569 /* reset flag for exact dirichlet trick */ 1570 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 1571 } 1572 PetscFunctionReturn(0); 1573 } 1574 1575 /* 1576 PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner 1577 by setting data structures and options. 1578 1579 Input Parameter: 1580 + pc - the preconditioner context 1581 1582 Application Interface Routine: PCSetUp() 1583 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; 1602 1603 PetscFunctionBegin; 1604 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&ismatis);CHKERRQ(ierr); 1605 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PCBDDC preconditioner requires matrix of type MATIS"); 1606 ierr = MatGetSize(pc->pmat,&nrows,&ncols);CHKERRQ(ierr); 1607 if (nrows != ncols) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCBDDC preconditioner requires a square preconditioning matrix"); 1608 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRMPI(ierr); 1609 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 ierr = MPIU_Allreduce(&rl,&pcbddc->recompute_topography,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)pc));CHKERRMPI(ierr); 1617 if (pcbddc->recompute_topography) { 1618 pcbddc->graphanalyzed = PETSC_FALSE; 1619 computeconstraintsmatrix = PETSC_TRUE; 1620 } else { 1621 computeconstraintsmatrix = PETSC_FALSE; 1622 } 1623 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; 1631 1632 computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling); 1633 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 } 1640 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 ierr = PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer);CHKERRQ(ierr); 1647 ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1648 } 1649 1650 /* process topology information */ 1651 ierr = PetscLogEventBegin(PC_BDDC_Topology[pcbddc->current_level],pc,0,0,0);CHKERRQ(ierr); 1652 if (pcbddc->recompute_topography) { 1653 ierr = PCBDDCComputeLocalTopologyInfo(pc);CHKERRQ(ierr); 1654 if (pcbddc->discretegradient) { 1655 ierr = PCBDDCNedelecSupport(pc);CHKERRQ(ierr); 1656 } 1657 } 1658 if (pcbddc->corner_selected) pcbddc->use_vertices = PETSC_TRUE; 1659 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 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 1665 } else { 1666 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1667 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1668 pcbddc->local_mat = matis->A; 1669 } 1670 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 ierr = PCBDDCBenignShellMat(pc,PETSC_TRUE);CHKERRQ(ierr); 1676 if (pcbddc->benign_saddle_point) { 1677 PC_IS* pcis = (PC_IS*)pc->data; 1678 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 ierr = PCBDDCBenignDetectSaddlePoint(pc,(PetscBool)(!pcbddc->recompute_topography),&zerodiag);CHKERRQ(ierr); 1682 /* pop B0 mat from local mat */ 1683 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 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 } 1695 1696 /* propagate relevant information */ 1697 if (matis->A->symmetric_set) { 1698 ierr = MatSetOption(pcbddc->local_mat,MAT_SYMMETRIC,matis->A->symmetric);CHKERRQ(ierr); 1699 } 1700 if (matis->A->spd_set) { 1701 ierr = MatSetOption(pcbddc->local_mat,MAT_SPD,matis->A->spd);CHKERRQ(ierr); 1702 } 1703 1704 /* Set up all the "iterative substructuring" common block without computing solvers */ 1705 { 1706 Mat temp_mat; 1707 1708 temp_mat = matis->A; 1709 matis->A = pcbddc->local_mat; 1710 ierr = PCISSetUp(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1711 pcbddc->local_mat = matis->A; 1712 matis->A = temp_mat; 1713 } 1714 1715 /* Analyze interface */ 1716 if (!pcbddc->graphanalyzed) { 1717 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 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; 1724 1725 if (!pcbddc->divudotp) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Missing divudotp operator"); 1726 ierr = PCBDDCComputeNoNetFlux(pc->pmat,pcbddc->divudotp,pcbddc->divudotp_trans,pcbddc->divudotp_vl2l,pcbddc->mat_graph,&nnfnnsp);CHKERRQ(ierr); 1727 /* TODO what if a nearnullspace is already attached? */ 1728 if (nnfnnsp) { 1729 ierr = MatSetNearNullSpace(pc->pmat,nnfnnsp);CHKERRQ(ierr); 1730 ierr = MatNullSpaceDestroy(&nnfnnsp);CHKERRQ(ierr); 1731 } 1732 } 1733 } 1734 ierr = PetscLogEventEnd(PC_BDDC_Topology[pcbddc->current_level],pc,0,0,0);CHKERRQ(ierr); 1735 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 ierr = PCBDDCBenignCheck(pc,zerodiag);CHKERRQ(ierr); 1741 } 1742 ierr = ISDestroy(&zerodiag);CHKERRQ(ierr); 1743 1744 /* Setup local dirichlet solver ksp_D and sub_schurs solvers */ 1745 if (computesubschurs && pcbddc->recompute_topography) { 1746 ierr = PCBDDCInitSubSchurs(pc);CHKERRQ(ierr); 1747 } 1748 /* SetUp Scaling operator (scaling matrices could be needed in SubSchursSetUp)*/ 1749 if (!pcbddc->use_deluxe_scaling) { 1750 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1751 } 1752 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 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1758 } 1759 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1760 } else { 1761 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1762 if (computesubschurs) { 1763 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1764 } 1765 } 1766 if (pcbddc->adaptive_selection) { 1767 ierr = PCBDDCAdaptiveSelection(pc);CHKERRQ(ierr); 1768 computeconstraintsmatrix = PETSC_TRUE; 1769 } 1770 1771 /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */ 1772 new_nearnullspace_provided = PETSC_FALSE; 1773 ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr); 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 ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 1787 for (i=0;i<nnsp_size;i++) { 1788 ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr); 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 } 1803 1804 /* Setup constraints and related work vectors */ 1805 /* reset primal space flags */ 1806 ierr = PetscLogEventBegin(PC_BDDC_LocalWork[pcbddc->current_level],pc,0,0,0);CHKERRQ(ierr); 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 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 1812 } 1813 /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */ 1814 ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr); 1815 1816 if (pcbddc->use_change_of_basis) { 1817 PC_IS *pcis = (PC_IS*)(pc->data); 1818 1819 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 1820 if (pcbddc->benign_change) { 1821 ierr = MatDestroy(&pcbddc->benign_B0);CHKERRQ(ierr); 1822 /* pop B0 from pcbddc->local_mat */ 1823 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1824 } 1825 /* get submatrices */ 1826 ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr); 1827 ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr); 1828 ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr); 1829 ierr = MatCreateSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);CHKERRQ(ierr); 1830 ierr = MatCreateSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);CHKERRQ(ierr); 1831 ierr = MatCreateSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);CHKERRQ(ierr); 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 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1836 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1837 pcbddc->local_mat = matis->A; 1838 } 1839 1840 /* interface pressure block row for B_C */ 1841 ierr = PetscObjectQuery((PetscObject)pc,"__KSPFETIDP_lP" ,(PetscObject*)&lP);CHKERRQ(ierr); 1842 ierr = PetscObjectQuery((PetscObject)pc,"__KSPFETIDP_lA" ,(PetscObject*)&lA);CHKERRQ(ierr); 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 ierr = MatIsSymmetric(lA,PETSC_SMALL,&issym);CHKERRQ(ierr); 1848 if (issym) { 1849 ierr = MatCreateSubMatrix(lA,lP,pcis->is_I_local,MAT_INITIAL_MATRIX,&B_BI);CHKERRQ(ierr); 1850 ierr = MatCreateSubMatrix(lA,lP,pcis->is_B_local,MAT_INITIAL_MATRIX,&B_BB);CHKERRQ(ierr); 1851 ierr = MatCreateTranspose(B_BI,&Bt_BI);CHKERRQ(ierr); 1852 ierr = MatCreateTranspose(B_BB,&Bt_BB);CHKERRQ(ierr); 1853 } else { 1854 ierr = MatCreateSubMatrix(lA,lP,pcis->is_I_local,MAT_INITIAL_MATRIX,&B_BI);CHKERRQ(ierr); 1855 ierr = MatCreateSubMatrix(lA,lP,pcis->is_B_local,MAT_INITIAL_MATRIX,&B_BB);CHKERRQ(ierr); 1856 ierr = MatCreateSubMatrix(lA,pcis->is_I_local,lP,MAT_INITIAL_MATRIX,&Bt_BI);CHKERRQ(ierr); 1857 ierr = MatCreateSubMatrix(lA,pcis->is_B_local,lP,MAT_INITIAL_MATRIX,&Bt_BB);CHKERRQ(ierr); 1858 } 1859 ierr = PetscObjectCompose((PetscObject)pc,"__KSPFETIDP_B_BI",(PetscObject)B_BI);CHKERRQ(ierr); 1860 ierr = PetscObjectCompose((PetscObject)pc,"__KSPFETIDP_B_BB",(PetscObject)B_BB);CHKERRQ(ierr); 1861 ierr = PetscObjectCompose((PetscObject)pc,"__KSPFETIDP_Bt_BI",(PetscObject)Bt_BI);CHKERRQ(ierr); 1862 ierr = PetscObjectCompose((PetscObject)pc,"__KSPFETIDP_Bt_BB",(PetscObject)Bt_BB);CHKERRQ(ierr); 1863 ierr = MatDestroy(&B_BI);CHKERRQ(ierr); 1864 ierr = MatDestroy(&B_BB);CHKERRQ(ierr); 1865 ierr = MatDestroy(&Bt_BI);CHKERRQ(ierr); 1866 ierr = MatDestroy(&Bt_BB);CHKERRQ(ierr); 1867 } 1868 ierr = PetscLogEventEnd(PC_BDDC_LocalWork[pcbddc->current_level],pc,0,0,0);CHKERRQ(ierr); 1869 1870 /* SetUp coarse and local Neumann solvers */ 1871 ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr); 1872 /* SetUp Scaling operator */ 1873 if (pcbddc->use_deluxe_scaling) { 1874 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1875 } 1876 1877 /* mark topography as done */ 1878 pcbddc->recompute_topography = PETSC_FALSE; 1879 1880 /* wrap pcis->A_IB and pcis->A_BI if we did not change explicitly the variables on the pressures */ 1881 ierr = PCBDDCBenignShellMat(pc,PETSC_FALSE);CHKERRQ(ierr); 1882 1883 if (pcbddc->dbg_flag) { 1884 ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1885 ierr = PetscViewerASCIIPopSynchronized(pcbddc->dbg_viewer);CHKERRQ(ierr); 1886 } 1887 PetscFunctionReturn(0); 1888 } 1889 1890 /* 1891 PCApply_BDDC - Applies the BDDC operator to a vector. 1892 1893 Input Parameters: 1894 + pc - the preconditioner context 1895 - r - input vector (global) 1896 1897 Output Parameter: 1898 . z - output vector (global) 1899 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) */ 1915 1916 PetscFunctionBegin; 1917 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 1918 if (pcbddc->switch_static) { 1919 ierr = MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat,&lA);CHKERRQ(ierr); 1920 } 1921 1922 if (pcbddc->ChangeOfBasisMatrix) { 1923 Vec swap; 1924 1925 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 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 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 1932 ierr = VecLockReadPush(pcis->vec1_global);CHKERRQ(ierr); 1933 } 1934 } 1935 if (pcbddc->benign_have_null) { /* get p0 from r */ 1936 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 1937 } 1938 if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_DIRICHLET && !pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 1939 ierr = VecCopy(r,z);CHKERRQ(ierr); 1940 /* First Dirichlet solve */ 1941 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1942 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 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 ierr = PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 1950 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1951 ierr = PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 1952 ierr = KSPCheckSolve(pcbddc->ksp_D,pc,pcis->vec2_D);CHKERRQ(ierr); 1953 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1954 if (pcbddc->switch_static) { 1955 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 1956 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1957 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1958 if (!pcbddc->switch_static_change) { 1959 ierr = MatMult(lA,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1960 } else { 1961 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1962 ierr = MatMult(lA,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1963 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1964 } 1965 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1966 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1967 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1968 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1969 } else { 1970 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1971 } 1972 } else { 1973 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1974 } 1975 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1976 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1977 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1978 } else { 1979 if (!pcbddc->benign_apply_coarse_only) { 1980 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 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 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1986 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1987 } 1988 1989 /* Apply interface preconditioner 1990 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1991 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr); 1992 1993 /* Apply transpose of partition of unity operator */ 1994 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1995 if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_LUMP) { 1996 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1997 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1998 PetscFunctionReturn(0); 1999 } 2000 /* Second Dirichlet solve and assembling of output */ 2001 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2002 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2003 if (n_B) { 2004 if (pcbddc->switch_static) { 2005 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2006 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2007 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2008 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2009 if (!pcbddc->switch_static_change) { 2010 ierr = MatMult(lA,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2011 } else { 2012 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2013 ierr = MatMult(lA,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 2014 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2015 } 2016 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2017 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2018 } else { 2019 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 2020 } 2021 } else if (pcbddc->switch_static) { /* n_B is zero */ 2022 if (!pcbddc->switch_static_change) { 2023 ierr = MatMult(lA,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 2024 } else { 2025 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 2026 ierr = MatMult(lA,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2027 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 2028 } 2029 } 2030 ierr = PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 2031 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 2032 ierr = PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 2033 ierr = KSPCheckSolve(pcbddc->ksp_D,pc,pcis->vec4_D);CHKERRQ(ierr); 2034 2035 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 2036 if (pcbddc->switch_static) { 2037 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 2038 } else { 2039 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 2040 } 2041 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2042 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2043 } else { 2044 if (pcbddc->switch_static) { 2045 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 2046 } else { 2047 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 2048 } 2049 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2050 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2051 } 2052 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 2053 if (pcbddc->benign_apply_coarse_only) { 2054 ierr = PetscArrayzero(pcbddc->benign_p0,pcbddc->benign_n);CHKERRQ(ierr); 2055 } 2056 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 2057 } 2058 2059 if (pcbddc->ChangeOfBasisMatrix) { 2060 pcbddc->work_change = r; 2061 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 2062 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 2063 } 2064 PetscFunctionReturn(0); 2065 } 2066 2067 /* 2068 PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector. 2069 2070 Input Parameters: 2071 + pc - the preconditioner context 2072 - r - input vector (global) 2073 2074 Output Parameter: 2075 . z - output vector (global) 2076 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; 2089 2090 PetscFunctionBegin; 2091 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 2092 if (pcbddc->switch_static) { 2093 ierr = MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat,&lA);CHKERRQ(ierr); 2094 } 2095 if (pcbddc->ChangeOfBasisMatrix) { 2096 Vec swap; 2097 2098 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 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 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 2105 ierr = VecLockReadPush(pcis->vec1_global);CHKERRQ(ierr); 2106 } 2107 } 2108 if (pcbddc->benign_have_null) { /* get p0 from r */ 2109 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 2110 } 2111 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 2112 ierr = VecCopy(r,z);CHKERRQ(ierr); 2113 /* First Dirichlet solve */ 2114 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2115 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 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 ierr = PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 2123 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 2124 ierr = PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 2125 ierr = KSPCheckSolve(pcbddc->ksp_D,pc,pcis->vec2_D);CHKERRQ(ierr); 2126 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 2127 if (pcbddc->switch_static) { 2128 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 2129 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2130 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2131 if (!pcbddc->switch_static_change) { 2132 ierr = MatMultTranspose(lA,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2133 } else { 2134 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2135 ierr = MatMultTranspose(lA,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 2136 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2137 } 2138 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2139 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2140 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2141 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2142 } else { 2143 ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 2144 } 2145 } else { 2146 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 2147 } 2148 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2149 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2150 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 2151 } else { 2152 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 2153 } 2154 2155 /* Apply interface preconditioner 2156 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 2157 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr); 2158 2159 /* Apply transpose of partition of unity operator */ 2160 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 2161 2162 /* Second Dirichlet solve and assembling of output */ 2163 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2164 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2165 if (n_B) { 2166 if (pcbddc->switch_static) { 2167 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2168 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2169 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2170 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2171 if (!pcbddc->switch_static_change) { 2172 ierr = MatMultTranspose(lA,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2173 } else { 2174 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2175 ierr = MatMultTranspose(lA,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 2176 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2177 } 2178 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2179 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2180 } else { 2181 ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 2182 } 2183 } else if (pcbddc->switch_static) { /* n_B is zero */ 2184 if (!pcbddc->switch_static_change) { 2185 ierr = MatMultTranspose(lA,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 2186 } else { 2187 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 2188 ierr = MatMultTranspose(lA,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2189 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 2190 } 2191 } 2192 ierr = PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 2193 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 2194 ierr = PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 2195 ierr = KSPCheckSolve(pcbddc->ksp_D,pc,pcis->vec4_D);CHKERRQ(ierr); 2196 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 2197 if (pcbddc->switch_static) { 2198 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 2199 } else { 2200 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 2201 } 2202 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2203 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2204 } else { 2205 if (pcbddc->switch_static) { 2206 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 2207 } else { 2208 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 2209 } 2210 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2211 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2212 } 2213 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 2214 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 2215 } 2216 if (pcbddc->ChangeOfBasisMatrix) { 2217 pcbddc->work_change = r; 2218 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 2219 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 2220 } 2221 PetscFunctionReturn(0); 2222 } 2223 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; 2229 PetscErrorCode ierr; 2230 2231 PetscFunctionBegin; 2232 /* free BDDC custom data */ 2233 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 2234 /* destroy objects related to topography */ 2235 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 2236 /* destroy objects for scaling operator */ 2237 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 2238 /* free solvers stuff */ 2239 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 2240 /* free global vectors needed in presolve */ 2241 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 2242 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 2243 /* free data created by PCIS */ 2244 ierr = PCISDestroy(pc);CHKERRQ(ierr); 2245 2246 /* restore defaults */ 2247 kspD = pcbddc->ksp_D; 2248 kspR = pcbddc->ksp_R; 2249 kspC = pcbddc->coarse_ksp; 2250 ierr = PetscMemzero(pc->data,sizeof(*pcbddc));CHKERRQ(ierr); 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 PetscFunctionReturn(0); 2273 } 2274 2275 PetscErrorCode PCDestroy_BDDC(PC pc) 2276 { 2277 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 2278 PetscErrorCode ierr; 2279 2280 PetscFunctionBegin; 2281 ierr = PCReset_BDDC(pc);CHKERRQ(ierr); 2282 ierr = KSPDestroy(&pcbddc->ksp_D);CHKERRQ(ierr); 2283 ierr = KSPDestroy(&pcbddc->ksp_R);CHKERRQ(ierr); 2284 ierr = KSPDestroy(&pcbddc->coarse_ksp);CHKERRQ(ierr); 2285 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDiscreteGradient_C",NULL);CHKERRQ(ierr); 2286 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",NULL);CHKERRQ(ierr); 2287 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);CHKERRQ(ierr); 2288 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 2289 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",NULL);CHKERRQ(ierr); 2290 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 2291 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr); 2292 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr); 2293 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr); 2294 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2295 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2296 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2297 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2298 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2299 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2300 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2301 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2302 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 2303 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr); 2304 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 2305 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 2306 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 2307 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 2308 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",NULL);CHKERRQ(ierr); 2309 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL);CHKERRQ(ierr); 2310 ierr = PetscFree(pc->data);CHKERRQ(ierr); 2311 PetscFunctionReturn(0); 2312 } 2313 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; 2318 PetscErrorCode ierr; 2319 2320 PetscFunctionBegin; 2321 ierr = PetscFree(mat_graph->coords);CHKERRQ(ierr); 2322 ierr = PetscMalloc1(nloc*dim,&mat_graph->coords);CHKERRQ(ierr); 2323 ierr = PetscArraycpy(mat_graph->coords,coords,nloc*dim);CHKERRQ(ierr); 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 PetscFunctionReturn(0); 2331 } 2332 2333 static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool* change) 2334 { 2335 PetscFunctionBegin; 2336 *change = PETSC_TRUE; 2337 PetscFunctionReturn(0); 2338 } 2339 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; 2346 PetscErrorCode ierr; 2347 2348 PetscFunctionBegin; 2349 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2350 pcis = (PC_IS*)mat_ctx->pc->data; 2351 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2352 2353 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 2354 /* copy rhs since we may change it during PCPreSolve_BDDC */ 2355 if (!pcbddc->original_rhs) { 2356 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 2357 } 2358 if (mat_ctx->rhs_flip) { 2359 ierr = VecPointwiseMult(pcbddc->original_rhs,standard_rhs,mat_ctx->rhs_flip);CHKERRQ(ierr); 2360 } else { 2361 ierr = VecCopy(standard_rhs,pcbddc->original_rhs);CHKERRQ(ierr); 2362 } 2363 if (mat_ctx->g2g_p) { 2364 /* interface pressure rhs */ 2365 ierr = VecScatterBegin(mat_ctx->g2g_p,fetidp_flux_rhs,pcbddc->original_rhs,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2366 ierr = VecScatterEnd(mat_ctx->g2g_p,fetidp_flux_rhs,pcbddc->original_rhs,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2367 ierr = VecScatterBegin(mat_ctx->g2g_p,standard_rhs,fetidp_flux_rhs,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2368 ierr = VecScatterEnd(mat_ctx->g2g_p,standard_rhs,fetidp_flux_rhs,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2369 if (!mat_ctx->rhs_flip) { 2370 ierr = VecScale(fetidp_flux_rhs,-1.);CHKERRQ(ierr); 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 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,pcbddc->original_rhs,NULL);CHKERRQ(ierr); 2378 if (pcbddc->ChangeOfBasisMatrix) { 2379 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,pcbddc->original_rhs,pcbddc->work_change);CHKERRQ(ierr); 2380 work = pcbddc->work_change; 2381 } else { 2382 work = pcbddc->original_rhs; 2383 } 2384 /* store vectors for computation of fetidp final solution */ 2385 ierr = VecScatterBegin(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2386 ierr = VecScatterEnd(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2387 /* scale rhs since it should be unassembled */ 2388 /* TODO use counter scaling? (also below) */ 2389 ierr = VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2390 ierr = VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2391 /* Apply partition of unity */ 2392 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2393 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2394 if (!pcbddc->switch_static) { 2395 /* compute partially subassembled Schur complement right-hand side */ 2396 ierr = PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],mat_ctx->pc,0,0,0);CHKERRQ(ierr); 2397 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2398 ierr = PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],mat_ctx->pc,0,0,0);CHKERRQ(ierr); 2399 /* Cannot propagate up error in KSPSolve() because there is no access to the PC */ 2400 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 2401 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 2402 ierr = VecSet(work,0.0);CHKERRQ(ierr); 2403 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2404 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2405 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2406 ierr = VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2407 ierr = VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2408 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2409 } 2410 /* BDDC rhs */ 2411 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 2412 if (pcbddc->switch_static) { 2413 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2414 } 2415 /* apply BDDC */ 2416 ierr = PetscArrayzero(pcbddc->benign_p0,pcbddc->benign_n);CHKERRQ(ierr); 2417 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2418 ierr = PetscArrayzero(pcbddc->benign_p0,pcbddc->benign_n);CHKERRQ(ierr); 2419 2420 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 2421 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 2422 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2423 ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2424 /* Add contribution to interface pressures */ 2425 if (mat_ctx->l2g_p) { 2426 ierr = MatMult(mat_ctx->B_BB,pcis->vec1_B,mat_ctx->vP);CHKERRQ(ierr); 2427 if (pcbddc->switch_static) { 2428 ierr = MatMultAdd(mat_ctx->B_BI,pcis->vec1_D,mat_ctx->vP,mat_ctx->vP);CHKERRQ(ierr); 2429 } 2430 ierr = VecScatterBegin(mat_ctx->l2g_p,mat_ctx->vP,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2431 ierr = VecScatterEnd(mat_ctx->l2g_p,mat_ctx->vP,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2432 } 2433 PetscFunctionReturn(0); 2434 } 2435 2436 /*@ 2437 PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETI-DP linear system using the physical right-hand side 2438 2439 Collective 2440 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 2444 2445 Output Parameters: 2446 . fetidp_flux_rhs - the right-hand side for the FETI-DP linear system 2447 2448 Level: developer 2449 2450 Notes: 2451 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; 2457 PetscErrorCode ierr; 2458 2459 PetscFunctionBegin; 2460 PetscValidHeaderSpecific(fetidp_mat,MAT_CLASSID,1); 2461 PetscValidHeaderSpecific(standard_rhs,VEC_CLASSID,2); 2462 PetscValidHeaderSpecific(fetidp_flux_rhs,VEC_CLASSID,3); 2463 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2464 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 2465 PetscFunctionReturn(0); 2466 } 2467 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; 2473 PetscErrorCode ierr; 2474 Vec work; 2475 2476 PetscFunctionBegin; 2477 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2478 pcis = (PC_IS*)mat_ctx->pc->data; 2479 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2480 2481 /* apply B_delta^T */ 2482 ierr = VecSet(pcis->vec1_B,0.);CHKERRQ(ierr); 2483 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2484 ierr = VecScatterEnd(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2485 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 2486 if (mat_ctx->l2g_p) { 2487 ierr = VecScatterBegin(mat_ctx->l2g_p,fetidp_flux_sol,mat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2488 ierr = VecScatterEnd(mat_ctx->l2g_p,fetidp_flux_sol,mat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2489 ierr = MatMultAdd(mat_ctx->Bt_BB,mat_ctx->vP,pcis->vec1_B,pcis->vec1_B);CHKERRQ(ierr); 2490 } 2491 2492 /* compute rhs for BDDC application */ 2493 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2494 if (pcbddc->switch_static) { 2495 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2496 if (mat_ctx->l2g_p) { 2497 ierr = VecScale(mat_ctx->vP,-1.);CHKERRQ(ierr); 2498 ierr = MatMultAdd(mat_ctx->Bt_BI,mat_ctx->vP,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); 2499 } 2500 } 2501 2502 /* apply BDDC */ 2503 ierr = PetscArrayzero(pcbddc->benign_p0,pcbddc->benign_n);CHKERRQ(ierr); 2504 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2505 2506 /* put values into global vector */ 2507 if (pcbddc->ChangeOfBasisMatrix) work = pcbddc->work_change; 2508 else work = standard_sol; 2509 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2510 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2511 if (!pcbddc->switch_static) { 2512 /* compute values into the interior if solved for the partially subassembled Schur complement */ 2513 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 2514 ierr = VecAYPX(pcis->vec1_D,-1.0,mat_ctx->temp_solution_D);CHKERRQ(ierr); 2515 ierr = PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],mat_ctx->pc,0,0,0);CHKERRQ(ierr); 2516 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); 2517 ierr = PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],mat_ctx->pc,0,0,0);CHKERRQ(ierr); 2518 /* Cannot propagate up error in KSPSolve() because there is no access to the PC */ 2519 } 2520 2521 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2522 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_D,work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2523 /* add p0 solution to final solution */ 2524 ierr = PCBDDCBenignGetOrSetP0(mat_ctx->pc,work,PETSC_FALSE);CHKERRQ(ierr); 2525 if (pcbddc->ChangeOfBasisMatrix) { 2526 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,work,standard_sol);CHKERRQ(ierr); 2527 } 2528 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 2529 if (mat_ctx->g2g_p) { 2530 ierr = VecScatterBegin(mat_ctx->g2g_p,fetidp_flux_sol,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2531 ierr = VecScatterEnd(mat_ctx->g2g_p,fetidp_flux_sol,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2532 } 2533 PetscFunctionReturn(0); 2534 } 2535 2536 static PetscErrorCode PCView_BDDCIPC(PC pc, PetscViewer viewer) 2537 { 2538 PetscErrorCode ierr; 2539 BDDCIPC_ctx bddcipc_ctx; 2540 PetscBool isascii; 2541 2542 PetscFunctionBegin; 2543 ierr = PCShellGetContext(pc,&bddcipc_ctx);CHKERRQ(ierr); 2544 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 2545 if (isascii) { 2546 ierr = PetscViewerASCIIPrintf(viewer,"BDDC interface preconditioner\n");CHKERRQ(ierr); 2547 } 2548 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 2549 ierr = PCView(bddcipc_ctx->bddc,viewer);CHKERRQ(ierr); 2550 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 2551 PetscFunctionReturn(0); 2552 } 2553 2554 static PetscErrorCode PCSetUp_BDDCIPC(PC pc) 2555 { 2556 PetscErrorCode ierr; 2557 BDDCIPC_ctx bddcipc_ctx; 2558 PetscBool isbddc; 2559 Vec vv; 2560 IS is; 2561 PC_IS *pcis; 2562 2563 PetscFunctionBegin; 2564 ierr = PCShellGetContext(pc,&bddcipc_ctx);CHKERRQ(ierr); 2565 ierr = PetscObjectTypeCompare((PetscObject)bddcipc_ctx->bddc,PCBDDC,&isbddc);CHKERRQ(ierr); 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 ierr = PCSetUp(bddcipc_ctx->bddc);CHKERRQ(ierr); 2568 2569 /* create interface scatter */ 2570 pcis = (PC_IS*)(bddcipc_ctx->bddc->data); 2571 ierr = VecScatterDestroy(&bddcipc_ctx->g2l);CHKERRQ(ierr); 2572 ierr = MatCreateVecs(pc->pmat,&vv,NULL);CHKERRQ(ierr); 2573 ierr = ISRenumber(pcis->is_B_global,NULL,NULL,&is);CHKERRQ(ierr); 2574 ierr = VecScatterCreate(vv,is,pcis->vec1_B,NULL,&bddcipc_ctx->g2l);CHKERRQ(ierr); 2575 ierr = ISDestroy(&is);CHKERRQ(ierr); 2576 ierr = VecDestroy(&vv);CHKERRQ(ierr); 2577 PetscFunctionReturn(0); 2578 } 2579 2580 static PetscErrorCode PCApply_BDDCIPC(PC pc, Vec r, Vec x) 2581 { 2582 PetscErrorCode ierr; 2583 BDDCIPC_ctx bddcipc_ctx; 2584 PC_IS *pcis; 2585 VecScatter tmps; 2586 2587 PetscFunctionBegin; 2588 ierr = PCShellGetContext(pc,&bddcipc_ctx);CHKERRQ(ierr); 2589 pcis = (PC_IS*)(bddcipc_ctx->bddc->data); 2590 tmps = pcis->global_to_B; 2591 pcis->global_to_B = bddcipc_ctx->g2l; 2592 ierr = PCBDDCScalingRestriction(bddcipc_ctx->bddc,r,pcis->vec1_B);CHKERRQ(ierr); 2593 ierr = PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc,PETSC_FALSE);CHKERRQ(ierr); 2594 ierr = PCBDDCScalingExtension(bddcipc_ctx->bddc,pcis->vec1_B,x);CHKERRQ(ierr); 2595 pcis->global_to_B = tmps; 2596 PetscFunctionReturn(0); 2597 } 2598 2599 static PetscErrorCode PCApplyTranspose_BDDCIPC(PC pc, Vec r, Vec x) 2600 { 2601 PetscErrorCode ierr; 2602 BDDCIPC_ctx bddcipc_ctx; 2603 PC_IS *pcis; 2604 VecScatter tmps; 2605 2606 PetscFunctionBegin; 2607 ierr = PCShellGetContext(pc,&bddcipc_ctx);CHKERRQ(ierr); 2608 pcis = (PC_IS*)(bddcipc_ctx->bddc->data); 2609 tmps = pcis->global_to_B; 2610 pcis->global_to_B = bddcipc_ctx->g2l; 2611 ierr = PCBDDCScalingRestriction(bddcipc_ctx->bddc,r,pcis->vec1_B);CHKERRQ(ierr); 2612 ierr = PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc,PETSC_TRUE);CHKERRQ(ierr); 2613 ierr = PCBDDCScalingExtension(bddcipc_ctx->bddc,pcis->vec1_B,x);CHKERRQ(ierr); 2614 pcis->global_to_B = tmps; 2615 PetscFunctionReturn(0); 2616 } 2617 2618 static PetscErrorCode PCDestroy_BDDCIPC(PC pc) 2619 { 2620 PetscErrorCode ierr; 2621 BDDCIPC_ctx bddcipc_ctx; 2622 2623 PetscFunctionBegin; 2624 ierr = PCShellGetContext(pc,&bddcipc_ctx);CHKERRQ(ierr); 2625 ierr = PCDestroy(&bddcipc_ctx->bddc);CHKERRQ(ierr); 2626 ierr = VecScatterDestroy(&bddcipc_ctx->g2l);CHKERRQ(ierr); 2627 ierr = PetscFree(bddcipc_ctx);CHKERRQ(ierr); 2628 PetscFunctionReturn(0); 2629 } 2630 2631 /*@ 2632 PCBDDCMatFETIDPGetSolution - Compute the physical solution using the solution of the FETI-DP linear system 2633 2634 Collective 2635 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 2639 2640 Output Parameters: 2641 . standard_sol - the solution defined on the physical domain 2642 2643 Level: developer 2644 2645 Notes: 2646 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; 2652 PetscErrorCode ierr; 2653 2654 PetscFunctionBegin; 2655 PetscValidHeaderSpecific(fetidp_mat,MAT_CLASSID,1); 2656 PetscValidHeaderSpecific(fetidp_flux_sol,VEC_CLASSID,2); 2657 PetscValidHeaderSpecific(standard_sol,VEC_CLASSID,3); 2658 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2659 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 2660 PetscFunctionReturn(0); 2661 } 2662 2663 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, PetscBool fully_redundant, const char* prefix, Mat *fetidp_mat, PC *fetidp_pc) 2664 { 2665 2666 FETIDPMat_ctx fetidpmat_ctx; 2667 Mat newmat; 2668 FETIDPPC_ctx fetidppc_ctx; 2669 PC newpc; 2670 MPI_Comm comm; 2671 PetscErrorCode ierr; 2672 2673 PetscFunctionBegin; 2674 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 2675 /* FETI-DP matrix */ 2676 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 2677 fetidpmat_ctx->fully_redundant = fully_redundant; 2678 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 2679 ierr = MatCreateShell(comm,fetidpmat_ctx->n,fetidpmat_ctx->n,fetidpmat_ctx->N,fetidpmat_ctx->N,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 2680 ierr = PetscObjectSetName((PetscObject)newmat,!fetidpmat_ctx->l2g_lambda_only ? "F" : "G");CHKERRQ(ierr); 2681 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 2682 ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr); 2683 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 2684 /* propagate MatOptions */ 2685 { 2686 PC_BDDC *pcbddc = (PC_BDDC*)fetidpmat_ctx->pc->data; 2687 PetscBool issym; 2688 2689 ierr = MatGetOption(pc->mat,MAT_SYMMETRIC,&issym);CHKERRQ(ierr); 2690 if (issym || pcbddc->symmetric_primal) { 2691 ierr = MatSetOption(newmat,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 2692 } 2693 } 2694 ierr = MatSetOptionsPrefix(newmat,prefix);CHKERRQ(ierr); 2695 ierr = MatAppendOptionsPrefix(newmat,"fetidp_");CHKERRQ(ierr); 2696 ierr = MatSetUp(newmat);CHKERRQ(ierr); 2697 /* FETI-DP preconditioner */ 2698 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 2699 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 2700 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 2701 ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr); 2702 ierr = PCSetOptionsPrefix(newpc,prefix);CHKERRQ(ierr); 2703 ierr = PCAppendOptionsPrefix(newpc,"fetidp_");CHKERRQ(ierr); 2704 ierr = PCSetErrorIfFailure(newpc,pc->erroriffailure);CHKERRQ(ierr); 2705 if (!fetidpmat_ctx->l2g_lambda_only) { /* standard FETI-DP */ 2706 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 2707 ierr = PCShellSetName(newpc,"FETI-DP multipliers");CHKERRQ(ierr); 2708 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 2709 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 2710 ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 2711 ierr = PCShellSetView(newpc,FETIDPPCView);CHKERRQ(ierr); 2712 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 2713 } else { /* saddle-point FETI-DP */ 2714 Mat M; 2715 PetscInt psize; 2716 PetscBool fake = PETSC_FALSE, isfieldsplit; 2717 2718 ierr = ISViewFromOptions(fetidpmat_ctx->lagrange,NULL,"-lag_view");CHKERRQ(ierr); 2719 ierr = ISViewFromOptions(fetidpmat_ctx->pressure,NULL,"-press_view");CHKERRQ(ierr); 2720 ierr = PetscObjectQuery((PetscObject)pc,"__KSPFETIDP_PPmat",(PetscObject*)&M);CHKERRQ(ierr); 2721 ierr = PCSetType(newpc,PCFIELDSPLIT);CHKERRQ(ierr); 2722 ierr = PCFieldSplitSetIS(newpc,"lag",fetidpmat_ctx->lagrange);CHKERRQ(ierr); 2723 ierr = PCFieldSplitSetIS(newpc,"p",fetidpmat_ctx->pressure);CHKERRQ(ierr); 2724 ierr = PCFieldSplitSetType(newpc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 2725 ierr = PCFieldSplitSetSchurFactType(newpc,PC_FIELDSPLIT_SCHUR_FACT_DIAG);CHKERRQ(ierr); 2726 ierr = ISGetSize(fetidpmat_ctx->pressure,&psize);CHKERRQ(ierr); 2727 if (psize != M->rmap->N) { 2728 Mat M2; 2729 PetscInt lpsize; 2730 2731 fake = PETSC_TRUE; 2732 ierr = ISGetLocalSize(fetidpmat_ctx->pressure,&lpsize);CHKERRQ(ierr); 2733 ierr = MatCreate(comm,&M2);CHKERRQ(ierr); 2734 ierr = MatSetType(M2,MATAIJ);CHKERRQ(ierr); 2735 ierr = MatSetSizes(M2,lpsize,lpsize,psize,psize);CHKERRQ(ierr); 2736 ierr = MatSetUp(M2);CHKERRQ(ierr); 2737 ierr = MatAssemblyBegin(M2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2738 ierr = MatAssemblyEnd(M2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2739 ierr = PCFieldSplitSetSchurPre(newpc,PC_FIELDSPLIT_SCHUR_PRE_USER,M2);CHKERRQ(ierr); 2740 ierr = MatDestroy(&M2);CHKERRQ(ierr); 2741 } else { 2742 ierr = PCFieldSplitSetSchurPre(newpc,PC_FIELDSPLIT_SCHUR_PRE_USER,M);CHKERRQ(ierr); 2743 } 2744 ierr = PCFieldSplitSetSchurScale(newpc,1.0);CHKERRQ(ierr); 2745 2746 /* we need to setfromoptions and setup here to access the blocks */ 2747 ierr = PCSetFromOptions(newpc);CHKERRQ(ierr); 2748 ierr = PCSetUp(newpc);CHKERRQ(ierr); 2749 2750 /* user may have changed the type (e.g. -fetidp_pc_type none) */ 2751 ierr = PetscObjectTypeCompare((PetscObject)newpc,PCFIELDSPLIT,&isfieldsplit);CHKERRQ(ierr); 2752 if (isfieldsplit) { 2753 KSP *ksps; 2754 PC ppc,lagpc; 2755 PetscInt nn; 2756 PetscBool ismatis,matisok = PETSC_FALSE,check = PETSC_FALSE; 2757 2758 /* set the solver for the (0,0) block */ 2759 ierr = PCFieldSplitSchurGetSubKSP(newpc,&nn,&ksps);CHKERRQ(ierr); 2760 if (!nn) { /* not of type PC_COMPOSITE_SCHUR */ 2761 ierr = PCFieldSplitGetSubKSP(newpc,&nn,&ksps);CHKERRQ(ierr); 2762 if (!fake) { /* pass pmat to the pressure solver */ 2763 Mat F; 2764 2765 ierr = KSPGetOperators(ksps[1],&F,NULL);CHKERRQ(ierr); 2766 ierr = KSPSetOperators(ksps[1],F,M);CHKERRQ(ierr); 2767 } 2768 } else { 2769 PetscBool issym; 2770 Mat S; 2771 2772 ierr = PCFieldSplitSchurGetS(newpc,&S);CHKERRQ(ierr); 2773 2774 ierr = MatGetOption(newmat,MAT_SYMMETRIC,&issym);CHKERRQ(ierr); 2775 if (issym) { 2776 ierr = MatSetOption(S,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 2777 } 2778 } 2779 ierr = KSPGetPC(ksps[0],&lagpc);CHKERRQ(ierr); 2780 ierr = PCSetType(lagpc,PCSHELL);CHKERRQ(ierr); 2781 ierr = PCShellSetName(lagpc,"FETI-DP multipliers");CHKERRQ(ierr); 2782 ierr = PCShellSetContext(lagpc,fetidppc_ctx);CHKERRQ(ierr); 2783 ierr = PCShellSetApply(lagpc,FETIDPPCApply);CHKERRQ(ierr); 2784 ierr = PCShellSetApplyTranspose(lagpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 2785 ierr = PCShellSetView(lagpc,FETIDPPCView);CHKERRQ(ierr); 2786 ierr = PCShellSetDestroy(lagpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 2787 2788 /* Olof's idea: interface Schur complement preconditioner for the mass matrix */ 2789 ierr = KSPGetPC(ksps[1],&ppc);CHKERRQ(ierr); 2790 if (fake) { 2791 BDDCIPC_ctx bddcipc_ctx; 2792 PetscContainer c; 2793 2794 matisok = PETSC_TRUE; 2795 2796 /* create inner BDDC solver */ 2797 ierr = PetscNew(&bddcipc_ctx);CHKERRQ(ierr); 2798 ierr = PCCreate(comm,&bddcipc_ctx->bddc);CHKERRQ(ierr); 2799 ierr = PCSetType(bddcipc_ctx->bddc,PCBDDC);CHKERRQ(ierr); 2800 ierr = PCSetOperators(bddcipc_ctx->bddc,M,M);CHKERRQ(ierr); 2801 ierr = PetscObjectQuery((PetscObject)pc,"__KSPFETIDP_pCSR",(PetscObject*)&c);CHKERRQ(ierr); 2802 ierr = PetscObjectTypeCompare((PetscObject)M,MATIS,&ismatis);CHKERRQ(ierr); 2803 if (c && ismatis) { 2804 Mat lM; 2805 PetscInt *csr,n; 2806 2807 ierr = MatISGetLocalMat(M,&lM);CHKERRQ(ierr); 2808 ierr = MatGetSize(lM,&n,NULL);CHKERRQ(ierr); 2809 ierr = PetscContainerGetPointer(c,(void**)&csr);CHKERRQ(ierr); 2810 ierr = PCBDDCSetLocalAdjacencyGraph(bddcipc_ctx->bddc,n,csr,csr + (n + 1),PETSC_COPY_VALUES);CHKERRQ(ierr); 2811 ierr = MatISRestoreLocalMat(M,&lM);CHKERRQ(ierr); 2812 } 2813 ierr = PCSetOptionsPrefix(bddcipc_ctx->bddc,((PetscObject)ksps[1])->prefix);CHKERRQ(ierr); 2814 ierr = PCSetErrorIfFailure(bddcipc_ctx->bddc,pc->erroriffailure);CHKERRQ(ierr); 2815 ierr = PCSetFromOptions(bddcipc_ctx->bddc);CHKERRQ(ierr); 2816 2817 /* wrap the interface application */ 2818 ierr = PCSetType(ppc,PCSHELL);CHKERRQ(ierr); 2819 ierr = PCShellSetName(ppc,"FETI-DP pressure");CHKERRQ(ierr); 2820 ierr = PCShellSetContext(ppc,bddcipc_ctx);CHKERRQ(ierr); 2821 ierr = PCShellSetSetUp(ppc,PCSetUp_BDDCIPC);CHKERRQ(ierr); 2822 ierr = PCShellSetApply(ppc,PCApply_BDDCIPC);CHKERRQ(ierr); 2823 ierr = PCShellSetApplyTranspose(ppc,PCApplyTranspose_BDDCIPC);CHKERRQ(ierr); 2824 ierr = PCShellSetView(ppc,PCView_BDDCIPC);CHKERRQ(ierr); 2825 ierr = PCShellSetDestroy(ppc,PCDestroy_BDDCIPC);CHKERRQ(ierr); 2826 } 2827 2828 /* determine if we need to assemble M to construct a preconditioner */ 2829 if (!matisok) { 2830 ierr = PetscObjectTypeCompare((PetscObject)M,MATIS,&ismatis);CHKERRQ(ierr); 2831 ierr = PetscObjectTypeCompareAny((PetscObject)ppc,&matisok,PCBDDC,PCJACOBI,PCNONE,PCMG,"");CHKERRQ(ierr); 2832 if (ismatis && !matisok) { 2833 ierr = MatConvert(M,MATAIJ,MAT_INPLACE_MATRIX,&M);CHKERRQ(ierr); 2834 } 2835 } 2836 2837 /* run the subproblems to check convergence */ 2838 ierr = PetscOptionsGetBool(NULL,((PetscObject)newmat)->prefix,"-check_saddlepoint",&check,NULL);CHKERRQ(ierr); 2839 if (check) { 2840 PetscInt i; 2841 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; 2848 2849 ierr = KSPCreate(PetscObjectComm((PetscObject)ksps[i]),&kspC);CHKERRQ(ierr); 2850 ierr = KSPSetOptionsPrefix(kspC,((PetscObject)ksps[i])->prefix);CHKERRQ(ierr); 2851 ierr = KSPAppendOptionsPrefix(kspC,"check_");CHKERRQ(ierr); 2852 ierr = KSPGetOperators(ksps[i],&F,&pF);CHKERRQ(ierr); 2853 ierr = PetscObjectTypeCompare((PetscObject)F,MATSCHURCOMPLEMENT,&isschur);CHKERRQ(ierr); 2854 if (isschur) { 2855 KSP kspS,kspS2; 2856 Mat A00,pA00,A10,A01,A11; 2857 char prefix[256]; 2858 2859 ierr = MatSchurComplementGetKSP(F,&kspS);CHKERRQ(ierr); 2860 ierr = MatSchurComplementGetSubMatrices(F,&A00,&pA00,&A01,&A10,&A11);CHKERRQ(ierr); 2861 ierr = MatCreateSchurComplement(A00,pA00,A01,A10,A11,&F);CHKERRQ(ierr); 2862 ierr = MatSchurComplementGetKSP(F,&kspS2);CHKERRQ(ierr); 2863 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sschur_",((PetscObject)kspC)->prefix);CHKERRQ(ierr); 2864 ierr = KSPSetOptionsPrefix(kspS2,prefix);CHKERRQ(ierr); 2865 ierr = KSPGetPC(kspS2,&pc);CHKERRQ(ierr); 2866 ierr = PCSetType(pc,PCKSP);CHKERRQ(ierr); 2867 ierr = PCKSPSetKSP(pc,kspS);CHKERRQ(ierr); 2868 ierr = KSPSetFromOptions(kspS2);CHKERRQ(ierr); 2869 ierr = KSPGetPC(kspS2,&pc);CHKERRQ(ierr); 2870 ierr = PCSetUseAmat(pc,PETSC_TRUE);CHKERRQ(ierr); 2871 } else { 2872 ierr = PetscObjectReference((PetscObject)F);CHKERRQ(ierr); 2873 } 2874 ierr = KSPSetFromOptions(kspC);CHKERRQ(ierr); 2875 ierr = PetscOptionsGetBool(NULL,((PetscObject)kspC)->prefix,"-preconditioned",&prec,NULL);CHKERRQ(ierr); 2876 if (prec) { 2877 ierr = KSPGetPC(ksps[i],&pc);CHKERRQ(ierr); 2878 ierr = KSPSetPC(kspC,pc);CHKERRQ(ierr); 2879 } 2880 ierr = KSPSetOperators(kspC,F,pF);CHKERRQ(ierr); 2881 ierr = MatCreateVecs(F,&x,&y);CHKERRQ(ierr); 2882 ierr = VecSetRandom(x,NULL);CHKERRQ(ierr); 2883 ierr = MatMult(F,x,y);CHKERRQ(ierr); 2884 ierr = KSPSolve(kspC,y,x);CHKERRQ(ierr); 2885 ierr = KSPCheckSolve(kspC,pc,x);CHKERRQ(ierr); 2886 ierr = KSPDestroy(&kspC);CHKERRQ(ierr); 2887 ierr = MatDestroy(&F);CHKERRQ(ierr); 2888 ierr = VecDestroy(&x);CHKERRQ(ierr); 2889 ierr = VecDestroy(&y);CHKERRQ(ierr); 2890 } 2891 } 2892 ierr = PetscFree(ksps);CHKERRQ(ierr); 2893 } 2894 } 2895 /* return pointers for objects created */ 2896 *fetidp_mat = newmat; 2897 *fetidp_pc = newpc; 2898 PetscFunctionReturn(0); 2899 } 2900 2901 /*@C 2902 PCBDDCCreateFETIDPOperators - Create FETI-DP operators 2903 2904 Collective 2905 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) 2910 2911 Output Parameters: 2912 + fetidp_mat - shell FETI-DP matrix object 2913 - fetidp_pc - shell Dirichlet preconditioner for FETI-DP matrix 2914 2915 Level: developer 2916 2917 Notes: 2918 Currently the only operations provided for FETI-DP matrix are MatMult and MatMultTranspose 2919 2920 .seealso: PCBDDC, PCBDDCMatFETIDPGetRHS, PCBDDCMatFETIDPGetSolution 2921 @*/ 2922 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc) 2923 { 2924 PetscErrorCode ierr; 2925 2926 PetscFunctionBegin; 2927 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2928 if (pc->setupcalled) { 2929 ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,PetscBool,const char*,Mat*,PC*),(pc,fully_redundant,prefix,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 2930 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"You must call PCSetup_BDDC() first"); 2931 PetscFunctionReturn(0); 2932 } 2933 /* -------------------------------------------------------------------------- */ 2934 /*MC 2935 PCBDDC - Balancing Domain Decomposition by Constraints. 2936 2937 An implementation of the BDDC preconditioner based on the bibliography found below. 2938 2939 The matrix to be preconditioned (Pmat) must be of type MATIS. 2940 2941 Currently works with MATIS matrices with local matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers. 2942 2943 It also works with unsymmetric and indefinite problems. 2944 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. 2946 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). 2948 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. 2951 2952 Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). Non-singular modes are retained via SVD. 2953 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() 2956 2957 The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using a MatPartitioning object. 2958 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. 2960 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. 2962 2963 Options Database Keys (some of them, run with -help for a complete list): 2964 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 2979 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. 2987 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 3000 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 3006 3007 Level: intermediate 3008 3009 Developer Notes: 3010 3011 Contributed by Stefano Zampini 3012 3013 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 3014 M*/ 3015 3016 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 3017 { 3018 PetscErrorCode ierr; 3019 PC_BDDC *pcbddc; 3020 3021 PetscFunctionBegin; 3022 ierr = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr); 3023 pc->data = pcbddc; 3024 3025 /* create PCIS data structure */ 3026 ierr = PCISCreate(pc);CHKERRQ(ierr); 3027 3028 /* create local graph structure */ 3029 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 3030 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; 3050 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; 3064 3065 /* composing function */ 3066 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDiscreteGradient_C",PCBDDCSetDiscreteGradient_BDDC);CHKERRQ(ierr); 3067 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",PCBDDCSetDivergenceMat_BDDC);CHKERRQ(ierr); 3068 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);CHKERRQ(ierr); 3069 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 3070 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",PCBDDCSetPrimalVerticesIS_BDDC);CHKERRQ(ierr); 3071 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetPrimalVerticesLocalIS_C",PCBDDCGetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 3072 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetPrimalVerticesIS_C",PCBDDCGetPrimalVerticesIS_BDDC);CHKERRQ(ierr); 3073 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 3074 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr); 3075 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr); 3076 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr); 3077 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 3078 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 3079 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 3080 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 3081 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 3082 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 3083 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 3084 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 3085 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 3086 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr); 3087 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 3088 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 3089 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 3090 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 3091 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_BDDC);CHKERRQ(ierr); 3092 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_BDDC);CHKERRQ(ierr); 3093 PetscFunctionReturn(0); 3094 } 3095 3096 /*@C 3097 PCBDDCInitializePackage - This function initializes everything in the PCBDDC package. It is called 3098 from PCInitializePackage(). 3099 3100 Level: developer 3101 3102 .seealso: PetscInitialize() 3103 @*/ 3104 PetscErrorCode PCBDDCInitializePackage(void) 3105 { 3106 PetscErrorCode ierr; 3107 int i; 3108 3109 PetscFunctionBegin; 3110 if (PCBDDCPackageInitialized) PetscFunctionReturn(0); 3111 PCBDDCPackageInitialized = PETSC_TRUE; 3112 ierr = PetscRegisterFinalize(PCBDDCFinalizePackage);CHKERRQ(ierr); 3113 3114 /* general events */ 3115 ierr = PetscLogEventRegister("PCBDDCTopo",PC_CLASSID,&PC_BDDC_Topology[0]);CHKERRQ(ierr); 3116 ierr = PetscLogEventRegister("PCBDDCLKSP",PC_CLASSID,&PC_BDDC_LocalSolvers[0]);CHKERRQ(ierr); 3117 ierr = PetscLogEventRegister("PCBDDCLWor",PC_CLASSID,&PC_BDDC_LocalWork[0]);CHKERRQ(ierr); 3118 ierr = PetscLogEventRegister("PCBDDCCorr",PC_CLASSID,&PC_BDDC_CorrectionSetUp[0]);CHKERRQ(ierr); 3119 ierr = PetscLogEventRegister("PCBDDCASet",PC_CLASSID,&PC_BDDC_ApproxSetUp[0]);CHKERRQ(ierr); 3120 ierr = PetscLogEventRegister("PCBDDCAApp",PC_CLASSID,&PC_BDDC_ApproxApply[0]);CHKERRQ(ierr); 3121 ierr = PetscLogEventRegister("PCBDDCCSet",PC_CLASSID,&PC_BDDC_CoarseSetUp[0]);CHKERRQ(ierr); 3122 ierr = PetscLogEventRegister("PCBDDCCKSP",PC_CLASSID,&PC_BDDC_CoarseSolver[0]);CHKERRQ(ierr); 3123 ierr = PetscLogEventRegister("PCBDDCAdap",PC_CLASSID,&PC_BDDC_AdaptiveSetUp[0]);CHKERRQ(ierr); 3124 ierr = PetscLogEventRegister("PCBDDCScal",PC_CLASSID,&PC_BDDC_Scaling[0]);CHKERRQ(ierr); 3125 ierr = PetscLogEventRegister("PCBDDCSchr",PC_CLASSID,&PC_BDDC_Schurs[0]);CHKERRQ(ierr); 3126 ierr = PetscLogEventRegister("PCBDDCDirS",PC_CLASSID,&PC_BDDC_Solves[0][0]);CHKERRQ(ierr); 3127 ierr = PetscLogEventRegister("PCBDDCNeuS",PC_CLASSID,&PC_BDDC_Solves[0][1]);CHKERRQ(ierr); 3128 ierr = PetscLogEventRegister("PCBDDCCoaS",PC_CLASSID,&PC_BDDC_Solves[0][2]);CHKERRQ(ierr); 3129 for (i=1;i<PETSC_PCBDDC_MAXLEVELS;i++) { 3130 char ename[32]; 3131 3132 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCTopo l%02d",i);CHKERRQ(ierr); 3133 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Topology[i]);CHKERRQ(ierr); 3134 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCLKSP l%02d",i);CHKERRQ(ierr); 3135 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_LocalSolvers[i]);CHKERRQ(ierr); 3136 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCLWor l%02d",i);CHKERRQ(ierr); 3137 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_LocalWork[i]);CHKERRQ(ierr); 3138 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCCorr l%02d",i);CHKERRQ(ierr); 3139 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_CorrectionSetUp[i]);CHKERRQ(ierr); 3140 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCASet l%02d",i);CHKERRQ(ierr); 3141 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_ApproxSetUp[i]);CHKERRQ(ierr); 3142 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCAApp l%02d",i);CHKERRQ(ierr); 3143 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_ApproxApply[i]);CHKERRQ(ierr); 3144 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCCSet l%02d",i);CHKERRQ(ierr); 3145 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_CoarseSetUp[i]);CHKERRQ(ierr); 3146 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCCKSP l%02d",i);CHKERRQ(ierr); 3147 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_CoarseSolver[i]);CHKERRQ(ierr); 3148 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCAdap l%02d",i);CHKERRQ(ierr); 3149 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_AdaptiveSetUp[i]);CHKERRQ(ierr); 3150 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCScal l%02d",i);CHKERRQ(ierr); 3151 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Scaling[i]);CHKERRQ(ierr); 3152 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCSchr l%02d",i);CHKERRQ(ierr); 3153 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Schurs[i]);CHKERRQ(ierr); 3154 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCDirS l%02d",i);CHKERRQ(ierr); 3155 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Solves[i][0]);CHKERRQ(ierr); 3156 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCNeuS l%02d",i);CHKERRQ(ierr); 3157 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Solves[i][1]);CHKERRQ(ierr); 3158 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCCoaS l%02d",i);CHKERRQ(ierr); 3159 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Solves[i][2]);CHKERRQ(ierr); 3160 } 3161 PetscFunctionReturn(0); 3162 } 3163 3164 /*@C 3165 PCBDDCFinalizePackage - This function frees everything from the PCBDDC package. It is 3166 called from PetscFinalize() automatically. 3167 3168 Level: developer 3169 3170 .seealso: PetscFinalize() 3171 @*/ 3172 PetscErrorCode PCBDDCFinalizePackage(void) 3173 { 3174 PetscFunctionBegin; 3175 PCBDDCPackageInitialized = PETSC_FALSE; 3176 PetscFunctionReturn(0); 3177 } 3178