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