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