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