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 (change explicit %d)\n",pcbddc->benign_saddle_point,pcbddc->benign_change_explicit);CHKERRQ(ierr); 151 ierr = PetscViewerASCIIPrintf(viewer," BDDC: Benign subspace trick is active: %d (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 if (pcbddc->benign_apply_coarse_only || pcbddc->switch_static) { 1159 ierr = PCBDDCSetUseExactDirichlet(pc,PETSC_FALSE);CHKERRQ(ierr); 1160 } 1161 1162 /* Creates parallel work vectors used in presolve */ 1163 if (!pcbddc->original_rhs) { 1164 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 1165 } 1166 if (!pcbddc->temp_solution) { 1167 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);CHKERRQ(ierr); 1168 } 1169 1170 pcbddc->temp_solution_used = PETSC_FALSE; 1171 if (x) { 1172 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr); 1173 used_vec = x; 1174 } else { /* it can only happen when calling PCBDDCMatFETIDPGetRHS */ 1175 ierr = PetscObjectReference((PetscObject)pcbddc->temp_solution);CHKERRQ(ierr); 1176 used_vec = pcbddc->temp_solution; 1177 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 1178 pcbddc->temp_solution_used = PETSC_TRUE; 1179 } 1180 1181 /* hack into ksp data structure since PCPreSolve comes earlier than setting to zero the guess in src/ksp/ksp/interface/itfunc.c */ 1182 if (ksp) { 1183 /* store the flag for the initial guess since it will be restored back during PCPostSolve_BDDC */ 1184 ierr = KSPGetInitialGuessNonzero(ksp,&pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 1185 if (!pcbddc->ksp_guess_nonzero) { 1186 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 1187 } 1188 } 1189 1190 pcbddc->rhs_change = PETSC_FALSE; 1191 /* Take into account zeroed rows -> change rhs and store solution removed */ 1192 if (rhs && pcbddc->eliminate_dirdofs) { 1193 IS dirIS = NULL; 1194 1195 /* DirichletBoundariesLocal may not be consistent among neighbours; gets a dirichlet dofs IS from graph (may be cached) */ 1196 ierr = PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph,&dirIS);CHKERRQ(ierr); 1197 if (dirIS) { 1198 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 1199 PetscInt dirsize,i,*is_indices; 1200 PetscScalar *array_x; 1201 const PetscScalar *array_diagonal; 1202 1203 ierr = MatGetDiagonal(pc->pmat,pcis->vec1_global);CHKERRQ(ierr); 1204 ierr = VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);CHKERRQ(ierr); 1205 ierr = VecScatterBegin(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1206 ierr = VecScatterEnd(matis->rctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1207 ierr = VecScatterBegin(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1208 ierr = VecScatterEnd(matis->rctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1209 ierr = ISGetLocalSize(dirIS,&dirsize);CHKERRQ(ierr); 1210 ierr = VecGetArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 1211 ierr = VecGetArrayRead(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 1212 ierr = ISGetIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1213 for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]]; 1214 ierr = ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1215 ierr = VecRestoreArrayRead(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 1216 ierr = VecRestoreArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 1217 ierr = VecScatterBegin(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1218 ierr = VecScatterEnd(matis->rctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1219 pcbddc->rhs_change = PETSC_TRUE; 1220 ierr = ISDestroy(&dirIS);CHKERRQ(ierr); 1221 } 1222 } 1223 1224 /* remove the computed solution or the initial guess from the rhs */ 1225 if (pcbddc->rhs_change || (ksp && pcbddc->ksp_guess_nonzero) ) { 1226 /* save the original rhs */ 1227 if (save_rhs) { 1228 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1229 save_rhs = PETSC_FALSE; 1230 } 1231 pcbddc->rhs_change = PETSC_TRUE; 1232 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 1233 ierr = MatMultAdd(pc->mat,used_vec,pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1234 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 1235 ierr = VecCopy(used_vec,pcbddc->temp_solution);CHKERRQ(ierr); 1236 pcbddc->temp_solution_used = PETSC_TRUE; 1237 if (ksp) { 1238 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_FALSE);CHKERRQ(ierr); 1239 } 1240 } 1241 ierr = VecDestroy(&used_vec);CHKERRQ(ierr); 1242 1243 /* compute initial vector in benign space if needed 1244 and remove non-benign solution from the rhs */ 1245 benign_correction_computed = PETSC_FALSE; 1246 if (rhs && pcbddc->benign_compute_correction && pcbddc->benign_have_null) { 1247 /* compute u^*_h using ideas similar to those in Xuemin Tu's PhD thesis (see Section 4.8.1) 1248 Recursively apply BDDC in the multilevel case */ 1249 if (!pcbddc->benign_vec) { 1250 ierr = VecDuplicate(rhs,&pcbddc->benign_vec);CHKERRQ(ierr); 1251 } 1252 pcbddc->benign_apply_coarse_only = PETSC_TRUE; 1253 if (!pcbddc->benign_skip_correction) { 1254 ierr = PCApply_BDDC(pc,rhs,pcbddc->benign_vec);CHKERRQ(ierr); 1255 benign_correction_computed = PETSC_TRUE; 1256 if (pcbddc->temp_solution_used) { 1257 ierr = VecAXPY(pcbddc->temp_solution,1.0,pcbddc->benign_vec);CHKERRQ(ierr); 1258 } 1259 ierr = VecScale(pcbddc->benign_vec,-1.0);CHKERRQ(ierr); 1260 /* store the original rhs if not done earlier */ 1261 if (save_rhs) { 1262 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1263 save_rhs = PETSC_FALSE; 1264 } 1265 if (pcbddc->rhs_change) { 1266 ierr = MatMultAdd(pc->mat,pcbddc->benign_vec,rhs,rhs);CHKERRQ(ierr); 1267 } else { 1268 ierr = MatMultAdd(pc->mat,pcbddc->benign_vec,pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1269 } 1270 pcbddc->rhs_change = PETSC_TRUE; 1271 } 1272 pcbddc->benign_apply_coarse_only = PETSC_FALSE; 1273 } 1274 #if 0 1275 if (pcbddc->dbg_flag && benign_correction_computed) { 1276 Vec v; 1277 ierr = VecDuplicate(pcis->vec1_global,&v);CHKERRQ(ierr); 1278 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,rhs,v);CHKERRQ(ierr); 1279 ierr = PCBDDCBenignGetOrSetP0(pc,v,PETSC_TRUE);CHKERRQ(ierr); 1280 PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)),"LEVEL %d: IS CORRECTION BENIGN?\n",pcbddc->current_level); 1281 PetscScalarView(pcbddc->benign_n,pcbddc->benign_p0,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc))); 1282 ierr = VecDestroy(&v);CHKERRQ(ierr); 1283 } 1284 #endif 1285 1286 /* set initial guess if using PCG */ 1287 if (x && pcbddc->use_exact_dirichlet_trick) { 1288 ierr = VecSet(x,0.0);CHKERRQ(ierr); 1289 if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) { 1290 if (benign_correction_computed) { /* we have already saved the changed rhs */ 1291 ierr = VecLockPop(pcis->vec1_global);CHKERRQ(ierr); 1292 } else { 1293 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,rhs,pcis->vec1_global);CHKERRQ(ierr); 1294 } 1295 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_global,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1296 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_global,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1297 } else { 1298 ierr = VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1299 ierr = VecScatterEnd(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1300 } 1301 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1302 if (pcbddc->ChangeOfBasisMatrix && pcbddc->change_interior) { 1303 ierr = VecSet(pcis->vec1_global,0.);CHKERRQ(ierr); 1304 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1305 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,pcis->vec1_global,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1306 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcis->vec1_global,x);CHKERRQ(ierr); 1307 } else { 1308 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1309 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1310 } 1311 if (ksp) { 1312 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 1313 } 1314 } 1315 PetscFunctionReturn(0); 1316 } 1317 1318 /* -------------------------------------------------------------------------- */ 1319 #undef __FUNCT__ 1320 #define __FUNCT__ "PCPostSolve_BDDC" 1321 /* -------------------------------------------------------------------------- */ 1322 /* 1323 PCPostSolve_BDDC - Changes the computed solution if a transformation of basis 1324 approach has been selected. Also, restores rhs to its original state. 1325 1326 Input Parameter: 1327 + pc - the preconditioner contex 1328 1329 Application Interface Routine: PCPostSolve() 1330 1331 Notes: 1332 The interface routine PCPostSolve() is not usually called directly by 1333 the user, but instead is called by KSPSolve(). 1334 */ 1335 static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 1336 { 1337 PetscErrorCode ierr; 1338 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1339 1340 PetscFunctionBegin; 1341 /* add solution removed in presolve */ 1342 if (x && pcbddc->rhs_change) { 1343 if (pcbddc->temp_solution_used) { 1344 ierr = VecAXPY(x,1.0,pcbddc->temp_solution);CHKERRQ(ierr); 1345 } else if (pcbddc->benign_compute_correction) { 1346 ierr = VecAXPY(x,-1.0,pcbddc->benign_vec);CHKERRQ(ierr); 1347 } 1348 } 1349 pcbddc->temp_solution_used = PETSC_FALSE; 1350 1351 /* restore rhs to its original state */ 1352 if (rhs && pcbddc->rhs_change) { 1353 ierr = VecSwap(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 1354 } 1355 pcbddc->rhs_change = PETSC_FALSE; 1356 /* restore ksp guess state */ 1357 if (ksp) { 1358 ierr = KSPSetInitialGuessNonzero(ksp,pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 1359 } 1360 PetscFunctionReturn(0); 1361 } 1362 /* -------------------------------------------------------------------------- */ 1363 #undef __FUNCT__ 1364 #define __FUNCT__ "PCSetUp_BDDC" 1365 /* -------------------------------------------------------------------------- */ 1366 /* 1367 PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner 1368 by setting data structures and options. 1369 1370 Input Parameter: 1371 + pc - the preconditioner context 1372 1373 Application Interface Routine: PCSetUp() 1374 1375 Notes: 1376 The interface routine PCSetUp() is not usually called directly by 1377 the user, but instead is called by PCApply() if necessary. 1378 */ 1379 PetscErrorCode PCSetUp_BDDC(PC pc) 1380 { 1381 PetscErrorCode ierr; 1382 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 1383 Mat_IS* matis; 1384 MatNullSpace nearnullspace; 1385 IS zerodiag = NULL; 1386 PetscInt nrows,ncols; 1387 PetscBool computetopography,computesolvers,computesubschurs; 1388 PetscBool computeconstraintsmatrix; 1389 PetscBool new_nearnullspace_provided,ismatis; 1390 1391 PetscFunctionBegin; 1392 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&ismatis);CHKERRQ(ierr); 1393 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PCBDDC preconditioner requires matrix of type MATIS"); 1394 ierr = MatGetSize(pc->pmat,&nrows,&ncols);CHKERRQ(ierr); 1395 if (nrows != ncols) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCBDDC preconditioner requires a square preconditioning matrix"); 1396 matis = (Mat_IS*)pc->pmat->data; 1397 /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */ 1398 /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup 1399 Also, BDDC directly build the Dirichlet problem */ 1400 /* split work */ 1401 if (pc->setupcalled) { 1402 if (pc->flag == SAME_NONZERO_PATTERN) { 1403 computetopography = PETSC_FALSE; 1404 computesolvers = PETSC_TRUE; 1405 } else { /* DIFFERENT_NONZERO_PATTERN */ 1406 computetopography = PETSC_TRUE; 1407 computesolvers = PETSC_TRUE; 1408 } 1409 } else { 1410 computetopography = PETSC_TRUE; 1411 computesolvers = PETSC_TRUE; 1412 } 1413 if (pcbddc->recompute_topography) { 1414 computetopography = PETSC_TRUE; 1415 } 1416 pcbddc->recompute_topography = computetopography; 1417 computeconstraintsmatrix = PETSC_FALSE; 1418 1419 /* check parameters' compatibility */ 1420 if (!pcbddc->use_deluxe_scaling) pcbddc->deluxe_zerorows = PETSC_FALSE; 1421 pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_threshold > 0.0); 1422 pcbddc->adaptive_userdefined = (PetscBool)(pcbddc->adaptive_selection && pcbddc->adaptive_userdefined); 1423 if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE; 1424 1425 computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling); 1426 if (pcbddc->switch_static) { 1427 PetscBool ismatis; 1428 ierr = PetscObjectTypeCompare((PetscObject)pc->mat,MATIS,&ismatis);CHKERRQ(ierr); 1429 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"When the static switch is one, the iteration matrix should be of type MATIS"); 1430 } 1431 1432 /* Get stdout for dbg */ 1433 if (pcbddc->dbg_flag) { 1434 if (!pcbddc->dbg_viewer) { 1435 pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)); 1436 ierr = PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer);CHKERRQ(ierr); 1437 } 1438 ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1439 } 1440 1441 if (pcbddc->user_ChangeOfBasisMatrix) { 1442 /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */ 1443 pcbddc->use_change_of_basis = PETSC_FALSE; 1444 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 1445 } else { 1446 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1447 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1448 pcbddc->local_mat = matis->A; 1449 } 1450 1451 /* detect local disconnected subdomains if requested and not done before */ 1452 if (pcbddc->detect_disconnected && !pcbddc->n_local_subs) { 1453 ierr = MatDetectDisconnectedComponents(pcbddc->local_mat,PETSC_FALSE,&pcbddc->n_local_subs,&pcbddc->local_subs);CHKERRQ(ierr); 1454 } 1455 1456 /* compute topology info in local ordering */ 1457 if (pcbddc->recompute_topography) { 1458 ierr = PCBDDCComputeLocalTopologyInfo(pc);CHKERRQ(ierr); 1459 } 1460 1461 /* 1462 Compute change of basis on local pressures (aka zerodiag dofs) 1463 This should come earlier then PCISSetUp for extracting the correct subdomain matrices 1464 */ 1465 ierr = PCBDDCBenignShellMat(pc,PETSC_TRUE);CHKERRQ(ierr); 1466 if (pcbddc->benign_saddle_point) { 1467 PC_IS* pcis = (PC_IS*)pc->data; 1468 1469 if (pcbddc->user_ChangeOfBasisMatrix || pcbddc->use_change_of_basis || !computesubschurs) pcbddc->benign_change_explicit = PETSC_TRUE; 1470 /* detect local saddle point and change the basis in pcbddc->local_mat */ 1471 ierr = PCBDDCBenignDetectSaddlePoint(pc,&zerodiag);CHKERRQ(ierr); 1472 /* pop B0 mat from local mat */ 1473 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1474 /* give pcis a hint to not reuse submatrices during PCISCreate */ 1475 if (pc->flag == SAME_NONZERO_PATTERN && pcis->reusesubmatrices == PETSC_TRUE) { 1476 if (pcbddc->benign_n && (pcbddc->benign_change_explicit || pcbddc->dbg_flag)) { 1477 pcis->reusesubmatrices = PETSC_FALSE; 1478 } else { 1479 pcis->reusesubmatrices = PETSC_TRUE; 1480 } 1481 } else { 1482 pcis->reusesubmatrices = PETSC_FALSE; 1483 } 1484 } 1485 1486 /* propagate relevant information */ 1487 #if !defined(PETSC_USE_COMPLEX) /* workaround for reals */ 1488 if (matis->A->symmetric_set) { 1489 ierr = MatSetOption(pcbddc->local_mat,MAT_HERMITIAN,matis->A->symmetric);CHKERRQ(ierr); 1490 } 1491 #endif 1492 if (matis->A->symmetric_set) { 1493 ierr = MatSetOption(pcbddc->local_mat,MAT_SYMMETRIC,matis->A->symmetric);CHKERRQ(ierr); 1494 } 1495 if (matis->A->spd_set) { 1496 ierr = MatSetOption(pcbddc->local_mat,MAT_SPD,matis->A->spd);CHKERRQ(ierr); 1497 } 1498 1499 /* Set up all the "iterative substructuring" common block without computing solvers */ 1500 { 1501 Mat temp_mat; 1502 1503 temp_mat = matis->A; 1504 matis->A = pcbddc->local_mat; 1505 ierr = PCISSetUp(pc,PETSC_FALSE);CHKERRQ(ierr); 1506 pcbddc->local_mat = matis->A; 1507 matis->A = temp_mat; 1508 } 1509 1510 /* Analyze interface */ 1511 if (computetopography) { 1512 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 1513 computeconstraintsmatrix = PETSC_TRUE; 1514 if (pcbddc->adaptive_selection && !pcbddc->use_deluxe_scaling && !pcbddc->mat_graph->twodim) { 1515 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot compute the adaptive primal space for a problem with 3D edges without deluxe scaling"); 1516 } 1517 } 1518 1519 /* check existence of a divergence free extension, i.e. 1520 b(v_I,p_0) = 0 for all v_I (raise error if not). 1521 Also, check that PCBDDCBenignGetOrSetP0 works */ 1522 #if defined(PETSC_USE_DEBUG) 1523 if (pcbddc->benign_saddle_point) { 1524 ierr = PCBDDCBenignCheck(pc,zerodiag);CHKERRQ(ierr); 1525 } 1526 #endif 1527 ierr = ISDestroy(&zerodiag);CHKERRQ(ierr); 1528 1529 /* Setup local dirichlet solver ksp_D and sub_schurs solvers */ 1530 if (computesolvers) { 1531 PCBDDCSubSchurs sub_schurs=pcbddc->sub_schurs; 1532 1533 if (computesubschurs && computetopography) { 1534 ierr = PCBDDCInitSubSchurs(pc);CHKERRQ(ierr); 1535 } 1536 /* SetUp Scaling operator (scaling matrices could be needed in SubSchursSetUp)*/ 1537 if (!pcbddc->use_deluxe_scaling) { 1538 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1539 } 1540 if (sub_schurs->schur_explicit) { 1541 if (computesubschurs) { 1542 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1543 } 1544 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1545 } else { 1546 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1547 if (computesubschurs) { 1548 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1549 } 1550 } 1551 if (pcbddc->adaptive_selection) { 1552 ierr = PCBDDCAdaptiveSelection(pc);CHKERRQ(ierr); 1553 computeconstraintsmatrix = PETSC_TRUE; 1554 } 1555 } 1556 1557 /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */ 1558 new_nearnullspace_provided = PETSC_FALSE; 1559 ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr); 1560 if (pcbddc->onearnullspace) { /* already used nearnullspace */ 1561 if (!nearnullspace) { /* near null space attached to mat has been destroyed */ 1562 new_nearnullspace_provided = PETSC_TRUE; 1563 } else { 1564 /* determine if the two nullspaces are different (should be lightweight) */ 1565 if (nearnullspace != pcbddc->onearnullspace) { 1566 new_nearnullspace_provided = PETSC_TRUE; 1567 } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */ 1568 PetscInt i; 1569 const Vec *nearnullvecs; 1570 PetscObjectState state; 1571 PetscInt nnsp_size; 1572 ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 1573 for (i=0;i<nnsp_size;i++) { 1574 ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr); 1575 if (pcbddc->onearnullvecs_state[i] != state) { 1576 new_nearnullspace_provided = PETSC_TRUE; 1577 break; 1578 } 1579 } 1580 } 1581 } 1582 } else { 1583 if (!nearnullspace) { /* both nearnullspaces are null */ 1584 new_nearnullspace_provided = PETSC_FALSE; 1585 } else { /* nearnullspace attached later */ 1586 new_nearnullspace_provided = PETSC_TRUE; 1587 } 1588 } 1589 1590 /* Setup constraints and related work vectors */ 1591 /* reset primal space flags */ 1592 pcbddc->new_primal_space = PETSC_FALSE; 1593 pcbddc->new_primal_space_local = PETSC_FALSE; 1594 if (computeconstraintsmatrix || new_nearnullspace_provided) { 1595 /* It also sets the primal space flags */ 1596 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 1597 /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */ 1598 ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr); 1599 } 1600 1601 if (computesolvers || pcbddc->new_primal_space) { 1602 if (pcbddc->use_change_of_basis) { 1603 PC_IS *pcis = (PC_IS*)(pc->data); 1604 1605 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 1606 if (pcbddc->benign_change) { 1607 ierr = MatDestroy(&pcbddc->benign_B0);CHKERRQ(ierr); 1608 /* pop B0 from pcbddc->local_mat */ 1609 ierr = PCBDDCBenignPopOrPushB0(pc,PETSC_TRUE);CHKERRQ(ierr); 1610 } 1611 /* get submatrices */ 1612 ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr); 1613 ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr); 1614 ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr); 1615 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);CHKERRQ(ierr); 1616 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);CHKERRQ(ierr); 1617 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);CHKERRQ(ierr); 1618 /* set flag in pcis to not reuse submatrices during PCISCreate */ 1619 pcis->reusesubmatrices = PETSC_FALSE; 1620 } else if (!pcbddc->user_ChangeOfBasisMatrix && !pcbddc->benign_change) { 1621 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1622 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1623 pcbddc->local_mat = matis->A; 1624 } 1625 /* SetUp coarse and local Neumann solvers */ 1626 ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr); 1627 /* SetUp Scaling operator */ 1628 if (pcbddc->use_deluxe_scaling) { 1629 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1630 } 1631 } 1632 /* mark topography as done */ 1633 pcbddc->recompute_topography = PETSC_FALSE; 1634 1635 /* wrap pcis->A_IB and pcis->A_BI if we did not change explicitly the variables on the pressures */ 1636 ierr = PCBDDCBenignShellMat(pc,PETSC_FALSE);CHKERRQ(ierr); 1637 1638 if (pcbddc->dbg_flag) { 1639 ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1640 } 1641 PetscFunctionReturn(0); 1642 } 1643 1644 /* -------------------------------------------------------------------------- */ 1645 /* 1646 PCApply_BDDC - Applies the BDDC operator to a vector. 1647 1648 Input Parameters: 1649 + pc - the preconditioner context 1650 - r - input vector (global) 1651 1652 Output Parameter: 1653 . z - output vector (global) 1654 1655 Application Interface Routine: PCApply() 1656 */ 1657 #undef __FUNCT__ 1658 #define __FUNCT__ "PCApply_BDDC" 1659 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 1660 { 1661 PC_IS *pcis = (PC_IS*)(pc->data); 1662 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1663 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1664 PetscErrorCode ierr; 1665 const PetscScalar one = 1.0; 1666 const PetscScalar m_one = -1.0; 1667 const PetscScalar zero = 0.0; 1668 1669 /* This code is similar to that provided in nn.c for PCNN 1670 NN interface preconditioner changed to BDDC 1671 Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */ 1672 1673 PetscFunctionBegin; 1674 if (pcbddc->ChangeOfBasisMatrix) { 1675 Vec swap; 1676 1677 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 1678 swap = pcbddc->work_change; 1679 pcbddc->work_change = r; 1680 r = swap; 1681 /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */ 1682 if (pcbddc->benign_apply_coarse_only && pcbddc->use_exact_dirichlet_trick && pcbddc->change_interior) { 1683 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 1684 ierr = VecLockPush(pcis->vec1_global);CHKERRQ(ierr); 1685 } 1686 } 1687 if (pcbddc->benign_have_null) { /* get p0 from r */ 1688 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 1689 } 1690 if (!pcbddc->use_exact_dirichlet_trick && !pcbddc->benign_apply_coarse_only) { 1691 ierr = VecCopy(r,z);CHKERRQ(ierr); 1692 /* First Dirichlet solve */ 1693 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1694 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1695 /* 1696 Assembling right hand side for BDDC operator 1697 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1698 - pcis->vec1_B the interface part of the global vector z 1699 */ 1700 if (n_D) { 1701 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1702 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1703 if (pcbddc->switch_static) { 1704 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1705 1706 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 1707 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1708 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1709 if (!pcbddc->switch_static_change) { 1710 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1711 } else { 1712 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1713 ierr = MatMult(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1714 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1715 } 1716 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1717 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1718 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1719 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1720 } else { 1721 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1722 } 1723 } else { 1724 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1725 } 1726 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1727 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1728 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1729 } else { 1730 if (!pcbddc->benign_apply_coarse_only) { 1731 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1732 } 1733 } 1734 1735 /* Apply interface preconditioner 1736 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1737 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr); 1738 1739 /* Apply transpose of partition of unity operator */ 1740 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1741 1742 /* Second Dirichlet solve and assembling of output */ 1743 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1744 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1745 if (n_B) { 1746 if (pcbddc->switch_static) { 1747 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1748 1749 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1750 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1751 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1752 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1753 if (!pcbddc->switch_static_change) { 1754 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1755 } else { 1756 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1757 ierr = MatMult(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1758 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1759 } 1760 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1761 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1762 } else { 1763 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1764 } 1765 } else if (pcbddc->switch_static) { /* n_B is zero */ 1766 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1767 1768 if (!pcbddc->switch_static_change) { 1769 ierr = MatMult(matis->A,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 1770 } else { 1771 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 1772 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1773 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 1774 } 1775 } 1776 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1777 1778 if (!pcbddc->use_exact_dirichlet_trick && !pcbddc->benign_apply_coarse_only) { 1779 if (pcbddc->switch_static) { 1780 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 1781 } else { 1782 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 1783 } 1784 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1785 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1786 } else { 1787 if (pcbddc->switch_static) { 1788 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 1789 } else { 1790 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1791 } 1792 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1793 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1794 } 1795 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 1796 if (pcbddc->benign_apply_coarse_only) { 1797 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 1798 } 1799 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 1800 } 1801 1802 if (pcbddc->ChangeOfBasisMatrix) { 1803 Vec swap; 1804 1805 swap = r; 1806 r = pcbddc->work_change; 1807 pcbddc->work_change = swap; 1808 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 1809 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 1810 } 1811 PetscFunctionReturn(0); 1812 } 1813 1814 /* -------------------------------------------------------------------------- */ 1815 /* 1816 PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector. 1817 1818 Input Parameters: 1819 + pc - the preconditioner context 1820 - r - input vector (global) 1821 1822 Output Parameter: 1823 . z - output vector (global) 1824 1825 Application Interface Routine: PCApplyTranspose() 1826 */ 1827 #undef __FUNCT__ 1828 #define __FUNCT__ "PCApplyTranspose_BDDC" 1829 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z) 1830 { 1831 PC_IS *pcis = (PC_IS*)(pc->data); 1832 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1833 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1834 PetscErrorCode ierr; 1835 const PetscScalar one = 1.0; 1836 const PetscScalar m_one = -1.0; 1837 const PetscScalar zero = 0.0; 1838 1839 PetscFunctionBegin; 1840 if (pcbddc->ChangeOfBasisMatrix) { 1841 Vec swap; 1842 1843 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,r,pcbddc->work_change);CHKERRQ(ierr); 1844 swap = pcbddc->work_change; 1845 pcbddc->work_change = r; 1846 r = swap; 1847 /* save rhs so that we don't need to apply the change of basis for the exact dirichlet trick in PreSolve */ 1848 if (pcbddc->benign_apply_coarse_only && pcbddc->use_exact_dirichlet_trick && pcbddc->change_interior) { 1849 ierr = VecCopy(r,pcis->vec1_global);CHKERRQ(ierr); 1850 ierr = VecLockPush(pcis->vec1_global);CHKERRQ(ierr); 1851 } 1852 } 1853 if (pcbddc->benign_have_null) { /* get p0 from r */ 1854 ierr = PCBDDCBenignGetOrSetP0(pc,r,PETSC_TRUE);CHKERRQ(ierr); 1855 } 1856 if (!pcbddc->use_exact_dirichlet_trick && !pcbddc->benign_apply_coarse_only) { 1857 ierr = VecCopy(r,z);CHKERRQ(ierr); 1858 /* First Dirichlet solve */ 1859 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1860 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1861 /* 1862 Assembling right hand side for BDDC operator 1863 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1864 - pcis->vec1_B the interface part of the global vector z 1865 */ 1866 if (n_D) { 1867 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1868 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1869 if (pcbddc->switch_static) { 1870 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1871 1872 ierr = VecSet(pcis->vec1_N,0.);CHKERRQ(ierr); 1873 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1874 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1875 if (!pcbddc->switch_static_change) { 1876 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1877 } else { 1878 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1879 ierr = MatMultTranspose(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1880 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1881 } 1882 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1883 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec1_D,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1884 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1885 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec2_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1886 } else { 1887 ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1888 } 1889 } else { 1890 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1891 } 1892 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1893 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1894 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1895 } else { 1896 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1897 } 1898 1899 /* Apply interface preconditioner 1900 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1901 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr); 1902 1903 /* Apply transpose of partition of unity operator */ 1904 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1905 1906 /* Second Dirichlet solve and assembling of output */ 1907 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1908 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1909 if (n_B) { 1910 if (pcbddc->switch_static) { 1911 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1912 1913 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1914 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec1_D,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1915 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1916 ierr = VecScatterEnd(pcis->N_to_B,pcis->vec1_B,pcis->vec1_N,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1917 if (!pcbddc->switch_static_change) { 1918 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1919 } else { 1920 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1921 ierr = MatMultTranspose(matis->A,pcis->vec2_N,pcis->vec1_N);CHKERRQ(ierr); 1922 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1923 } 1924 ierr = VecScatterBegin(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1925 ierr = VecScatterEnd(pcis->N_to_D,pcis->vec2_N,pcis->vec3_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1926 } else { 1927 ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1928 } 1929 } else if (pcbddc->switch_static) { /* n_B is zero */ 1930 Mat_IS *matis = (Mat_IS*)(pc->mat->data); 1931 1932 if (!pcbddc->switch_static_change) { 1933 ierr = MatMultTranspose(matis->A,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 1934 } else { 1935 ierr = MatMult(pcbddc->switch_static_change,pcis->vec1_D,pcis->vec1_N);CHKERRQ(ierr); 1936 ierr = MatMultTranspose(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1937 ierr = MatMultTranspose(pcbddc->switch_static_change,pcis->vec2_N,pcis->vec3_D);CHKERRQ(ierr); 1938 } 1939 } 1940 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1941 if (!pcbddc->use_exact_dirichlet_trick && !pcbddc->benign_apply_coarse_only) { 1942 if (pcbddc->switch_static) { 1943 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 1944 } else { 1945 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 1946 } 1947 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1948 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1949 } else { 1950 if (pcbddc->switch_static) { 1951 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 1952 } else { 1953 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1954 } 1955 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1956 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1957 } 1958 if (pcbddc->benign_have_null) { /* set p0 (computed in PCBDDCApplyInterface) */ 1959 ierr = PCBDDCBenignGetOrSetP0(pc,z,PETSC_FALSE);CHKERRQ(ierr); 1960 } 1961 if (pcbddc->ChangeOfBasisMatrix) { 1962 Vec swap; 1963 1964 swap = r; 1965 r = pcbddc->work_change; 1966 pcbddc->work_change = swap; 1967 ierr = VecCopy(z,pcbddc->work_change);CHKERRQ(ierr); 1968 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcbddc->work_change,z);CHKERRQ(ierr); 1969 } 1970 PetscFunctionReturn(0); 1971 } 1972 /* -------------------------------------------------------------------------- */ 1973 1974 #undef __FUNCT__ 1975 #define __FUNCT__ "PCDestroy_BDDC" 1976 PetscErrorCode PCDestroy_BDDC(PC pc) 1977 { 1978 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1979 PetscErrorCode ierr; 1980 1981 PetscFunctionBegin; 1982 /* free BDDC custom data */ 1983 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 1984 /* destroy objects related to topography */ 1985 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 1986 /* free allocated graph structure */ 1987 ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr); 1988 /* free allocated sub schurs structure */ 1989 ierr = PetscFree(pcbddc->sub_schurs);CHKERRQ(ierr); 1990 /* destroy objects for scaling operator */ 1991 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 1992 ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr); 1993 /* free solvers stuff */ 1994 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 1995 /* free global vectors needed in presolve */ 1996 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 1997 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 1998 /* free data created by PCIS */ 1999 ierr = PCISDestroy(pc);CHKERRQ(ierr); 2000 /* remove functions */ 2001 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);CHKERRQ(ierr); 2002 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 2003 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",NULL);CHKERRQ(ierr); 2004 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 2005 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr); 2006 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr); 2007 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr); 2008 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2009 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2010 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2011 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2012 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 2013 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 2014 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 2015 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 2016 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 2017 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr); 2018 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 2019 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 2020 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 2021 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 2022 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",NULL);CHKERRQ(ierr); 2023 /* Free the private data structure */ 2024 ierr = PetscFree(pc->data);CHKERRQ(ierr); 2025 PetscFunctionReturn(0); 2026 } 2027 /* -------------------------------------------------------------------------- */ 2028 2029 #undef __FUNCT__ 2030 #define __FUNCT__ "PCPreSolveChangeRHS_BDDC" 2031 static PetscErrorCode PCPreSolveChangeRHS_BDDC(PC pc, PetscBool* change) 2032 { 2033 PetscFunctionBegin; 2034 *change = PETSC_TRUE; 2035 PetscFunctionReturn(0); 2036 } 2037 2038 #undef __FUNCT__ 2039 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC" 2040 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2041 { 2042 FETIDPMat_ctx mat_ctx; 2043 Vec copy_standard_rhs; 2044 PC_IS* pcis; 2045 PC_BDDC* pcbddc; 2046 PetscErrorCode ierr; 2047 2048 PetscFunctionBegin; 2049 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2050 pcis = (PC_IS*)mat_ctx->pc->data; 2051 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2052 2053 /* 2054 change of basis for physical rhs if needed 2055 It also changes the rhs in case of dirichlet boundaries 2056 TODO: better management when FETIDP will have its own class 2057 */ 2058 ierr = VecDuplicate(standard_rhs,©_standard_rhs);CHKERRQ(ierr); 2059 if (pcbddc->ChangeOfBasisMatrix) { 2060 Vec v2; 2061 ierr = VecDuplicate(standard_rhs,&v2);CHKERRQ(ierr); 2062 ierr = VecCopy(standard_rhs,v2);CHKERRQ(ierr); 2063 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,v2,NULL);CHKERRQ(ierr); 2064 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,v2,copy_standard_rhs);CHKERRQ(ierr); 2065 ierr = VecDestroy(&v2);CHKERRQ(ierr); 2066 } else { 2067 ierr = VecCopy(standard_rhs,copy_standard_rhs);CHKERRQ(ierr); 2068 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,copy_standard_rhs,NULL);CHKERRQ(ierr); 2069 } 2070 /* store vectors for computation of fetidp final solution */ 2071 ierr = VecScatterBegin(pcis->global_to_D,copy_standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2072 ierr = VecScatterEnd(pcis->global_to_D,copy_standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2073 /* scale rhs since it should be unassembled */ 2074 /* TODO use counter scaling? (also below) */ 2075 ierr = VecScatterBegin(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2076 ierr = VecScatterEnd(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2077 /* Apply partition of unity */ 2078 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2079 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,copy_standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2080 if (!pcbddc->switch_static) { 2081 /* compute partially subassembled Schur complement right-hand side */ 2082 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2083 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 2084 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 2085 ierr = VecSet(copy_standard_rhs,0.0);CHKERRQ(ierr); 2086 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,copy_standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2087 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,copy_standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2088 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,copy_standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 2089 ierr = VecScatterBegin(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2090 ierr = VecScatterEnd(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2091 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2092 } 2093 ierr = VecDestroy(©_standard_rhs);CHKERRQ(ierr); 2094 /* BDDC rhs */ 2095 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 2096 if (pcbddc->switch_static) { 2097 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2098 } 2099 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2100 /* apply BDDC */ 2101 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2102 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 2103 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 2104 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 2105 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2106 ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2107 PetscFunctionReturn(0); 2108 } 2109 2110 #undef __FUNCT__ 2111 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS" 2112 /*@ 2113 PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETI-DP linear system using the physical right-hand side 2114 2115 Collective 2116 2117 Input Parameters: 2118 + fetidp_mat - the FETI-DP matrix object obtained by a call to PCBDDCCreateFETIDPOperators 2119 - standard_rhs - the right-hand side of the original linear system 2120 2121 Output Parameters: 2122 . fetidp_flux_rhs - the right-hand side for the FETI-DP linear system 2123 2124 Level: developer 2125 2126 Notes: 2127 2128 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetSolution 2129 @*/ 2130 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 2131 { 2132 FETIDPMat_ctx mat_ctx; 2133 PetscErrorCode ierr; 2134 2135 PetscFunctionBegin; 2136 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2137 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 2138 PetscFunctionReturn(0); 2139 } 2140 /* -------------------------------------------------------------------------- */ 2141 2142 #undef __FUNCT__ 2143 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC" 2144 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2145 { 2146 FETIDPMat_ctx mat_ctx; 2147 PC_IS* pcis; 2148 PC_BDDC* pcbddc; 2149 PetscErrorCode ierr; 2150 2151 PetscFunctionBegin; 2152 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2153 pcis = (PC_IS*)mat_ctx->pc->data; 2154 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 2155 2156 /* apply B_delta^T */ 2157 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2158 ierr = VecScatterEnd(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2159 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 2160 /* compute rhs for BDDC application */ 2161 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 2162 if (pcbddc->switch_static) { 2163 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2164 } 2165 ierr = PetscMemzero(pcbddc->benign_p0,pcbddc->benign_n*sizeof(PetscScalar));CHKERRQ(ierr); 2166 /* apply BDDC */ 2167 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 2168 /* put values into standard global vector */ 2169 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2170 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2171 if (!pcbddc->switch_static) { 2172 /* compute values into the interior if solved for the partially subassembled Schur complement */ 2173 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 2174 ierr = VecAXPY(mat_ctx->temp_solution_D,-1.0,pcis->vec1_D);CHKERRQ(ierr); 2175 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 2176 } 2177 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2178 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2179 if (pcbddc->ChangeOfBasisMatrix) { 2180 Vec v2; 2181 ierr = VecDuplicate(standard_sol,&v2);CHKERRQ(ierr); 2182 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,standard_sol,v2);CHKERRQ(ierr); 2183 ierr = VecCopy(v2,standard_sol);CHKERRQ(ierr); 2184 ierr = VecDestroy(&v2);CHKERRQ(ierr); 2185 } 2186 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 2187 PetscFunctionReturn(0); 2188 } 2189 2190 #undef __FUNCT__ 2191 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution" 2192 /*@ 2193 PCBDDCMatFETIDPGetSolution - Compute the physical solution using the solution of the FETI-DP linear system 2194 2195 Collective 2196 2197 Input Parameters: 2198 + fetidp_mat - the FETI-DP matrix obtained by a call to PCBDDCCreateFETIDPOperators 2199 - fetidp_flux_sol - the solution of the FETI-DP linear system 2200 2201 Output Parameters: 2202 . standard_sol - the solution defined on the physical domain 2203 2204 Level: developer 2205 2206 Notes: 2207 2208 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetRHS 2209 @*/ 2210 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 2211 { 2212 FETIDPMat_ctx mat_ctx; 2213 PetscErrorCode ierr; 2214 2215 PetscFunctionBegin; 2216 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 2217 ierr = PetscUseMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 2218 PetscFunctionReturn(0); 2219 } 2220 /* -------------------------------------------------------------------------- */ 2221 2222 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec); 2223 extern PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec); 2224 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat); 2225 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec); 2226 extern PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec); 2227 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC); 2228 2229 #undef __FUNCT__ 2230 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC" 2231 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 2232 { 2233 2234 FETIDPMat_ctx fetidpmat_ctx; 2235 Mat newmat; 2236 FETIDPPC_ctx fetidppc_ctx; 2237 PC newpc; 2238 MPI_Comm comm; 2239 PetscErrorCode ierr; 2240 2241 PetscFunctionBegin; 2242 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 2243 /* FETIDP linear matrix */ 2244 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 2245 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 2246 ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 2247 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 2248 ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr); 2249 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 2250 ierr = MatSetUp(newmat);CHKERRQ(ierr); 2251 /* FETIDP preconditioner */ 2252 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 2253 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 2254 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 2255 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 2256 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 2257 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 2258 ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 2259 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 2260 ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr); 2261 ierr = PCSetUp(newpc);CHKERRQ(ierr); 2262 /* return pointers for objects created */ 2263 *fetidp_mat=newmat; 2264 *fetidp_pc=newpc; 2265 PetscFunctionReturn(0); 2266 } 2267 2268 #undef __FUNCT__ 2269 #define __FUNCT__ "PCBDDCCreateFETIDPOperators" 2270 /*@ 2271 PCBDDCCreateFETIDPOperators - Create FETI-DP operators 2272 2273 Collective 2274 2275 Input Parameters: 2276 . pc - the BDDC preconditioning context (setup should have been called before) 2277 2278 Output Parameters: 2279 + fetidp_mat - shell FETI-DP matrix object 2280 - fetidp_pc - shell Dirichlet preconditioner for FETI-DP matrix 2281 2282 Options Database Keys: 2283 . -fetidp_fullyredundant <false> - use or not a fully redundant set of Lagrange multipliers 2284 2285 Level: developer 2286 2287 Notes: 2288 Currently the only operations provided for FETI-DP matrix are MatMult and MatMultTranspose 2289 2290 .seealso: PCBDDC, PCBDDCMatFETIDPGetRHS, PCBDDCMatFETIDPGetSolution 2291 @*/ 2292 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 2293 { 2294 PetscErrorCode ierr; 2295 2296 PetscFunctionBegin; 2297 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2298 if (pc->setupcalled) { 2299 ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,Mat*,PC*),(pc,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 2300 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n"); 2301 PetscFunctionReturn(0); 2302 } 2303 /* -------------------------------------------------------------------------- */ 2304 /*MC 2305 PCBDDC - Balancing Domain Decomposition by Constraints. 2306 2307 An implementation of the BDDC preconditioner based on 2308 2309 .vb 2310 [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 2311 [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 2312 [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977 2313 [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 2314 .ve 2315 2316 The matrix to be preconditioned (Pmat) must be of type MATIS. 2317 2318 Currently works with MATIS matrices with local matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers. 2319 2320 It also works with unsymmetric and indefinite problems. 2321 2322 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. 2323 2324 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). 2325 2326 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() 2327 Additional information on dofs can be provided by using PCBDDCSetDofsSplitting(), PCBDDCSetDirichletBoundaries(), PCBDDCSetNeumannBoundaries(), and PCBDDCSetPrimalVerticesIS() and their local counterparts. 2328 2329 Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). Non-singular modes are retained via SVD. 2330 2331 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. 2332 User defined change of basis can be passed to PCBDDC by using PCBDDCSetChangeOfBasisMat() 2333 2334 The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using a MatPartitioning object. 2335 2336 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. 2337 2338 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. 2339 Deluxe scaling is not supported yet for FETI-DP. 2340 2341 Options Database Keys (some of them, run with -h for a complete list): 2342 2343 . -pc_bddc_use_vertices <true> - use or not vertices in primal space 2344 . -pc_bddc_use_edges <true> - use or not edges in primal space 2345 . -pc_bddc_use_faces <false> - use or not faces in primal space 2346 . -pc_bddc_symmetric <true> - symmetric computation of primal basis functions. Specify false for unsymmetric problems 2347 . -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only) 2348 . -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested 2349 . -pc_bddc_switch_static <false> - switches from M_2 (default) to M_3 operator (see reference article [1]) 2350 . -pc_bddc_levels <0> - maximum number of levels for multilevel 2351 . -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) 2352 . -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) 2353 . -pc_bddc_use_deluxe_scaling <false> - use deluxe scaling 2354 . -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) 2355 . -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) 2356 - -pc_bddc_check_level <0> - set verbosity level of debugging output 2357 2358 Options for Dirichlet, Neumann or coarse solver can be set with 2359 .vb 2360 -pc_bddc_dirichlet_ 2361 -pc_bddc_neumann_ 2362 -pc_bddc_coarse_ 2363 .ve 2364 e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. PCBDDC uses by default KPSPREONLY and PCLU. 2365 2366 When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified as 2367 .vb 2368 -pc_bddc_dirichlet_lN_ 2369 -pc_bddc_neumann_lN_ 2370 -pc_bddc_coarse_lN_ 2371 .ve 2372 Note that level number ranges from the finest (0) to the coarsest (N). 2373 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. 2374 .vb 2375 -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3 2376 .ve 2377 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 2378 2379 Level: intermediate 2380 2381 Developer notes: 2382 2383 Contributed by Stefano Zampini 2384 2385 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 2386 M*/ 2387 2388 #undef __FUNCT__ 2389 #define __FUNCT__ "PCCreate_BDDC" 2390 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 2391 { 2392 PetscErrorCode ierr; 2393 PC_BDDC *pcbddc; 2394 2395 PetscFunctionBegin; 2396 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 2397 ierr = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr); 2398 pc->data = (void*)pcbddc; 2399 2400 /* create PCIS data structure */ 2401 ierr = PCISCreate(pc);CHKERRQ(ierr); 2402 2403 /* BDDC customization */ 2404 pcbddc->use_local_adj = PETSC_TRUE; 2405 pcbddc->use_vertices = PETSC_TRUE; 2406 pcbddc->use_edges = PETSC_TRUE; 2407 pcbddc->use_faces = PETSC_FALSE; 2408 pcbddc->use_change_of_basis = PETSC_FALSE; 2409 pcbddc->use_change_on_faces = PETSC_FALSE; 2410 pcbddc->switch_static = PETSC_FALSE; 2411 pcbddc->use_nnsp_true = PETSC_FALSE; 2412 pcbddc->use_qr_single = PETSC_FALSE; 2413 pcbddc->symmetric_primal = PETSC_TRUE; 2414 pcbddc->benign_saddle_point = PETSC_FALSE; 2415 pcbddc->benign_have_null = PETSC_FALSE; 2416 pcbddc->vertex_size = 1; 2417 pcbddc->dbg_flag = 0; 2418 /* private */ 2419 pcbddc->local_primal_size = 0; 2420 pcbddc->local_primal_size_cc = 0; 2421 pcbddc->local_primal_ref_node = 0; 2422 pcbddc->local_primal_ref_mult = 0; 2423 pcbddc->n_vertices = 0; 2424 pcbddc->primal_indices_local_idxs = 0; 2425 pcbddc->recompute_topography = PETSC_FALSE; 2426 pcbddc->coarse_size = -1; 2427 pcbddc->new_primal_space = PETSC_FALSE; 2428 pcbddc->new_primal_space_local = PETSC_FALSE; 2429 pcbddc->global_primal_indices = 0; 2430 pcbddc->onearnullspace = 0; 2431 pcbddc->onearnullvecs_state = 0; 2432 pcbddc->user_primal_vertices = 0; 2433 pcbddc->user_primal_vertices_local = 0; 2434 pcbddc->temp_solution = 0; 2435 pcbddc->original_rhs = 0; 2436 pcbddc->local_mat = 0; 2437 pcbddc->ChangeOfBasisMatrix = 0; 2438 pcbddc->user_ChangeOfBasisMatrix = 0; 2439 pcbddc->coarse_vec = 0; 2440 pcbddc->coarse_ksp = 0; 2441 pcbddc->coarse_phi_B = 0; 2442 pcbddc->coarse_phi_D = 0; 2443 pcbddc->coarse_psi_B = 0; 2444 pcbddc->coarse_psi_D = 0; 2445 pcbddc->vec1_P = 0; 2446 pcbddc->vec1_R = 0; 2447 pcbddc->vec2_R = 0; 2448 pcbddc->local_auxmat1 = 0; 2449 pcbddc->local_auxmat2 = 0; 2450 pcbddc->R_to_B = 0; 2451 pcbddc->R_to_D = 0; 2452 pcbddc->ksp_D = 0; 2453 pcbddc->ksp_R = 0; 2454 pcbddc->NeumannBoundaries = 0; 2455 pcbddc->NeumannBoundariesLocal = 0; 2456 pcbddc->DirichletBoundaries = 0; 2457 pcbddc->DirichletBoundariesLocal = 0; 2458 pcbddc->user_provided_isfordofs = PETSC_FALSE; 2459 pcbddc->n_ISForDofs = 0; 2460 pcbddc->n_ISForDofsLocal = 0; 2461 pcbddc->ISForDofs = 0; 2462 pcbddc->ISForDofsLocal = 0; 2463 pcbddc->ConstraintMatrix = 0; 2464 pcbddc->use_exact_dirichlet_trick = PETSC_TRUE; 2465 pcbddc->coarse_loc_to_glob = 0; 2466 pcbddc->coarsening_ratio = 8; 2467 pcbddc->coarse_adj_red = 0; 2468 pcbddc->current_level = 0; 2469 pcbddc->max_levels = 0; 2470 pcbddc->use_coarse_estimates = PETSC_FALSE; 2471 pcbddc->coarse_eqs_per_proc = 1; 2472 pcbddc->coarse_subassembling = 0; 2473 pcbddc->detect_disconnected = PETSC_FALSE; 2474 pcbddc->n_local_subs = 0; 2475 pcbddc->local_subs = NULL; 2476 2477 /* benign subspace trick */ 2478 pcbddc->benign_change = 0; 2479 pcbddc->benign_compute_correction = PETSC_TRUE; 2480 pcbddc->benign_vec = 0; 2481 pcbddc->benign_original_mat = 0; 2482 pcbddc->benign_sf = 0; 2483 pcbddc->benign_B0 = 0; 2484 pcbddc->benign_n = 0; 2485 pcbddc->benign_p0 = NULL; 2486 pcbddc->benign_p0_lidx = NULL; 2487 pcbddc->benign_p0_gidx = NULL; 2488 pcbddc->benign_null = PETSC_FALSE; 2489 2490 /* create local graph structure */ 2491 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 2492 2493 /* scaling */ 2494 pcbddc->work_scaling = 0; 2495 pcbddc->use_deluxe_scaling = PETSC_FALSE; 2496 2497 /* create sub schurs structure */ 2498 ierr = PCBDDCSubSchursCreate(&pcbddc->sub_schurs);CHKERRQ(ierr); 2499 pcbddc->sub_schurs_rebuild = PETSC_FALSE; 2500 pcbddc->sub_schurs_layers = -1; 2501 pcbddc->sub_schurs_use_useradj = PETSC_FALSE; 2502 2503 pcbddc->computed_rowadj = PETSC_FALSE; 2504 2505 /* adaptivity */ 2506 pcbddc->adaptive_threshold = 0.0; 2507 pcbddc->adaptive_nmax = 0; 2508 pcbddc->adaptive_nmin = 0; 2509 2510 /* function pointers */ 2511 pc->ops->apply = PCApply_BDDC; 2512 pc->ops->applytranspose = PCApplyTranspose_BDDC; 2513 pc->ops->setup = PCSetUp_BDDC; 2514 pc->ops->destroy = PCDestroy_BDDC; 2515 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 2516 pc->ops->view = PCView_BDDC; 2517 pc->ops->applyrichardson = 0; 2518 pc->ops->applysymmetricleft = 0; 2519 pc->ops->applysymmetricright = 0; 2520 pc->ops->presolve = PCPreSolve_BDDC; 2521 pc->ops->postsolve = PCPostSolve_BDDC; 2522 2523 /* composing function */ 2524 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);CHKERRQ(ierr); 2525 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 2526 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesIS_C",PCBDDCSetPrimalVerticesIS_BDDC);CHKERRQ(ierr); 2527 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 2528 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr); 2529 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr); 2530 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr); 2531 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2532 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2533 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2534 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2535 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2536 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2537 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2538 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2539 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 2540 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr); 2541 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 2542 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 2543 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 2544 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 2545 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_BDDC);CHKERRQ(ierr); 2546 PetscFunctionReturn(0); 2547 } 2548 2549