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