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