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