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," Corner selection: %d (selected %d)\n",pcbddc->corner_selection,pcbddc->corner_selected);CHKERRQ(ierr); 169 ierr = PetscViewerASCIIPrintf(viewer," Use vertices: %d (vertex size %D)\n",pcbddc->use_vertices,pcbddc->vertex_size);CHKERRQ(ierr); 170 ierr = PetscViewerASCIIPrintf(viewer," Use edges: %d\n",pcbddc->use_edges);CHKERRQ(ierr); 171 ierr = PetscViewerASCIIPrintf(viewer," Use faces: %d\n",pcbddc->use_faces);CHKERRQ(ierr); 172 ierr = PetscViewerASCIIPrintf(viewer," Use true near null space: %d\n",pcbddc->use_nnsp_true);CHKERRQ(ierr); 173 ierr = PetscViewerASCIIPrintf(viewer," Use QR for single constraints on cc: %d\n",pcbddc->use_qr_single);CHKERRQ(ierr); 174 ierr = PetscViewerASCIIPrintf(viewer," Use change of basis on local edge nodes: %d\n",pcbddc->use_change_of_basis);CHKERRQ(ierr); 175 ierr = PetscViewerASCIIPrintf(viewer," Use change of basis on local face nodes: %d\n",pcbddc->use_change_on_faces);CHKERRQ(ierr); 176 ierr = PetscViewerASCIIPrintf(viewer," User defined change of basis matrix: %d\n",!!pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 177 ierr = PetscViewerASCIIPrintf(viewer," Has change of basis matrix: %d\n",!!pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 178 ierr = PetscViewerASCIIPrintf(viewer," Eliminate dirichlet boundary dofs: %d\n",pcbddc->eliminate_dirdofs);CHKERRQ(ierr); 179 ierr = PetscViewerASCIIPrintf(viewer," Switch on static condensation ops around the interface preconditioner: %d\n",pcbddc->switch_static);CHKERRQ(ierr); 180 ierr = PetscViewerASCIIPrintf(viewer," Use exact dirichlet trick: %d\n",pcbddc->use_exact_dirichlet_trick);CHKERRQ(ierr); 181 ierr = PetscViewerASCIIPrintf(viewer," Interface extension: %s\n",PCBDDCInterfaceExtTypes[pcbddc->interface_extension]);CHKERRQ(ierr); 182 ierr = PetscViewerASCIIPrintf(viewer," Multilevel max levels: %D\n",pcbddc->max_levels);CHKERRQ(ierr); 183 ierr = PetscViewerASCIIPrintf(viewer," Multilevel coarsening ratio: %D\n",pcbddc->coarsening_ratio);CHKERRQ(ierr); 184 ierr = PetscViewerASCIIPrintf(viewer," Use estimated eigs for coarse problem: %d\n",pcbddc->use_coarse_estimates);CHKERRQ(ierr); 185 ierr = PetscViewerASCIIPrintf(viewer," Use deluxe scaling: %d\n",pcbddc->use_deluxe_scaling);CHKERRQ(ierr); 186 ierr = PetscViewerASCIIPrintf(viewer," Use deluxe zerorows: %d\n",pcbddc->deluxe_zerorows);CHKERRQ(ierr); 187 ierr = PetscViewerASCIIPrintf(viewer," Use deluxe singlemat: %d\n",pcbddc->deluxe_singlemat);CHKERRQ(ierr); 188 ierr = PetscViewerASCIIPrintf(viewer," Rebuild interface graph for Schur principal minors: %d\n",pcbddc->sub_schurs_rebuild);CHKERRQ(ierr); 189 ierr = PetscViewerASCIIPrintf(viewer," Number of dofs' layers for the computation of principal minors: %D\n",pcbddc->sub_schurs_layers);CHKERRQ(ierr); 190 ierr = PetscViewerASCIIPrintf(viewer," Use user CSR graph to compute successive layers: %d\n",pcbddc->sub_schurs_use_useradj);CHKERRQ(ierr); 191 if (pcbddc->adaptive_threshold[1] != pcbddc->adaptive_threshold[0]) { 192 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); 193 } else { 194 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); 195 } 196 ierr = PetscViewerASCIIPrintf(viewer," Min constraints / connected component: %D\n",pcbddc->adaptive_nmin);CHKERRQ(ierr); 197 ierr = PetscViewerASCIIPrintf(viewer," Max constraints / connected component: %D\n",pcbddc->adaptive_nmax);CHKERRQ(ierr); 198 ierr = PetscViewerASCIIPrintf(viewer," Invert exact Schur complement for adaptive selection: %d\n",pcbddc->sub_schurs_exact_schur);CHKERRQ(ierr); 199 ierr = PetscViewerASCIIPrintf(viewer," Symmetric computation of primal basis functions: %d\n",pcbddc->symmetric_primal);CHKERRQ(ierr); 200 ierr = PetscViewerASCIIPrintf(viewer," Num. Procs. to map coarse adjacency list: %D\n",pcbddc->coarse_adj_red);CHKERRQ(ierr); 201 ierr = PetscViewerASCIIPrintf(viewer," Coarse eqs per proc (significant at the coarsest level): %D\n",pcbddc->coarse_eqs_per_proc);CHKERRQ(ierr); 202 ierr = PetscViewerASCIIPrintf(viewer," Detect disconnected: %d (filter %d)\n",pcbddc->detect_disconnected,pcbddc->detect_disconnected_filter);CHKERRQ(ierr); 203 ierr = PetscViewerASCIIPrintf(viewer," Benign subspace trick: %d (change explicit %d)\n",pcbddc->benign_saddle_point,pcbddc->benign_change_explicit);CHKERRQ(ierr); 204 ierr = PetscViewerASCIIPrintf(viewer," Benign subspace trick is active: %d\n",pcbddc->benign_have_null);CHKERRQ(ierr); 205 ierr = PetscViewerASCIIPrintf(viewer," Algebraic computation of no-net-flux: %d\n",pcbddc->compute_nonetflux);CHKERRQ(ierr); 206 if (!pc->setupcalled) PetscFunctionReturn(0); 207 208 /* compute interface size */ 209 ierr = VecSet(pcis->vec1_B,1.0);CHKERRQ(ierr); 210 ierr = MatCreateVecs(pc->pmat,&counter,NULL);CHKERRQ(ierr); 211 ierr = VecSet(counter,0.0);CHKERRQ(ierr); 212 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,counter,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 213 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,counter,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 214 ierr = VecSum(counter,&interface_size);CHKERRQ(ierr); 215 ierr = VecDestroy(&counter);CHKERRQ(ierr); 216 217 /* compute some statistics on the domain decomposition */ 218 gsum[0] = 1; 219 gsum[1] = gsum[2] = gsum[3] = gsum[4] = gsum[5] = 0; 220 loc[0] = !!pcis->n; 221 loc[1] = pcis->n - pcis->n_B; 222 loc[2] = pcis->n_B; 223 loc[3] = pcbddc->local_primal_size; 224 loc[4] = pcis->n; 225 loc[5] = pcbddc->n_local_subs > 0 ? pcbddc->n_local_subs : (pcis->n ? 1 : 0); 226 loc[6] = pcbddc->benign_n; 227 ierr = MPI_Reduce(loc,gsum,6,MPIU_INT64,MPI_SUM,0,PetscObjectComm((PetscObject)pc));CHKERRMPI(ierr); 228 if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = loc[5] = -1; 229 ierr = MPI_Reduce(loc,gmax,6,MPIU_INT64,MPI_MAX,0,PetscObjectComm((PetscObject)pc));CHKERRMPI(ierr); 230 if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = loc[5] = PETSC_MAX_INT; 231 ierr = MPI_Reduce(loc,gmin,6,MPIU_INT64,MPI_MIN,0,PetscObjectComm((PetscObject)pc));CHKERRMPI(ierr); 232 ierr = MPI_Reduce(&loc[6],&totbenign,1,MPIU_INT64,MPI_SUM,0,PetscObjectComm((PetscObject)pc));CHKERRMPI(ierr); 233 if (pcbddc->coarse_size) { 234 ratio1 = pc->pmat->rmap->N/(1.*pcbddc->coarse_size); 235 ratio2 = PetscRealPart(interface_size)/pcbddc->coarse_size; 236 } 237 ierr = PetscViewerASCIIPrintf(viewer,"********************************** STATISTICS AT LEVEL %d **********************************\n",pcbddc->current_level);CHKERRQ(ierr); 238 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); 239 ierr = PetscViewerASCIIPrintf(viewer," Coarsening ratios: all/coarse %D interface/coarse %D\n",(PetscInt)ratio1,(PetscInt)ratio2);CHKERRQ(ierr); 240 ierr = PetscViewerASCIIPrintf(viewer," Active processes : %D\n",(PetscInt)gsum[0]);CHKERRQ(ierr); 241 ierr = PetscViewerASCIIPrintf(viewer," Total subdomains : %D\n",(PetscInt)gsum[5]);CHKERRQ(ierr); 242 if (pcbddc->benign_have_null) { 243 ierr = PetscViewerASCIIPrintf(viewer," Benign subs : %D\n",(PetscInt)totbenign);CHKERRQ(ierr); 244 } 245 ierr = PetscViewerASCIIPrintf(viewer," Dofs type :\tMIN\tMAX\tMEAN\n");CHKERRQ(ierr); 246 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); 247 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); 248 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); 249 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); 250 ierr = PetscViewerASCIIPrintf(viewer," Local subs :\t%D\t%D\n" ,(PetscInt)gmin[5],(PetscInt)gmax[5]);CHKERRQ(ierr); 251 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 252 253 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);CHKERRMPI(ierr); 254 255 /* local solvers */ 256 ierr = PetscViewerGetSubViewer(viewer,PetscObjectComm((PetscObject)pcbddc->ksp_D),&subviewer);CHKERRQ(ierr); 257 if (rank == 0) { 258 ierr = PetscViewerASCIIPrintf(subviewer,"--- Interior solver (rank 0)\n");CHKERRQ(ierr); 259 ierr = PetscViewerASCIIPushTab(subviewer);CHKERRQ(ierr); 260 ierr = KSPView(pcbddc->ksp_D,subviewer);CHKERRQ(ierr); 261 ierr = PetscViewerASCIIPopTab(subviewer);CHKERRQ(ierr); 262 ierr = PetscViewerASCIIPrintf(subviewer,"--- Correction solver (rank 0)\n");CHKERRQ(ierr); 263 ierr = PetscViewerASCIIPushTab(subviewer);CHKERRQ(ierr); 264 ierr = KSPView(pcbddc->ksp_R,subviewer);CHKERRQ(ierr); 265 ierr = PetscViewerASCIIPopTab(subviewer);CHKERRQ(ierr); 266 ierr = PetscViewerFlush(subviewer);CHKERRQ(ierr); 267 } 268 ierr = PetscViewerRestoreSubViewer(viewer,PetscObjectComm((PetscObject)pcbddc->ksp_D),&subviewer);CHKERRQ(ierr); 269 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 270 271 /* the coarse problem can be handled by a different communicator */ 272 if (pcbddc->coarse_ksp) color = 1; 273 else color = 0; 274 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRMPI(ierr); 275 ierr = PetscSubcommCreate(PetscObjectComm((PetscObject)pc),&subcomm);CHKERRQ(ierr); 276 ierr = PetscSubcommSetNumber(subcomm,PetscMin(size,2));CHKERRQ(ierr); 277 ierr = PetscSubcommSetTypeGeneral(subcomm,color,rank);CHKERRQ(ierr); 278 ierr = PetscViewerGetSubViewer(viewer,PetscSubcommChild(subcomm),&subviewer);CHKERRQ(ierr); 279 if (color == 1) { 280 ierr = PetscViewerASCIIPrintf(subviewer,"--- Coarse solver\n");CHKERRQ(ierr); 281 ierr = PetscViewerASCIIPushTab(subviewer);CHKERRQ(ierr); 282 ierr = KSPView(pcbddc->coarse_ksp,subviewer);CHKERRQ(ierr); 283 ierr = PetscViewerASCIIPopTab(subviewer);CHKERRQ(ierr); 284 ierr = PetscViewerFlush(subviewer);CHKERRQ(ierr); 285 } 286 ierr = PetscViewerRestoreSubViewer(viewer,PetscSubcommChild(subcomm),&subviewer);CHKERRQ(ierr); 287 ierr = PetscSubcommDestroy(&subcomm);CHKERRQ(ierr); 288 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 289 } 290 PetscFunctionReturn(0); 291 } 292 293 static PetscErrorCode PCBDDCSetDiscreteGradient_BDDC(PC pc, Mat G, PetscInt order, PetscInt field, PetscBool global, PetscBool conforming) 294 { 295 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 296 PetscErrorCode ierr; 297 298 PetscFunctionBegin; 299 ierr = PetscObjectReference((PetscObject)G);CHKERRQ(ierr); 300 ierr = MatDestroy(&pcbddc->discretegradient);CHKERRQ(ierr); 301 pcbddc->discretegradient = G; 302 pcbddc->nedorder = order > 0 ? order : -order; 303 pcbddc->nedfield = field; 304 pcbddc->nedglobal = global; 305 pcbddc->conforming = conforming; 306 PetscFunctionReturn(0); 307 } 308 309 /*@ 310 PCBDDCSetDiscreteGradient - Sets the discrete gradient 311 312 Collective on PC 313 314 Input Parameters: 315 + pc - the preconditioning context 316 . G - the discrete gradient matrix (should be in AIJ format) 317 . order - the order of the Nedelec space (1 for the lowest order) 318 . field - the field id of the Nedelec dofs (not used if the fields have not been specified) 319 . global - the type of global ordering for the rows of G 320 - conforming - whether the mesh is conforming or not 321 322 Level: advanced 323 324 Notes: 325 The discrete gradient matrix G is used to analyze the subdomain edges, and it should not contain any zero entry. 326 For variable order spaces, the order should be set to zero. 327 If global is true, the rows of G should be given in global ordering for the whole dofs; 328 if false, the ordering should be global for the Nedelec field. 329 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 330 and geid the one for the Nedelec field. 331 332 .seealso: PCBDDC,PCBDDCSetDofsSplitting(),PCBDDCSetDofsSplittingLocal() 333 @*/ 334 PetscErrorCode PCBDDCSetDiscreteGradient(PC pc, Mat G, PetscInt order, PetscInt field, PetscBool global, PetscBool conforming) 335 { 336 PetscErrorCode ierr; 337 338 PetscFunctionBegin; 339 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 340 PetscValidHeaderSpecific(G,MAT_CLASSID,2); 341 PetscValidLogicalCollectiveInt(pc,order,3); 342 PetscValidLogicalCollectiveInt(pc,field,4); 343 PetscValidLogicalCollectiveBool(pc,global,5); 344 PetscValidLogicalCollectiveBool(pc,conforming,6); 345 PetscCheckSameComm(pc,1,G,2); 346 ierr = PetscTryMethod(pc,"PCBDDCSetDiscreteGradient_C",(PC,Mat,PetscInt,PetscInt,PetscBool,PetscBool),(pc,G,order,field,global,conforming));CHKERRQ(ierr); 347 PetscFunctionReturn(0); 348 } 349 350 static PetscErrorCode PCBDDCSetDivergenceMat_BDDC(PC pc, Mat divudotp, PetscBool trans, IS vl2l) 351 { 352 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 353 PetscErrorCode ierr; 354 355 PetscFunctionBegin; 356 ierr = PetscObjectReference((PetscObject)divudotp);CHKERRQ(ierr); 357 ierr = MatDestroy(&pcbddc->divudotp);CHKERRQ(ierr); 358 pcbddc->divudotp = divudotp; 359 pcbddc->divudotp_trans = trans; 360 pcbddc->compute_nonetflux = PETSC_TRUE; 361 if (vl2l) { 362 ierr = PetscObjectReference((PetscObject)vl2l);CHKERRQ(ierr); 363 ierr = ISDestroy(&pcbddc->divudotp_vl2l);CHKERRQ(ierr); 364 pcbddc->divudotp_vl2l = vl2l; 365 } 366 PetscFunctionReturn(0); 367 } 368 369 /*@ 370 PCBDDCSetDivergenceMat - Sets the linear operator representing \int_\Omega \div {\bf u} \cdot p dx 371 372 Collective on PC 373 374 Input Parameters: 375 + pc - the preconditioning context 376 . divudotp - the matrix (must be of type MATIS) 377 . trans - if trans if false (resp. true), then pressures are in the test (trial) space and velocities are in the trial (test) space. 378 - vl2l - optional index set describing the local (wrt the local matrix in divudotp) to local (wrt the local matrix in the preconditioning matrix) map for the velocities 379 380 Level: advanced 381 382 Notes: 383 This auxiliary matrix is used to compute quadrature weights representing the net-flux across subdomain boundaries 384 If vl2l is NULL, the local ordering for velocities in divudotp should match that of the preconditioning matrix 385 386 .seealso: PCBDDC 387 @*/ 388 PetscErrorCode PCBDDCSetDivergenceMat(PC pc, Mat divudotp, PetscBool trans, IS vl2l) 389 { 390 PetscBool ismatis; 391 PetscErrorCode ierr; 392 393 PetscFunctionBegin; 394 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 395 PetscValidHeaderSpecific(divudotp,MAT_CLASSID,2); 396 PetscCheckSameComm(pc,1,divudotp,2); 397 PetscValidLogicalCollectiveBool(pc,trans,3); 398 if (vl2l) PetscValidHeaderSpecific(vl2l,IS_CLASSID,4); 399 ierr = PetscObjectTypeCompare((PetscObject)divudotp,MATIS,&ismatis);CHKERRQ(ierr); 400 PetscCheckFalse(!ismatis,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Divergence matrix needs to be of type MATIS"); 401 ierr = PetscTryMethod(pc,"PCBDDCSetDivergenceMat_C",(PC,Mat,PetscBool,IS),(pc,divudotp,trans,vl2l));CHKERRQ(ierr); 402 PetscFunctionReturn(0); 403 } 404 405 static PetscErrorCode PCBDDCSetChangeOfBasisMat_BDDC(PC pc, Mat change, PetscBool interior) 406 { 407 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 408 PetscErrorCode ierr; 409 410 PetscFunctionBegin; 411 ierr = PetscObjectReference((PetscObject)change);CHKERRQ(ierr); 412 ierr = MatDestroy(&pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 413 pcbddc->user_ChangeOfBasisMatrix = change; 414 pcbddc->change_interior = interior; 415 PetscFunctionReturn(0); 416 } 417 /*@ 418 PCBDDCSetChangeOfBasisMat - Set user defined change of basis for dofs 419 420 Collective on PC 421 422 Input Parameters: 423 + pc - the preconditioning context 424 . change - the change of basis matrix 425 - interior - whether or not the change of basis modifies interior dofs 426 427 Level: intermediate 428 429 Notes: 430 431 .seealso: PCBDDC 432 @*/ 433 PetscErrorCode PCBDDCSetChangeOfBasisMat(PC pc, Mat change, PetscBool interior) 434 { 435 PetscErrorCode ierr; 436 437 PetscFunctionBegin; 438 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 439 PetscValidHeaderSpecific(change,MAT_CLASSID,2); 440 PetscCheckSameComm(pc,1,change,2); 441 if (pc->mat) { 442 PetscInt rows_c,cols_c,rows,cols; 443 ierr = MatGetSize(pc->mat,&rows,&cols);CHKERRQ(ierr); 444 ierr = MatGetSize(change,&rows_c,&cols_c);CHKERRQ(ierr); 445 PetscCheckFalse(rows_c != rows,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of rows for change of basis matrix! %D != %D",rows_c,rows); 446 PetscCheckFalse(cols_c != cols,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of columns for change of basis matrix! %D != %D",cols_c,cols); 447 ierr = MatGetLocalSize(pc->mat,&rows,&cols);CHKERRQ(ierr); 448 ierr = MatGetLocalSize(change,&rows_c,&cols_c);CHKERRQ(ierr); 449 PetscCheckFalse(rows_c != rows,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of local rows for change of basis matrix! %D != %D",rows_c,rows); 450 PetscCheckFalse(cols_c != cols,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of local columns for change of basis matrix! %D != %D",cols_c,cols); 451 } 452 ierr = PetscTryMethod(pc,"PCBDDCSetChangeOfBasisMat_C",(PC,Mat,PetscBool),(pc,change,interior));CHKERRQ(ierr); 453 PetscFunctionReturn(0); 454 } 455 456 static PetscErrorCode PCBDDCSetPrimalVerticesIS_BDDC(PC pc, IS PrimalVertices) 457 { 458 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 459 PetscBool isequal = PETSC_FALSE; 460 PetscErrorCode ierr; 461 462 PetscFunctionBegin; 463 ierr = PetscObjectReference((PetscObject)PrimalVertices);CHKERRQ(ierr); 464 if (pcbddc->user_primal_vertices) { 465 ierr = ISEqual(PrimalVertices,pcbddc->user_primal_vertices,&isequal);CHKERRQ(ierr); 466 } 467 ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr); 468 ierr = ISDestroy(&pcbddc->user_primal_vertices_local);CHKERRQ(ierr); 469 pcbddc->user_primal_vertices = PrimalVertices; 470 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 471 PetscFunctionReturn(0); 472 } 473 474 /*@ 475 PCBDDCSetPrimalVerticesIS - Set additional user defined primal vertices in PCBDDC 476 477 Collective 478 479 Input Parameters: 480 + pc - the preconditioning context 481 - PrimalVertices - index set of primal vertices in global numbering (can be empty) 482 483 Level: intermediate 484 485 Notes: 486 Any process can list any global node 487 488 .seealso: PCBDDC, PCBDDCGetPrimalVerticesIS(), PCBDDCSetPrimalVerticesLocalIS(), PCBDDCGetPrimalVerticesLocalIS() 489 @*/ 490 PetscErrorCode PCBDDCSetPrimalVerticesIS(PC pc, IS PrimalVertices) 491 { 492 PetscErrorCode ierr; 493 494 PetscFunctionBegin; 495 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 496 PetscValidHeaderSpecific(PrimalVertices,IS_CLASSID,2); 497 PetscCheckSameComm(pc,1,PrimalVertices,2); 498 ierr = PetscTryMethod(pc,"PCBDDCSetPrimalVerticesIS_C",(PC,IS),(pc,PrimalVertices));CHKERRQ(ierr); 499 PetscFunctionReturn(0); 500 } 501 502 static PetscErrorCode PCBDDCGetPrimalVerticesIS_BDDC(PC pc, IS *is) 503 { 504 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 505 506 PetscFunctionBegin; 507 *is = pcbddc->user_primal_vertices; 508 PetscFunctionReturn(0); 509 } 510 511 /*@ 512 PCBDDCGetPrimalVerticesIS - Get user defined primal vertices set with PCBDDCSetPrimalVerticesIS() 513 514 Collective 515 516 Input Parameters: 517 . pc - the preconditioning context 518 519 Output Parameters: 520 . is - index set of primal vertices in global numbering (NULL if not set) 521 522 Level: intermediate 523 524 Notes: 525 526 .seealso: PCBDDC, PCBDDCSetPrimalVerticesIS(), PCBDDCSetPrimalVerticesLocalIS(), PCBDDCGetPrimalVerticesLocalIS() 527 @*/ 528 PetscErrorCode PCBDDCGetPrimalVerticesIS(PC pc, IS *is) 529 { 530 PetscErrorCode ierr; 531 532 PetscFunctionBegin; 533 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 534 PetscValidPointer(is,2); 535 ierr = PetscUseMethod(pc,"PCBDDCGetPrimalVerticesIS_C",(PC,IS*),(pc,is));CHKERRQ(ierr); 536 PetscFunctionReturn(0); 537 } 538 539 static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices) 540 { 541 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 542 PetscBool isequal = PETSC_FALSE; 543 PetscErrorCode ierr; 544 545 PetscFunctionBegin; 546 ierr = PetscObjectReference((PetscObject)PrimalVertices);CHKERRQ(ierr); 547 if (pcbddc->user_primal_vertices_local) { 548 ierr = ISEqual(PrimalVertices,pcbddc->user_primal_vertices_local,&isequal);CHKERRQ(ierr); 549 } 550 ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr); 551 ierr = ISDestroy(&pcbddc->user_primal_vertices_local);CHKERRQ(ierr); 552 pcbddc->user_primal_vertices_local = PrimalVertices; 553 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 554 PetscFunctionReturn(0); 555 } 556 557 /*@ 558 PCBDDCSetPrimalVerticesLocalIS - Set additional user defined primal vertices in PCBDDC 559 560 Collective 561 562 Input Parameters: 563 + pc - the preconditioning context 564 - PrimalVertices - index set of primal vertices in local numbering (can be empty) 565 566 Level: intermediate 567 568 Notes: 569 570 .seealso: PCBDDC, PCBDDCSetPrimalVerticesIS(), PCBDDCGetPrimalVerticesIS(), PCBDDCGetPrimalVerticesLocalIS() 571 @*/ 572 PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices) 573 { 574 PetscErrorCode ierr; 575 576 PetscFunctionBegin; 577 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 578 PetscValidHeaderSpecific(PrimalVertices,IS_CLASSID,2); 579 PetscCheckSameComm(pc,1,PrimalVertices,2); 580 ierr = PetscTryMethod(pc,"PCBDDCSetPrimalVerticesLocalIS_C",(PC,IS),(pc,PrimalVertices));CHKERRQ(ierr); 581 PetscFunctionReturn(0); 582 } 583 584 static PetscErrorCode PCBDDCGetPrimalVerticesLocalIS_BDDC(PC pc, IS *is) 585 { 586 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 587 588 PetscFunctionBegin; 589 *is = pcbddc->user_primal_vertices_local; 590 PetscFunctionReturn(0); 591 } 592 593 /*@ 594 PCBDDCGetPrimalVerticesLocalIS - Get user defined primal vertices set with PCBDDCSetPrimalVerticesLocalIS() 595 596 Collective 597 598 Input Parameters: 599 . pc - the preconditioning context 600 601 Output Parameters: 602 . is - index set of primal vertices in local numbering (NULL if not set) 603 604 Level: intermediate 605 606 Notes: 607 608 .seealso: PCBDDC, PCBDDCSetPrimalVerticesIS(), PCBDDCGetPrimalVerticesIS(), PCBDDCSetPrimalVerticesLocalIS() 609 @*/ 610 PetscErrorCode PCBDDCGetPrimalVerticesLocalIS(PC pc, IS *is) 611 { 612 PetscErrorCode ierr; 613 614 PetscFunctionBegin; 615 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 616 PetscValidPointer(is,2); 617 ierr = PetscUseMethod(pc,"PCBDDCGetPrimalVerticesLocalIS_C",(PC,IS*),(pc,is));CHKERRQ(ierr); 618 PetscFunctionReturn(0); 619 } 620 621 static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc,PetscInt k) 622 { 623 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 624 625 PetscFunctionBegin; 626 pcbddc->coarsening_ratio = k; 627 PetscFunctionReturn(0); 628 } 629 630 /*@ 631 PCBDDCSetCoarseningRatio - Set coarsening ratio used in multilevel 632 633 Logically collective on PC 634 635 Input Parameters: 636 + pc - the preconditioning context 637 - k - coarsening ratio (H/h at the coarser level) 638 639 Options Database Keys: 640 . -pc_bddc_coarsening_ratio 641 642 Level: intermediate 643 644 Notes: 645 Approximatively k subdomains at the finer level will be aggregated into a single subdomain at the coarser level 646 647 .seealso: PCBDDC, PCBDDCSetLevels() 648 @*/ 649 PetscErrorCode PCBDDCSetCoarseningRatio(PC pc,PetscInt k) 650 { 651 PetscErrorCode ierr; 652 653 PetscFunctionBegin; 654 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 655 PetscValidLogicalCollectiveInt(pc,k,2); 656 ierr = PetscTryMethod(pc,"PCBDDCSetCoarseningRatio_C",(PC,PetscInt),(pc,k));CHKERRQ(ierr); 657 PetscFunctionReturn(0); 658 } 659 660 /* The following functions (PCBDDCSetUseExactDirichlet PCBDDCSetLevel) are not public */ 661 static PetscErrorCode PCBDDCSetUseExactDirichlet_BDDC(PC pc,PetscBool flg) 662 { 663 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 664 665 PetscFunctionBegin; 666 pcbddc->use_exact_dirichlet_trick = flg; 667 PetscFunctionReturn(0); 668 } 669 670 PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool flg) 671 { 672 PetscErrorCode ierr; 673 674 PetscFunctionBegin; 675 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 676 PetscValidLogicalCollectiveBool(pc,flg,2); 677 ierr = PetscTryMethod(pc,"PCBDDCSetUseExactDirichlet_C",(PC,PetscBool),(pc,flg));CHKERRQ(ierr); 678 PetscFunctionReturn(0); 679 } 680 681 static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc,PetscInt level) 682 { 683 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 684 685 PetscFunctionBegin; 686 pcbddc->current_level = level; 687 PetscFunctionReturn(0); 688 } 689 690 PetscErrorCode PCBDDCSetLevel(PC pc,PetscInt level) 691 { 692 PetscErrorCode ierr; 693 694 PetscFunctionBegin; 695 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 696 PetscValidLogicalCollectiveInt(pc,level,2); 697 ierr = PetscTryMethod(pc,"PCBDDCSetLevel_C",(PC,PetscInt),(pc,level));CHKERRQ(ierr); 698 PetscFunctionReturn(0); 699 } 700 701 static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc,PetscInt levels) 702 { 703 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 704 705 PetscFunctionBegin; 706 PetscCheckFalse(levels > PETSC_PCBDDC_MAXLEVELS-1,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Maximum number of additional levels for BDDC is %d",PETSC_PCBDDC_MAXLEVELS-1); 707 pcbddc->max_levels = levels; 708 PetscFunctionReturn(0); 709 } 710 711 /*@ 712 PCBDDCSetLevels - Sets the maximum number of additional levels allowed for multilevel BDDC 713 714 Logically collective on PC 715 716 Input Parameters: 717 + pc - the preconditioning context 718 - levels - the maximum number of levels 719 720 Options Database Keys: 721 . -pc_bddc_levels 722 723 Level: intermediate 724 725 Notes: 726 The default value is 0, that gives the classical two-levels BDDC 727 728 .seealso: PCBDDC, PCBDDCSetCoarseningRatio() 729 @*/ 730 PetscErrorCode PCBDDCSetLevels(PC pc,PetscInt levels) 731 { 732 PetscErrorCode ierr; 733 734 PetscFunctionBegin; 735 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 736 PetscValidLogicalCollectiveInt(pc,levels,2); 737 ierr = PetscTryMethod(pc,"PCBDDCSetLevels_C",(PC,PetscInt),(pc,levels));CHKERRQ(ierr); 738 PetscFunctionReturn(0); 739 } 740 741 static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc,IS DirichletBoundaries) 742 { 743 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 744 PetscBool isequal = PETSC_FALSE; 745 PetscErrorCode ierr; 746 747 PetscFunctionBegin; 748 ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr); 749 if (pcbddc->DirichletBoundaries) { 750 ierr = ISEqual(DirichletBoundaries,pcbddc->DirichletBoundaries,&isequal);CHKERRQ(ierr); 751 } 752 /* last user setting takes precedence -> destroy any other customization */ 753 ierr = ISDestroy(&pcbddc->DirichletBoundariesLocal);CHKERRQ(ierr); 754 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 755 pcbddc->DirichletBoundaries = DirichletBoundaries; 756 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 757 PetscFunctionReturn(0); 758 } 759 760 /*@ 761 PCBDDCSetDirichletBoundaries - Set IS defining Dirichlet boundaries for the global problem. 762 763 Collective 764 765 Input Parameters: 766 + pc - the preconditioning context 767 - DirichletBoundaries - parallel IS defining the Dirichlet boundaries 768 769 Level: intermediate 770 771 Notes: 772 Provide the information if you used MatZeroRows/Columns routines. Any process can list any global node 773 774 .seealso: PCBDDC, PCBDDCSetDirichletBoundariesLocal(), MatZeroRows(), MatZeroRowsColumns() 775 @*/ 776 PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc,IS DirichletBoundaries) 777 { 778 PetscErrorCode ierr; 779 780 PetscFunctionBegin; 781 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 782 PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2); 783 PetscCheckSameComm(pc,1,DirichletBoundaries,2); 784 ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundaries_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr); 785 PetscFunctionReturn(0); 786 } 787 788 static PetscErrorCode PCBDDCSetDirichletBoundariesLocal_BDDC(PC pc,IS DirichletBoundaries) 789 { 790 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 791 PetscBool isequal = PETSC_FALSE; 792 PetscErrorCode ierr; 793 794 PetscFunctionBegin; 795 ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr); 796 if (pcbddc->DirichletBoundariesLocal) { 797 ierr = ISEqual(DirichletBoundaries,pcbddc->DirichletBoundariesLocal,&isequal);CHKERRQ(ierr); 798 } 799 /* last user setting takes precedence -> destroy any other customization */ 800 ierr = ISDestroy(&pcbddc->DirichletBoundariesLocal);CHKERRQ(ierr); 801 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 802 pcbddc->DirichletBoundariesLocal = DirichletBoundaries; 803 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 804 PetscFunctionReturn(0); 805 } 806 807 /*@ 808 PCBDDCSetDirichletBoundariesLocal - Set IS defining Dirichlet boundaries for the global problem in local ordering. 809 810 Collective 811 812 Input Parameters: 813 + pc - the preconditioning context 814 - DirichletBoundaries - parallel IS defining the Dirichlet boundaries (in local ordering) 815 816 Level: intermediate 817 818 Notes: 819 820 .seealso: PCBDDC, PCBDDCSetDirichletBoundaries(), MatZeroRows(), MatZeroRowsColumns() 821 @*/ 822 PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc,IS DirichletBoundaries) 823 { 824 PetscErrorCode ierr; 825 826 PetscFunctionBegin; 827 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 828 PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2); 829 PetscCheckSameComm(pc,1,DirichletBoundaries,2); 830 ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundariesLocal_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr); 831 PetscFunctionReturn(0); 832 } 833 834 static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc,IS NeumannBoundaries) 835 { 836 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 837 PetscBool isequal = PETSC_FALSE; 838 PetscErrorCode ierr; 839 840 PetscFunctionBegin; 841 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 842 if (pcbddc->NeumannBoundaries) { 843 ierr = ISEqual(NeumannBoundaries,pcbddc->NeumannBoundaries,&isequal);CHKERRQ(ierr); 844 } 845 /* last user setting takes precedence -> destroy any other customization */ 846 ierr = ISDestroy(&pcbddc->NeumannBoundariesLocal);CHKERRQ(ierr); 847 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 848 pcbddc->NeumannBoundaries = NeumannBoundaries; 849 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 850 PetscFunctionReturn(0); 851 } 852 853 /*@ 854 PCBDDCSetNeumannBoundaries - Set IS defining Neumann boundaries for the global problem. 855 856 Collective 857 858 Input Parameters: 859 + pc - the preconditioning context 860 - NeumannBoundaries - parallel IS defining the Neumann boundaries 861 862 Level: intermediate 863 864 Notes: 865 Any process can list any global node 866 867 .seealso: PCBDDC, PCBDDCSetNeumannBoundariesLocal() 868 @*/ 869 PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc,IS NeumannBoundaries) 870 { 871 PetscErrorCode ierr; 872 873 PetscFunctionBegin; 874 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 875 PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2); 876 PetscCheckSameComm(pc,1,NeumannBoundaries,2); 877 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundaries_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 878 PetscFunctionReturn(0); 879 } 880 881 static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc,IS NeumannBoundaries) 882 { 883 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 884 PetscBool isequal = PETSC_FALSE; 885 PetscErrorCode ierr; 886 887 PetscFunctionBegin; 888 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 889 if (pcbddc->NeumannBoundariesLocal) { 890 ierr = ISEqual(NeumannBoundaries,pcbddc->NeumannBoundariesLocal,&isequal);CHKERRQ(ierr); 891 } 892 /* last user setting takes precedence -> destroy any other customization */ 893 ierr = ISDestroy(&pcbddc->NeumannBoundariesLocal);CHKERRQ(ierr); 894 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 895 pcbddc->NeumannBoundariesLocal = NeumannBoundaries; 896 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 897 PetscFunctionReturn(0); 898 } 899 900 /*@ 901 PCBDDCSetNeumannBoundariesLocal - Set IS defining Neumann boundaries for the global problem in local ordering. 902 903 Collective 904 905 Input Parameters: 906 + pc - the preconditioning context 907 - NeumannBoundaries - parallel IS defining the subdomain part of Neumann boundaries (in local ordering) 908 909 Level: intermediate 910 911 Notes: 912 913 .seealso: PCBDDC, PCBDDCSetNeumannBoundaries() 914 @*/ 915 PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc,IS NeumannBoundaries) 916 { 917 PetscErrorCode ierr; 918 919 PetscFunctionBegin; 920 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 921 PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2); 922 PetscCheckSameComm(pc,1,NeumannBoundaries,2); 923 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundariesLocal_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 924 PetscFunctionReturn(0); 925 } 926 927 static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc,IS *DirichletBoundaries) 928 { 929 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 930 931 PetscFunctionBegin; 932 *DirichletBoundaries = pcbddc->DirichletBoundaries; 933 PetscFunctionReturn(0); 934 } 935 936 /*@ 937 PCBDDCGetDirichletBoundaries - Get parallel IS for Dirichlet boundaries 938 939 Collective 940 941 Input Parameters: 942 . pc - the preconditioning context 943 944 Output Parameters: 945 . DirichletBoundaries - index set defining the Dirichlet boundaries 946 947 Level: intermediate 948 949 Notes: 950 The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetDirichletBoundaries 951 952 .seealso: PCBDDC 953 @*/ 954 PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc,IS *DirichletBoundaries) 955 { 956 PetscErrorCode ierr; 957 958 PetscFunctionBegin; 959 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 960 ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundaries_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr); 961 PetscFunctionReturn(0); 962 } 963 964 static PetscErrorCode PCBDDCGetDirichletBoundariesLocal_BDDC(PC pc,IS *DirichletBoundaries) 965 { 966 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 967 968 PetscFunctionBegin; 969 *DirichletBoundaries = pcbddc->DirichletBoundariesLocal; 970 PetscFunctionReturn(0); 971 } 972 973 /*@ 974 PCBDDCGetDirichletBoundariesLocal - Get parallel IS for Dirichlet boundaries (in local ordering) 975 976 Collective 977 978 Input Parameters: 979 . pc - the preconditioning context 980 981 Output Parameters: 982 . DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries 983 984 Level: intermediate 985 986 Notes: 987 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). 988 In the latter case, the IS will be available after PCSetUp. 989 990 .seealso: PCBDDC 991 @*/ 992 PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc,IS *DirichletBoundaries) 993 { 994 PetscErrorCode ierr; 995 996 PetscFunctionBegin; 997 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 998 ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundariesLocal_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr); 999 PetscFunctionReturn(0); 1000 } 1001 1002 static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc,IS *NeumannBoundaries) 1003 { 1004 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1005 1006 PetscFunctionBegin; 1007 *NeumannBoundaries = pcbddc->NeumannBoundaries; 1008 PetscFunctionReturn(0); 1009 } 1010 1011 /*@ 1012 PCBDDCGetNeumannBoundaries - Get parallel IS for Neumann boundaries 1013 1014 Collective 1015 1016 Input Parameters: 1017 . pc - the preconditioning context 1018 1019 Output Parameters: 1020 . NeumannBoundaries - index set defining the Neumann boundaries 1021 1022 Level: intermediate 1023 1024 Notes: 1025 The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetNeumannBoundaries 1026 1027 .seealso: PCBDDC 1028 @*/ 1029 PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc,IS *NeumannBoundaries) 1030 { 1031 PetscErrorCode ierr; 1032 1033 PetscFunctionBegin; 1034 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1035 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundaries_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 1036 PetscFunctionReturn(0); 1037 } 1038 1039 static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc,IS *NeumannBoundaries) 1040 { 1041 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1042 1043 PetscFunctionBegin; 1044 *NeumannBoundaries = pcbddc->NeumannBoundariesLocal; 1045 PetscFunctionReturn(0); 1046 } 1047 1048 /*@ 1049 PCBDDCGetNeumannBoundariesLocal - Get parallel IS for Neumann boundaries (in local ordering) 1050 1051 Collective 1052 1053 Input Parameters: 1054 . pc - the preconditioning context 1055 1056 Output Parameters: 1057 . NeumannBoundaries - index set defining the subdomain part of Neumann boundaries 1058 1059 Level: intermediate 1060 1061 Notes: 1062 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). 1063 In the latter case, the IS will be available after PCSetUp. 1064 1065 .seealso: PCBDDC 1066 @*/ 1067 PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc,IS *NeumannBoundaries) 1068 { 1069 PetscErrorCode ierr; 1070 1071 PetscFunctionBegin; 1072 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1073 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundariesLocal_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 1074 PetscFunctionReturn(0); 1075 } 1076 1077 static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 1078 { 1079 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1080 PCBDDCGraph mat_graph = pcbddc->mat_graph; 1081 PetscBool same_data = PETSC_FALSE; 1082 PetscErrorCode ierr; 1083 1084 PetscFunctionBegin; 1085 if (!nvtxs) { 1086 if (copymode == PETSC_OWN_POINTER) { 1087 ierr = PetscFree(xadj);CHKERRQ(ierr); 1088 ierr = PetscFree(adjncy);CHKERRQ(ierr); 1089 } 1090 ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr); 1091 PetscFunctionReturn(0); 1092 } 1093 if (mat_graph->nvtxs == nvtxs && mat_graph->freecsr) { /* we own the data */ 1094 if (mat_graph->xadj == xadj && mat_graph->adjncy == adjncy) same_data = PETSC_TRUE; 1095 if (!same_data && mat_graph->xadj[nvtxs] == xadj[nvtxs]) { 1096 ierr = PetscArraycmp(xadj,mat_graph->xadj,nvtxs+1,&same_data);CHKERRQ(ierr); 1097 if (same_data) { 1098 ierr = PetscArraycmp(adjncy,mat_graph->adjncy,xadj[nvtxs],&same_data);CHKERRQ(ierr); 1099 } 1100 } 1101 } 1102 if (!same_data) { 1103 /* free old CSR */ 1104 ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr); 1105 /* get CSR into graph structure */ 1106 if (copymode == PETSC_COPY_VALUES) { 1107 ierr = PetscMalloc1(nvtxs+1,&mat_graph->xadj);CHKERRQ(ierr); 1108 ierr = PetscMalloc1(xadj[nvtxs],&mat_graph->adjncy);CHKERRQ(ierr); 1109 ierr = PetscArraycpy(mat_graph->xadj,xadj,nvtxs+1);CHKERRQ(ierr); 1110 ierr = PetscArraycpy(mat_graph->adjncy,adjncy,xadj[nvtxs]);CHKERRQ(ierr); 1111 mat_graph->freecsr = PETSC_TRUE; 1112 } else if (copymode == PETSC_OWN_POINTER) { 1113 mat_graph->xadj = (PetscInt*)xadj; 1114 mat_graph->adjncy = (PetscInt*)adjncy; 1115 mat_graph->freecsr = PETSC_TRUE; 1116 } else if (copymode == PETSC_USE_POINTER) { 1117 mat_graph->xadj = (PetscInt*)xadj; 1118 mat_graph->adjncy = (PetscInt*)adjncy; 1119 mat_graph->freecsr = PETSC_FALSE; 1120 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported copy mode %D",copymode); 1121 mat_graph->nvtxs_csr = nvtxs; 1122 pcbddc->recompute_topography = PETSC_TRUE; 1123 } 1124 PetscFunctionReturn(0); 1125 } 1126 1127 /*@ 1128 PCBDDCSetLocalAdjacencyGraph - Set adjacency structure (CSR graph) of the local degrees of freedom. 1129 1130 Not collective 1131 1132 Input Parameters: 1133 + pc - the preconditioning context. 1134 . nvtxs - number of local vertices of the graph (i.e., the number of local dofs). 1135 . xadj, adjncy - the connectivity of the dofs in CSR format. 1136 - copymode - supported modes are PETSC_COPY_VALUES, PETSC_USE_POINTER or PETSC_OWN_POINTER. 1137 1138 Level: intermediate 1139 1140 Notes: 1141 A dof is considered connected with all local dofs if xadj[dof+1]-xadj[dof] == 1 and adjncy[xadj[dof]] is negative. 1142 1143 .seealso: PCBDDC,PetscCopyMode 1144 @*/ 1145 PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc,PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 1146 { 1147 void (*f)(void) = NULL; 1148 PetscErrorCode ierr; 1149 1150 PetscFunctionBegin; 1151 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1152 if (nvtxs) { 1153 PetscValidIntPointer(xadj,3); 1154 if (xadj[nvtxs]) PetscValidIntPointer(adjncy,4); 1155 } 1156 ierr = PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));CHKERRQ(ierr); 1157 /* free arrays if PCBDDC is not the PC type */ 1158 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);CHKERRQ(ierr); 1159 if (!f && copymode == PETSC_OWN_POINTER) { 1160 ierr = PetscFree(xadj);CHKERRQ(ierr); 1161 ierr = PetscFree(adjncy);CHKERRQ(ierr); 1162 } 1163 PetscFunctionReturn(0); 1164 } 1165 1166 static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 1167 { 1168 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1169 PetscInt i; 1170 PetscBool isequal = PETSC_FALSE; 1171 PetscErrorCode ierr; 1172 1173 PetscFunctionBegin; 1174 if (pcbddc->n_ISForDofsLocal == n_is) { 1175 for (i=0;i<n_is;i++) { 1176 PetscBool isequalt; 1177 ierr = ISEqual(ISForDofs[i],pcbddc->ISForDofsLocal[i],&isequalt);CHKERRQ(ierr); 1178 if (!isequalt) break; 1179 } 1180 if (i == n_is) isequal = PETSC_TRUE; 1181 } 1182 for (i=0;i<n_is;i++) { 1183 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 1184 } 1185 /* Destroy ISes if they were already set */ 1186 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 1187 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 1188 } 1189 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 1190 /* last user setting takes precedence -> destroy any other customization */ 1191 for (i=0;i<pcbddc->n_ISForDofs;i++) { 1192 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 1193 } 1194 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 1195 pcbddc->n_ISForDofs = 0; 1196 /* allocate space then set */ 1197 if (n_is) { 1198 ierr = PetscMalloc1(n_is,&pcbddc->ISForDofsLocal);CHKERRQ(ierr); 1199 } 1200 for (i=0;i<n_is;i++) { 1201 pcbddc->ISForDofsLocal[i] = ISForDofs[i]; 1202 } 1203 pcbddc->n_ISForDofsLocal = n_is; 1204 if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE; 1205 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 1206 PetscFunctionReturn(0); 1207 } 1208 1209 /*@ 1210 PCBDDCSetDofsSplittingLocal - Set index sets defining fields of the local subdomain matrix 1211 1212 Collective 1213 1214 Input Parameters: 1215 + pc - the preconditioning context 1216 . n_is - number of index sets defining the fields 1217 - ISForDofs - array of IS describing the fields in local ordering 1218 1219 Level: intermediate 1220 1221 Notes: 1222 n_is should be the same among processes. Not all nodes need to be listed: unlisted nodes will belong to the complement field. 1223 1224 .seealso: PCBDDC 1225 @*/ 1226 PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc,PetscInt n_is, IS ISForDofs[]) 1227 { 1228 PetscInt i; 1229 PetscErrorCode ierr; 1230 1231 PetscFunctionBegin; 1232 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1233 PetscValidLogicalCollectiveInt(pc,n_is,2); 1234 for (i=0;i<n_is;i++) { 1235 PetscCheckSameComm(pc,1,ISForDofs[i],3); 1236 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 1237 } 1238 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplittingLocal_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 1239 PetscFunctionReturn(0); 1240 } 1241 1242 static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 1243 { 1244 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1245 PetscInt i; 1246 PetscBool isequal = PETSC_FALSE; 1247 PetscErrorCode ierr; 1248 1249 PetscFunctionBegin; 1250 if (pcbddc->n_ISForDofs == n_is) { 1251 for (i=0;i<n_is;i++) { 1252 PetscBool isequalt; 1253 ierr = ISEqual(ISForDofs[i],pcbddc->ISForDofs[i],&isequalt);CHKERRQ(ierr); 1254 if (!isequalt) break; 1255 } 1256 if (i == n_is) isequal = PETSC_TRUE; 1257 } 1258 for (i=0;i<n_is;i++) { 1259 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 1260 } 1261 /* Destroy ISes if they were already set */ 1262 for (i=0;i<pcbddc->n_ISForDofs;i++) { 1263 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 1264 } 1265 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 1266 /* last user setting takes precedence -> destroy any other customization */ 1267 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 1268 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 1269 } 1270 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 1271 pcbddc->n_ISForDofsLocal = 0; 1272 /* allocate space then set */ 1273 if (n_is) { 1274 ierr = PetscMalloc1(n_is,&pcbddc->ISForDofs);CHKERRQ(ierr); 1275 } 1276 for (i=0;i<n_is;i++) { 1277 pcbddc->ISForDofs[i] = ISForDofs[i]; 1278 } 1279 pcbddc->n_ISForDofs = n_is; 1280 if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE; 1281 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 1282 PetscFunctionReturn(0); 1283 } 1284 1285 /*@ 1286 PCBDDCSetDofsSplitting - Set index sets defining fields of the global matrix 1287 1288 Collective 1289 1290 Input Parameters: 1291 + pc - the preconditioning context 1292 . n_is - number of index sets defining the fields 1293 - ISForDofs - array of IS describing the fields in global ordering 1294 1295 Level: intermediate 1296 1297 Notes: 1298 Any process can list any global node. Not all nodes need to be listed: unlisted nodes will belong to the complement field. 1299 1300 .seealso: PCBDDC 1301 @*/ 1302 PetscErrorCode PCBDDCSetDofsSplitting(PC pc,PetscInt n_is, IS ISForDofs[]) 1303 { 1304 PetscInt i; 1305 PetscErrorCode ierr; 1306 1307 PetscFunctionBegin; 1308 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1309 PetscValidLogicalCollectiveInt(pc,n_is,2); 1310 for (i=0;i<n_is;i++) { 1311 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 1312 PetscCheckSameComm(pc,1,ISForDofs[i],3); 1313 } 1314 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 1315 PetscFunctionReturn(0); 1316 } 1317 1318 /* 1319 PCPreSolve_BDDC - Changes the right hand side and (if necessary) the initial 1320 guess if a transformation of basis approach has been selected. 1321 1322 Input Parameter: 1323 + pc - the preconditioner context 1324 1325 Application Interface Routine: PCPreSolve() 1326 1327 Notes: 1328 The interface routine PCPreSolve() is not usually called directly by 1329 the user, but instead is called by KSPSolve(). 1330 */ 1331 static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 1332 { 1333 PetscErrorCode ierr; 1334 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1335 PC_IS *pcis = (PC_IS*)(pc->data); 1336 Vec used_vec; 1337 PetscBool iscg = PETSC_FALSE, save_rhs = PETSC_TRUE, benign_correction_computed; 1338 1339 PetscFunctionBegin; 1340 /* if we are working with CG, one dirichlet solve can be avoided during Krylov iterations */ 1341 if (ksp) { 1342 PetscBool isgroppcg, ispipecg, ispipelcg, ispipecgrr; 1343 1344 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPCG,&iscg);CHKERRQ(ierr); 1345 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPGROPPCG,&isgroppcg);CHKERRQ(ierr); 1346 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECG,&ispipecg);CHKERRQ(ierr); 1347 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECG,&ispipelcg);CHKERRQ(ierr); 1348 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECGRR,&ispipecgrr);CHKERRQ(ierr); 1349 iscg = (PetscBool)(iscg || isgroppcg || ispipecg || ispipelcg || ispipecgrr); 1350 if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || !iscg || pc->mat != pc->pmat) { 1351 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 1352 } 1353 } 1354 if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || pc->mat != pc->pmat) { 1355 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 1356 } 1357 1358 /* Creates parallel work vectors used in presolve */ 1359 if (!pcbddc->original_rhs) { 1360 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 1361 } 1362 if (!pcbddc->temp_solution) { 1363 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);CHKERRQ(ierr); 1364 } 1365 1366 pcbddc->temp_solution_used = PETSC_FALSE; 1367 if (x) { 1368 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr); 1369 used_vec = x; 1370 } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */ 1371 ierr = PetscObjectReference((PetscObject)pcbddc->temp_solution);CHKERRQ(ierr); 1372 used_vec = pcbddc->temp_solution; 1373 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 1374 pcbddc->temp_solution_used = PETSC_TRUE; 1375 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1376 save_rhs = PETSC_FALSE; 1377 pcbddc->eliminate_dirdofs = PETSC_TRUE; 1378 } 1379 1380 /* hack into ksp data structure since PCPreSolve comes earlier than setting to zero the guess in src/ksp/ksp/interface/itfunc.c */ 1381 if (ksp) { 1382 /* store the flag for the initial guess since it will be restored back during PCPostSolve_BDDC */ 1383 ierr = KSPGetInitialGuessNonzero(ksp,&pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 1384 if (!pcbddc->ksp_guess_nonzero) { 1385 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 1386 } 1387 } 1388 1389 pcbddc->rhs_change = PETSC_FALSE; 1390 /* Take into account zeroed rows -> change rhs and store solution removed */ 1391 if (rhs && pcbddc->eliminate_dirdofs) { 1392 IS dirIS = NULL; 1393 1394 /* DirichletBoundariesLocal may not be consistent among neighbours; gets a dirichlet dofs IS from graph (may be cached) */ 1395 ierr = PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph,&dirIS);CHKERRQ(ierr); 1396 if (dirIS) { 1397 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 1398 PetscInt dirsize,i,*is_indices; 1399 PetscScalar *array_x; 1400 const PetscScalar *array_diagonal; 1401 1402 ierr = MatGetDiagonal(pc->pmat,pcis->vec1_global);CHKERRQ(ierr); 1403 ierr = VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);CHKERRQ(ierr); 1404 ierr = VecScatterBegin(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1405 ierr = VecScatterEnd(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1406 ierr = VecScatterBegin(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1407 ierr = VecScatterEnd(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1408 ierr = ISGetLocalSize(dirIS,&dirsize);CHKERRQ(ierr); 1409 ierr = VecGetArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 1410 ierr = VecGetArrayRead(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 1411 ierr = ISGetIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1412 for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]]; 1413 ierr = ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1414 ierr = VecRestoreArrayRead(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 1415 ierr = VecRestoreArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 1416 ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1417 ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1418 pcbddc->rhs_change = PETSC_TRUE; 1419 ierr = ISDestroy(&dirIS);CHKERRQ(ierr); 1420 } 1421 } 1422 1423 /* remove the computed solution or the initial guess from the rhs */ 1424 if (pcbddc->rhs_change || (ksp && pcbddc->ksp_guess_nonzero)) { 1425 /* save the original rhs */ 1426 if (save_rhs) { 1427 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1428 save_rhs = PETSC_FALSE; 1429 } 1430 pcbddc->rhs_change = PETSC_TRUE; 1431 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 1432 ierr = MatMultAdd(pc->mat,used_vec,pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1433 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 1434 ierr = VecCopy(used_vec,pcbddc->temp_solution);CHKERRQ(ierr); 1435 pcbddc->temp_solution_used = PETSC_TRUE; 1436 if (ksp) { 1437 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_FALSE);CHKERRQ(ierr); 1438 } 1439 } 1440 ierr = VecDestroy(&used_vec);CHKERRQ(ierr); 1441 1442 /* compute initial vector in benign space if needed 1443 and remove non-benign solution from the rhs */ 1444 benign_correction_computed = PETSC_FALSE; 1445 if (rhs && pcbddc->benign_compute_correction && (pcbddc->benign_have_null || pcbddc->benign_apply_coarse_only)) { 1446 /* compute u^*_h using ideas similar to those in Xuemin Tu's PhD thesis (see Section 4.8.1) 1447 Recursively apply BDDC in the multilevel case */ 1448 if (!pcbddc->benign_vec) { 1449 ierr = VecDuplicate(rhs,&pcbddc->benign_vec);CHKERRQ(ierr); 1450 } 1451 /* keep applying coarse solver unless we no longer have benign subdomains */ 1452 pcbddc->benign_apply_coarse_only = pcbddc->benign_have_null ? PETSC_TRUE : PETSC_FALSE; 1453 if (!pcbddc->benign_skip_correction) { 1454 ierr = PCApply_BDDC(pc,rhs,pcbddc->benign_vec);CHKERRQ(ierr); 1455 benign_correction_computed = PETSC_TRUE; 1456 if (pcbddc->temp_solution_used) { 1457 ierr = VecAXPY(pcbddc->temp_solution,1.0,pcbddc->benign_vec);CHKERRQ(ierr); 1458 } 1459 ierr = VecScale(pcbddc->benign_vec,-1.0);CHKERRQ(ierr); 1460 /* store the original rhs if not done earlier */ 1461 if (save_rhs) { 1462 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1463 } 1464 if (pcbddc->rhs_change) { 1465 ierr = MatMultAdd(pc->mat,pcbddc->benign_vec,rhs,rhs);CHKERRQ(ierr); 1466 } else { 1467 ierr = MatMultAdd(pc->mat,pcbddc->benign_vec,pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1468 } 1469 pcbddc->rhs_change = PETSC_TRUE; 1470 } 1471 pcbddc->benign_apply_coarse_only = PETSC_FALSE; 1472 } else { 1473 ierr = VecDestroy(&pcbddc->benign_vec);CHKERRQ(ierr); 1474 } 1475 1476 /* dbg output */ 1477 if (pcbddc->dbg_flag && benign_correction_computed) { 1478 Vec v; 1479 1480 ierr = VecDuplicate(pcis->vec1_global,&v);CHKERRQ(ierr); 1481 if (pcbddc->ChangeOfBasisMatrix) { 1482 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,rhs,v);CHKERRQ(ierr); 1483 } else { 1484 ierr = VecCopy(rhs,v);CHKERRQ(ierr); 1485 } 1486 ierr = PCBDDCBenignGetOrSetP0(pc,v,PETSC_TRUE);CHKERRQ(ierr); 1487 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"LEVEL %D: is the correction benign?\n",pcbddc->current_level);CHKERRQ(ierr); 1488 ierr = PetscScalarView(pcbddc->benign_n,pcbddc->benign_p0,pcbddc->dbg_viewer);CHKERRQ(ierr); 1489 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 1490 ierr = VecDestroy(&v);CHKERRQ(ierr); 1491 } 1492 1493 /* set initial guess if using PCG */ 1494 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 1495 if (x && pcbddc->use_exact_dirichlet_trick) { 1496 ierr = VecSet(x,0.0);CHKERRQ(ierr); 1497 if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) { 1498 if (benign_correction_computed) { /* we have already saved the changed rhs */ 1499 ierr = VecLockReadPop(pcis->vec1_global);CHKERRQ(ierr); 1500 } else { 1501 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,rhs,pcis->vec1_global);CHKERRQ(ierr); 1502 } 1503 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_global,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1504 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_global,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1505 } else { 1506 ierr = VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1507 ierr = VecScatterEnd(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1508 } 1509 ierr = PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 1510 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1511 ierr = PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 1512 ierr = KSPCheckSolve(pcbddc->ksp_D,pc,pcis->vec2_D);CHKERRQ(ierr); 1513 if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) { 1514 ierr = VecSet(pcis->vec1_global,0.);CHKERRQ(ierr); 1515 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1516 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1517 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcis->vec1_global,x);CHKERRQ(ierr); 1518 } else { 1519 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1520 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1521 } 1522 if (ksp) { 1523 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 1524 } 1525 pcbddc->exact_dirichlet_trick_app = PETSC_TRUE; 1526 } else if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior && benign_correction_computed && pcbddc->use_exact_dirichlet_trick) { 1527 ierr = VecLockReadPop(pcis->vec1_global);CHKERRQ(ierr); 1528 } 1529 PetscFunctionReturn(0); 1530 } 1531 1532 /* 1533 PCPostSolve_BDDC - Changes the computed solution if a transformation of basis 1534 approach has been selected. Also, restores rhs to its original state. 1535 1536 Input Parameter: 1537 + pc - the preconditioner context 1538 1539 Application Interface Routine: PCPostSolve() 1540 1541 Notes: 1542 The interface routine PCPostSolve() is not usually called directly by 1543 the user, but instead is called by KSPSolve(). 1544 */ 1545 static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 1546 { 1547 PetscErrorCode ierr; 1548 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1549 1550 PetscFunctionBegin; 1551 /* add solution removed in presolve */ 1552 if (x && pcbddc->rhs_change) { 1553 if (pcbddc->temp_solution_used) { 1554 ierr = VecAXPY(x,1.0,pcbddc->temp_solution);CHKERRQ(ierr); 1555 } else if (pcbddc->benign_compute_correction && pcbddc->benign_vec) { 1556 ierr = VecAXPY(x,-1.0,pcbddc->benign_vec);CHKERRQ(ierr); 1557 } 1558 /* restore to original state (not for FETI-DP) */ 1559 if (ksp) pcbddc->temp_solution_used = PETSC_FALSE; 1560 } 1561 1562 /* restore rhs to its original state (not needed for FETI-DP) */ 1563 if (rhs && pcbddc->rhs_change) { 1564 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1565 pcbddc->rhs_change = PETSC_FALSE; 1566 } 1567 /* restore ksp guess state */ 1568 if (ksp) { 1569 ierr = KSPSetInitialGuessNonzero(ksp,pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 1570 /* reset flag for exact dirichlet trick */ 1571 pcbddc->exact_dirichlet_trick_app = PETSC_FALSE; 1572 } 1573 PetscFunctionReturn(0); 1574 } 1575 1576 /* 1577 PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner 1578 by setting data structures and options. 1579 1580 Input Parameter: 1581 + pc - the preconditioner context 1582 1583 Application Interface Routine: PCSetUp() 1584 1585 Notes: 1586 The interface routine PCSetUp() is not usually called directly by 1587 the user, but instead is called by PCApply() if necessary. 1588 */ 1589 PetscErrorCode PCSetUp_BDDC(PC pc) 1590 { 1591 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 1592 PCBDDCSubSchurs sub_schurs; 1593 Mat_IS* matis; 1594 MatNullSpace nearnullspace; 1595 Mat lA; 1596 IS lP,zerodiag = NULL; 1597 PetscInt nrows,ncols; 1598 PetscMPIInt size; 1599 PetscBool computesubschurs; 1600 PetscBool computeconstraintsmatrix; 1601 PetscBool new_nearnullspace_provided,ismatis,rl; 1602 PetscErrorCode ierr; 1603 1604 PetscFunctionBegin; 1605 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&ismatis);CHKERRQ(ierr); 1606 PetscCheckFalse(!ismatis,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PCBDDC preconditioner requires matrix of type MATIS"); 1607 ierr = MatGetSize(pc->pmat,&nrows,&ncols);CHKERRQ(ierr); 1608 PetscCheckFalse(nrows != ncols,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCBDDC preconditioner requires a square preconditioning matrix"); 1609 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRMPI(ierr); 1610 1611 matis = (Mat_IS*)pc->pmat->data; 1612 /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */ 1613 /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup 1614 Also, BDDC builds its own KSP for the Dirichlet problem */ 1615 rl = pcbddc->recompute_topography; 1616 if (!pc->setupcalled || pc->flag == DIFFERENT_NONZERO_PATTERN) rl = PETSC_TRUE; 1617 ierr = MPIU_Allreduce(&rl,&pcbddc->recompute_topography,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)pc));CHKERRMPI(ierr); 1618 if (pcbddc->recompute_topography) { 1619 pcbddc->graphanalyzed = PETSC_FALSE; 1620 computeconstraintsmatrix = PETSC_TRUE; 1621 } else { 1622 computeconstraintsmatrix = PETSC_FALSE; 1623 } 1624 1625 /* check parameters' compatibility */ 1626 if (!pcbddc->use_deluxe_scaling) pcbddc->deluxe_zerorows = PETSC_FALSE; 1627 pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_threshold[0] != 0.0 || pcbddc->adaptive_threshold[1] != 0.0); 1628 pcbddc->use_deluxe_scaling = (PetscBool)(pcbddc->use_deluxe_scaling && size > 1); 1629 pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_selection && size > 1); 1630 pcbddc->adaptive_userdefined = (PetscBool)(pcbddc->adaptive_selection && pcbddc->adaptive_userdefined); 1631 if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE; 1632 1633 computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling); 1634 1635 /* activate all connected components if the netflux has been requested */ 1636 if (pcbddc->compute_nonetflux) { 1637 pcbddc->use_vertices = PETSC_TRUE; 1638 pcbddc->use_edges = PETSC_TRUE; 1639 pcbddc->use_faces = PETSC_TRUE; 1640 } 1641 1642 /* Get stdout for dbg */ 1643 if (pcbddc->dbg_flag) { 1644 if (!pcbddc->dbg_viewer) { 1645 pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)); 1646 } 1647 ierr = PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer);CHKERRQ(ierr); 1648 ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1649 } 1650 1651 /* process topology information */ 1652 ierr = PetscLogEventBegin(PC_BDDC_Topology[pcbddc->current_level],pc,0,0,0);CHKERRQ(ierr); 1653 if (pcbddc->recompute_topography) { 1654 ierr = PCBDDCComputeLocalTopologyInfo(pc);CHKERRQ(ierr); 1655 if (pcbddc->discretegradient) { 1656 ierr = PCBDDCNedelecSupport(pc);CHKERRQ(ierr); 1657 } 1658 } 1659 if (pcbddc->corner_selected) pcbddc->use_vertices = PETSC_TRUE; 1660 1661 /* change basis if requested by the user */ 1662 if (pcbddc->user_ChangeOfBasisMatrix) { 1663 /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */ 1664 pcbddc->use_change_of_basis = PETSC_FALSE; 1665 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 1666 } else { 1667 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1668 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1669 pcbddc->local_mat = matis->A; 1670 } 1671 1672 /* 1673 Compute change of basis on local pressures (aka zerodiag dofs) with the benign trick 1674 This should come earlier then PCISSetUp for extracting the correct subdomain matrices 1675 */ 1676 ierr = PCBDDCBenignShellMat(pc,PETSC_TRUE);CHKERRQ(ierr); 1677 if (pcbddc->benign_saddle_point) { 1678 PC_IS* pcis = (PC_IS*)pc->data; 1679 1680 if (pcbddc->user_ChangeOfBasisMatrix || pcbddc->use_change_of_basis || !computesubschurs) pcbddc->benign_change_explicit = PETSC_TRUE; 1681 /* detect local saddle point and change the basis in pcbddc->local_mat */ 1682 ierr = PCBDDCBenignDetectSaddlePoint(pc,(PetscBool)(!pcbddc->recompute_topography),&zerodiag);CHKERRQ(ierr); 1683 /* pop B0 mat from local mat */ 1684 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1685 /* give pcis a hint to not reuse submatrices during PCISCreate */ 1686 if (pc->flag == SAME_NONZERO_PATTERN && pcis->reusesubmatrices == PETSC_TRUE) { 1687 if (pcbddc->benign_n && (pcbddc->benign_change_explicit || pcbddc->dbg_flag)) { 1688 pcis->reusesubmatrices = PETSC_FALSE; 1689 } else { 1690 pcis->reusesubmatrices = PETSC_TRUE; 1691 } 1692 } else { 1693 pcis->reusesubmatrices = PETSC_FALSE; 1694 } 1695 } 1696 1697 /* propagate relevant information */ 1698 if (matis->A->symmetric_set) { 1699 ierr = MatSetOption(pcbddc->local_mat,MAT_SYMMETRIC,matis->A->symmetric);CHKERRQ(ierr); 1700 } 1701 if (matis->A->spd_set) { 1702 ierr = MatSetOption(pcbddc->local_mat,MAT_SPD,matis->A->spd);CHKERRQ(ierr); 1703 } 1704 1705 /* Set up all the "iterative substructuring" common block without computing solvers */ 1706 { 1707 Mat temp_mat; 1708 1709 temp_mat = matis->A; 1710 matis->A = pcbddc->local_mat; 1711 ierr = PCISSetUp(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1712 pcbddc->local_mat = matis->A; 1713 matis->A = temp_mat; 1714 } 1715 1716 /* Analyze interface */ 1717 if (!pcbddc->graphanalyzed) { 1718 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 1719 computeconstraintsmatrix = PETSC_TRUE; 1720 if (pcbddc->adaptive_selection && !pcbddc->use_deluxe_scaling && !pcbddc->mat_graph->twodim) { 1721 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot compute the adaptive primal space for a problem with 3D edges without deluxe scaling"); 1722 } 1723 if (pcbddc->compute_nonetflux) { 1724 MatNullSpace nnfnnsp; 1725 1726 PetscCheckFalse(!pcbddc->divudotp,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Missing divudotp operator"); 1727 ierr = PCBDDCComputeNoNetFlux(pc->pmat,pcbddc->divudotp,pcbddc->divudotp_trans,pcbddc->divudotp_vl2l,pcbddc->mat_graph,&nnfnnsp);CHKERRQ(ierr); 1728 /* TODO what if a nearnullspace is already attached? */ 1729 if (nnfnnsp) { 1730 ierr = MatSetNearNullSpace(pc->pmat,nnfnnsp);CHKERRQ(ierr); 1731 ierr = MatNullSpaceDestroy(&nnfnnsp);CHKERRQ(ierr); 1732 } 1733 } 1734 } 1735 ierr = PetscLogEventEnd(PC_BDDC_Topology[pcbddc->current_level],pc,0,0,0);CHKERRQ(ierr); 1736 1737 /* check existence of a divergence free extension, i.e. 1738 b(v_I,p_0) = 0 for all v_I (raise error if not). 1739 Also, check that PCBDDCBenignGetOrSetP0 works */ 1740 if (pcbddc->benign_saddle_point && pcbddc->dbg_flag > 1) { 1741 ierr = PCBDDCBenignCheck(pc,zerodiag);CHKERRQ(ierr); 1742 } 1743 ierr = ISDestroy(&zerodiag);CHKERRQ(ierr); 1744 1745 /* Setup local dirichlet solver ksp_D and sub_schurs solvers */ 1746 if (computesubschurs && pcbddc->recompute_topography) { 1747 ierr = PCBDDCInitSubSchurs(pc);CHKERRQ(ierr); 1748 } 1749 /* SetUp Scaling operator (scaling matrices could be needed in SubSchursSetUp)*/ 1750 if (!pcbddc->use_deluxe_scaling) { 1751 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1752 } 1753 1754 /* finish setup solvers and do adaptive selection of constraints */ 1755 sub_schurs = pcbddc->sub_schurs; 1756 if (sub_schurs && sub_schurs->schur_explicit) { 1757 if (computesubschurs) { 1758 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1759 } 1760 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1761 } else { 1762 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1763 if (computesubschurs) { 1764 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1765 } 1766 } 1767 if (pcbddc->adaptive_selection) { 1768 ierr = PCBDDCAdaptiveSelection(pc);CHKERRQ(ierr); 1769 computeconstraintsmatrix = PETSC_TRUE; 1770 } 1771 1772 /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */ 1773 new_nearnullspace_provided = PETSC_FALSE; 1774 ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr); 1775 if (pcbddc->onearnullspace) { /* already used nearnullspace */ 1776 if (!nearnullspace) { /* near null space attached to mat has been destroyed */ 1777 new_nearnullspace_provided = PETSC_TRUE; 1778 } else { 1779 /* determine if the two nullspaces are different (should be lightweight) */ 1780 if (nearnullspace != pcbddc->onearnullspace) { 1781 new_nearnullspace_provided = PETSC_TRUE; 1782 } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */ 1783 PetscInt i; 1784 const Vec *nearnullvecs; 1785 PetscObjectState state; 1786 PetscInt nnsp_size; 1787 ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 1788 for (i=0;i<nnsp_size;i++) { 1789 ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr); 1790 if (pcbddc->onearnullvecs_state[i] != state) { 1791 new_nearnullspace_provided = PETSC_TRUE; 1792 break; 1793 } 1794 } 1795 } 1796 } 1797 } else { 1798 if (!nearnullspace) { /* both nearnullspaces are null */ 1799 new_nearnullspace_provided = PETSC_FALSE; 1800 } else { /* nearnullspace attached later */ 1801 new_nearnullspace_provided = PETSC_TRUE; 1802 } 1803 } 1804 1805 /* Setup constraints and related work vectors */ 1806 /* reset primal space flags */ 1807 ierr = PetscLogEventBegin(PC_BDDC_LocalWork[pcbddc->current_level],pc,0,0,0);CHKERRQ(ierr); 1808 pcbddc->new_primal_space = PETSC_FALSE; 1809 pcbddc->new_primal_space_local = PETSC_FALSE; 1810 if (computeconstraintsmatrix || new_nearnullspace_provided) { 1811 /* It also sets the primal space flags */ 1812 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 1813 } 1814 /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */ 1815 ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr); 1816 1817 if (pcbddc->use_change_of_basis) { 1818 PC_IS *pcis = (PC_IS*)(pc->data); 1819 1820 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 1821 if (pcbddc->benign_change) { 1822 ierr = MatDestroy(&pcbddc->benign_B0);CHKERRQ(ierr); 1823 /* pop B0 from pcbddc->local_mat */ 1824 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1825 } 1826 /* get submatrices */ 1827 ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr); 1828 ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr); 1829 ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr); 1830 ierr = MatCreateSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);CHKERRQ(ierr); 1831 ierr = MatCreateSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);CHKERRQ(ierr); 1832 ierr = MatCreateSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);CHKERRQ(ierr); 1833 /* set flag in pcis to not reuse submatrices during PCISCreate */ 1834 pcis->reusesubmatrices = PETSC_FALSE; 1835 } else if (!pcbddc->user_ChangeOfBasisMatrix && !pcbddc->benign_change) { 1836 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1837 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1838 pcbddc->local_mat = matis->A; 1839 } 1840 1841 /* interface pressure block row for B_C */ 1842 ierr = PetscObjectQuery((PetscObject)pc,"__KSPFETIDP_lP" ,(PetscObject*)&lP);CHKERRQ(ierr); 1843 ierr = PetscObjectQuery((PetscObject)pc,"__KSPFETIDP_lA" ,(PetscObject*)&lA);CHKERRQ(ierr); 1844 if (lA && lP) { 1845 PC_IS* pcis = (PC_IS*)pc->data; 1846 Mat B_BI,B_BB,Bt_BI,Bt_BB; 1847 PetscBool issym; 1848 ierr = MatIsSymmetric(lA,PETSC_SMALL,&issym);CHKERRQ(ierr); 1849 if (issym) { 1850 ierr = MatCreateSubMatrix(lA,lP,pcis->is_I_local,MAT_INITIAL_MATRIX,&B_BI);CHKERRQ(ierr); 1851 ierr = MatCreateSubMatrix(lA,lP,pcis->is_B_local,MAT_INITIAL_MATRIX,&B_BB);CHKERRQ(ierr); 1852 ierr = MatCreateTranspose(B_BI,&Bt_BI);CHKERRQ(ierr); 1853 ierr = MatCreateTranspose(B_BB,&Bt_BB);CHKERRQ(ierr); 1854 } else { 1855 ierr = MatCreateSubMatrix(lA,lP,pcis->is_I_local,MAT_INITIAL_MATRIX,&B_BI);CHKERRQ(ierr); 1856 ierr = MatCreateSubMatrix(lA,lP,pcis->is_B_local,MAT_INITIAL_MATRIX,&B_BB);CHKERRQ(ierr); 1857 ierr = MatCreateSubMatrix(lA,pcis->is_I_local,lP,MAT_INITIAL_MATRIX,&Bt_BI);CHKERRQ(ierr); 1858 ierr = MatCreateSubMatrix(lA,pcis->is_B_local,lP,MAT_INITIAL_MATRIX,&Bt_BB);CHKERRQ(ierr); 1859 } 1860 ierr = PetscObjectCompose((PetscObject)pc,"__KSPFETIDP_B_BI",(PetscObject)B_BI);CHKERRQ(ierr); 1861 ierr = PetscObjectCompose((PetscObject)pc,"__KSPFETIDP_B_BB",(PetscObject)B_BB);CHKERRQ(ierr); 1862 ierr = PetscObjectCompose((PetscObject)pc,"__KSPFETIDP_Bt_BI",(PetscObject)Bt_BI);CHKERRQ(ierr); 1863 ierr = PetscObjectCompose((PetscObject)pc,"__KSPFETIDP_Bt_BB",(PetscObject)Bt_BB);CHKERRQ(ierr); 1864 ierr = MatDestroy(&B_BI);CHKERRQ(ierr); 1865 ierr = MatDestroy(&B_BB);CHKERRQ(ierr); 1866 ierr = MatDestroy(&Bt_BI);CHKERRQ(ierr); 1867 ierr = MatDestroy(&Bt_BB);CHKERRQ(ierr); 1868 } 1869 ierr = PetscLogEventEnd(PC_BDDC_LocalWork[pcbddc->current_level],pc,0,0,0);CHKERRQ(ierr); 1870 1871 /* SetUp coarse and local Neumann solvers */ 1872 ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr); 1873 /* SetUp Scaling operator */ 1874 if (pcbddc->use_deluxe_scaling) { 1875 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1876 } 1877 1878 /* mark topography as done */ 1879 pcbddc->recompute_topography = PETSC_FALSE; 1880 1881 /* wrap pcis->A_IB and pcis->A_BI if we did not change explicitly the variables on the pressures */ 1882 ierr = PCBDDCBenignShellMat(pc,PETSC_FALSE);CHKERRQ(ierr); 1883 1884 if (pcbddc->dbg_flag) { 1885 ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1886 ierr = PetscViewerASCIIPopSynchronized(pcbddc->dbg_viewer);CHKERRQ(ierr); 1887 } 1888 PetscFunctionReturn(0); 1889 } 1890 1891 /* 1892 PCApply_BDDC - Applies the BDDC operator to a vector. 1893 1894 Input Parameters: 1895 + pc - the preconditioner context 1896 - r - input vector (global) 1897 1898 Output Parameter: 1899 . z - output vector (global) 1900 1901 Application Interface Routine: PCApply() 1902 */ 1903 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 1904 { 1905 PC_IS *pcis = (PC_IS*)(pc->data); 1906 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1907 Mat lA = NULL; 1908 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1909 PetscErrorCode ierr; 1910 const PetscScalar one = 1.0; 1911 const PetscScalar m_one = -1.0; 1912 const PetscScalar zero = 0.0; 1913 /* This code is similar to that provided in nn.c for PCNN 1914 NN interface preconditioner changed to BDDC 1915 Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */ 1916 1917 PetscFunctionBegin; 1918 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 1919 if (pcbddc->switch_static) { 1920 ierr = MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat,&lA);CHKERRQ(ierr); 1921 } 1922 1923 if (pcbddc->ChangeOfBasisMatrix) { 1924 Vec swap; 1925 1926 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 1927 swap = pcbddc->work_change; 1928 pcbddc->work_change = r; 1929 r = swap; 1930 /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */ 1931 if (pcbddc->benign_apply_coarse_only && pcbddc->use_exact_dirichlet_trick && pcbddc->change_interior) { 1932 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 1933 ierr = VecLockReadPush(pcis->vec1_global);CHKERRQ(ierr); 1934 } 1935 } 1936 if (pcbddc->benign_have_null) { /* get p0 from r */ 1937 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 1938 } 1939 if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_DIRICHLET && !pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 1940 ierr = VecCopy(r,z);CHKERRQ(ierr); 1941 /* First Dirichlet solve */ 1942 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1943 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1944 /* 1945 Assembling right hand side for BDDC operator 1946 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1947 - pcis->vec1_B the interface part of the global vector z 1948 */ 1949 if (n_D) { 1950 ierr = PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 1951 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1952 ierr = PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 1953 ierr = KSPCheckSolve(pcbddc->ksp_D,pc,pcis->vec2_D);CHKERRQ(ierr); 1954 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1955 if (pcbddc->switch_static) { 1956 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 1957 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1958 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1959 if (!pcbddc->switch_static_change) { 1960 ierr = MatMult(lA,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1961 } else { 1962 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1963 ierr = MatMult(lA,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1964 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1965 } 1966 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1967 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1968 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1969 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1970 } else { 1971 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1972 } 1973 } else { 1974 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1975 } 1976 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1977 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1978 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1979 } else { 1980 if (!pcbddc->benign_apply_coarse_only) { 1981 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1982 } 1983 } 1984 if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_LUMP) { 1985 PetscCheckFalse(!pcbddc->switch_static,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"You forgot to pass -pc_bddc_switch_static"); 1986 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1987 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1988 } 1989 1990 /* Apply interface preconditioner 1991 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1992 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr); 1993 1994 /* Apply transpose of partition of unity operator */ 1995 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1996 if (pcbddc->interface_extension == PC_BDDC_INTERFACE_EXT_LUMP) { 1997 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1998 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1999 PetscFunctionReturn(0); 2000 } 2001 /* Second Dirichlet solve and assembling of output */ 2002 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2003 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2004 if (n_B) { 2005 if (pcbddc->switch_static) { 2006 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2007 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2008 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2009 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2010 if (!pcbddc->switch_static_change) { 2011 ierr = MatMult(lA,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2012 } else { 2013 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2014 ierr = MatMult(lA,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 2015 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2016 } 2017 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2018 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2019 } else { 2020 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 2021 } 2022 } else if (pcbddc->switch_static) { /* n_B is zero */ 2023 if (!pcbddc->switch_static_change) { 2024 ierr = MatMult(lA,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 2025 } else { 2026 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 2027 ierr = MatMult(lA,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2028 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 2029 } 2030 } 2031 ierr = PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 2032 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 2033 ierr = PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 2034 ierr = KSPCheckSolve(pcbddc->ksp_D,pc,pcis->vec4_D);CHKERRQ(ierr); 2035 2036 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 2037 if (pcbddc->switch_static) { 2038 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 2039 } else { 2040 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 2041 } 2042 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2043 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2044 } else { 2045 if (pcbddc->switch_static) { 2046 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 2047 } else { 2048 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 2049 } 2050 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2051 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2052 } 2053 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 2054 if (pcbddc->benign_apply_coarse_only) { 2055 ierr = PetscArrayzero(pcbddc->benign_p0,pcbddc->benign_n);CHKERRQ(ierr); 2056 } 2057 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 2058 } 2059 2060 if (pcbddc->ChangeOfBasisMatrix) { 2061 pcbddc->work_change = r; 2062 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 2063 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 2064 } 2065 PetscFunctionReturn(0); 2066 } 2067 2068 /* 2069 PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector. 2070 2071 Input Parameters: 2072 + pc - the preconditioner context 2073 - r - input vector (global) 2074 2075 Output Parameter: 2076 . z - output vector (global) 2077 2078 Application Interface Routine: PCApplyTranspose() 2079 */ 2080 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z) 2081 { 2082 PC_IS *pcis = (PC_IS*)(pc->data); 2083 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 2084 Mat lA = NULL; 2085 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 2086 PetscErrorCode ierr; 2087 const PetscScalar one = 1.0; 2088 const PetscScalar m_one = -1.0; 2089 const PetscScalar zero = 0.0; 2090 2091 PetscFunctionBegin; 2092 ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr); 2093 if (pcbddc->switch_static) { 2094 ierr = MatISGetLocalMat(pc->useAmat ? pc->mat : pc->pmat,&lA);CHKERRQ(ierr); 2095 } 2096 if (pcbddc->ChangeOfBasisMatrix) { 2097 Vec swap; 2098 2099 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 2100 swap = pcbddc->work_change; 2101 pcbddc->work_change = r; 2102 r = swap; 2103 /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */ 2104 if (pcbddc->benign_apply_coarse_only && pcbddc->exact_dirichlet_trick_app && pcbddc->change_interior) { 2105 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 2106 ierr = VecLockReadPush(pcis->vec1_global);CHKERRQ(ierr); 2107 } 2108 } 2109 if (pcbddc->benign_have_null) { /* get p0 from r */ 2110 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 2111 } 2112 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 2113 ierr = VecCopy(r,z);CHKERRQ(ierr); 2114 /* First Dirichlet solve */ 2115 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2116 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2117 /* 2118 Assembling right hand side for BDDC operator 2119 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 2120 - pcis->vec1_B the interface part of the global vector z 2121 */ 2122 if (n_D) { 2123 ierr = PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 2124 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 2125 ierr = PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 2126 ierr = KSPCheckSolve(pcbddc->ksp_D,pc,pcis->vec2_D);CHKERRQ(ierr); 2127 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 2128 if (pcbddc->switch_static) { 2129 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 2130 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2131 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2132 if (!pcbddc->switch_static_change) { 2133 ierr = MatMultTranspose(lA,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2134 } else { 2135 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2136 ierr = MatMultTranspose(lA,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 2137 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2138 } 2139 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2140 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2141 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2142 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2143 } else { 2144 ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 2145 } 2146 } else { 2147 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 2148 } 2149 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2150 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2151 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 2152 } else { 2153 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 2154 } 2155 2156 /* Apply interface preconditioner 2157 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 2158 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr); 2159 2160 /* Apply transpose of partition of unity operator */ 2161 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 2162 2163 /* Second Dirichlet solve and assembling of output */ 2164 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2165 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2166 if (n_B) { 2167 if (pcbddc->switch_static) { 2168 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2169 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2170 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2171 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2172 if (!pcbddc->switch_static_change) { 2173 ierr = MatMultTranspose(lA,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2174 } else { 2175 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2176 ierr = MatMultTranspose(lA,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 2177 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2178 } 2179 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2180 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2181 } else { 2182 ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 2183 } 2184 } else if (pcbddc->switch_static) { /* n_B is zero */ 2185 if (!pcbddc->switch_static_change) { 2186 ierr = MatMultTranspose(lA,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 2187 } else { 2188 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 2189 ierr = MatMultTranspose(lA,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2190 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 2191 } 2192 } 2193 ierr = PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 2194 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 2195 ierr = PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],pc,0,0,0);CHKERRQ(ierr); 2196 ierr = KSPCheckSolve(pcbddc->ksp_D,pc,pcis->vec4_D);CHKERRQ(ierr); 2197 if (!pcbddc->exact_dirichlet_trick_app && !pcbddc->benign_apply_coarse_only) { 2198 if (pcbddc->switch_static) { 2199 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 2200 } else { 2201 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 2202 } 2203 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2204 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2205 } else { 2206 if (pcbddc->switch_static) { 2207 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 2208 } else { 2209 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 2210 } 2211 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2212 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2213 } 2214 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 2215 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 2216 } 2217 if (pcbddc->ChangeOfBasisMatrix) { 2218 pcbddc->work_change = r; 2219 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 2220 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 2221 } 2222 PetscFunctionReturn(0); 2223 } 2224 2225 PetscErrorCode PCReset_BDDC(PC pc) 2226 { 2227 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 2228 PC_IS *pcis = (PC_IS*)pc->data; 2229 KSP kspD,kspR,kspC; 2230 PetscErrorCode ierr; 2231 2232 PetscFunctionBegin; 2233 /* free BDDC custom data */ 2234 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 2235 /* destroy objects related to topography */ 2236 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 2237 /* destroy objects for scaling operator */ 2238 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 2239 /* free solvers stuff */ 2240 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 2241 /* free global vectors needed in presolve */ 2242 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 2243 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 2244 /* free data created by PCIS */ 2245 ierr = PCISDestroy(pc);CHKERRQ(ierr); 2246 2247 /* restore defaults */ 2248 kspD = pcbddc->ksp_D; 2249 kspR = pcbddc->ksp_R; 2250 kspC = pcbddc->coarse_ksp; 2251 ierr = PetscMemzero(pc->data,sizeof(*pcbddc));CHKERRQ(ierr); 2252 pcis->n_neigh = -1; 2253 pcis->scaling_factor = 1.0; 2254 pcis->reusesubmatrices = PETSC_TRUE; 2255 pcbddc->use_local_adj = PETSC_TRUE; 2256 pcbddc->use_vertices = PETSC_TRUE; 2257 pcbddc->use_edges = PETSC_TRUE; 2258 pcbddc->symmetric_primal = PETSC_TRUE; 2259 pcbddc->vertex_size = 1; 2260 pcbddc->recompute_topography = PETSC_TRUE; 2261 pcbddc->coarse_size = -1; 2262 pcbddc->use_exact_dirichlet_trick = PETSC_TRUE; 2263 pcbddc->coarsening_ratio = 8; 2264 pcbddc->coarse_eqs_per_proc = 1; 2265 pcbddc->benign_compute_correction = PETSC_TRUE; 2266 pcbddc->nedfield = -1; 2267 pcbddc->nedglobal = PETSC_TRUE; 2268 pcbddc->graphmaxcount = PETSC_MAX_INT; 2269 pcbddc->sub_schurs_layers = -1; 2270 pcbddc->ksp_D = kspD; 2271 pcbddc->ksp_R = kspR; 2272 pcbddc->coarse_ksp = kspC; 2273 PetscFunctionReturn(0); 2274 } 2275 2276 PetscErrorCode PCDestroy_BDDC(PC pc) 2277 { 2278 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 2279 PetscErrorCode ierr; 2280 2281 PetscFunctionBegin; 2282 ierr = PCReset_BDDC(pc);CHKERRQ(ierr); 2283 ierr = KSPDestroy(&pcbddc->ksp_D);CHKERRQ(ierr); 2284 ierr = KSPDestroy(&pcbddc->ksp_R);CHKERRQ(ierr); 2285 ierr = KSPDestroy(&pcbddc->coarse_ksp);CHKERRQ(ierr); 2286 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDiscreteGradient_C",NULL);CHKERRQ(ierr); 2287 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",NULL);CHKERRQ(ierr); 2288 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);CHKERRQ(ierr); 2289 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 2290 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",NULL);CHKERRQ(ierr); 2291 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 2292 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr); 2293 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr); 2294 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr); 2295 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2296 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2297 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2298 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2299 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2300 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2301 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2302 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2303 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 2304 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr); 2305 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 2306 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 2307 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 2308 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 2309 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",NULL);CHKERRQ(ierr); 2310 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL);CHKERRQ(ierr); 2311 ierr = PetscFree(pc->data);CHKERRQ(ierr); 2312 PetscFunctionReturn(0); 2313 } 2314 2315 static PetscErrorCode PCSetCoordinates_BDDC(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords) 2316 { 2317 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 2318 PCBDDCGraph mat_graph = pcbddc->mat_graph; 2319 PetscErrorCode ierr; 2320 2321 PetscFunctionBegin; 2322 ierr = PetscFree(mat_graph->coords);CHKERRQ(ierr); 2323 ierr = PetscMalloc1(nloc*dim,&mat_graph->coords);CHKERRQ(ierr); 2324 ierr = PetscArraycpy(mat_graph->coords,coords,nloc*dim);CHKERRQ(ierr); 2325 mat_graph->cnloc = nloc; 2326 mat_graph->cdim = dim; 2327 mat_graph->cloc = PETSC_FALSE; 2328 /* flg setup */ 2329 pcbddc->recompute_topography = PETSC_TRUE; 2330 pcbddc->corner_selected = PETSC_FALSE; 2331 PetscFunctionReturn(0); 2332 } 2333 2334 static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool* change) 2335 { 2336 PetscFunctionBegin; 2337 *change = PETSC_TRUE; 2338 PetscFunctionReturn(0); 2339 } 2340 2341 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2342 { 2343 FETIDPMat_ctx mat_ctx; 2344 Vec work; 2345 PC_IS* pcis; 2346 PC_BDDC* pcbddc; 2347 PetscErrorCode ierr; 2348 2349 PetscFunctionBegin; 2350 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2351 pcis = (PC_IS*)mat_ctx->pc->data; 2352 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2353 2354 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 2355 /* copy rhs since we may change it during PCPreSolve_BDDC */ 2356 if (!pcbddc->original_rhs) { 2357 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 2358 } 2359 if (mat_ctx->rhs_flip) { 2360 ierr = VecPointwiseMult(pcbddc->original_rhs,standard_rhs,mat_ctx->rhs_flip);CHKERRQ(ierr); 2361 } else { 2362 ierr = VecCopy(standard_rhs,pcbddc->original_rhs);CHKERRQ(ierr); 2363 } 2364 if (mat_ctx->g2g_p) { 2365 /* interface pressure rhs */ 2366 ierr = VecScatterBegin(mat_ctx->g2g_p,fetidp_flux_rhs,pcbddc->original_rhs,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2367 ierr = VecScatterEnd(mat_ctx->g2g_p,fetidp_flux_rhs,pcbddc->original_rhs,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2368 ierr = VecScatterBegin(mat_ctx->g2g_p,standard_rhs,fetidp_flux_rhs,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2369 ierr = VecScatterEnd(mat_ctx->g2g_p,standard_rhs,fetidp_flux_rhs,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2370 if (!mat_ctx->rhs_flip) { 2371 ierr = VecScale(fetidp_flux_rhs,-1.);CHKERRQ(ierr); 2372 } 2373 } 2374 /* 2375 change of basis for physical rhs if needed 2376 It also changes the rhs in case of dirichlet boundaries 2377 */ 2378 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,pcbddc->original_rhs,NULL);CHKERRQ(ierr); 2379 if (pcbddc->ChangeOfBasisMatrix) { 2380 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,pcbddc->original_rhs,pcbddc->work_change);CHKERRQ(ierr); 2381 work = pcbddc->work_change; 2382 } else { 2383 work = pcbddc->original_rhs; 2384 } 2385 /* store vectors for computation of fetidp final solution */ 2386 ierr = VecScatterBegin(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2387 ierr = VecScatterEnd(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2388 /* scale rhs since it should be unassembled */ 2389 /* TODO use counter scaling? (also below) */ 2390 ierr = VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2391 ierr = VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2392 /* Apply partition of unity */ 2393 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2394 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2395 if (!pcbddc->switch_static) { 2396 /* compute partially subassembled Schur complement right-hand side */ 2397 ierr = PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],mat_ctx->pc,0,0,0);CHKERRQ(ierr); 2398 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2399 ierr = PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],mat_ctx->pc,0,0,0);CHKERRQ(ierr); 2400 /* Cannot propagate up error in KSPSolve() because there is no access to the PC */ 2401 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 2402 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 2403 ierr = VecSet(work,0.0);CHKERRQ(ierr); 2404 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2405 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2406 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2407 ierr = VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2408 ierr = VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2409 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2410 } 2411 /* BDDC rhs */ 2412 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 2413 if (pcbddc->switch_static) { 2414 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2415 } 2416 /* apply BDDC */ 2417 ierr = PetscArrayzero(pcbddc->benign_p0,pcbddc->benign_n);CHKERRQ(ierr); 2418 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2419 ierr = PetscArrayzero(pcbddc->benign_p0,pcbddc->benign_n);CHKERRQ(ierr); 2420 2421 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 2422 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 2423 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2424 ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2425 /* Add contribution to interface pressures */ 2426 if (mat_ctx->l2g_p) { 2427 ierr = MatMult(mat_ctx->B_BB,pcis->vec1_B,mat_ctx->vP);CHKERRQ(ierr); 2428 if (pcbddc->switch_static) { 2429 ierr = MatMultAdd(mat_ctx->B_BI,pcis->vec1_D,mat_ctx->vP,mat_ctx->vP);CHKERRQ(ierr); 2430 } 2431 ierr = VecScatterBegin(mat_ctx->l2g_p,mat_ctx->vP,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2432 ierr = VecScatterEnd(mat_ctx->l2g_p,mat_ctx->vP,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2433 } 2434 PetscFunctionReturn(0); 2435 } 2436 2437 /*@ 2438 PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETI-DP linear system using the physical right-hand side 2439 2440 Collective 2441 2442 Input Parameters: 2443 + fetidp_mat - the FETI-DP matrix object obtained by a call to PCBDDCCreateFETIDPOperators 2444 - standard_rhs - the right-hand side of the original linear system 2445 2446 Output Parameters: 2447 . fetidp_flux_rhs - the right-hand side for the FETI-DP linear system 2448 2449 Level: developer 2450 2451 Notes: 2452 2453 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetSolution 2454 @*/ 2455 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2456 { 2457 FETIDPMat_ctx mat_ctx; 2458 PetscErrorCode ierr; 2459 2460 PetscFunctionBegin; 2461 PetscValidHeaderSpecific(fetidp_mat,MAT_CLASSID,1); 2462 PetscValidHeaderSpecific(standard_rhs,VEC_CLASSID,2); 2463 PetscValidHeaderSpecific(fetidp_flux_rhs,VEC_CLASSID,3); 2464 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2465 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 2466 PetscFunctionReturn(0); 2467 } 2468 2469 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2470 { 2471 FETIDPMat_ctx mat_ctx; 2472 PC_IS* pcis; 2473 PC_BDDC* pcbddc; 2474 PetscErrorCode ierr; 2475 Vec work; 2476 2477 PetscFunctionBegin; 2478 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2479 pcis = (PC_IS*)mat_ctx->pc->data; 2480 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2481 2482 /* apply B_delta^T */ 2483 ierr = VecSet(pcis->vec1_B,0.);CHKERRQ(ierr); 2484 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2485 ierr = VecScatterEnd(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2486 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 2487 if (mat_ctx->l2g_p) { 2488 ierr = VecScatterBegin(mat_ctx->l2g_p,fetidp_flux_sol,mat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2489 ierr = VecScatterEnd(mat_ctx->l2g_p,fetidp_flux_sol,mat_ctx->vP,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2490 ierr = MatMultAdd(mat_ctx->Bt_BB,mat_ctx->vP,pcis->vec1_B,pcis->vec1_B);CHKERRQ(ierr); 2491 } 2492 2493 /* compute rhs for BDDC application */ 2494 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2495 if (pcbddc->switch_static) { 2496 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2497 if (mat_ctx->l2g_p) { 2498 ierr = VecScale(mat_ctx->vP,-1.);CHKERRQ(ierr); 2499 ierr = MatMultAdd(mat_ctx->Bt_BI,mat_ctx->vP,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); 2500 } 2501 } 2502 2503 /* apply BDDC */ 2504 ierr = PetscArrayzero(pcbddc->benign_p0,pcbddc->benign_n);CHKERRQ(ierr); 2505 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2506 2507 /* put values into global vector */ 2508 if (pcbddc->ChangeOfBasisMatrix) work = pcbddc->work_change; 2509 else work = standard_sol; 2510 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2511 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2512 if (!pcbddc->switch_static) { 2513 /* compute values into the interior if solved for the partially subassembled Schur complement */ 2514 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 2515 ierr = VecAYPX(pcis->vec1_D,-1.0,mat_ctx->temp_solution_D);CHKERRQ(ierr); 2516 ierr = PetscLogEventBegin(PC_BDDC_Solves[pcbddc->current_level][0],mat_ctx->pc,0,0,0);CHKERRQ(ierr); 2517 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); 2518 ierr = PetscLogEventEnd(PC_BDDC_Solves[pcbddc->current_level][0],mat_ctx->pc,0,0,0);CHKERRQ(ierr); 2519 /* Cannot propagate up error in KSPSolve() because there is no access to the PC */ 2520 } 2521 2522 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2523 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_D,work,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2524 /* add p0 solution to final solution */ 2525 ierr = PCBDDCBenignGetOrSetP0(mat_ctx->pc,work,PETSC_FALSE);CHKERRQ(ierr); 2526 if (pcbddc->ChangeOfBasisMatrix) { 2527 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,work,standard_sol);CHKERRQ(ierr); 2528 } 2529 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 2530 if (mat_ctx->g2g_p) { 2531 ierr = VecScatterBegin(mat_ctx->g2g_p,fetidp_flux_sol,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2532 ierr = VecScatterEnd(mat_ctx->g2g_p,fetidp_flux_sol,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2533 } 2534 PetscFunctionReturn(0); 2535 } 2536 2537 static PetscErrorCode PCView_BDDCIPC(PC pc, PetscViewer viewer) 2538 { 2539 PetscErrorCode ierr; 2540 BDDCIPC_ctx bddcipc_ctx; 2541 PetscBool isascii; 2542 2543 PetscFunctionBegin; 2544 ierr = PCShellGetContext(pc,&bddcipc_ctx);CHKERRQ(ierr); 2545 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 2546 if (isascii) { 2547 ierr = PetscViewerASCIIPrintf(viewer,"BDDC interface preconditioner\n");CHKERRQ(ierr); 2548 } 2549 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 2550 ierr = PCView(bddcipc_ctx->bddc,viewer);CHKERRQ(ierr); 2551 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 2552 PetscFunctionReturn(0); 2553 } 2554 2555 static PetscErrorCode PCSetUp_BDDCIPC(PC pc) 2556 { 2557 PetscErrorCode ierr; 2558 BDDCIPC_ctx bddcipc_ctx; 2559 PetscBool isbddc; 2560 Vec vv; 2561 IS is; 2562 PC_IS *pcis; 2563 2564 PetscFunctionBegin; 2565 ierr = PCShellGetContext(pc,&bddcipc_ctx);CHKERRQ(ierr); 2566 ierr = PetscObjectTypeCompare((PetscObject)bddcipc_ctx->bddc,PCBDDC,&isbddc);CHKERRQ(ierr); 2567 PetscCheckFalse(!isbddc,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid type %s. Must be of type bddc",((PetscObject)bddcipc_ctx->bddc)->type_name); 2568 ierr = PCSetUp(bddcipc_ctx->bddc);CHKERRQ(ierr); 2569 2570 /* create interface scatter */ 2571 pcis = (PC_IS*)(bddcipc_ctx->bddc->data); 2572 ierr = VecScatterDestroy(&bddcipc_ctx->g2l);CHKERRQ(ierr); 2573 ierr = MatCreateVecs(pc->pmat,&vv,NULL);CHKERRQ(ierr); 2574 ierr = ISRenumber(pcis->is_B_global,NULL,NULL,&is);CHKERRQ(ierr); 2575 ierr = VecScatterCreate(vv,is,pcis->vec1_B,NULL,&bddcipc_ctx->g2l);CHKERRQ(ierr); 2576 ierr = ISDestroy(&is);CHKERRQ(ierr); 2577 ierr = VecDestroy(&vv);CHKERRQ(ierr); 2578 PetscFunctionReturn(0); 2579 } 2580 2581 static PetscErrorCode PCApply_BDDCIPC(PC pc, Vec r, Vec x) 2582 { 2583 PetscErrorCode ierr; 2584 BDDCIPC_ctx bddcipc_ctx; 2585 PC_IS *pcis; 2586 VecScatter tmps; 2587 2588 PetscFunctionBegin; 2589 ierr = PCShellGetContext(pc,&bddcipc_ctx);CHKERRQ(ierr); 2590 pcis = (PC_IS*)(bddcipc_ctx->bddc->data); 2591 tmps = pcis->global_to_B; 2592 pcis->global_to_B = bddcipc_ctx->g2l; 2593 ierr = PCBDDCScalingRestriction(bddcipc_ctx->bddc,r,pcis->vec1_B);CHKERRQ(ierr); 2594 ierr = PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc,PETSC_FALSE);CHKERRQ(ierr); 2595 ierr = PCBDDCScalingExtension(bddcipc_ctx->bddc,pcis->vec1_B,x);CHKERRQ(ierr); 2596 pcis->global_to_B = tmps; 2597 PetscFunctionReturn(0); 2598 } 2599 2600 static PetscErrorCode PCApplyTranspose_BDDCIPC(PC pc, Vec r, Vec x) 2601 { 2602 PetscErrorCode ierr; 2603 BDDCIPC_ctx bddcipc_ctx; 2604 PC_IS *pcis; 2605 VecScatter tmps; 2606 2607 PetscFunctionBegin; 2608 ierr = PCShellGetContext(pc,&bddcipc_ctx);CHKERRQ(ierr); 2609 pcis = (PC_IS*)(bddcipc_ctx->bddc->data); 2610 tmps = pcis->global_to_B; 2611 pcis->global_to_B = bddcipc_ctx->g2l; 2612 ierr = PCBDDCScalingRestriction(bddcipc_ctx->bddc,r,pcis->vec1_B);CHKERRQ(ierr); 2613 ierr = PCBDDCApplyInterfacePreconditioner(bddcipc_ctx->bddc,PETSC_TRUE);CHKERRQ(ierr); 2614 ierr = PCBDDCScalingExtension(bddcipc_ctx->bddc,pcis->vec1_B,x);CHKERRQ(ierr); 2615 pcis->global_to_B = tmps; 2616 PetscFunctionReturn(0); 2617 } 2618 2619 static PetscErrorCode PCDestroy_BDDCIPC(PC pc) 2620 { 2621 PetscErrorCode ierr; 2622 BDDCIPC_ctx bddcipc_ctx; 2623 2624 PetscFunctionBegin; 2625 ierr = PCShellGetContext(pc,&bddcipc_ctx);CHKERRQ(ierr); 2626 ierr = PCDestroy(&bddcipc_ctx->bddc);CHKERRQ(ierr); 2627 ierr = VecScatterDestroy(&bddcipc_ctx->g2l);CHKERRQ(ierr); 2628 ierr = PetscFree(bddcipc_ctx);CHKERRQ(ierr); 2629 PetscFunctionReturn(0); 2630 } 2631 2632 /*@ 2633 PCBDDCMatFETIDPGetSolution - Compute the physical solution using the solution of the FETI-DP linear system 2634 2635 Collective 2636 2637 Input Parameters: 2638 + fetidp_mat - the FETI-DP matrix obtained by a call to PCBDDCCreateFETIDPOperators 2639 - fetidp_flux_sol - the solution of the FETI-DP linear system 2640 2641 Output Parameters: 2642 . standard_sol - the solution defined on the physical domain 2643 2644 Level: developer 2645 2646 Notes: 2647 2648 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetRHS 2649 @*/ 2650 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2651 { 2652 FETIDPMat_ctx mat_ctx; 2653 PetscErrorCode ierr; 2654 2655 PetscFunctionBegin; 2656 PetscValidHeaderSpecific(fetidp_mat,MAT_CLASSID,1); 2657 PetscValidHeaderSpecific(fetidp_flux_sol,VEC_CLASSID,2); 2658 PetscValidHeaderSpecific(standard_sol,VEC_CLASSID,3); 2659 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2660 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 2661 PetscFunctionReturn(0); 2662 } 2663 2664 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, PetscBool fully_redundant, const char* prefix, Mat *fetidp_mat, PC *fetidp_pc) 2665 { 2666 2667 FETIDPMat_ctx fetidpmat_ctx; 2668 Mat newmat; 2669 FETIDPPC_ctx fetidppc_ctx; 2670 PC newpc; 2671 MPI_Comm comm; 2672 PetscErrorCode ierr; 2673 2674 PetscFunctionBegin; 2675 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 2676 /* FETI-DP matrix */ 2677 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 2678 fetidpmat_ctx->fully_redundant = fully_redundant; 2679 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 2680 ierr = MatCreateShell(comm,fetidpmat_ctx->n,fetidpmat_ctx->n,fetidpmat_ctx->N,fetidpmat_ctx->N,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 2681 ierr = PetscObjectSetName((PetscObject)newmat,!fetidpmat_ctx->l2g_lambda_only ? "F" : "G");CHKERRQ(ierr); 2682 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 2683 ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr); 2684 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 2685 /* propagate MatOptions */ 2686 { 2687 PC_BDDC *pcbddc = (PC_BDDC*)fetidpmat_ctx->pc->data; 2688 PetscBool issym; 2689 2690 ierr = MatGetOption(pc->mat,MAT_SYMMETRIC,&issym);CHKERRQ(ierr); 2691 if (issym || pcbddc->symmetric_primal) { 2692 ierr = MatSetOption(newmat,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 2693 } 2694 } 2695 ierr = MatSetOptionsPrefix(newmat,prefix);CHKERRQ(ierr); 2696 ierr = MatAppendOptionsPrefix(newmat,"fetidp_");CHKERRQ(ierr); 2697 ierr = MatSetUp(newmat);CHKERRQ(ierr); 2698 /* FETI-DP preconditioner */ 2699 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 2700 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 2701 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 2702 ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr); 2703 ierr = PCSetOptionsPrefix(newpc,prefix);CHKERRQ(ierr); 2704 ierr = PCAppendOptionsPrefix(newpc,"fetidp_");CHKERRQ(ierr); 2705 ierr = PCSetErrorIfFailure(newpc,pc->erroriffailure);CHKERRQ(ierr); 2706 if (!fetidpmat_ctx->l2g_lambda_only) { /* standard FETI-DP */ 2707 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 2708 ierr = PCShellSetName(newpc,"FETI-DP multipliers");CHKERRQ(ierr); 2709 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 2710 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 2711 ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 2712 ierr = PCShellSetView(newpc,FETIDPPCView);CHKERRQ(ierr); 2713 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 2714 } else { /* saddle-point FETI-DP */ 2715 Mat M; 2716 PetscInt psize; 2717 PetscBool fake = PETSC_FALSE, isfieldsplit; 2718 2719 ierr = ISViewFromOptions(fetidpmat_ctx->lagrange,NULL,"-lag_view");CHKERRQ(ierr); 2720 ierr = ISViewFromOptions(fetidpmat_ctx->pressure,NULL,"-press_view");CHKERRQ(ierr); 2721 ierr = PetscObjectQuery((PetscObject)pc,"__KSPFETIDP_PPmat",(PetscObject*)&M);CHKERRQ(ierr); 2722 ierr = PCSetType(newpc,PCFIELDSPLIT);CHKERRQ(ierr); 2723 ierr = PCFieldSplitSetIS(newpc,"lag",fetidpmat_ctx->lagrange);CHKERRQ(ierr); 2724 ierr = PCFieldSplitSetIS(newpc,"p",fetidpmat_ctx->pressure);CHKERRQ(ierr); 2725 ierr = PCFieldSplitSetType(newpc,PC_COMPOSITE_SCHUR);CHKERRQ(ierr); 2726 ierr = PCFieldSplitSetSchurFactType(newpc,PC_FIELDSPLIT_SCHUR_FACT_DIAG);CHKERRQ(ierr); 2727 ierr = ISGetSize(fetidpmat_ctx->pressure,&psize);CHKERRQ(ierr); 2728 if (psize != M->rmap->N) { 2729 Mat M2; 2730 PetscInt lpsize; 2731 2732 fake = PETSC_TRUE; 2733 ierr = ISGetLocalSize(fetidpmat_ctx->pressure,&lpsize);CHKERRQ(ierr); 2734 ierr = MatCreate(comm,&M2);CHKERRQ(ierr); 2735 ierr = MatSetType(M2,MATAIJ);CHKERRQ(ierr); 2736 ierr = MatSetSizes(M2,lpsize,lpsize,psize,psize);CHKERRQ(ierr); 2737 ierr = MatSetUp(M2);CHKERRQ(ierr); 2738 ierr = MatAssemblyBegin(M2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2739 ierr = MatAssemblyEnd(M2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2740 ierr = PCFieldSplitSetSchurPre(newpc,PC_FIELDSPLIT_SCHUR_PRE_USER,M2);CHKERRQ(ierr); 2741 ierr = MatDestroy(&M2);CHKERRQ(ierr); 2742 } else { 2743 ierr = PCFieldSplitSetSchurPre(newpc,PC_FIELDSPLIT_SCHUR_PRE_USER,M);CHKERRQ(ierr); 2744 } 2745 ierr = PCFieldSplitSetSchurScale(newpc,1.0);CHKERRQ(ierr); 2746 2747 /* we need to setfromoptions and setup here to access the blocks */ 2748 ierr = PCSetFromOptions(newpc);CHKERRQ(ierr); 2749 ierr = PCSetUp(newpc);CHKERRQ(ierr); 2750 2751 /* user may have changed the type (e.g. -fetidp_pc_type none) */ 2752 ierr = PetscObjectTypeCompare((PetscObject)newpc,PCFIELDSPLIT,&isfieldsplit);CHKERRQ(ierr); 2753 if (isfieldsplit) { 2754 KSP *ksps; 2755 PC ppc,lagpc; 2756 PetscInt nn; 2757 PetscBool ismatis,matisok = PETSC_FALSE,check = PETSC_FALSE; 2758 2759 /* set the solver for the (0,0) block */ 2760 ierr = PCFieldSplitSchurGetSubKSP(newpc,&nn,&ksps);CHKERRQ(ierr); 2761 if (!nn) { /* not of type PC_COMPOSITE_SCHUR */ 2762 ierr = PCFieldSplitGetSubKSP(newpc,&nn,&ksps);CHKERRQ(ierr); 2763 if (!fake) { /* pass pmat to the pressure solver */ 2764 Mat F; 2765 2766 ierr = KSPGetOperators(ksps[1],&F,NULL);CHKERRQ(ierr); 2767 ierr = KSPSetOperators(ksps[1],F,M);CHKERRQ(ierr); 2768 } 2769 } else { 2770 PetscBool issym; 2771 Mat S; 2772 2773 ierr = PCFieldSplitSchurGetS(newpc,&S);CHKERRQ(ierr); 2774 2775 ierr = MatGetOption(newmat,MAT_SYMMETRIC,&issym);CHKERRQ(ierr); 2776 if (issym) { 2777 ierr = MatSetOption(S,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 2778 } 2779 } 2780 ierr = KSPGetPC(ksps[0],&lagpc);CHKERRQ(ierr); 2781 ierr = PCSetType(lagpc,PCSHELL);CHKERRQ(ierr); 2782 ierr = PCShellSetName(lagpc,"FETI-DP multipliers");CHKERRQ(ierr); 2783 ierr = PCShellSetContext(lagpc,fetidppc_ctx);CHKERRQ(ierr); 2784 ierr = PCShellSetApply(lagpc,FETIDPPCApply);CHKERRQ(ierr); 2785 ierr = PCShellSetApplyTranspose(lagpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 2786 ierr = PCShellSetView(lagpc,FETIDPPCView);CHKERRQ(ierr); 2787 ierr = PCShellSetDestroy(lagpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 2788 2789 /* Olof's idea: interface Schur complement preconditioner for the mass matrix */ 2790 ierr = KSPGetPC(ksps[1],&ppc);CHKERRQ(ierr); 2791 if (fake) { 2792 BDDCIPC_ctx bddcipc_ctx; 2793 PetscContainer c; 2794 2795 matisok = PETSC_TRUE; 2796 2797 /* create inner BDDC solver */ 2798 ierr = PetscNew(&bddcipc_ctx);CHKERRQ(ierr); 2799 ierr = PCCreate(comm,&bddcipc_ctx->bddc);CHKERRQ(ierr); 2800 ierr = PCSetType(bddcipc_ctx->bddc,PCBDDC);CHKERRQ(ierr); 2801 ierr = PCSetOperators(bddcipc_ctx->bddc,M,M);CHKERRQ(ierr); 2802 ierr = PetscObjectQuery((PetscObject)pc,"__KSPFETIDP_pCSR",(PetscObject*)&c);CHKERRQ(ierr); 2803 ierr = PetscObjectTypeCompare((PetscObject)M,MATIS,&ismatis);CHKERRQ(ierr); 2804 if (c && ismatis) { 2805 Mat lM; 2806 PetscInt *csr,n; 2807 2808 ierr = MatISGetLocalMat(M,&lM);CHKERRQ(ierr); 2809 ierr = MatGetSize(lM,&n,NULL);CHKERRQ(ierr); 2810 ierr = PetscContainerGetPointer(c,(void**)&csr);CHKERRQ(ierr); 2811 ierr = PCBDDCSetLocalAdjacencyGraph(bddcipc_ctx->bddc,n,csr,csr + (n + 1),PETSC_COPY_VALUES);CHKERRQ(ierr); 2812 ierr = MatISRestoreLocalMat(M,&lM);CHKERRQ(ierr); 2813 } 2814 ierr = PCSetOptionsPrefix(bddcipc_ctx->bddc,((PetscObject)ksps[1])->prefix);CHKERRQ(ierr); 2815 ierr = PCSetErrorIfFailure(bddcipc_ctx->bddc,pc->erroriffailure);CHKERRQ(ierr); 2816 ierr = PCSetFromOptions(bddcipc_ctx->bddc);CHKERRQ(ierr); 2817 2818 /* wrap the interface application */ 2819 ierr = PCSetType(ppc,PCSHELL);CHKERRQ(ierr); 2820 ierr = PCShellSetName(ppc,"FETI-DP pressure");CHKERRQ(ierr); 2821 ierr = PCShellSetContext(ppc,bddcipc_ctx);CHKERRQ(ierr); 2822 ierr = PCShellSetSetUp(ppc,PCSetUp_BDDCIPC);CHKERRQ(ierr); 2823 ierr = PCShellSetApply(ppc,PCApply_BDDCIPC);CHKERRQ(ierr); 2824 ierr = PCShellSetApplyTranspose(ppc,PCApplyTranspose_BDDCIPC);CHKERRQ(ierr); 2825 ierr = PCShellSetView(ppc,PCView_BDDCIPC);CHKERRQ(ierr); 2826 ierr = PCShellSetDestroy(ppc,PCDestroy_BDDCIPC);CHKERRQ(ierr); 2827 } 2828 2829 /* determine if we need to assemble M to construct a preconditioner */ 2830 if (!matisok) { 2831 ierr = PetscObjectTypeCompare((PetscObject)M,MATIS,&ismatis);CHKERRQ(ierr); 2832 ierr = PetscObjectTypeCompareAny((PetscObject)ppc,&matisok,PCBDDC,PCJACOBI,PCNONE,PCMG,"");CHKERRQ(ierr); 2833 if (ismatis && !matisok) { 2834 ierr = MatConvert(M,MATAIJ,MAT_INPLACE_MATRIX,&M);CHKERRQ(ierr); 2835 } 2836 } 2837 2838 /* run the subproblems to check convergence */ 2839 ierr = PetscOptionsGetBool(NULL,((PetscObject)newmat)->prefix,"-check_saddlepoint",&check,NULL);CHKERRQ(ierr); 2840 if (check) { 2841 PetscInt i; 2842 2843 for (i=0;i<nn;i++) { 2844 KSP kspC; 2845 PC pc; 2846 Mat F,pF; 2847 Vec x,y; 2848 PetscBool isschur,prec = PETSC_TRUE; 2849 2850 ierr = KSPCreate(PetscObjectComm((PetscObject)ksps[i]),&kspC);CHKERRQ(ierr); 2851 ierr = KSPSetOptionsPrefix(kspC,((PetscObject)ksps[i])->prefix);CHKERRQ(ierr); 2852 ierr = KSPAppendOptionsPrefix(kspC,"check_");CHKERRQ(ierr); 2853 ierr = KSPGetOperators(ksps[i],&F,&pF);CHKERRQ(ierr); 2854 ierr = PetscObjectTypeCompare((PetscObject)F,MATSCHURCOMPLEMENT,&isschur);CHKERRQ(ierr); 2855 if (isschur) { 2856 KSP kspS,kspS2; 2857 Mat A00,pA00,A10,A01,A11; 2858 char prefix[256]; 2859 2860 ierr = MatSchurComplementGetKSP(F,&kspS);CHKERRQ(ierr); 2861 ierr = MatSchurComplementGetSubMatrices(F,&A00,&pA00,&A01,&A10,&A11);CHKERRQ(ierr); 2862 ierr = MatCreateSchurComplement(A00,pA00,A01,A10,A11,&F);CHKERRQ(ierr); 2863 ierr = MatSchurComplementGetKSP(F,&kspS2);CHKERRQ(ierr); 2864 ierr = PetscSNPrintf(prefix,sizeof(prefix),"%sschur_",((PetscObject)kspC)->prefix);CHKERRQ(ierr); 2865 ierr = KSPSetOptionsPrefix(kspS2,prefix);CHKERRQ(ierr); 2866 ierr = KSPGetPC(kspS2,&pc);CHKERRQ(ierr); 2867 ierr = PCSetType(pc,PCKSP);CHKERRQ(ierr); 2868 ierr = PCKSPSetKSP(pc,kspS);CHKERRQ(ierr); 2869 ierr = KSPSetFromOptions(kspS2);CHKERRQ(ierr); 2870 ierr = KSPGetPC(kspS2,&pc);CHKERRQ(ierr); 2871 ierr = PCSetUseAmat(pc,PETSC_TRUE);CHKERRQ(ierr); 2872 } else { 2873 ierr = PetscObjectReference((PetscObject)F);CHKERRQ(ierr); 2874 } 2875 ierr = KSPSetFromOptions(kspC);CHKERRQ(ierr); 2876 ierr = PetscOptionsGetBool(NULL,((PetscObject)kspC)->prefix,"-preconditioned",&prec,NULL);CHKERRQ(ierr); 2877 if (prec) { 2878 ierr = KSPGetPC(ksps[i],&pc);CHKERRQ(ierr); 2879 ierr = KSPSetPC(kspC,pc);CHKERRQ(ierr); 2880 } 2881 ierr = KSPSetOperators(kspC,F,pF);CHKERRQ(ierr); 2882 ierr = MatCreateVecs(F,&x,&y);CHKERRQ(ierr); 2883 ierr = VecSetRandom(x,NULL);CHKERRQ(ierr); 2884 ierr = MatMult(F,x,y);CHKERRQ(ierr); 2885 ierr = KSPSolve(kspC,y,x);CHKERRQ(ierr); 2886 ierr = KSPCheckSolve(kspC,pc,x);CHKERRQ(ierr); 2887 ierr = KSPDestroy(&kspC);CHKERRQ(ierr); 2888 ierr = MatDestroy(&F);CHKERRQ(ierr); 2889 ierr = VecDestroy(&x);CHKERRQ(ierr); 2890 ierr = VecDestroy(&y);CHKERRQ(ierr); 2891 } 2892 } 2893 ierr = PetscFree(ksps);CHKERRQ(ierr); 2894 } 2895 } 2896 /* return pointers for objects created */ 2897 *fetidp_mat = newmat; 2898 *fetidp_pc = newpc; 2899 PetscFunctionReturn(0); 2900 } 2901 2902 /*@C 2903 PCBDDCCreateFETIDPOperators - Create FETI-DP operators 2904 2905 Collective 2906 2907 Input Parameters: 2908 + pc - the BDDC preconditioning context (setup should have been called before) 2909 . fully_redundant - true for a fully redundant set of Lagrange multipliers 2910 - prefix - optional options database prefix for the objects to be created (can be NULL) 2911 2912 Output Parameters: 2913 + fetidp_mat - shell FETI-DP matrix object 2914 - fetidp_pc - shell Dirichlet preconditioner for FETI-DP matrix 2915 2916 Level: developer 2917 2918 Notes: 2919 Currently the only operations provided for FETI-DP matrix are MatMult and MatMultTranspose 2920 2921 .seealso: PCBDDC, PCBDDCMatFETIDPGetRHS, PCBDDCMatFETIDPGetSolution 2922 @*/ 2923 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, PetscBool fully_redundant, const char *prefix, Mat *fetidp_mat, PC *fetidp_pc) 2924 { 2925 PetscErrorCode ierr; 2926 2927 PetscFunctionBegin; 2928 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2929 if (pc->setupcalled) { 2930 ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,PetscBool,const char*,Mat*,PC*),(pc,fully_redundant,prefix,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 2931 } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"You must call PCSetup_BDDC() first"); 2932 PetscFunctionReturn(0); 2933 } 2934 /* -------------------------------------------------------------------------- */ 2935 /*MC 2936 PCBDDC - Balancing Domain Decomposition by Constraints. 2937 2938 An implementation of the BDDC preconditioner based on the bibliography found below. 2939 2940 The matrix to be preconditioned (Pmat) must be of type MATIS. 2941 2942 Currently works with MATIS matrices with local matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers. 2943 2944 It also works with unsymmetric and indefinite problems. 2945 2946 Unlike 'conventional' interface preconditioners, PCBDDC iterates over all degrees of freedom, not just those on the interface. This allows the use of approximate solvers on the subdomains. 2947 2948 Approximate local solvers are automatically adapted (see [1]) if the user has attached a nullspace object to the subdomain matrices, and informed BDDC of using approximate solvers (via the command line). 2949 2950 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() 2951 Additional information on dofs can be provided by using PCBDDCSetDofsSplitting(), PCBDDCSetDirichletBoundaries(), PCBDDCSetNeumannBoundaries(), and PCBDDCSetPrimalVerticesIS() and their local counterparts. 2952 2953 Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). Non-singular modes are retained via SVD. 2954 2955 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. 2956 User defined change of basis can be passed to PCBDDC by using PCBDDCSetChangeOfBasisMat() 2957 2958 The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using a MatPartitioning object. 2959 2960 Adaptive selection of primal constraints [4] is supported for SPD systems with high-contrast in the coefficients if MUMPS or MKL_PARDISO are present. Future versions of the code will also consider using PASTIX. 2961 2962 An experimental interface to the FETI-DP method is available. FETI-DP operators could be created using PCBDDCCreateFETIDPOperators(). A stand-alone class for the FETI-DP method will be provided in the next releases. 2963 2964 Options Database Keys (some of them, run with -help for a complete list): 2965 2966 + -pc_bddc_use_vertices <true> - use or not vertices in primal space 2967 . -pc_bddc_use_edges <true> - use or not edges in primal space 2968 . -pc_bddc_use_faces <false> - use or not faces in primal space 2969 . -pc_bddc_symmetric <true> - symmetric computation of primal basis functions. Specify false for unsymmetric problems 2970 . -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only) 2971 . -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested 2972 . -pc_bddc_switch_static <false> - switches from M_2 (default) to M_3 operator (see reference article [1]) 2973 . -pc_bddc_levels <0> - maximum number of levels for multilevel 2974 . -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) 2975 . -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) 2976 . -pc_bddc_use_deluxe_scaling <false> - use deluxe scaling 2977 . -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) 2978 . -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) 2979 - -pc_bddc_check_level <0> - set verbosity level of debugging output 2980 2981 Options for Dirichlet, Neumann or coarse solver can be set with 2982 .vb 2983 -pc_bddc_dirichlet_ 2984 -pc_bddc_neumann_ 2985 -pc_bddc_coarse_ 2986 .ve 2987 e.g. -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. PCBDDC uses by default KSPPREONLY and PCLU. 2988 2989 When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified as 2990 .vb 2991 -pc_bddc_dirichlet_lN_ 2992 -pc_bddc_neumann_lN_ 2993 -pc_bddc_coarse_lN_ 2994 .ve 2995 Note that level number ranges from the finest (0) to the coarsest (N). 2996 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. 2997 .vb 2998 -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3 2999 .ve 3000 will use a threshold of 5 for constraints' selection at the first coarse level and will redistribute the coarse problem of the first coarse level on 3 processors 3001 3002 References: 3003 + [1] - C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149--168, March 2007 3004 . [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 3005 . [3] - J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", Computing Volume 83, Issue 2--3, pages 55--85, November 2008 3006 - [4] - C. Pechstein and C. R. Dohrmann. "Modern domain decomposition methods BDDC, deluxe scaling, and an algebraic approach", Seminar talk, Linz, December 2013, http://people.ricam.oeaw.ac.at/c.pechstein/pechstein-bddc2013.pdf 3007 3008 Level: intermediate 3009 3010 Developer Notes: 3011 3012 Contributed by Stefano Zampini 3013 3014 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 3015 M*/ 3016 3017 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 3018 { 3019 PetscErrorCode ierr; 3020 PC_BDDC *pcbddc; 3021 3022 PetscFunctionBegin; 3023 ierr = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr); 3024 pc->data = pcbddc; 3025 3026 /* create PCIS data structure */ 3027 ierr = PCISCreate(pc);CHKERRQ(ierr); 3028 3029 /* create local graph structure */ 3030 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 3031 3032 /* BDDC nonzero defaults */ 3033 pcbddc->use_nnsp = PETSC_TRUE; 3034 pcbddc->use_local_adj = PETSC_TRUE; 3035 pcbddc->use_vertices = PETSC_TRUE; 3036 pcbddc->use_edges = PETSC_TRUE; 3037 pcbddc->symmetric_primal = PETSC_TRUE; 3038 pcbddc->vertex_size = 1; 3039 pcbddc->recompute_topography = PETSC_TRUE; 3040 pcbddc->coarse_size = -1; 3041 pcbddc->use_exact_dirichlet_trick = PETSC_TRUE; 3042 pcbddc->coarsening_ratio = 8; 3043 pcbddc->coarse_eqs_per_proc = 1; 3044 pcbddc->benign_compute_correction = PETSC_TRUE; 3045 pcbddc->nedfield = -1; 3046 pcbddc->nedglobal = PETSC_TRUE; 3047 pcbddc->graphmaxcount = PETSC_MAX_INT; 3048 pcbddc->sub_schurs_layers = -1; 3049 pcbddc->adaptive_threshold[0] = 0.0; 3050 pcbddc->adaptive_threshold[1] = 0.0; 3051 3052 /* function pointers */ 3053 pc->ops->apply = PCApply_BDDC; 3054 pc->ops->applytranspose = PCApplyTranspose_BDDC; 3055 pc->ops->setup = PCSetUp_BDDC; 3056 pc->ops->destroy = PCDestroy_BDDC; 3057 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 3058 pc->ops->view = PCView_BDDC; 3059 pc->ops->applyrichardson = NULL; 3060 pc->ops->applysymmetricleft = NULL; 3061 pc->ops->applysymmetricright = NULL; 3062 pc->ops->presolve = PCPreSolve_BDDC; 3063 pc->ops->postsolve = PCPostSolve_BDDC; 3064 pc->ops->reset = PCReset_BDDC; 3065 3066 /* composing function */ 3067 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDiscreteGradient_C",PCBDDCSetDiscreteGradient_BDDC);CHKERRQ(ierr); 3068 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",PCBDDCSetDivergenceMat_BDDC);CHKERRQ(ierr); 3069 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);CHKERRQ(ierr); 3070 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 3071 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",PCBDDCSetPrimalVerticesIS_BDDC);CHKERRQ(ierr); 3072 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetPrimalVerticesLocalIS_C",PCBDDCGetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 3073 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetPrimalVerticesIS_C",PCBDDCGetPrimalVerticesIS_BDDC);CHKERRQ(ierr); 3074 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 3075 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr); 3076 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr); 3077 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr); 3078 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 3079 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 3080 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 3081 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 3082 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 3083 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 3084 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 3085 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 3086 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 3087 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr); 3088 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 3089 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 3090 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 3091 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 3092 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_BDDC);CHKERRQ(ierr); 3093 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_BDDC);CHKERRQ(ierr); 3094 PetscFunctionReturn(0); 3095 } 3096 3097 /*@C 3098 PCBDDCInitializePackage - This function initializes everything in the PCBDDC package. It is called 3099 from PCInitializePackage(). 3100 3101 Level: developer 3102 3103 .seealso: PetscInitialize() 3104 @*/ 3105 PetscErrorCode PCBDDCInitializePackage(void) 3106 { 3107 PetscErrorCode ierr; 3108 int i; 3109 3110 PetscFunctionBegin; 3111 if (PCBDDCPackageInitialized) PetscFunctionReturn(0); 3112 PCBDDCPackageInitialized = PETSC_TRUE; 3113 ierr = PetscRegisterFinalize(PCBDDCFinalizePackage);CHKERRQ(ierr); 3114 3115 /* general events */ 3116 ierr = PetscLogEventRegister("PCBDDCTopo",PC_CLASSID,&PC_BDDC_Topology[0]);CHKERRQ(ierr); 3117 ierr = PetscLogEventRegister("PCBDDCLKSP",PC_CLASSID,&PC_BDDC_LocalSolvers[0]);CHKERRQ(ierr); 3118 ierr = PetscLogEventRegister("PCBDDCLWor",PC_CLASSID,&PC_BDDC_LocalWork[0]);CHKERRQ(ierr); 3119 ierr = PetscLogEventRegister("PCBDDCCorr",PC_CLASSID,&PC_BDDC_CorrectionSetUp[0]);CHKERRQ(ierr); 3120 ierr = PetscLogEventRegister("PCBDDCASet",PC_CLASSID,&PC_BDDC_ApproxSetUp[0]);CHKERRQ(ierr); 3121 ierr = PetscLogEventRegister("PCBDDCAApp",PC_CLASSID,&PC_BDDC_ApproxApply[0]);CHKERRQ(ierr); 3122 ierr = PetscLogEventRegister("PCBDDCCSet",PC_CLASSID,&PC_BDDC_CoarseSetUp[0]);CHKERRQ(ierr); 3123 ierr = PetscLogEventRegister("PCBDDCCKSP",PC_CLASSID,&PC_BDDC_CoarseSolver[0]);CHKERRQ(ierr); 3124 ierr = PetscLogEventRegister("PCBDDCAdap",PC_CLASSID,&PC_BDDC_AdaptiveSetUp[0]);CHKERRQ(ierr); 3125 ierr = PetscLogEventRegister("PCBDDCScal",PC_CLASSID,&PC_BDDC_Scaling[0]);CHKERRQ(ierr); 3126 ierr = PetscLogEventRegister("PCBDDCSchr",PC_CLASSID,&PC_BDDC_Schurs[0]);CHKERRQ(ierr); 3127 ierr = PetscLogEventRegister("PCBDDCDirS",PC_CLASSID,&PC_BDDC_Solves[0][0]);CHKERRQ(ierr); 3128 ierr = PetscLogEventRegister("PCBDDCNeuS",PC_CLASSID,&PC_BDDC_Solves[0][1]);CHKERRQ(ierr); 3129 ierr = PetscLogEventRegister("PCBDDCCoaS",PC_CLASSID,&PC_BDDC_Solves[0][2]);CHKERRQ(ierr); 3130 for (i=1;i<PETSC_PCBDDC_MAXLEVELS;i++) { 3131 char ename[32]; 3132 3133 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCTopo l%02d",i);CHKERRQ(ierr); 3134 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Topology[i]);CHKERRQ(ierr); 3135 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCLKSP l%02d",i);CHKERRQ(ierr); 3136 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_LocalSolvers[i]);CHKERRQ(ierr); 3137 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCLWor l%02d",i);CHKERRQ(ierr); 3138 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_LocalWork[i]);CHKERRQ(ierr); 3139 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCCorr l%02d",i);CHKERRQ(ierr); 3140 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_CorrectionSetUp[i]);CHKERRQ(ierr); 3141 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCASet l%02d",i);CHKERRQ(ierr); 3142 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_ApproxSetUp[i]);CHKERRQ(ierr); 3143 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCAApp l%02d",i);CHKERRQ(ierr); 3144 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_ApproxApply[i]);CHKERRQ(ierr); 3145 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCCSet l%02d",i);CHKERRQ(ierr); 3146 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_CoarseSetUp[i]);CHKERRQ(ierr); 3147 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCCKSP l%02d",i);CHKERRQ(ierr); 3148 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_CoarseSolver[i]);CHKERRQ(ierr); 3149 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCAdap l%02d",i);CHKERRQ(ierr); 3150 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_AdaptiveSetUp[i]);CHKERRQ(ierr); 3151 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCScal l%02d",i);CHKERRQ(ierr); 3152 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Scaling[i]);CHKERRQ(ierr); 3153 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCSchr l%02d",i);CHKERRQ(ierr); 3154 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Schurs[i]);CHKERRQ(ierr); 3155 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCDirS l%02d",i);CHKERRQ(ierr); 3156 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Solves[i][0]);CHKERRQ(ierr); 3157 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCNeuS l%02d",i);CHKERRQ(ierr); 3158 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Solves[i][1]);CHKERRQ(ierr); 3159 ierr = PetscSNPrintf(ename,sizeof(ename),"PCBDDCCoaS l%02d",i);CHKERRQ(ierr); 3160 ierr = PetscLogEventRegister(ename,PC_CLASSID,&PC_BDDC_Solves[i][2]);CHKERRQ(ierr); 3161 } 3162 PetscFunctionReturn(0); 3163 } 3164 3165 /*@C 3166 PCBDDCFinalizePackage - This function frees everything from the PCBDDC package. It is 3167 called from PetscFinalize() automatically. 3168 3169 Level: developer 3170 3171 .seealso: PetscFinalize() 3172 @*/ 3173 PetscErrorCode PCBDDCFinalizePackage(void) 3174 { 3175 PetscFunctionBegin; 3176 PCBDDCPackageInitialized = PETSC_FALSE; 3177 PetscFunctionReturn(0); 3178 } 3179