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