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