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 FETIDP 18 - Move FETIDP code to its own classes 19 20 MATIS related operations contained in BDDC code 21 - Provide general case for subassembling 22 23 */ 24 25 #include <../src/ksp/pc/impls/bddc/bddc.h> /*I "petscpc.h" I*/ /* includes for fortran wrappers */ 26 #include <../src/ksp/pc/impls/bddc/bddcprivate.h> 27 #include <petscblaslapack.h> 28 29 /* temporarily declare it */ 30 PetscErrorCode PCApply_BDDC(PC,Vec,Vec); 31 32 /* -------------------------------------------------------------------------- */ 33 #undef __FUNCT__ 34 #define __FUNCT__ "PCSetFromOptions_BDDC" 35 PetscErrorCode PCSetFromOptions_BDDC(PetscOptionItems *PetscOptionsObject,PC pc) 36 { 37 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 ierr = PetscOptionsHead(PetscOptionsObject,"BDDC options");CHKERRQ(ierr); 42 /* Verbose debugging */ 43 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); 44 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); 45 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); 46 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); 47 ierr = PetscOptionsInt("-pc_bddc_check_level","Verbose output for PCBDDC (intended for debug)","none",pcbddc->dbg_flag,&pcbddc->dbg_flag,NULL);CHKERRQ(ierr); 48 /* Primal space customization */ 49 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); 50 ierr = PetscOptionsBool("-pc_bddc_use_vertices","Use or not corner dofs in coarse space","none",pcbddc->use_vertices,&pcbddc->use_vertices,NULL);CHKERRQ(ierr); 51 ierr = PetscOptionsBool("-pc_bddc_use_edges","Use or not edge constraints in coarse space","none",pcbddc->use_edges,&pcbddc->use_edges,NULL);CHKERRQ(ierr); 52 ierr = PetscOptionsBool("-pc_bddc_use_faces","Use or not face constraints in coarse space","none",pcbddc->use_faces,&pcbddc->use_faces,NULL);CHKERRQ(ierr); 53 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); 54 ierr = PetscOptionsBool("-pc_bddc_use_true_nnsp","Use near null space attached to the matrix without modifications","none",pcbddc->use_nnsp_true,&pcbddc->use_nnsp_true,NULL);CHKERRQ(ierr); 55 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); 56 /* Change of basis */ 57 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); 58 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); 59 if (!pcbddc->use_change_of_basis) { 60 pcbddc->use_change_on_faces = PETSC_FALSE; 61 } 62 /* Switch between M_2 (default) and M_3 preconditioners (as defined by C. Dohrmann in the ref. article) */ 63 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); 64 ierr = PetscOptionsInt("-pc_bddc_coarse_eqs_per_proc","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); 65 ierr = PetscOptionsInt("-pc_bddc_coarsening_ratio","Set coarsening ratio used in multilevel coarsening","none",pcbddc->coarsening_ratio,&pcbddc->coarsening_ratio,NULL);CHKERRQ(ierr); 66 ierr = PetscOptionsInt("-pc_bddc_levels","Set maximum number of levels for multilevel","none",pcbddc->max_levels,&pcbddc->max_levels,NULL);CHKERRQ(ierr); 67 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); 68 ierr = PetscOptionsBool("-pc_bddc_use_deluxe_scaling","Use deluxe scaling for BDDC","none",pcbddc->use_deluxe_scaling,&pcbddc->use_deluxe_scaling,NULL);CHKERRQ(ierr); 69 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); 70 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); 71 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); 72 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); 73 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); 74 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); 75 ierr = PetscOptionsReal("-pc_bddc_adaptive_threshold","Threshold to be used for adaptive selection of constraints","none",pcbddc->adaptive_threshold,&pcbddc->adaptive_threshold,NULL);CHKERRQ(ierr); 76 ierr = PetscOptionsInt("-pc_bddc_adaptive_nmin","Minimum number of constraints per connected components","none",pcbddc->adaptive_nmin,&pcbddc->adaptive_nmin,NULL);CHKERRQ(ierr); 77 ierr = PetscOptionsInt("-pc_bddc_adaptive_nmax","Maximum number of constraints per connected components","none",pcbddc->adaptive_nmax,&pcbddc->adaptive_nmax,NULL);CHKERRQ(ierr); 78 ierr = PetscOptionsBool("-pc_bddc_symmetric","Symmetric computation of primal basis functions","none",pcbddc->symmetric_primal,&pcbddc->symmetric_primal,NULL);CHKERRQ(ierr); 79 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); 80 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); 81 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); 82 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); 83 ierr = PetscOptionsBool("-pc_bddc_nonetflux","Automatic computation of no-net-flux quadrature weights","none",pcbddc->compute_nonetflux,&pcbddc->compute_nonetflux,NULL);CHKERRQ(ierr); 84 ierr = PetscOptionsBool("-pc_bddc_detect_disconnected","Detects disconnected subdomains","none",pcbddc->detect_disconnected,&pcbddc->detect_disconnected,NULL);CHKERRQ(ierr); 85 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); 86 ierr = PetscOptionsTail();CHKERRQ(ierr); 87 PetscFunctionReturn(0); 88 } 89 90 /* -------------------------------------------------------------------------- */ 91 #undef __FUNCT__ 92 #define __FUNCT__ "PCView_BDDC" 93 static PetscErrorCode PCView_BDDC(PC pc,PetscViewer viewer) 94 { 95 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 96 PC_IS *pcis = (PC_IS*)pc->data; 97 PetscErrorCode ierr; 98 PetscBool isascii,isstring; 99 PetscSubcomm subcomm; 100 PetscViewer subviewer; 101 102 PetscFunctionBegin; 103 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 104 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);CHKERRQ(ierr); 105 /* Nothing printed for the String viewer */ 106 /* ASCII viewer */ 107 if (isascii) { 108 PetscMPIInt color,rank,size; 109 Petsc64bitInt loc[6],gsum[5],gmax[5],gmin[5],totbenign; 110 PetscScalar interface_size; 111 PetscReal ratio1=0.,ratio2=0.; 112 Vec counter; 113 114 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use verbose output: %d\n",pcbddc->dbg_flag);CHKERRQ(ierr); 115 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use user-defined CSR: %d\n",!!pcbddc->mat_graph->nvtxs_csr);CHKERRQ(ierr); 116 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use local mat graph: %d\n",pcbddc->use_local_adj);CHKERRQ(ierr); 117 if (pcbddc->mat_graph->twodim) { 118 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Connectivity graph topological dimension: 2\n");CHKERRQ(ierr); 119 } else { 120 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Connectivity graph topological dimension: 3\n");CHKERRQ(ierr); 121 } 122 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use vertices: %d (vertex size %d)\n",pcbddc->use_vertices,pcbddc->vertex_size);CHKERRQ(ierr); 123 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use edges: %d\n",pcbddc->use_edges);CHKERRQ(ierr); 124 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use faces: %d\n",pcbddc->use_faces);CHKERRQ(ierr); 125 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use true near null space: %d\n",pcbddc->use_nnsp_true);CHKERRQ(ierr); 126 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use QR for single constraints on cc: %d\n",pcbddc->use_qr_single);CHKERRQ(ierr); 127 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use change of basis on local edge nodes: %d\n",pcbddc->use_change_of_basis);CHKERRQ(ierr); 128 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use change of basis on local face nodes: %d\n",pcbddc->use_change_on_faces);CHKERRQ(ierr); 129 ierr = PetscViewerASCIIPrintf(viewer," BDDC: User defined change of basis matrix: %d\n",!!pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 130 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Has change of basis matrix: %d\n",!!pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 131 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Eliminate dirichlet boundary dofs: %d\n",pcbddc->eliminate_dirdofs);CHKERRQ(ierr); 132 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Switch on static condensation ops around the interface preconditioner: %d\n",pcbddc->switch_static);CHKERRQ(ierr); 133 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use exact dirichlet trick: %d\n",pcbddc->use_exact_dirichlet_trick);CHKERRQ(ierr); 134 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Multilevel max levels: %d\n",pcbddc->max_levels);CHKERRQ(ierr); 135 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Multilevel coarsening ratio: %d\n",pcbddc->coarsening_ratio);CHKERRQ(ierr); 136 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use estimated eigs for coarse problem: %d\n",pcbddc->use_coarse_estimates);CHKERRQ(ierr); 137 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use deluxe scaling: %d\n",pcbddc->use_deluxe_scaling);CHKERRQ(ierr); 138 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use deluxe zerorows: %d\n",pcbddc->deluxe_zerorows);CHKERRQ(ierr); 139 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Rebuild interface graph for Schur principal minors: %d\n",pcbddc->sub_schurs_rebuild);CHKERRQ(ierr); 140 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Number of dofs' layers for the computation of principal minors: %d\n",pcbddc->sub_schurs_layers);CHKERRQ(ierr); 141 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Use user CSR graph to compute successive layers: %d\n",pcbddc->sub_schurs_use_useradj);CHKERRQ(ierr); 142 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Adaptive constraint selection threshold (active %d, userdefined %d): %g\n",pcbddc->adaptive_threshold,pcbddc->adaptive_selection,pcbddc->adaptive_userdefined);CHKERRQ(ierr); 143 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Min constraints / connected component: %d\n",pcbddc->adaptive_nmin);CHKERRQ(ierr); 144 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Max constraints / connected component: %d\n",pcbddc->adaptive_nmax);CHKERRQ(ierr); 145 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Invert exact Schur complement for adaptive selection: %d\n",pcbddc->sub_schurs_exact_schur);CHKERRQ(ierr); 146 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Symmetric computation of primal basis functions: %d\n",pcbddc->symmetric_primal);CHKERRQ(ierr); 147 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Num. Procs. to map coarse adjacency list: %d\n",pcbddc->coarse_adj_red);CHKERRQ(ierr); 148 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Coarse eqs per proc (significant at the coarsest level): %d\n",pcbddc->coarse_eqs_per_proc);CHKERRQ(ierr); 149 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Detect disconnected: %d\n",pcbddc->detect_disconnected);CHKERRQ(ierr); 150 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Benign subspace trick: %d (change explicit %d)\n",pcbddc->benign_saddle_point,pcbddc->benign_change_explicit);CHKERRQ(ierr); 151 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Benign subspace trick is active: %d\n",pcbddc->benign_have_null);CHKERRQ(ierr); 152 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Algebraic computation of no-net-flux %d\n",pcbddc->compute_nonetflux);CHKERRQ(ierr); 153 154 /* compute some numbers on the domain decomposition */ 155 ierr = VecSet(pcis->vec1_B,1.0);CHKERRQ(ierr); 156 ierr = MatCreateVecs(pc->pmat,&counter,0);CHKERRQ(ierr); 157 ierr = VecSet(counter,0.0);CHKERRQ(ierr); 158 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,counter,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 159 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,counter,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 160 ierr = VecSum(counter,&interface_size);CHKERRQ(ierr); 161 ierr = VecDestroy(&counter);CHKERRQ(ierr); 162 gsum[0] = 1; 163 gsum[1] = gsum[2] = gsum[3] = gsum[4] = 0; 164 loc[0] = !!pcis->n; 165 loc[1] = pcis->n - pcis->n_B; 166 loc[2] = pcis->n_B; 167 loc[3] = pcbddc->local_primal_size; 168 loc[4] = pcis->n; 169 loc[5] = pcbddc->benign_n; 170 ierr = MPI_Reduce(loc,gsum,5,MPIU_INT64,MPI_SUM,0,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 171 if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = -1; 172 ierr = MPI_Reduce(loc,gmax,5,MPIU_INT64,MPI_MAX,0,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 173 if (!loc[0]) loc[1] = loc[2] = loc[3] = loc[4] = PETSC_MAX_INT; 174 ierr = MPI_Reduce(loc,gmin,5,MPIU_INT64,MPI_MIN,0,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 175 ierr = MPI_Reduce(&loc[5],&totbenign,1,MPIU_INT64,MPI_SUM,0,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 176 if (pcbddc->coarse_size) { 177 ratio1 = pc->pmat->rmap->N/(1.*pcbddc->coarse_size); 178 ratio2 = PetscRealPart(interface_size)/pcbddc->coarse_size; 179 } 180 ierr = PetscViewerASCIIPrintf(viewer," BDDC: ********************************** STATISTICS AT LEVEL %d **********************************\n",pcbddc->current_level);CHKERRQ(ierr); 181 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Global dofs sizes: all %d interface %d coarse %d\n",pc->pmat->rmap->N,(PetscInt)PetscRealPart(interface_size),pcbddc->coarse_size);CHKERRQ(ierr); 182 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Coarsening ratios: all/coarse %d interface/coarse %d\n",(PetscInt)ratio1,(PetscInt)ratio2);CHKERRQ(ierr); 183 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Active processes : %d\n",(PetscInt)gsum[0]);CHKERRQ(ierr); 184 if (pcbddc->benign_have_null) { 185 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Benign subs : %d\n",(PetscInt)totbenign);CHKERRQ(ierr); 186 } 187 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Dofs type :\tMIN\tMAX\tMEAN\n",(PetscInt)gmin[1],(PetscInt)gmax[1],(PetscInt)(gsum[1]/gsum[0]));CHKERRQ(ierr); 188 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Interior dofs :\t%d\t%d\t%d\n",(PetscInt)gmin[1],(PetscInt)gmax[1],(PetscInt)(gsum[1]/gsum[0]));CHKERRQ(ierr); 189 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Interface dofs :\t%d\t%d\t%d\n",(PetscInt)gmin[2],(PetscInt)gmax[2],(PetscInt)(gsum[2]/gsum[0]));CHKERRQ(ierr); 190 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Primal dofs :\t%d\t%d\t%d\n",(PetscInt)gmin[3],(PetscInt)gmax[3],(PetscInt)(gsum[3]/gsum[0]));CHKERRQ(ierr); 191 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Local dofs :\t%d\t%d\t%d\n",(PetscInt)gmin[4],(PetscInt)gmax[4],(PetscInt)(gsum[4]/gsum[0]));CHKERRQ(ierr); 192 ierr = PetscViewerASCIIPrintf(viewer," BDDC: ********************************** COARSE PROBLEM DETAILS *********************************\n",pcbddc->current_level);CHKERRQ(ierr); 193 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 194 195 /* the coarse problem can be handled with a different communicator */ 196 if (pcbddc->coarse_ksp) color = 1; 197 else color = 0; 198 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);CHKERRQ(ierr); 199 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr); 200 ierr = PetscSubcommCreate(PetscObjectComm((PetscObject)pc),&subcomm);CHKERRQ(ierr); 201 ierr = PetscSubcommSetNumber(subcomm,PetscMin(size,2));CHKERRQ(ierr); 202 ierr = PetscSubcommSetTypeGeneral(subcomm,color,rank);CHKERRQ(ierr); 203 ierr = PetscViewerGetSubViewer(viewer,PetscSubcommChild(subcomm),&subviewer);CHKERRQ(ierr); 204 if (color == 1) { 205 ierr = KSPView(pcbddc->coarse_ksp,subviewer);CHKERRQ(ierr); 206 ierr = PetscViewerFlush(subviewer);CHKERRQ(ierr); 207 } 208 ierr = PetscViewerRestoreSubViewer(viewer,PetscSubcommChild(subcomm),&subviewer);CHKERRQ(ierr); 209 ierr = PetscSubcommDestroy(&subcomm);CHKERRQ(ierr); 210 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 211 } 212 PetscFunctionReturn(0); 213 } 214 /* -------------------------------------------------------------------------- */ 215 #undef __FUNCT__ 216 #define __FUNCT__ "PCBDDCSetDivergenceMat_BDDC" 217 static PetscErrorCode PCBDDCSetDivergenceMat_BDDC(PC pc, Mat divudotp, ISLocalToGlobalMapping vl2l) 218 { 219 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 220 PetscErrorCode ierr; 221 222 PetscFunctionBegin; 223 ierr = PetscObjectReference((PetscObject)divudotp);CHKERRQ(ierr); 224 ierr = MatDestroy(&pcbddc->divudotp);CHKERRQ(ierr); 225 pcbddc->divudotp = divudotp; 226 pcbddc->compute_nonetflux = PETSC_TRUE; 227 if (vl2l) { 228 ierr = PetscObjectReference((PetscObject)vl2l);CHKERRQ(ierr); 229 ierr = ISLocalToGlobalMappingDestroy(&pcbddc->divudotp_vl2l);CHKERRQ(ierr); 230 pcbddc->divudotp_vl2l = vl2l; 231 } 232 PetscFunctionReturn(0); 233 } 234 #undef __FUNCT__ 235 #define __FUNCT__ "PCBDDCSetDivergenceMat" 236 /*@ 237 PCBDDCSetDivergenceMat - Sets the linear operator representing \int_\Omega \div {\bf u} \cdot p dx 238 239 Collective on PC 240 241 Input Parameters: 242 + pc - the preconditioning context 243 . divudotp - the matrix (must be of type MATIS) 244 - vl2l - optional ISLocalToGlobalMapping describing the local (to divudotp) to local (to pc->pmat) map for the velocities 245 246 Level: advanced 247 248 Notes: It triggers the algebraic computation of no-net-flux condition 249 250 .seealso: PCBDDC 251 @*/ 252 PetscErrorCode PCBDDCSetDivergenceMat(PC pc, Mat divudotp, ISLocalToGlobalMapping vl2l) 253 { 254 PetscBool ismatis; 255 PetscErrorCode ierr; 256 257 PetscFunctionBegin; 258 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 259 PetscValidHeaderSpecific(divudotp,MAT_CLASSID,2); 260 PetscCheckSameComm(pc,1,divudotp,2); 261 if (vl2l) PetscValidHeaderSpecific(divudotp,IS_LTOGM_CLASSID,3); 262 ierr = PetscObjectTypeCompare((PetscObject)divudotp,MATIS,&ismatis);CHKERRQ(ierr); 263 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Divergence matrix needs to be of type MATIS"); 264 ierr = PetscTryMethod(pc,"PCBDDCSetDivergenceMat_C",(PC,Mat,ISLocalToGlobalMapping),(pc,divudotp,vl2l));CHKERRQ(ierr); 265 PetscFunctionReturn(0); 266 } 267 /* -------------------------------------------------------------------------- */ 268 #undef __FUNCT__ 269 #define __FUNCT__ "PCBDDCSetChangeOfBasisMat_BDDC" 270 static PetscErrorCode PCBDDCSetChangeOfBasisMat_BDDC(PC pc, Mat change, PetscBool interior) 271 { 272 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 273 PetscErrorCode ierr; 274 275 PetscFunctionBegin; 276 ierr = PetscObjectReference((PetscObject)change);CHKERRQ(ierr); 277 ierr = MatDestroy(&pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 278 pcbddc->user_ChangeOfBasisMatrix = change; 279 pcbddc->change_interior = interior; 280 PetscFunctionReturn(0); 281 } 282 #undef __FUNCT__ 283 #define __FUNCT__ "PCBDDCSetChangeOfBasisMat" 284 /*@ 285 PCBDDCSetChangeOfBasisMat - Set user defined change of basis for dofs 286 287 Collective on PC 288 289 Input Parameters: 290 + pc - the preconditioning context 291 . change - the change of basis matrix 292 - interior - whether or not the change of basis modifies interior dofs 293 294 Level: intermediate 295 296 Notes: 297 298 .seealso: PCBDDC 299 @*/ 300 PetscErrorCode PCBDDCSetChangeOfBasisMat(PC pc, Mat change, PetscBool interior) 301 { 302 PetscErrorCode ierr; 303 304 PetscFunctionBegin; 305 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 306 PetscValidHeaderSpecific(change,MAT_CLASSID,2); 307 PetscCheckSameComm(pc,1,change,2); 308 if (pc->mat) { 309 PetscInt rows_c,cols_c,rows,cols; 310 ierr = MatGetSize(pc->mat,&rows,&cols);CHKERRQ(ierr); 311 ierr = MatGetSize(change,&rows_c,&cols_c);CHKERRQ(ierr); 312 if (rows_c != rows) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of rows for change of basis matrix! %d != %d",rows_c,rows); 313 if (cols_c != cols) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of columns for change of basis matrix! %d != %d",cols_c,cols); 314 ierr = MatGetLocalSize(pc->mat,&rows,&cols);CHKERRQ(ierr); 315 ierr = MatGetLocalSize(change,&rows_c,&cols_c);CHKERRQ(ierr); 316 if (rows_c != rows) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of local rows for change of basis matrix! %d != %d",rows_c,rows); 317 if (cols_c != cols) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Invalid number of local columns for change of basis matrix! %d != %d",cols_c,cols); 318 } 319 ierr = PetscTryMethod(pc,"PCBDDCSetChangeOfBasisMat_C",(PC,Mat,PetscBool),(pc,change,interior));CHKERRQ(ierr); 320 PetscFunctionReturn(0); 321 } 322 /* -------------------------------------------------------------------------- */ 323 #undef __FUNCT__ 324 #define __FUNCT__ "PCBDDCSetPrimalVerticesIS_BDDC" 325 static PetscErrorCode PCBDDCSetPrimalVerticesIS_BDDC(PC pc, IS PrimalVertices) 326 { 327 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 328 PetscBool isequal = PETSC_FALSE; 329 PetscErrorCode ierr; 330 331 PetscFunctionBegin; 332 ierr = PetscObjectReference((PetscObject)PrimalVertices);CHKERRQ(ierr); 333 if (pcbddc->user_primal_vertices) { 334 ierr = ISEqual(PrimalVertices,pcbddc->user_primal_vertices,&isequal);CHKERRQ(ierr); 335 } 336 ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr); 337 ierr = ISDestroy(&pcbddc->user_primal_vertices_local);CHKERRQ(ierr); 338 pcbddc->user_primal_vertices = PrimalVertices; 339 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 340 PetscFunctionReturn(0); 341 } 342 #undef __FUNCT__ 343 #define __FUNCT__ "PCBDDCSetPrimalVerticesIS" 344 /*@ 345 PCBDDCSetPrimalVerticesIS - Set additional user defined primal vertices in PCBDDC 346 347 Collective 348 349 Input Parameters: 350 + pc - the preconditioning context 351 - PrimalVertices - index set of primal vertices in global numbering (can be empty) 352 353 Level: intermediate 354 355 Notes: 356 Any process can list any global node 357 358 .seealso: PCBDDC 359 @*/ 360 PetscErrorCode PCBDDCSetPrimalVerticesIS(PC pc, IS PrimalVertices) 361 { 362 PetscErrorCode ierr; 363 364 PetscFunctionBegin; 365 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 366 PetscValidHeaderSpecific(PrimalVertices,IS_CLASSID,2); 367 PetscCheckSameComm(pc,1,PrimalVertices,2); 368 ierr = PetscTryMethod(pc,"PCBDDCSetPrimalVerticesIS_C",(PC,IS),(pc,PrimalVertices));CHKERRQ(ierr); 369 PetscFunctionReturn(0); 370 } 371 /* -------------------------------------------------------------------------- */ 372 #undef __FUNCT__ 373 #define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS_BDDC" 374 static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices) 375 { 376 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 377 PetscBool isequal = PETSC_FALSE; 378 PetscErrorCode ierr; 379 380 PetscFunctionBegin; 381 ierr = PetscObjectReference((PetscObject)PrimalVertices);CHKERRQ(ierr); 382 if (pcbddc->user_primal_vertices_local) { 383 ierr = ISEqual(PrimalVertices,pcbddc->user_primal_vertices_local,&isequal);CHKERRQ(ierr); 384 } 385 ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr); 386 ierr = ISDestroy(&pcbddc->user_primal_vertices_local);CHKERRQ(ierr); 387 pcbddc->user_primal_vertices_local = PrimalVertices; 388 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 389 PetscFunctionReturn(0); 390 } 391 #undef __FUNCT__ 392 #define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS" 393 /*@ 394 PCBDDCSetPrimalVerticesLocalIS - Set additional user defined primal vertices in PCBDDC 395 396 Collective 397 398 Input Parameters: 399 + pc - the preconditioning context 400 - PrimalVertices - index set of primal vertices in local numbering (can be empty) 401 402 Level: intermediate 403 404 Notes: 405 406 .seealso: PCBDDC 407 @*/ 408 PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices) 409 { 410 PetscErrorCode ierr; 411 412 PetscFunctionBegin; 413 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 414 PetscValidHeaderSpecific(PrimalVertices,IS_CLASSID,2); 415 PetscCheckSameComm(pc,1,PrimalVertices,2); 416 ierr = PetscTryMethod(pc,"PCBDDCSetPrimalVerticesLocalIS_C",(PC,IS),(pc,PrimalVertices));CHKERRQ(ierr); 417 PetscFunctionReturn(0); 418 } 419 /* -------------------------------------------------------------------------- */ 420 #undef __FUNCT__ 421 #define __FUNCT__ "PCBDDCSetCoarseningRatio_BDDC" 422 static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc,PetscInt k) 423 { 424 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 425 426 PetscFunctionBegin; 427 pcbddc->coarsening_ratio = k; 428 PetscFunctionReturn(0); 429 } 430 431 #undef __FUNCT__ 432 #define __FUNCT__ "PCBDDCSetCoarseningRatio" 433 /*@ 434 PCBDDCSetCoarseningRatio - Set coarsening ratio used in multilevel 435 436 Logically collective on PC 437 438 Input Parameters: 439 + pc - the preconditioning context 440 - k - coarsening ratio (H/h at the coarser level) 441 442 Options Database Keys: 443 . -pc_bddc_coarsening_ratio 444 445 Level: intermediate 446 447 Notes: 448 Approximatively k subdomains at the finer level will be aggregated into a single subdomain at the coarser level 449 450 .seealso: PCBDDC, PCBDDCSetLevels() 451 @*/ 452 PetscErrorCode PCBDDCSetCoarseningRatio(PC pc,PetscInt k) 453 { 454 PetscErrorCode ierr; 455 456 PetscFunctionBegin; 457 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 458 PetscValidLogicalCollectiveInt(pc,k,2); 459 ierr = PetscTryMethod(pc,"PCBDDCSetCoarseningRatio_C",(PC,PetscInt),(pc,k));CHKERRQ(ierr); 460 PetscFunctionReturn(0); 461 } 462 463 /* The following functions (PCBDDCSetUseExactDirichlet PCBDDCSetLevel) are not public */ 464 #undef __FUNCT__ 465 #define __FUNCT__ "PCBDDCSetUseExactDirichlet_BDDC" 466 static PetscErrorCode PCBDDCSetUseExactDirichlet_BDDC(PC pc,PetscBool flg) 467 { 468 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 469 470 PetscFunctionBegin; 471 pcbddc->use_exact_dirichlet_trick = flg; 472 PetscFunctionReturn(0); 473 } 474 475 #undef __FUNCT__ 476 #define __FUNCT__ "PCBDDCSetUseExactDirichlet" 477 PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool flg) 478 { 479 PetscErrorCode ierr; 480 481 PetscFunctionBegin; 482 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 483 PetscValidLogicalCollectiveBool(pc,flg,2); 484 ierr = PetscTryMethod(pc,"PCBDDCSetUseExactDirichlet_C",(PC,PetscBool),(pc,flg));CHKERRQ(ierr); 485 PetscFunctionReturn(0); 486 } 487 488 #undef __FUNCT__ 489 #define __FUNCT__ "PCBDDCSetLevel_BDDC" 490 static PetscErrorCode PCBDDCSetLevel_BDDC(PC pc,PetscInt level) 491 { 492 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 493 494 PetscFunctionBegin; 495 pcbddc->current_level = level; 496 PetscFunctionReturn(0); 497 } 498 499 #undef __FUNCT__ 500 #define __FUNCT__ "PCBDDCSetLevel" 501 PetscErrorCode PCBDDCSetLevel(PC pc,PetscInt level) 502 { 503 PetscErrorCode ierr; 504 505 PetscFunctionBegin; 506 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 507 PetscValidLogicalCollectiveInt(pc,level,2); 508 ierr = PetscTryMethod(pc,"PCBDDCSetLevel_C",(PC,PetscInt),(pc,level));CHKERRQ(ierr); 509 PetscFunctionReturn(0); 510 } 511 512 #undef __FUNCT__ 513 #define __FUNCT__ "PCBDDCSetLevels_BDDC" 514 static PetscErrorCode PCBDDCSetLevels_BDDC(PC pc,PetscInt levels) 515 { 516 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 517 518 PetscFunctionBegin; 519 pcbddc->max_levels = levels; 520 PetscFunctionReturn(0); 521 } 522 523 #undef __FUNCT__ 524 #define __FUNCT__ "PCBDDCSetLevels" 525 /*@ 526 PCBDDCSetLevels - Sets the maximum number of levels for multilevel 527 528 Logically collective on PC 529 530 Input Parameters: 531 + pc - the preconditioning context 532 - levels - the maximum number of levels (max 9) 533 534 Options Database Keys: 535 . -pc_bddc_levels 536 537 Level: intermediate 538 539 Notes: 540 Default value is 0, i.e. traditional one-level BDDC 541 542 .seealso: PCBDDC, PCBDDCSetCoarseningRatio() 543 @*/ 544 PetscErrorCode PCBDDCSetLevels(PC pc,PetscInt levels) 545 { 546 PetscErrorCode ierr; 547 548 PetscFunctionBegin; 549 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 550 PetscValidLogicalCollectiveInt(pc,levels,2); 551 if (levels > 99) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Maximum number of levels for bddc is 99\n"); 552 ierr = PetscTryMethod(pc,"PCBDDCSetLevels_C",(PC,PetscInt),(pc,levels));CHKERRQ(ierr); 553 PetscFunctionReturn(0); 554 } 555 /* -------------------------------------------------------------------------- */ 556 557 #undef __FUNCT__ 558 #define __FUNCT__ "PCBDDCSetDirichletBoundaries_BDDC" 559 static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc,IS DirichletBoundaries) 560 { 561 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 562 PetscBool isequal = PETSC_FALSE; 563 PetscErrorCode ierr; 564 565 PetscFunctionBegin; 566 ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr); 567 if (pcbddc->DirichletBoundaries) { 568 ierr = ISEqual(DirichletBoundaries,pcbddc->DirichletBoundaries,&isequal);CHKERRQ(ierr); 569 } 570 /* last user setting takes precendence -> destroy any other customization */ 571 ierr = ISDestroy(&pcbddc->DirichletBoundariesLocal);CHKERRQ(ierr); 572 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 573 pcbddc->DirichletBoundaries = DirichletBoundaries; 574 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 575 PetscFunctionReturn(0); 576 } 577 578 #undef __FUNCT__ 579 #define __FUNCT__ "PCBDDCSetDirichletBoundaries" 580 /*@ 581 PCBDDCSetDirichletBoundaries - Set IS defining Dirichlet boundaries for the global problem. 582 583 Collective 584 585 Input Parameters: 586 + pc - the preconditioning context 587 - DirichletBoundaries - parallel IS defining the Dirichlet boundaries 588 589 Level: intermediate 590 591 Notes: 592 Provide the information if you used MatZeroRows/Columns routines. Any process can list any global node 593 594 .seealso: PCBDDC, PCBDDCSetDirichletBoundariesLocal(), MatZeroRows(), MatZeroRowsColumns() 595 @*/ 596 PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc,IS DirichletBoundaries) 597 { 598 PetscErrorCode ierr; 599 600 PetscFunctionBegin; 601 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 602 PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2); 603 PetscCheckSameComm(pc,1,DirichletBoundaries,2); 604 ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundaries_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr); 605 PetscFunctionReturn(0); 606 } 607 /* -------------------------------------------------------------------------- */ 608 609 #undef __FUNCT__ 610 #define __FUNCT__ "PCBDDCSetDirichletBoundariesLocal_BDDC" 611 static PetscErrorCode PCBDDCSetDirichletBoundariesLocal_BDDC(PC pc,IS DirichletBoundaries) 612 { 613 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 614 PetscBool isequal = PETSC_FALSE; 615 PetscErrorCode ierr; 616 617 PetscFunctionBegin; 618 ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr); 619 if (pcbddc->DirichletBoundariesLocal) { 620 ierr = ISEqual(DirichletBoundaries,pcbddc->DirichletBoundariesLocal,&isequal);CHKERRQ(ierr); 621 } 622 /* last user setting takes precendence -> destroy any other customization */ 623 ierr = ISDestroy(&pcbddc->DirichletBoundariesLocal);CHKERRQ(ierr); 624 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 625 pcbddc->DirichletBoundariesLocal = DirichletBoundaries; 626 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 627 PetscFunctionReturn(0); 628 } 629 630 #undef __FUNCT__ 631 #define __FUNCT__ "PCBDDCSetDirichletBoundariesLocal" 632 /*@ 633 PCBDDCSetDirichletBoundariesLocal - Set IS defining Dirichlet boundaries for the global problem in local ordering. 634 635 Collective 636 637 Input Parameters: 638 + pc - the preconditioning context 639 - DirichletBoundaries - parallel IS defining the Dirichlet boundaries (in local ordering) 640 641 Level: intermediate 642 643 Notes: 644 645 .seealso: PCBDDC, PCBDDCSetDirichletBoundaries(), MatZeroRows(), MatZeroRowsColumns() 646 @*/ 647 PetscErrorCode PCBDDCSetDirichletBoundariesLocal(PC pc,IS DirichletBoundaries) 648 { 649 PetscErrorCode ierr; 650 651 PetscFunctionBegin; 652 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 653 PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2); 654 PetscCheckSameComm(pc,1,DirichletBoundaries,2); 655 ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundariesLocal_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr); 656 PetscFunctionReturn(0); 657 } 658 /* -------------------------------------------------------------------------- */ 659 660 #undef __FUNCT__ 661 #define __FUNCT__ "PCBDDCSetNeumannBoundaries_BDDC" 662 static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc,IS NeumannBoundaries) 663 { 664 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 665 PetscBool isequal = PETSC_FALSE; 666 PetscErrorCode ierr; 667 668 PetscFunctionBegin; 669 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 670 if (pcbddc->NeumannBoundaries) { 671 ierr = ISEqual(NeumannBoundaries,pcbddc->NeumannBoundaries,&isequal);CHKERRQ(ierr); 672 } 673 /* last user setting takes precendence -> destroy any other customization */ 674 ierr = ISDestroy(&pcbddc->NeumannBoundariesLocal);CHKERRQ(ierr); 675 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 676 pcbddc->NeumannBoundaries = NeumannBoundaries; 677 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 678 PetscFunctionReturn(0); 679 } 680 681 #undef __FUNCT__ 682 #define __FUNCT__ "PCBDDCSetNeumannBoundaries" 683 /*@ 684 PCBDDCSetNeumannBoundaries - Set IS defining Neumann boundaries for the global problem. 685 686 Collective 687 688 Input Parameters: 689 + pc - the preconditioning context 690 - NeumannBoundaries - parallel IS defining the Neumann boundaries 691 692 Level: intermediate 693 694 Notes: 695 Any process can list any global node 696 697 .seealso: PCBDDC, PCBDDCSetNeumannBoundariesLocal() 698 @*/ 699 PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc,IS NeumannBoundaries) 700 { 701 PetscErrorCode ierr; 702 703 PetscFunctionBegin; 704 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 705 PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2); 706 PetscCheckSameComm(pc,1,NeumannBoundaries,2); 707 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundaries_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 708 PetscFunctionReturn(0); 709 } 710 /* -------------------------------------------------------------------------- */ 711 712 #undef __FUNCT__ 713 #define __FUNCT__ "PCBDDCSetNeumannBoundariesLocal_BDDC" 714 static PetscErrorCode PCBDDCSetNeumannBoundariesLocal_BDDC(PC pc,IS NeumannBoundaries) 715 { 716 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 717 PetscBool isequal = PETSC_FALSE; 718 PetscErrorCode ierr; 719 720 PetscFunctionBegin; 721 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 722 if (pcbddc->NeumannBoundariesLocal) { 723 ierr = ISEqual(NeumannBoundaries,pcbddc->NeumannBoundariesLocal,&isequal);CHKERRQ(ierr); 724 } 725 /* last user setting takes precendence -> destroy any other customization */ 726 ierr = ISDestroy(&pcbddc->NeumannBoundariesLocal);CHKERRQ(ierr); 727 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 728 pcbddc->NeumannBoundariesLocal = NeumannBoundaries; 729 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 730 PetscFunctionReturn(0); 731 } 732 733 #undef __FUNCT__ 734 #define __FUNCT__ "PCBDDCSetNeumannBoundariesLocal" 735 /*@ 736 PCBDDCSetNeumannBoundariesLocal - Set IS defining Neumann boundaries for the global problem in local ordering. 737 738 Collective 739 740 Input Parameters: 741 + pc - the preconditioning context 742 - NeumannBoundaries - parallel IS defining the subdomain part of Neumann boundaries (in local ordering) 743 744 Level: intermediate 745 746 Notes: 747 748 .seealso: PCBDDC, PCBDDCSetNeumannBoundaries() 749 @*/ 750 PetscErrorCode PCBDDCSetNeumannBoundariesLocal(PC pc,IS NeumannBoundaries) 751 { 752 PetscErrorCode ierr; 753 754 PetscFunctionBegin; 755 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 756 PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2); 757 PetscCheckSameComm(pc,1,NeumannBoundaries,2); 758 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundariesLocal_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 759 PetscFunctionReturn(0); 760 } 761 /* -------------------------------------------------------------------------- */ 762 763 #undef __FUNCT__ 764 #define __FUNCT__ "PCBDDCGetDirichletBoundaries_BDDC" 765 static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc,IS *DirichletBoundaries) 766 { 767 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 768 769 PetscFunctionBegin; 770 *DirichletBoundaries = pcbddc->DirichletBoundaries; 771 PetscFunctionReturn(0); 772 } 773 774 #undef __FUNCT__ 775 #define __FUNCT__ "PCBDDCGetDirichletBoundaries" 776 /*@ 777 PCBDDCGetDirichletBoundaries - Get parallel IS for Dirichlet boundaries 778 779 Collective 780 781 Input Parameters: 782 . pc - the preconditioning context 783 784 Output Parameters: 785 . DirichletBoundaries - index set defining the Dirichlet boundaries 786 787 Level: intermediate 788 789 Notes: 790 The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetDirichletBoundaries 791 792 .seealso: PCBDDC 793 @*/ 794 PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc,IS *DirichletBoundaries) 795 { 796 PetscErrorCode ierr; 797 798 PetscFunctionBegin; 799 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 800 ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundaries_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr); 801 PetscFunctionReturn(0); 802 } 803 /* -------------------------------------------------------------------------- */ 804 805 #undef __FUNCT__ 806 #define __FUNCT__ "PCBDDCGetDirichletBoundariesLocal_BDDC" 807 static PetscErrorCode PCBDDCGetDirichletBoundariesLocal_BDDC(PC pc,IS *DirichletBoundaries) 808 { 809 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 810 811 PetscFunctionBegin; 812 *DirichletBoundaries = pcbddc->DirichletBoundariesLocal; 813 PetscFunctionReturn(0); 814 } 815 816 #undef __FUNCT__ 817 #define __FUNCT__ "PCBDDCGetDirichletBoundariesLocal" 818 /*@ 819 PCBDDCGetDirichletBoundariesLocal - Get parallel IS for Dirichlet boundaries (in local ordering) 820 821 Collective 822 823 Input Parameters: 824 . pc - the preconditioning context 825 826 Output Parameters: 827 . DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries 828 829 Level: intermediate 830 831 Notes: 832 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). 833 In the latter case, the IS will be available after PCSetUp. 834 835 .seealso: PCBDDC 836 @*/ 837 PetscErrorCode PCBDDCGetDirichletBoundariesLocal(PC pc,IS *DirichletBoundaries) 838 { 839 PetscErrorCode ierr; 840 841 PetscFunctionBegin; 842 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 843 ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundariesLocal_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr); 844 PetscFunctionReturn(0); 845 } 846 /* -------------------------------------------------------------------------- */ 847 848 #undef __FUNCT__ 849 #define __FUNCT__ "PCBDDCGetNeumannBoundaries_BDDC" 850 static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc,IS *NeumannBoundaries) 851 { 852 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 853 854 PetscFunctionBegin; 855 *NeumannBoundaries = pcbddc->NeumannBoundaries; 856 PetscFunctionReturn(0); 857 } 858 859 #undef __FUNCT__ 860 #define __FUNCT__ "PCBDDCGetNeumannBoundaries" 861 /*@ 862 PCBDDCGetNeumannBoundaries - Get parallel IS for Neumann boundaries 863 864 Collective 865 866 Input Parameters: 867 . pc - the preconditioning context 868 869 Output Parameters: 870 . NeumannBoundaries - index set defining the Neumann boundaries 871 872 Level: intermediate 873 874 Notes: 875 The IS returned (if any) is the same passed in earlier by the user with PCBDDCSetNeumannBoundaries 876 877 .seealso: PCBDDC 878 @*/ 879 PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc,IS *NeumannBoundaries) 880 { 881 PetscErrorCode ierr; 882 883 PetscFunctionBegin; 884 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 885 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundaries_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 886 PetscFunctionReturn(0); 887 } 888 /* -------------------------------------------------------------------------- */ 889 890 #undef __FUNCT__ 891 #define __FUNCT__ "PCBDDCGetNeumannBoundariesLocal_BDDC" 892 static PetscErrorCode PCBDDCGetNeumannBoundariesLocal_BDDC(PC pc,IS *NeumannBoundaries) 893 { 894 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 895 896 PetscFunctionBegin; 897 *NeumannBoundaries = pcbddc->NeumannBoundariesLocal; 898 PetscFunctionReturn(0); 899 } 900 901 #undef __FUNCT__ 902 #define __FUNCT__ "PCBDDCGetNeumannBoundariesLocal" 903 /*@ 904 PCBDDCGetNeumannBoundariesLocal - Get parallel IS for Neumann boundaries (in local ordering) 905 906 Collective 907 908 Input Parameters: 909 . pc - the preconditioning context 910 911 Output Parameters: 912 . NeumannBoundaries - index set defining the subdomain part of Neumann boundaries 913 914 Level: intermediate 915 916 Notes: 917 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). 918 In the latter case, the IS will be available after PCSetUp. 919 920 .seealso: PCBDDC 921 @*/ 922 PetscErrorCode PCBDDCGetNeumannBoundariesLocal(PC pc,IS *NeumannBoundaries) 923 { 924 PetscErrorCode ierr; 925 926 PetscFunctionBegin; 927 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 928 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundariesLocal_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 929 PetscFunctionReturn(0); 930 } 931 /* -------------------------------------------------------------------------- */ 932 933 #undef __FUNCT__ 934 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph_BDDC" 935 static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 936 { 937 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 938 PCBDDCGraph mat_graph = pcbddc->mat_graph; 939 PetscBool same_data = PETSC_FALSE; 940 PetscErrorCode ierr; 941 942 PetscFunctionBegin; 943 if (!nvtxs) { 944 ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr); 945 PetscFunctionReturn(0); 946 } 947 if (mat_graph->nvtxs_csr == nvtxs) { 948 if (mat_graph->xadj == xadj && mat_graph->adjncy == adjncy) same_data = PETSC_TRUE; 949 if (!same_data && mat_graph->xadj[nvtxs] == xadj[nvtxs]) { 950 ierr = PetscMemcmp(adjncy,mat_graph->adjncy,nvtxs*sizeof(PetscInt),&same_data);CHKERRQ(ierr); 951 } 952 } 953 if (!same_data) { 954 /* free old CSR */ 955 ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr); 956 /* get CSR into graph structure */ 957 if (copymode == PETSC_COPY_VALUES) { 958 ierr = PetscMalloc1(nvtxs+1,&mat_graph->xadj);CHKERRQ(ierr); 959 ierr = PetscMalloc1(xadj[nvtxs],&mat_graph->adjncy);CHKERRQ(ierr); 960 ierr = PetscMemcpy(mat_graph->xadj,xadj,(nvtxs+1)*sizeof(PetscInt));CHKERRQ(ierr); 961 ierr = PetscMemcpy(mat_graph->adjncy,adjncy,xadj[nvtxs]*sizeof(PetscInt));CHKERRQ(ierr); 962 } else if (copymode == PETSC_OWN_POINTER) { 963 mat_graph->xadj = (PetscInt*)xadj; 964 mat_graph->adjncy = (PetscInt*)adjncy; 965 } 966 mat_graph->nvtxs_csr = nvtxs; 967 pcbddc->recompute_topography = PETSC_TRUE; 968 } 969 PetscFunctionReturn(0); 970 } 971 972 #undef __FUNCT__ 973 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph" 974 /*@ 975 PCBDDCSetLocalAdjacencyGraph - Set adjacency structure (CSR graph) of the local matrix 976 977 Not collective 978 979 Input Parameters: 980 + pc - the preconditioning context 981 . nvtxs - number of local vertices of the graph (i.e., the size of the local problem) 982 . xadj, adjncy - the CSR graph 983 - copymode - either PETSC_COPY_VALUES or PETSC_OWN_POINTER. 984 985 Level: intermediate 986 987 Notes: 988 989 .seealso: PCBDDC,PetscCopyMode 990 @*/ 991 PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc,PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 992 { 993 void (*f)(void) = 0; 994 PetscErrorCode ierr; 995 996 PetscFunctionBegin; 997 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 998 if (nvtxs) { 999 PetscValidIntPointer(xadj,3); 1000 PetscValidIntPointer(adjncy,4); 1001 } 1002 if (copymode != PETSC_COPY_VALUES && copymode != PETSC_OWN_POINTER) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported copy mode %d",copymode); 1003 ierr = PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));CHKERRQ(ierr); 1004 /* free arrays if PCBDDC is not the PC type */ 1005 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);CHKERRQ(ierr); 1006 if (!f && copymode == PETSC_OWN_POINTER) { 1007 ierr = PetscFree(xadj);CHKERRQ(ierr); 1008 ierr = PetscFree(adjncy);CHKERRQ(ierr); 1009 } 1010 PetscFunctionReturn(0); 1011 } 1012 /* -------------------------------------------------------------------------- */ 1013 1014 #undef __FUNCT__ 1015 #define __FUNCT__ "PCBDDCSetDofsSplittingLocal_BDDC" 1016 static PetscErrorCode PCBDDCSetDofsSplittingLocal_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 1017 { 1018 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1019 PetscInt i; 1020 PetscBool isequal = PETSC_FALSE; 1021 PetscErrorCode ierr; 1022 1023 PetscFunctionBegin; 1024 if (pcbddc->n_ISForDofsLocal == n_is) { 1025 for (i=0;i<n_is;i++) { 1026 PetscBool isequalt; 1027 ierr = ISEqual(ISForDofs[i],pcbddc->ISForDofsLocal[i],&isequalt);CHKERRQ(ierr); 1028 if (!isequalt) break; 1029 } 1030 if (i == n_is) isequal = PETSC_TRUE; 1031 } 1032 for (i=0;i<n_is;i++) { 1033 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 1034 } 1035 /* Destroy ISes if they were already set */ 1036 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 1037 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 1038 } 1039 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 1040 /* last user setting takes precendence -> destroy any other customization */ 1041 for (i=0;i<pcbddc->n_ISForDofs;i++) { 1042 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 1043 } 1044 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 1045 pcbddc->n_ISForDofs = 0; 1046 /* allocate space then set */ 1047 if (n_is) { 1048 ierr = PetscMalloc1(n_is,&pcbddc->ISForDofsLocal);CHKERRQ(ierr); 1049 } 1050 for (i=0;i<n_is;i++) { 1051 pcbddc->ISForDofsLocal[i] = ISForDofs[i]; 1052 } 1053 pcbddc->n_ISForDofsLocal = n_is; 1054 if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE; 1055 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 1056 PetscFunctionReturn(0); 1057 } 1058 1059 #undef __FUNCT__ 1060 #define __FUNCT__ "PCBDDCSetDofsSplittingLocal" 1061 /*@ 1062 PCBDDCSetDofsSplittingLocal - Set index sets defining fields of the local subdomain matrix 1063 1064 Collective 1065 1066 Input Parameters: 1067 + pc - the preconditioning context 1068 . n_is - number of index sets defining the fields 1069 - ISForDofs - array of IS describing the fields in local ordering 1070 1071 Level: intermediate 1072 1073 Notes: 1074 n_is should be the same among processes. Not all nodes need to be listed: unlisted nodes will belong to the complement field. 1075 1076 .seealso: PCBDDC 1077 @*/ 1078 PetscErrorCode PCBDDCSetDofsSplittingLocal(PC pc,PetscInt n_is, IS ISForDofs[]) 1079 { 1080 PetscInt i; 1081 PetscErrorCode ierr; 1082 1083 PetscFunctionBegin; 1084 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1085 PetscValidLogicalCollectiveInt(pc,n_is,2); 1086 for (i=0;i<n_is;i++) { 1087 PetscCheckSameComm(pc,1,ISForDofs[i],3); 1088 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 1089 } 1090 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplittingLocal_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 1091 PetscFunctionReturn(0); 1092 } 1093 /* -------------------------------------------------------------------------- */ 1094 1095 #undef __FUNCT__ 1096 #define __FUNCT__ "PCBDDCSetDofsSplitting_BDDC" 1097 static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 1098 { 1099 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1100 PetscInt i; 1101 PetscBool isequal = PETSC_FALSE; 1102 PetscErrorCode ierr; 1103 1104 PetscFunctionBegin; 1105 if (pcbddc->n_ISForDofs == n_is) { 1106 for (i=0;i<n_is;i++) { 1107 PetscBool isequalt; 1108 ierr = ISEqual(ISForDofs[i],pcbddc->ISForDofs[i],&isequalt);CHKERRQ(ierr); 1109 if (!isequalt) break; 1110 } 1111 if (i == n_is) isequal = PETSC_TRUE; 1112 } 1113 for (i=0;i<n_is;i++) { 1114 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 1115 } 1116 /* Destroy ISes if they were already set */ 1117 for (i=0;i<pcbddc->n_ISForDofs;i++) { 1118 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 1119 } 1120 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 1121 /* last user setting takes precendence -> destroy any other customization */ 1122 for (i=0;i<pcbddc->n_ISForDofsLocal;i++) { 1123 ierr = ISDestroy(&pcbddc->ISForDofsLocal[i]);CHKERRQ(ierr); 1124 } 1125 ierr = PetscFree(pcbddc->ISForDofsLocal);CHKERRQ(ierr); 1126 pcbddc->n_ISForDofsLocal = 0; 1127 /* allocate space then set */ 1128 if (n_is) { 1129 ierr = PetscMalloc1(n_is,&pcbddc->ISForDofs);CHKERRQ(ierr); 1130 } 1131 for (i=0;i<n_is;i++) { 1132 pcbddc->ISForDofs[i] = ISForDofs[i]; 1133 } 1134 pcbddc->n_ISForDofs = n_is; 1135 if (n_is) pcbddc->user_provided_isfordofs = PETSC_TRUE; 1136 if (!isequal) pcbddc->recompute_topography = PETSC_TRUE; 1137 PetscFunctionReturn(0); 1138 } 1139 1140 #undef __FUNCT__ 1141 #define __FUNCT__ "PCBDDCSetDofsSplitting" 1142 /*@ 1143 PCBDDCSetDofsSplitting - Set index sets defining fields of the global matrix 1144 1145 Collective 1146 1147 Input Parameters: 1148 + pc - the preconditioning context 1149 . n_is - number of index sets defining the fields 1150 - ISForDofs - array of IS describing the fields in global ordering 1151 1152 Level: intermediate 1153 1154 Notes: 1155 Any process can list any global node. Not all nodes need to be listed: unlisted nodes will belong to the complement field. 1156 1157 .seealso: PCBDDC 1158 @*/ 1159 PetscErrorCode PCBDDCSetDofsSplitting(PC pc,PetscInt n_is, IS ISForDofs[]) 1160 { 1161 PetscInt i; 1162 PetscErrorCode ierr; 1163 1164 PetscFunctionBegin; 1165 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1166 PetscValidLogicalCollectiveInt(pc,n_is,2); 1167 for (i=0;i<n_is;i++) { 1168 PetscCheckSameComm(pc,1,ISForDofs[i],3); 1169 PetscValidHeaderSpecific(ISForDofs[i],IS_CLASSID,3); 1170 } 1171 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 1172 PetscFunctionReturn(0); 1173 } 1174 1175 /* -------------------------------------------------------------------------- */ 1176 #undef __FUNCT__ 1177 #define __FUNCT__ "PCPreSolve_BDDC" 1178 /* -------------------------------------------------------------------------- */ 1179 /* 1180 PCPreSolve_BDDC - Changes the right hand side and (if necessary) the initial 1181 guess if a transformation of basis approach has been selected. 1182 1183 Input Parameter: 1184 + pc - the preconditioner contex 1185 1186 Application Interface Routine: PCPreSolve() 1187 1188 Notes: 1189 The interface routine PCPreSolve() is not usually called directly by 1190 the user, but instead is called by KSPSolve(). 1191 */ 1192 static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 1193 { 1194 PetscErrorCode ierr; 1195 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1196 PC_IS *pcis = (PC_IS*)(pc->data); 1197 Vec used_vec; 1198 PetscBool save_rhs = PETSC_TRUE, benign_correction_computed; 1199 1200 PetscFunctionBegin; 1201 /* if we are working with CG, one dirichlet solve can be avoided during Krylov iterations */ 1202 if (ksp) { 1203 PetscBool iscg, isgroppcg, ispipecg, ispipecgrr; 1204 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPCG,&iscg);CHKERRQ(ierr); 1205 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPGROPPCG,&isgroppcg);CHKERRQ(ierr); 1206 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECG,&ispipecg);CHKERRQ(ierr); 1207 ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPIPECGRR,&ispipecgrr);CHKERRQ(ierr); 1208 if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static || (!iscg && !isgroppcg && !ispipecg && !ispipecgrr)) { 1209 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 1210 } 1211 } 1212 if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static) { 1213 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 1214 } 1215 1216 /* Creates parallel work vectors used in presolve */ 1217 if (!pcbddc->original_rhs) { 1218 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 1219 } 1220 if (!pcbddc->temp_solution) { 1221 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);CHKERRQ(ierr); 1222 } 1223 1224 pcbddc->temp_solution_used = PETSC_FALSE; 1225 if (x) { 1226 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr); 1227 used_vec = x; 1228 } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */ 1229 ierr = PetscObjectReference((PetscObject)pcbddc->temp_solution);CHKERRQ(ierr); 1230 used_vec = pcbddc->temp_solution; 1231 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 1232 pcbddc->temp_solution_used = PETSC_TRUE; 1233 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1234 save_rhs = PETSC_FALSE; 1235 pcbddc->eliminate_dirdofs = PETSC_TRUE; 1236 } 1237 1238 /* hack into ksp data structure since PCPreSolve comes earlier than setting to zero the guess in src/ksp/ksp/interface/itfunc.c */ 1239 if (ksp) { 1240 /* store the flag for the initial guess since it will be restored back during PCPostSolve_BDDC */ 1241 ierr = KSPGetInitialGuessNonzero(ksp,&pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 1242 if (!pcbddc->ksp_guess_nonzero) { 1243 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 1244 } 1245 } 1246 1247 pcbddc->rhs_change = PETSC_FALSE; 1248 /* Take into account zeroed rows -> change rhs and store solution removed */ 1249 if (rhs && pcbddc->eliminate_dirdofs) { 1250 IS dirIS = NULL; 1251 1252 /* DirichletBoundariesLocal may not be consistent among neighbours; gets a dirichlet dofs IS from graph (may be cached) */ 1253 ierr = PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph,&dirIS);CHKERRQ(ierr); 1254 if (dirIS) { 1255 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 1256 PetscInt dirsize,i,*is_indices; 1257 PetscScalar *array_x; 1258 const PetscScalar *array_diagonal; 1259 1260 ierr = MatGetDiagonal(pc->pmat,pcis->vec1_global);CHKERRQ(ierr); 1261 ierr = VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);CHKERRQ(ierr); 1262 ierr = VecScatterBegin(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1263 ierr = VecScatterEnd(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1264 ierr = VecScatterBegin(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1265 ierr = VecScatterEnd(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1266 ierr = ISGetLocalSize(dirIS,&dirsize);CHKERRQ(ierr); 1267 ierr = VecGetArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 1268 ierr = VecGetArrayRead(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 1269 ierr = ISGetIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1270 for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]]; 1271 ierr = ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1272 ierr = VecRestoreArrayRead(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 1273 ierr = VecRestoreArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 1274 ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1275 ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1276 pcbddc->rhs_change = PETSC_TRUE; 1277 ierr = ISDestroy(&dirIS);CHKERRQ(ierr); 1278 } 1279 } 1280 1281 /* remove the computed solution or the initial guess from the rhs */ 1282 if (pcbddc->rhs_change || (ksp && pcbddc->ksp_guess_nonzero) ) { 1283 /* save the original rhs */ 1284 if (save_rhs) { 1285 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1286 save_rhs = PETSC_FALSE; 1287 } 1288 pcbddc->rhs_change = PETSC_TRUE; 1289 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 1290 ierr = MatMultAdd(pc->mat,used_vec,pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1291 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 1292 ierr = VecCopy(used_vec,pcbddc->temp_solution);CHKERRQ(ierr); 1293 pcbddc->temp_solution_used = PETSC_TRUE; 1294 if (ksp) { 1295 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_FALSE);CHKERRQ(ierr); 1296 } 1297 } 1298 ierr = VecDestroy(&used_vec);CHKERRQ(ierr); 1299 1300 /* compute initial vector in benign space if needed 1301 and remove non-benign solution from the rhs */ 1302 benign_correction_computed = PETSC_FALSE; 1303 if (rhs && pcbddc->benign_compute_correction && pcbddc->benign_have_null) { 1304 /* compute u^*_h using ideas similar to those in Xuemin Tu's PhD thesis (see Section 4.8.1) 1305 Recursively apply BDDC in the multilevel case */ 1306 if (!pcbddc->benign_vec) { 1307 ierr = VecDuplicate(rhs,&pcbddc->benign_vec);CHKERRQ(ierr); 1308 } 1309 pcbddc->benign_apply_coarse_only = PETSC_TRUE; 1310 if (!pcbddc->benign_skip_correction) { 1311 ierr = PCApply_BDDC(pc,rhs,pcbddc->benign_vec);CHKERRQ(ierr); 1312 benign_correction_computed = PETSC_TRUE; 1313 if (pcbddc->temp_solution_used) { 1314 ierr = VecAXPY(pcbddc->temp_solution,1.0,pcbddc->benign_vec);CHKERRQ(ierr); 1315 } 1316 ierr = VecScale(pcbddc->benign_vec,-1.0);CHKERRQ(ierr); 1317 /* store the original rhs if not done earlier */ 1318 if (save_rhs) { 1319 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1320 save_rhs = PETSC_FALSE; 1321 } 1322 if (pcbddc->rhs_change) { 1323 ierr = MatMultAdd(pc->mat,pcbddc->benign_vec,rhs,rhs);CHKERRQ(ierr); 1324 } else { 1325 ierr = MatMultAdd(pc->mat,pcbddc->benign_vec,pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1326 } 1327 pcbddc->rhs_change = PETSC_TRUE; 1328 } 1329 pcbddc->benign_apply_coarse_only = PETSC_FALSE; 1330 } 1331 1332 /* dbg output */ 1333 if (pcbddc->dbg_flag && benign_correction_computed) { 1334 Vec v; 1335 ierr = VecDuplicate(pcis->vec1_global,&v);CHKERRQ(ierr); 1336 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,rhs,v);CHKERRQ(ierr); 1337 ierr = PCBDDCBenignGetOrSetP0(pc,v,PETSC_TRUE);CHKERRQ(ierr); 1338 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"LEVEL %d: is the correction benign?\n",pcbddc->current_level);CHKERRQ(ierr); 1339 ierr = PetscScalarView(pcbddc->benign_n,pcbddc->benign_p0,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)));CHKERRQ(ierr); 1340 ierr = VecDestroy(&v);CHKERRQ(ierr); 1341 } 1342 1343 /* set initial guess if using PCG */ 1344 if (x && pcbddc->use_exact_dirichlet_trick) { 1345 ierr = VecSet(x,0.0);CHKERRQ(ierr); 1346 if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) { 1347 if (benign_correction_computed) { /* we have already saved the changed rhs */ 1348 ierr = VecLockPop(pcis->vec1_global);CHKERRQ(ierr); 1349 } else { 1350 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,rhs,pcis->vec1_global);CHKERRQ(ierr); 1351 } 1352 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_global,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1353 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_global,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1354 } else { 1355 ierr = VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1356 ierr = VecScatterEnd(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1357 } 1358 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1359 if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) { 1360 ierr = VecSet(pcis->vec1_global,0.);CHKERRQ(ierr); 1361 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1362 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1363 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcis->vec1_global,x);CHKERRQ(ierr); 1364 } else { 1365 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1366 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1367 } 1368 if (ksp) { 1369 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 1370 } 1371 } else if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior && benign_correction_computed && pcbddc->use_exact_dirichlet_trick) { 1372 ierr = VecLockPop(pcis->vec1_global);CHKERRQ(ierr); 1373 } 1374 PetscFunctionReturn(0); 1375 } 1376 1377 /* -------------------------------------------------------------------------- */ 1378 #undef __FUNCT__ 1379 #define __FUNCT__ "PCPostSolve_BDDC" 1380 /* -------------------------------------------------------------------------- */ 1381 /* 1382 PCPostSolve_BDDC - Changes the computed solution if a transformation of basis 1383 approach has been selected. Also, restores rhs to its original state. 1384 1385 Input Parameter: 1386 + pc - the preconditioner contex 1387 1388 Application Interface Routine: PCPostSolve() 1389 1390 Notes: 1391 The interface routine PCPostSolve() is not usually called directly by 1392 the user, but instead is called by KSPSolve(). 1393 */ 1394 static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 1395 { 1396 PetscErrorCode ierr; 1397 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1398 1399 PetscFunctionBegin; 1400 /* add solution removed in presolve */ 1401 if (x && pcbddc->rhs_change) { 1402 if (pcbddc->temp_solution_used) { 1403 ierr = VecAXPY(x,1.0,pcbddc->temp_solution);CHKERRQ(ierr); 1404 } else if (pcbddc->benign_compute_correction) { 1405 ierr = VecAXPY(x,-1.0,pcbddc->benign_vec);CHKERRQ(ierr); 1406 } 1407 } 1408 pcbddc->temp_solution_used = PETSC_FALSE; 1409 1410 /* restore rhs to its original state (not needed for FETI-DP) */ 1411 if (rhs && pcbddc->rhs_change) { 1412 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1413 } 1414 pcbddc->rhs_change = PETSC_FALSE; 1415 /* restore ksp guess state */ 1416 if (ksp) { 1417 ierr = KSPSetInitialGuessNonzero(ksp,pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 1418 } 1419 PetscFunctionReturn(0); 1420 } 1421 /* -------------------------------------------------------------------------- */ 1422 #undef __FUNCT__ 1423 #define __FUNCT__ "PCSetUp_BDDC" 1424 /* -------------------------------------------------------------------------- */ 1425 /* 1426 PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner 1427 by setting data structures and options. 1428 1429 Input Parameter: 1430 + pc - the preconditioner context 1431 1432 Application Interface Routine: PCSetUp() 1433 1434 Notes: 1435 The interface routine PCSetUp() is not usually called directly by 1436 the user, but instead is called by PCApply() if necessary. 1437 */ 1438 PetscErrorCode PCSetUp_BDDC(PC pc) 1439 { 1440 PetscErrorCode ierr; 1441 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 1442 Mat_IS* matis; 1443 MatNullSpace nearnullspace; 1444 IS zerodiag = NULL; 1445 PetscInt nrows,ncols; 1446 PetscBool computetopography,computesolvers,computesubschurs; 1447 PetscBool computeconstraintsmatrix; 1448 PetscBool new_nearnullspace_provided,ismatis; 1449 1450 PetscFunctionBegin; 1451 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&ismatis);CHKERRQ(ierr); 1452 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PCBDDC preconditioner requires matrix of type MATIS"); 1453 ierr = MatGetSize(pc->pmat,&nrows,&ncols);CHKERRQ(ierr); 1454 if (nrows != ncols) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCBDDC preconditioner requires a square preconditioning matrix"); 1455 matis = (Mat_IS*)pc->pmat->data; 1456 /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */ 1457 /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup 1458 Also, BDDC directly build the Dirichlet problem */ 1459 /* split work */ 1460 if (pc->setupcalled) { 1461 if (pc->flag == SAME_NONZERO_PATTERN) { 1462 computetopography = PETSC_FALSE; 1463 computesolvers = PETSC_TRUE; 1464 } else { /* DIFFERENT_NONZERO_PATTERN */ 1465 computetopography = PETSC_TRUE; 1466 computesolvers = PETSC_TRUE; 1467 } 1468 } else { 1469 computetopography = PETSC_TRUE; 1470 computesolvers = PETSC_TRUE; 1471 } 1472 if (pcbddc->recompute_topography) { 1473 computetopography = PETSC_TRUE; 1474 } 1475 pcbddc->recompute_topography = computetopography; 1476 computeconstraintsmatrix = PETSC_FALSE; 1477 1478 /* check parameters' compatibility */ 1479 if (!pcbddc->use_deluxe_scaling) pcbddc->deluxe_zerorows = PETSC_FALSE; 1480 pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_threshold > 0.0); 1481 pcbddc->adaptive_userdefined = (PetscBool)(pcbddc->adaptive_selection && pcbddc->adaptive_userdefined); 1482 if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE; 1483 1484 computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling); 1485 if (pcbddc->switch_static) { 1486 PetscBool ismatis; 1487 ierr = PetscObjectTypeCompare((PetscObject)pc->mat,MATIS,&ismatis);CHKERRQ(ierr); 1488 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"When the static switch is one, the iteration matrix should be of type MATIS"); 1489 } 1490 1491 /* Get stdout for dbg */ 1492 if (pcbddc->dbg_flag) { 1493 if (!pcbddc->dbg_viewer) { 1494 pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)); 1495 ierr = PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer);CHKERRQ(ierr); 1496 } 1497 ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1498 } 1499 1500 if (pcbddc->user_ChangeOfBasisMatrix) { 1501 /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */ 1502 pcbddc->use_change_of_basis = PETSC_FALSE; 1503 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 1504 } else { 1505 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1506 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1507 pcbddc->local_mat = matis->A; 1508 } 1509 1510 /* detect local disconnected subdomains if requested and not done before */ 1511 if (pcbddc->detect_disconnected && !pcbddc->n_local_subs) { 1512 ierr = MatDetectDisconnectedComponents(pcbddc->local_mat,PETSC_FALSE,&pcbddc->n_local_subs,&pcbddc->local_subs);CHKERRQ(ierr); 1513 } 1514 1515 /* compute topology info in local ordering */ 1516 if (pcbddc->recompute_topography) { 1517 ierr = PCBDDCComputeLocalTopologyInfo(pc);CHKERRQ(ierr); 1518 } 1519 1520 /* 1521 Compute change of basis on local pressures (aka zerodiag dofs) 1522 This should come earlier then PCISSetUp for extracting the correct subdomain matrices 1523 */ 1524 ierr = PCBDDCBenignShellMat(pc,PETSC_TRUE);CHKERRQ(ierr); 1525 if (pcbddc->benign_saddle_point) { 1526 PC_IS* pcis = (PC_IS*)pc->data; 1527 1528 if (pcbddc->user_ChangeOfBasisMatrix || pcbddc->use_change_of_basis || !computesubschurs) pcbddc->benign_change_explicit = PETSC_TRUE; 1529 /* detect local saddle point and change the basis in pcbddc->local_mat (TODO: reuse case) */ 1530 ierr = PCBDDCBenignDetectSaddlePoint(pc,&zerodiag);CHKERRQ(ierr); 1531 /* pop B0 mat from local mat */ 1532 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1533 /* give pcis a hint to not reuse submatrices during PCISCreate */ 1534 if (pc->flag == SAME_NONZERO_PATTERN && pcis->reusesubmatrices == PETSC_TRUE) { 1535 if (pcbddc->benign_n && (pcbddc->benign_change_explicit || pcbddc->dbg_flag)) { 1536 pcis->reusesubmatrices = PETSC_FALSE; 1537 } else { 1538 pcis->reusesubmatrices = PETSC_TRUE; 1539 } 1540 } else { 1541 pcis->reusesubmatrices = PETSC_FALSE; 1542 } 1543 } 1544 1545 /* propagate relevant information */ 1546 #if !defined(PETSC_USE_COMPLEX) /* workaround for reals */ 1547 if (matis->A->symmetric_set) { 1548 ierr = MatSetOption(pcbddc->local_mat,MAT_HERMITIAN,matis->A->symmetric);CHKERRQ(ierr); 1549 } 1550 #endif 1551 if (matis->A->symmetric_set) { 1552 ierr = MatSetOption(pcbddc->local_mat,MAT_SYMMETRIC,matis->A->symmetric);CHKERRQ(ierr); 1553 } 1554 if (matis->A->spd_set) { 1555 ierr = MatSetOption(pcbddc->local_mat,MAT_SPD,matis->A->spd);CHKERRQ(ierr); 1556 } 1557 1558 /* Set up all the "iterative substructuring" common block without computing solvers */ 1559 { 1560 Mat temp_mat; 1561 1562 temp_mat = matis->A; 1563 matis->A = pcbddc->local_mat; 1564 ierr = PCISSetUp(pc,PETSC_FALSE);CHKERRQ(ierr); 1565 pcbddc->local_mat = matis->A; 1566 matis->A = temp_mat; 1567 } 1568 1569 /* Analyze interface */ 1570 if (computetopography) { 1571 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 1572 computeconstraintsmatrix = PETSC_TRUE; 1573 if (pcbddc->adaptive_selection && !pcbddc->use_deluxe_scaling && !pcbddc->mat_graph->twodim) { 1574 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot compute the adaptive primal space for a problem with 3D edges without deluxe scaling"); 1575 } 1576 if (pcbddc->compute_nonetflux) { 1577 MatNullSpace nnfnnsp; 1578 1579 ierr = PCBDDCComputeNoNetFlux(pc->pmat,pcbddc->divudotp,pcbddc->mat_graph,&nnfnnsp);CHKERRQ(ierr); 1580 ierr = MatSetNearNullSpace(pc->pmat,nnfnnsp);CHKERRQ(ierr); 1581 ierr = MatNullSpaceDestroy(&nnfnnsp);CHKERRQ(ierr); 1582 } 1583 } 1584 1585 /* check existence of a divergence free extension, i.e. 1586 b(v_I,p_0) = 0 for all v_I (raise error if not). 1587 Also, check that PCBDDCBenignGetOrSetP0 works */ 1588 #if defined(PETSC_USE_DEBUG) 1589 if (pcbddc->benign_saddle_point) { 1590 ierr = PCBDDCBenignCheck(pc,zerodiag);CHKERRQ(ierr); 1591 } 1592 #endif 1593 ierr = ISDestroy(&zerodiag);CHKERRQ(ierr); 1594 1595 /* Setup local dirichlet solver ksp_D and sub_schurs solvers */ 1596 if (computesolvers) { 1597 PCBDDCSubSchurs sub_schurs; 1598 1599 if (computesubschurs && computetopography) { 1600 ierr = PCBDDCInitSubSchurs(pc);CHKERRQ(ierr); 1601 } 1602 /* SetUp Scaling operator (scaling matrices could be needed in SubSchursSetUp)*/ 1603 if (!pcbddc->use_deluxe_scaling) { 1604 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1605 } 1606 sub_schurs = pcbddc->sub_schurs; 1607 if (sub_schurs && sub_schurs->schur_explicit) { 1608 if (computesubschurs) { 1609 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1610 } 1611 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1612 } else { 1613 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1614 if (computesubschurs) { 1615 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1616 } 1617 } 1618 if (pcbddc->adaptive_selection) { 1619 ierr = PCBDDCAdaptiveSelection(pc);CHKERRQ(ierr); 1620 computeconstraintsmatrix = PETSC_TRUE; 1621 } 1622 } 1623 1624 /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */ 1625 new_nearnullspace_provided = PETSC_FALSE; 1626 ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr); 1627 if (pcbddc->onearnullspace) { /* already used nearnullspace */ 1628 if (!nearnullspace) { /* near null space attached to mat has been destroyed */ 1629 new_nearnullspace_provided = PETSC_TRUE; 1630 } else { 1631 /* determine if the two nullspaces are different (should be lightweight) */ 1632 if (nearnullspace != pcbddc->onearnullspace) { 1633 new_nearnullspace_provided = PETSC_TRUE; 1634 } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */ 1635 PetscInt i; 1636 const Vec *nearnullvecs; 1637 PetscObjectState state; 1638 PetscInt nnsp_size; 1639 ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 1640 for (i=0;i<nnsp_size;i++) { 1641 ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr); 1642 if (pcbddc->onearnullvecs_state[i] != state) { 1643 new_nearnullspace_provided = PETSC_TRUE; 1644 break; 1645 } 1646 } 1647 } 1648 } 1649 } else { 1650 if (!nearnullspace) { /* both nearnullspaces are null */ 1651 new_nearnullspace_provided = PETSC_FALSE; 1652 } else { /* nearnullspace attached later */ 1653 new_nearnullspace_provided = PETSC_TRUE; 1654 } 1655 } 1656 1657 /* Setup constraints and related work vectors */ 1658 /* reset primal space flags */ 1659 pcbddc->new_primal_space = PETSC_FALSE; 1660 pcbddc->new_primal_space_local = PETSC_FALSE; 1661 if (computeconstraintsmatrix || new_nearnullspace_provided) { 1662 /* It also sets the primal space flags */ 1663 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 1664 /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */ 1665 ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr); 1666 } 1667 1668 if (computesolvers || pcbddc->new_primal_space) { 1669 if (pcbddc->use_change_of_basis) { 1670 PC_IS *pcis = (PC_IS*)(pc->data); 1671 1672 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 1673 if (pcbddc->benign_change) { 1674 ierr = MatDestroy(&pcbddc->benign_B0);CHKERRQ(ierr); 1675 /* pop B0 from pcbddc->local_mat */ 1676 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1677 } 1678 /* get submatrices */ 1679 ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr); 1680 ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr); 1681 ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr); 1682 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);CHKERRQ(ierr); 1683 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);CHKERRQ(ierr); 1684 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);CHKERRQ(ierr); 1685 /* set flag in pcis to not reuse submatrices during PCISCreate */ 1686 pcis->reusesubmatrices = PETSC_FALSE; 1687 } else if (!pcbddc->user_ChangeOfBasisMatrix && !pcbddc->benign_change) { 1688 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1689 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1690 pcbddc->local_mat = matis->A; 1691 } 1692 /* SetUp coarse and local Neumann solvers */ 1693 ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr); 1694 /* SetUp Scaling operator */ 1695 if (pcbddc->use_deluxe_scaling) { 1696 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1697 } 1698 } 1699 /* mark topography as done */ 1700 pcbddc->recompute_topography = PETSC_FALSE; 1701 1702 /* wrap pcis->A_IB and pcis->A_BI if we did not change explicitly the variables on the pressures */ 1703 ierr = PCBDDCBenignShellMat(pc,PETSC_FALSE);CHKERRQ(ierr); 1704 1705 if (pcbddc->dbg_flag) { 1706 ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1707 } 1708 PetscFunctionReturn(0); 1709 } 1710 1711 /* -------------------------------------------------------------------------- */ 1712 /* 1713 PCApply_BDDC - Applies the BDDC operator to a vector. 1714 1715 Input Parameters: 1716 + pc - the preconditioner context 1717 - r - input vector (global) 1718 1719 Output Parameter: 1720 . z - output vector (global) 1721 1722 Application Interface Routine: PCApply() 1723 */ 1724 #undef __FUNCT__ 1725 #define __FUNCT__ "PCApply_BDDC" 1726 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 1727 { 1728 PC_IS *pcis = (PC_IS*)(pc->data); 1729 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1730 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1731 PetscErrorCode ierr; 1732 const PetscScalar one = 1.0; 1733 const PetscScalar m_one = -1.0; 1734 const PetscScalar zero = 0.0; 1735 1736 /* This code is similar to that provided in nn.c for PCNN 1737 NN interface preconditioner changed to BDDC 1738 Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */ 1739 1740 PetscFunctionBegin; 1741 if (pcbddc->ChangeOfBasisMatrix) { 1742 Vec swap; 1743 1744 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 1745 swap = pcbddc->work_change; 1746 pcbddc->work_change = r; 1747 r = swap; 1748 /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */ 1749 if (pcbddc->benign_apply_coarse_only && pcbddc->use_exact_dirichlet_trick && pcbddc->change_interior) { 1750 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 1751 ierr = VecLockPush(pcis->vec1_global);CHKERRQ(ierr); 1752 } 1753 } 1754 if (pcbddc->benign_have_null) { /* get p0 from r */ 1755 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 1756 } 1757 if (!pcbddc->use_exact_dirichlet_trick && !pcbddc->benign_apply_coarse_only) { 1758 ierr = VecCopy(r,z);CHKERRQ(ierr); 1759 /* First Dirichlet solve */ 1760 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1761 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1762 /* 1763 Assembling right hand side for BDDC operator 1764 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1765 - pcis->vec1_B the interface part of the global vector z 1766 */ 1767 if (n_D) { 1768 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1769 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1770 if (pcbddc->switch_static) { 1771 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1772 1773 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 1774 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1775 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1776 if (!pcbddc->switch_static_change) { 1777 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1778 } else { 1779 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1780 ierr = MatMult(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1781 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1782 } 1783 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1784 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1785 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1786 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1787 } else { 1788 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1789 } 1790 } else { 1791 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1792 } 1793 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1794 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1795 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1796 } else { 1797 if (!pcbddc->benign_apply_coarse_only) { 1798 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1799 } 1800 } 1801 1802 /* Apply interface preconditioner 1803 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1804 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr); 1805 1806 /* Apply transpose of partition of unity operator */ 1807 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1808 1809 /* Second Dirichlet solve and assembling of output */ 1810 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1811 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1812 if (n_B) { 1813 if (pcbddc->switch_static) { 1814 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1815 1816 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1817 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1818 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1819 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1820 if (!pcbddc->switch_static_change) { 1821 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1822 } else { 1823 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1824 ierr = MatMult(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1825 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1826 } 1827 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1828 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1829 } else { 1830 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1831 } 1832 } else if (pcbddc->switch_static) { /* n_B is zero */ 1833 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1834 1835 if (!pcbddc->switch_static_change) { 1836 ierr = MatMult(matis->A,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 1837 } else { 1838 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 1839 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1840 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 1841 } 1842 } 1843 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1844 1845 if (!pcbddc->use_exact_dirichlet_trick && !pcbddc->benign_apply_coarse_only) { 1846 if (pcbddc->switch_static) { 1847 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 1848 } else { 1849 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 1850 } 1851 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1852 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1853 } else { 1854 if (pcbddc->switch_static) { 1855 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 1856 } else { 1857 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1858 } 1859 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1860 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1861 } 1862 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 1863 if (pcbddc->benign_apply_coarse_only) { 1864 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 1865 } 1866 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 1867 } 1868 1869 if (pcbddc->ChangeOfBasisMatrix) { 1870 Vec swap; 1871 1872 swap = r; 1873 r = pcbddc->work_change; 1874 pcbddc->work_change = swap; 1875 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 1876 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 1877 } 1878 PetscFunctionReturn(0); 1879 } 1880 1881 /* -------------------------------------------------------------------------- */ 1882 /* 1883 PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector. 1884 1885 Input Parameters: 1886 + pc - the preconditioner context 1887 - r - input vector (global) 1888 1889 Output Parameter: 1890 . z - output vector (global) 1891 1892 Application Interface Routine: PCApplyTranspose() 1893 */ 1894 #undef __FUNCT__ 1895 #define __FUNCT__ "PCApplyTranspose_BDDC" 1896 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z) 1897 { 1898 PC_IS *pcis = (PC_IS*)(pc->data); 1899 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1900 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1901 PetscErrorCode ierr; 1902 const PetscScalar one = 1.0; 1903 const PetscScalar m_one = -1.0; 1904 const PetscScalar zero = 0.0; 1905 1906 PetscFunctionBegin; 1907 if (pcbddc->ChangeOfBasisMatrix) { 1908 Vec swap; 1909 1910 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 1911 swap = pcbddc->work_change; 1912 pcbddc->work_change = r; 1913 r = swap; 1914 /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */ 1915 if (pcbddc->benign_apply_coarse_only && pcbddc->use_exact_dirichlet_trick && pcbddc->change_interior) { 1916 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 1917 ierr = VecLockPush(pcis->vec1_global);CHKERRQ(ierr); 1918 } 1919 } 1920 if (pcbddc->benign_have_null) { /* get p0 from r */ 1921 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 1922 } 1923 if (!pcbddc->use_exact_dirichlet_trick && !pcbddc->benign_apply_coarse_only) { 1924 ierr = VecCopy(r,z);CHKERRQ(ierr); 1925 /* First Dirichlet solve */ 1926 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1927 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1928 /* 1929 Assembling right hand side for BDDC operator 1930 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1931 - pcis->vec1_B the interface part of the global vector z 1932 */ 1933 if (n_D) { 1934 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1935 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1936 if (pcbddc->switch_static) { 1937 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1938 1939 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 1940 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1941 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1942 if (!pcbddc->switch_static_change) { 1943 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1944 } else { 1945 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1946 ierr = MatMultTranspose(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1947 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1948 } 1949 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1950 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1951 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1952 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1953 } else { 1954 ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1955 } 1956 } else { 1957 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1958 } 1959 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1960 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1961 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1962 } else { 1963 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1964 } 1965 1966 /* Apply interface preconditioner 1967 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1968 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr); 1969 1970 /* Apply transpose of partition of unity operator */ 1971 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1972 1973 /* Second Dirichlet solve and assembling of output */ 1974 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1975 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1976 if (n_B) { 1977 if (pcbddc->switch_static) { 1978 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1979 1980 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1981 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1982 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1983 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1984 if (!pcbddc->switch_static_change) { 1985 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1986 } else { 1987 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1988 ierr = MatMultTranspose(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1989 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1990 } 1991 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1992 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1993 } else { 1994 ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1995 } 1996 } else if (pcbddc->switch_static) { /* n_B is zero */ 1997 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1998 1999 if (!pcbddc->switch_static_change) { 2000 ierr = MatMultTranspose(matis->A,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 2001 } else { 2002 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 2003 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 2004 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 2005 } 2006 } 2007 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 2008 if (!pcbddc->use_exact_dirichlet_trick && !pcbddc->benign_apply_coarse_only) { 2009 if (pcbddc->switch_static) { 2010 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 2011 } else { 2012 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 2013 } 2014 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2015 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2016 } else { 2017 if (pcbddc->switch_static) { 2018 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 2019 } else { 2020 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 2021 } 2022 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2023 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2024 } 2025 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 2026 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 2027 } 2028 if (pcbddc->ChangeOfBasisMatrix) { 2029 Vec swap; 2030 2031 swap = r; 2032 r = pcbddc->work_change; 2033 pcbddc->work_change = swap; 2034 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 2035 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 2036 } 2037 PetscFunctionReturn(0); 2038 } 2039 /* -------------------------------------------------------------------------- */ 2040 2041 #undef __FUNCT__ 2042 #define __FUNCT__ "PCDestroy_BDDC" 2043 PetscErrorCode PCDestroy_BDDC(PC pc) 2044 { 2045 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 2046 PetscErrorCode ierr; 2047 2048 PetscFunctionBegin; 2049 /* free BDDC custom data */ 2050 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 2051 /* destroy objects related to topography */ 2052 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 2053 /* free allocated graph structure */ 2054 ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr); 2055 /* free allocated sub schurs structure */ 2056 ierr = PetscFree(pcbddc->sub_schurs);CHKERRQ(ierr); 2057 /* destroy objects for scaling operator */ 2058 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 2059 ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr); 2060 /* free solvers stuff */ 2061 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 2062 /* free global vectors needed in presolve */ 2063 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 2064 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 2065 /* free data created by PCIS */ 2066 ierr = PCISDestroy(pc);CHKERRQ(ierr); 2067 /* remove functions */ 2068 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",NULL);CHKERRQ(ierr); 2069 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);CHKERRQ(ierr); 2070 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 2071 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",NULL);CHKERRQ(ierr); 2072 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 2073 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr); 2074 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr); 2075 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr); 2076 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2077 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2078 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2079 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2080 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2081 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2082 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2083 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2084 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 2085 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr); 2086 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 2087 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 2088 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 2089 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 2090 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",NULL);CHKERRQ(ierr); 2091 /* Free the private data structure */ 2092 ierr = PetscFree(pc->data);CHKERRQ(ierr); 2093 PetscFunctionReturn(0); 2094 } 2095 /* -------------------------------------------------------------------------- */ 2096 2097 #undef __FUNCT__ 2098 #define __FUNCT__ "PCPreSolveChangeRHS_BDDC" 2099 static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool* change) 2100 { 2101 PetscFunctionBegin; 2102 *change = PETSC_TRUE; 2103 PetscFunctionReturn(0); 2104 } 2105 2106 #undef __FUNCT__ 2107 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC" 2108 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2109 { 2110 FETIDPMat_ctx mat_ctx; 2111 Vec work; 2112 PC_IS* pcis; 2113 PC_BDDC* pcbddc; 2114 PetscErrorCode ierr; 2115 2116 PetscFunctionBegin; 2117 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2118 pcis = (PC_IS*)mat_ctx->pc->data; 2119 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2120 2121 /* 2122 change of basis for physical rhs if needed 2123 It also changes the rhs in case of dirichlet boundaries 2124 */ 2125 if (!pcbddc->original_rhs) { 2126 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 2127 } 2128 ierr = VecCopy(standard_rhs,pcbddc->original_rhs);CHKERRQ(ierr); 2129 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,pcbddc->original_rhs,NULL);CHKERRQ(ierr); 2130 if (pcbddc->ChangeOfBasisMatrix) { 2131 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,pcbddc->original_rhs,pcbddc->work_change);CHKERRQ(ierr); 2132 work = pcbddc->work_change; 2133 } else { 2134 work = pcbddc->original_rhs; 2135 } 2136 /* store vectors for computation of fetidp final solution */ 2137 ierr = VecScatterBegin(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2138 ierr = VecScatterEnd(pcis->global_to_D,work,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2139 /* scale rhs since it should be unassembled */ 2140 /* TODO use counter scaling? (also below) */ 2141 ierr = VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2142 ierr = VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2143 /* Apply partition of unity */ 2144 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2145 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2146 if (!pcbddc->switch_static) { 2147 /* compute partially subassembled Schur complement right-hand side */ 2148 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2149 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 2150 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 2151 ierr = VecSet(work,0.0);CHKERRQ(ierr); 2152 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2153 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,work,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2154 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,work,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2155 ierr = VecScatterBegin(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2156 ierr = VecScatterEnd(pcis->global_to_B,work,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2157 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2158 } 2159 /* BDDC rhs */ 2160 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 2161 if (pcbddc->switch_static) { 2162 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2163 } 2164 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2165 /* apply BDDC */ 2166 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2167 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2168 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 2169 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 2170 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 2171 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2172 ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2173 PetscFunctionReturn(0); 2174 } 2175 2176 #undef __FUNCT__ 2177 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS" 2178 /*@ 2179 PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETI-DP linear system using the physical right-hand side 2180 2181 Collective 2182 2183 Input Parameters: 2184 + fetidp_mat - the FETI-DP matrix object obtained by a call to PCBDDCCreateFETIDPOperators 2185 - standard_rhs - the right-hand side of the original linear system 2186 2187 Output Parameters: 2188 . fetidp_flux_rhs - the right-hand side for the FETI-DP linear system 2189 2190 Level: developer 2191 2192 Notes: 2193 2194 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetSolution 2195 @*/ 2196 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2197 { 2198 FETIDPMat_ctx mat_ctx; 2199 PetscErrorCode ierr; 2200 2201 PetscFunctionBegin; 2202 PetscValidHeaderSpecific(fetidp_mat,MAT_CLASSID,1); 2203 PetscValidHeaderSpecific(standard_rhs,VEC_CLASSID,2); 2204 PetscValidHeaderSpecific(fetidp_flux_rhs,VEC_CLASSID,3); 2205 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2206 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 2207 PetscFunctionReturn(0); 2208 } 2209 /* -------------------------------------------------------------------------- */ 2210 2211 #undef __FUNCT__ 2212 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC" 2213 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2214 { 2215 FETIDPMat_ctx mat_ctx; 2216 PC_IS* pcis; 2217 PC_BDDC* pcbddc; 2218 PetscErrorCode ierr; 2219 2220 PetscFunctionBegin; 2221 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2222 pcis = (PC_IS*)mat_ctx->pc->data; 2223 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2224 2225 /* apply B_delta^T */ 2226 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2227 ierr = VecScatterEnd(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2228 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 2229 /* compute rhs for BDDC application */ 2230 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2231 if (pcbddc->switch_static) { 2232 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2233 } 2234 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2235 /* apply BDDC */ 2236 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2237 /* put values into standard global vector */ 2238 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2239 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2240 if (!pcbddc->switch_static) { 2241 /* compute values into the interior if solved for the partially subassembled Schur complement */ 2242 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 2243 ierr = VecAXPY(mat_ctx->temp_solution_D,-1.0,pcis->vec1_D);CHKERRQ(ierr); 2244 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2245 } 2246 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2247 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2248 /* add p0 solution to final solution */ 2249 ierr = PCBDDCBenignGetOrSetP0(mat_ctx->pc,standard_sol,PETSC_FALSE);CHKERRQ(ierr); 2250 if (pcbddc->ChangeOfBasisMatrix) { 2251 Vec v2; 2252 ierr = VecDuplicate(standard_sol,&v2);CHKERRQ(ierr); 2253 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,standard_sol,v2);CHKERRQ(ierr); 2254 ierr = VecCopy(v2,standard_sol);CHKERRQ(ierr); 2255 ierr = VecDestroy(&v2);CHKERRQ(ierr); 2256 } 2257 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 2258 PetscFunctionReturn(0); 2259 } 2260 2261 #undef __FUNCT__ 2262 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution" 2263 /*@ 2264 PCBDDCMatFETIDPGetSolution - Compute the physical solution using the solution of the FETI-DP linear system 2265 2266 Collective 2267 2268 Input Parameters: 2269 + fetidp_mat - the FETI-DP matrix obtained by a call to PCBDDCCreateFETIDPOperators 2270 - fetidp_flux_sol - the solution of the FETI-DP linear system 2271 2272 Output Parameters: 2273 . standard_sol - the solution defined on the physical domain 2274 2275 Level: developer 2276 2277 Notes: 2278 2279 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetRHS 2280 @*/ 2281 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2282 { 2283 FETIDPMat_ctx mat_ctx; 2284 PetscErrorCode ierr; 2285 2286 PetscFunctionBegin; 2287 PetscValidHeaderSpecific(fetidp_mat,MAT_CLASSID,1); 2288 PetscValidHeaderSpecific(fetidp_flux_sol,VEC_CLASSID,2); 2289 PetscValidHeaderSpecific(standard_sol,VEC_CLASSID,3); 2290 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2291 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 2292 PetscFunctionReturn(0); 2293 } 2294 /* -------------------------------------------------------------------------- */ 2295 2296 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec); 2297 extern PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec); 2298 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat); 2299 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec); 2300 extern PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec); 2301 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC); 2302 2303 #undef __FUNCT__ 2304 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC" 2305 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 2306 { 2307 2308 FETIDPMat_ctx fetidpmat_ctx; 2309 Mat newmat; 2310 FETIDPPC_ctx fetidppc_ctx; 2311 PC newpc; 2312 MPI_Comm comm; 2313 PetscErrorCode ierr; 2314 2315 PetscFunctionBegin; 2316 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 2317 /* FETIDP linear matrix */ 2318 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 2319 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 2320 ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 2321 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 2322 ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr); 2323 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 2324 ierr = MatSetUp(newmat);CHKERRQ(ierr); 2325 /* FETIDP preconditioner */ 2326 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 2327 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 2328 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 2329 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 2330 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 2331 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 2332 ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 2333 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 2334 ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr); 2335 ierr = PCSetUp(newpc);CHKERRQ(ierr); 2336 /* return pointers for objects created */ 2337 *fetidp_mat=newmat; 2338 *fetidp_pc=newpc; 2339 PetscFunctionReturn(0); 2340 } 2341 2342 #undef __FUNCT__ 2343 #define __FUNCT__ "PCBDDCCreateFETIDPOperators" 2344 /*@ 2345 PCBDDCCreateFETIDPOperators - Create FETI-DP operators 2346 2347 Collective 2348 2349 Input Parameters: 2350 . pc - the BDDC preconditioning context (setup should have been called before) 2351 2352 Output Parameters: 2353 + fetidp_mat - shell FETI-DP matrix object 2354 - fetidp_pc - shell Dirichlet preconditioner for FETI-DP matrix 2355 2356 Options Database Keys: 2357 . -fetidp_fullyredundant <false> - use or not a fully redundant set of Lagrange multipliers 2358 2359 Level: developer 2360 2361 Notes: 2362 Currently the only operations provided for FETI-DP matrix are MatMult and MatMultTranspose 2363 2364 .seealso: PCBDDC, PCBDDCMatFETIDPGetRHS, PCBDDCMatFETIDPGetSolution 2365 @*/ 2366 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 2367 { 2368 PetscErrorCode ierr; 2369 2370 PetscFunctionBegin; 2371 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2372 if (pc->setupcalled) { 2373 ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,Mat*,PC*),(pc,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 2374 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n"); 2375 PetscFunctionReturn(0); 2376 } 2377 /* -------------------------------------------------------------------------- */ 2378 /*MC 2379 PCBDDC - Balancing Domain Decomposition by Constraints. 2380 2381 An implementation of the BDDC preconditioner based on 2382 2383 .vb 2384 [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 2385 [2] A. Klawonn and O. B. Widlund. "Dual-Primal FETI Methods for Linear Elasticity", http://cs.nyu.edu/csweb/Research/TechReports/TR2004-855/TR2004-855.pdf 2386 [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977 2387 [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 2388 .ve 2389 2390 The matrix to be preconditioned (Pmat) must be of type MATIS. 2391 2392 Currently works with MATIS matrices with local matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers. 2393 2394 It also works with unsymmetric and indefinite problems. 2395 2396 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. 2397 2398 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). 2399 2400 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() 2401 Additional information on dofs can be provided by using PCBDDCSetDofsSplitting(), PCBDDCSetDirichletBoundaries(), PCBDDCSetNeumannBoundaries(), and PCBDDCSetPrimalVerticesIS() and their local counterparts. 2402 2403 Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). Non-singular modes are retained via SVD. 2404 2405 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. 2406 User defined change of basis can be passed to PCBDDC by using PCBDDCSetChangeOfBasisMat() 2407 2408 The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using a MatPartitioning object. 2409 2410 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. 2411 2412 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. 2413 Deluxe scaling is not supported yet for FETI-DP. 2414 2415 Options Database Keys (some of them, run with -h for a complete list): 2416 2417 . -pc_bddc_use_vertices <true> - use or not vertices in primal space 2418 . -pc_bddc_use_edges <true> - use or not edges in primal space 2419 . -pc_bddc_use_faces <false> - use or not faces in primal space 2420 . -pc_bddc_symmetric <true> - symmetric computation of primal basis functions. Specify false for unsymmetric problems 2421 . -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only) 2422 . -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested 2423 . -pc_bddc_switch_static <false> - switches from M_2 (default) to M_3 operator (see reference article [1]) 2424 . -pc_bddc_levels <0> - maximum number of levels for multilevel 2425 . -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) 2426 . -pc_bddc_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) 2427 . -pc_bddc_use_deluxe_scaling <false> - use deluxe scaling 2428 . -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) 2429 . -pc_bddc_adaptive_threshold <0.0> - when a value greater than one is specified, adaptive selection of constraints is performed on edges and faces (requires deluxe scaling and MUMPS or MKL_PARDISO installed) 2430 - -pc_bddc_check_level <0> - set verbosity level of debugging output 2431 2432 Options for Dirichlet, Neumann or coarse solver can be set with 2433 .vb 2434 -pc_bddc_dirichlet_ 2435 -pc_bddc_neumann_ 2436 -pc_bddc_coarse_ 2437 .ve 2438 e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. PCBDDC uses by default KPSPREONLY and PCLU. 2439 2440 When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified as 2441 .vb 2442 -pc_bddc_dirichlet_lN_ 2443 -pc_bddc_neumann_lN_ 2444 -pc_bddc_coarse_lN_ 2445 .ve 2446 Note that level number ranges from the finest (0) to the coarsest (N). 2447 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. 2448 .vb 2449 -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3 2450 .ve 2451 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 2452 2453 Level: intermediate 2454 2455 Developer notes: 2456 2457 Contributed by Stefano Zampini 2458 2459 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 2460 M*/ 2461 2462 #undef __FUNCT__ 2463 #define __FUNCT__ "PCCreate_BDDC" 2464 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 2465 { 2466 PetscErrorCode ierr; 2467 PC_BDDC *pcbddc; 2468 2469 PetscFunctionBegin; 2470 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 2471 ierr = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr); 2472 pc->data = (void*)pcbddc; 2473 2474 /* create PCIS data structure */ 2475 ierr = PCISCreate(pc);CHKERRQ(ierr); 2476 2477 /* BDDC customization */ 2478 pcbddc->use_local_adj = PETSC_TRUE; 2479 pcbddc->use_vertices = PETSC_TRUE; 2480 pcbddc->use_edges = PETSC_TRUE; 2481 pcbddc->use_faces = PETSC_FALSE; 2482 pcbddc->use_change_of_basis = PETSC_FALSE; 2483 pcbddc->use_change_on_faces = PETSC_FALSE; 2484 pcbddc->switch_static = PETSC_FALSE; 2485 pcbddc->use_nnsp_true = PETSC_FALSE; 2486 pcbddc->use_qr_single = PETSC_FALSE; 2487 pcbddc->symmetric_primal = PETSC_TRUE; 2488 pcbddc->benign_saddle_point = PETSC_FALSE; 2489 pcbddc->benign_have_null = PETSC_FALSE; 2490 pcbddc->vertex_size = 1; 2491 pcbddc->dbg_flag = 0; 2492 /* private */ 2493 pcbddc->local_primal_size = 0; 2494 pcbddc->local_primal_size_cc = 0; 2495 pcbddc->local_primal_ref_node = 0; 2496 pcbddc->local_primal_ref_mult = 0; 2497 pcbddc->n_vertices = 0; 2498 pcbddc->primal_indices_local_idxs = 0; 2499 pcbddc->recompute_topography = PETSC_FALSE; 2500 pcbddc->coarse_size = -1; 2501 pcbddc->new_primal_space = PETSC_FALSE; 2502 pcbddc->new_primal_space_local = PETSC_FALSE; 2503 pcbddc->global_primal_indices = 0; 2504 pcbddc->onearnullspace = 0; 2505 pcbddc->onearnullvecs_state = 0; 2506 pcbddc->user_primal_vertices = 0; 2507 pcbddc->user_primal_vertices_local = 0; 2508 pcbddc->temp_solution = 0; 2509 pcbddc->original_rhs = 0; 2510 pcbddc->local_mat = 0; 2511 pcbddc->ChangeOfBasisMatrix = 0; 2512 pcbddc->user_ChangeOfBasisMatrix = 0; 2513 pcbddc->coarse_vec = 0; 2514 pcbddc->coarse_ksp = 0; 2515 pcbddc->coarse_phi_B = 0; 2516 pcbddc->coarse_phi_D = 0; 2517 pcbddc->coarse_psi_B = 0; 2518 pcbddc->coarse_psi_D = 0; 2519 pcbddc->vec1_P = 0; 2520 pcbddc->vec1_R = 0; 2521 pcbddc->vec2_R = 0; 2522 pcbddc->local_auxmat1 = 0; 2523 pcbddc->local_auxmat2 = 0; 2524 pcbddc->R_to_B = 0; 2525 pcbddc->R_to_D = 0; 2526 pcbddc->ksp_D = 0; 2527 pcbddc->ksp_R = 0; 2528 pcbddc->NeumannBoundaries = 0; 2529 pcbddc->NeumannBoundariesLocal = 0; 2530 pcbddc->DirichletBoundaries = 0; 2531 pcbddc->DirichletBoundariesLocal = 0; 2532 pcbddc->user_provided_isfordofs = PETSC_FALSE; 2533 pcbddc->n_ISForDofs = 0; 2534 pcbddc->n_ISForDofsLocal = 0; 2535 pcbddc->ISForDofs = 0; 2536 pcbddc->ISForDofsLocal = 0; 2537 pcbddc->ConstraintMatrix = 0; 2538 pcbddc->use_exact_dirichlet_trick = PETSC_TRUE; 2539 pcbddc->coarse_loc_to_glob = 0; 2540 pcbddc->coarsening_ratio = 8; 2541 pcbddc->coarse_adj_red = 0; 2542 pcbddc->current_level = 0; 2543 pcbddc->max_levels = 0; 2544 pcbddc->use_coarse_estimates = PETSC_FALSE; 2545 pcbddc->coarse_eqs_per_proc = 1; 2546 pcbddc->coarse_subassembling = 0; 2547 pcbddc->detect_disconnected = PETSC_FALSE; 2548 pcbddc->n_local_subs = 0; 2549 pcbddc->local_subs = NULL; 2550 2551 /* benign subspace trick */ 2552 pcbddc->benign_change = 0; 2553 pcbddc->benign_compute_correction = PETSC_TRUE; 2554 pcbddc->benign_vec = 0; 2555 pcbddc->benign_original_mat = 0; 2556 pcbddc->benign_sf = 0; 2557 pcbddc->benign_B0 = 0; 2558 pcbddc->benign_n = 0; 2559 pcbddc->benign_p0 = NULL; 2560 pcbddc->benign_p0_lidx = NULL; 2561 pcbddc->benign_p0_gidx = NULL; 2562 pcbddc->benign_null = PETSC_FALSE; 2563 2564 /* create local graph structure */ 2565 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 2566 2567 /* scaling */ 2568 pcbddc->work_scaling = 0; 2569 pcbddc->use_deluxe_scaling = PETSC_FALSE; 2570 2571 /* sub schurs options */ 2572 pcbddc->sub_schurs_rebuild = PETSC_FALSE; 2573 pcbddc->sub_schurs_layers = -1; 2574 pcbddc->sub_schurs_use_useradj = PETSC_FALSE; 2575 2576 pcbddc->computed_rowadj = PETSC_FALSE; 2577 2578 /* adaptivity */ 2579 pcbddc->adaptive_threshold = 0.0; 2580 pcbddc->adaptive_nmax = 0; 2581 pcbddc->adaptive_nmin = 0; 2582 2583 /* function pointers */ 2584 pc->ops->apply = PCApply_BDDC; 2585 pc->ops->applytranspose = PCApplyTranspose_BDDC; 2586 pc->ops->setup = PCSetUp_BDDC; 2587 pc->ops->destroy = PCDestroy_BDDC; 2588 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 2589 pc->ops->view = PCView_BDDC; 2590 pc->ops->applyrichardson = 0; 2591 pc->ops->applysymmetricleft = 0; 2592 pc->ops->applysymmetricright = 0; 2593 pc->ops->presolve = PCPreSolve_BDDC; 2594 pc->ops->postsolve = PCPostSolve_BDDC; 2595 2596 /* composing function */ 2597 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDivergenceMat_C",PCBDDCSetDivergenceMat_BDDC);CHKERRQ(ierr); 2598 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);CHKERRQ(ierr); 2599 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 2600 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",PCBDDCSetPrimalVerticesIS_BDDC);CHKERRQ(ierr); 2601 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 2602 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr); 2603 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr); 2604 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr); 2605 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2606 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2607 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2608 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2609 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2610 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2611 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2612 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2613 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 2614 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr); 2615 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 2616 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 2617 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 2618 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 2619 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_BDDC);CHKERRQ(ierr); 2620 PetscFunctionReturn(0); 2621 } 2622 2623