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