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