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