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 copy_rhs = PETSC_FALSE; 1133 } 1134 pcbddc->rhs_change = PETSC_TRUE; 1135 ierr = MatNullSpaceRemove(pcbddc->NullSpace,rhs);CHKERRQ(ierr); 1136 } 1137 PetscFunctionReturn(0); 1138 } 1139 1140 /* -------------------------------------------------------------------------- */ 1141 #undef __FUNCT__ 1142 #define __FUNCT__ "PCPostSolve_BDDC" 1143 /* -------------------------------------------------------------------------- */ 1144 /* 1145 PCPostSolve_BDDC - Changes the computed solution if a transformation of basis 1146 approach has been selected. Also, restores rhs to its original state. 1147 1148 Input Parameter: 1149 + pc - the preconditioner contex 1150 1151 Application Interface Routine: PCPostSolve() 1152 1153 Notes: 1154 The interface routine PCPostSolve() is not usually called directly by 1155 the user, but instead is called by KSPSolve(). 1156 */ 1157 static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 1158 { 1159 PetscErrorCode ierr; 1160 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1161 1162 PetscFunctionBegin; 1163 if (pcbddc->ChangeOfBasisMatrix) { 1164 PCBDDCChange_ctx change_ctx; 1165 1166 /* get change ctx */ 1167 ierr = MatShellGetContext(pcbddc->new_global_mat,&change_ctx);CHKERRQ(ierr); 1168 1169 /* restore iteration matrix */ 1170 ierr = MatDestroy(&pc->mat);CHKERRQ(ierr); 1171 ierr = PetscObjectReference((PetscObject)change_ctx->original_mat);CHKERRQ(ierr); 1172 pc->mat = change_ctx->original_mat; 1173 1174 /* get solution in original basis */ 1175 if (x) { 1176 PC_IS *pcis = (PC_IS*)(pc->data); 1177 ierr = MatMult(change_ctx->global_change,x,pcis->vec1_global);CHKERRQ(ierr); 1178 ierr = VecCopy(pcis->vec1_global,x);CHKERRQ(ierr); 1179 } 1180 } 1181 1182 /* add solution removed in presolve */ 1183 if (x && pcbddc->rhs_change) { 1184 ierr = VecAXPY(x,1.0,pcbddc->temp_solution);CHKERRQ(ierr); 1185 } 1186 1187 /* restore rhs to its original state */ 1188 if (rhs && pcbddc->rhs_change) { 1189 ierr = VecCopy(pcbddc->original_rhs,rhs);CHKERRQ(ierr); 1190 } 1191 pcbddc->rhs_change = PETSC_FALSE; 1192 1193 /* restore ksp guess state */ 1194 if (ksp) { 1195 ierr = KSPSetInitialGuessNonzero(ksp,pcbddc->ksp_guess_nonzero);CHKERRQ(ierr); 1196 } 1197 PetscFunctionReturn(0); 1198 } 1199 /* -------------------------------------------------------------------------- */ 1200 #undef __FUNCT__ 1201 #define __FUNCT__ "PCSetUp_BDDC" 1202 /* -------------------------------------------------------------------------- */ 1203 /* 1204 PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner 1205 by setting data structures and options. 1206 1207 Input Parameter: 1208 + pc - the preconditioner context 1209 1210 Application Interface Routine: PCSetUp() 1211 1212 Notes: 1213 The interface routine PCSetUp() is not usually called directly by 1214 the user, but instead is called by PCApply() if necessary. 1215 */ 1216 PetscErrorCode PCSetUp_BDDC(PC pc) 1217 { 1218 PetscErrorCode ierr; 1219 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 1220 Mat_IS* matis; 1221 MatNullSpace nearnullspace; 1222 PetscInt nrows,ncols; 1223 PetscBool computetopography,computesolvers,computesubschurs; 1224 PetscBool computeconstraintsmatrix; 1225 PetscBool new_nearnullspace_provided,ismatis; 1226 1227 PetscFunctionBegin; 1228 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATIS,&ismatis);CHKERRQ(ierr); 1229 if (!ismatis) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PCBDDC preconditioner requires matrix of type MATIS"); 1230 ierr = MatGetSize(pc->pmat,&nrows,&ncols);CHKERRQ(ierr); 1231 if (nrows != ncols) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PCBDDC preconditioner requires a square preconditioning matrix"); 1232 matis = (Mat_IS*)pc->pmat->data; 1233 /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other future nonoverlapping preconditioners */ 1234 /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup 1235 Also, BDDC directly build the Dirichlet problem */ 1236 /* split work */ 1237 if (pc->setupcalled) { 1238 if (pc->flag == SAME_NONZERO_PATTERN) { 1239 computetopography = PETSC_FALSE; 1240 computesolvers = PETSC_TRUE; 1241 } else { /* DIFFERENT_NONZERO_PATTERN */ 1242 computetopography = PETSC_TRUE; 1243 computesolvers = PETSC_TRUE; 1244 } 1245 } else { 1246 computetopography = PETSC_TRUE; 1247 computesolvers = PETSC_TRUE; 1248 } 1249 if (pcbddc->recompute_topography) { 1250 computetopography = PETSC_TRUE; 1251 } 1252 computeconstraintsmatrix = PETSC_FALSE; 1253 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"); 1254 pcbddc->adaptive_selection = (PetscBool)(pcbddc->adaptive_threshold > 0.0 && pcbddc->use_deluxe_scaling); 1255 if (pcbddc->adaptive_selection) pcbddc->use_faces = PETSC_TRUE; 1256 1257 computesubschurs = (PetscBool)(pcbddc->adaptive_selection || pcbddc->use_deluxe_scaling); 1258 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"); 1259 1260 /* Get stdout for dbg */ 1261 if (pcbddc->dbg_flag) { 1262 if (!pcbddc->dbg_viewer) { 1263 pcbddc->dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)pc)); 1264 ierr = PetscViewerASCIIPushSynchronized(pcbddc->dbg_viewer);CHKERRQ(ierr); 1265 } 1266 ierr = PetscViewerASCIIAddTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1267 } 1268 1269 if (pcbddc->user_ChangeOfBasisMatrix) { 1270 /* use_change_of_basis flag is used to automatically compute a change of basis from constraints */ 1271 pcbddc->use_change_of_basis = PETSC_FALSE; 1272 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->user_ChangeOfBasisMatrix);CHKERRQ(ierr); 1273 } else { 1274 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1275 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1276 pcbddc->local_mat = matis->A; 1277 } 1278 1279 /* workaround for reals */ 1280 #if !defined(PETSC_USE_COMPLEX) 1281 if (matis->A->symmetric_set) { 1282 ierr = MatSetOption(pcbddc->local_mat,MAT_HERMITIAN,matis->A->symmetric);CHKERRQ(ierr); 1283 } 1284 #endif 1285 1286 /* Set up all the "iterative substructuring" common block without computing solvers */ 1287 { 1288 Mat temp_mat; 1289 1290 temp_mat = matis->A; 1291 matis->A = pcbddc->local_mat; 1292 ierr = PCISSetUp(pc,PETSC_FALSE);CHKERRQ(ierr); 1293 pcbddc->local_mat = matis->A; 1294 matis->A = temp_mat; 1295 } 1296 1297 /* Analyze interface and setup sub_schurs data */ 1298 if (computetopography) { 1299 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 1300 computeconstraintsmatrix = PETSC_TRUE; 1301 } 1302 1303 /* Setup local dirichlet solver ksp_D and sub_schurs solvers */ 1304 if (computesolvers) { 1305 PCBDDCSubSchurs sub_schurs=pcbddc->sub_schurs; 1306 1307 if (computesubschurs && computetopography) { 1308 ierr = PCBDDCInitSubSchurs(pc);CHKERRQ(ierr); 1309 } 1310 if (sub_schurs->use_mumps) { 1311 if (computesubschurs) { 1312 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1313 } 1314 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1315 } else { 1316 ierr = PCBDDCSetUpLocalSolvers(pc,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr); 1317 if (computesubschurs) { 1318 ierr = PCBDDCSetUpSubSchurs(pc);CHKERRQ(ierr); 1319 } 1320 } 1321 if (pcbddc->adaptive_selection) { 1322 ierr = PCBDDCAdaptiveSelection(pc);CHKERRQ(ierr); 1323 computeconstraintsmatrix = PETSC_TRUE; 1324 } 1325 } 1326 1327 /* infer if NullSpace object attached to Mat via MatSetNearNullSpace has changed */ 1328 new_nearnullspace_provided = PETSC_FALSE; 1329 ierr = MatGetNearNullSpace(pc->pmat,&nearnullspace);CHKERRQ(ierr); 1330 if (pcbddc->onearnullspace) { /* already used nearnullspace */ 1331 if (!nearnullspace) { /* near null space attached to mat has been destroyed */ 1332 new_nearnullspace_provided = PETSC_TRUE; 1333 } else { 1334 /* determine if the two nullspaces are different (should be lightweight) */ 1335 if (nearnullspace != pcbddc->onearnullspace) { 1336 new_nearnullspace_provided = PETSC_TRUE; 1337 } else { /* maybe the user has changed the content of the nearnullspace so check vectors ObjectStateId */ 1338 PetscInt i; 1339 const Vec *nearnullvecs; 1340 PetscObjectState state; 1341 PetscInt nnsp_size; 1342 ierr = MatNullSpaceGetVecs(nearnullspace,NULL,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 1343 for (i=0;i<nnsp_size;i++) { 1344 ierr = PetscObjectStateGet((PetscObject)nearnullvecs[i],&state);CHKERRQ(ierr); 1345 if (pcbddc->onearnullvecs_state[i] != state) { 1346 new_nearnullspace_provided = PETSC_TRUE; 1347 break; 1348 } 1349 } 1350 } 1351 } 1352 } else { 1353 if (!nearnullspace) { /* both nearnullspaces are null */ 1354 new_nearnullspace_provided = PETSC_FALSE; 1355 } else { /* nearnullspace attached later */ 1356 new_nearnullspace_provided = PETSC_TRUE; 1357 } 1358 } 1359 1360 /* Setup constraints and related work vectors */ 1361 /* reset primal space flags */ 1362 pcbddc->new_primal_space = PETSC_FALSE; 1363 pcbddc->new_primal_space_local = PETSC_FALSE; 1364 if (computeconstraintsmatrix || new_nearnullspace_provided) { 1365 /* It also sets the primal space flags */ 1366 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 1367 /* Allocate needed local vectors (which depends on quantities defined during ConstraintsSetUp) */ 1368 ierr = PCBDDCSetUpLocalWorkVectors(pc);CHKERRQ(ierr); 1369 } 1370 1371 if (computesolvers || pcbddc->new_primal_space) { 1372 if (pcbddc->use_change_of_basis) { 1373 PC_IS *pcis = (PC_IS*)(pc->data); 1374 1375 ierr = PCBDDCComputeLocalMatrix(pc,pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 1376 /* get submatrices */ 1377 ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr); 1378 ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr); 1379 ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr); 1380 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);CHKERRQ(ierr); 1381 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);CHKERRQ(ierr); 1382 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);CHKERRQ(ierr); 1383 /* set flag in pcis to not reuse submatrices during PCISCreate */ 1384 pcis->reusesubmatrices = PETSC_FALSE; 1385 } else if (!pcbddc->user_ChangeOfBasisMatrix) { 1386 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 1387 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1388 pcbddc->local_mat = matis->A; 1389 } 1390 /* SetUp coarse and local Neumann solvers */ 1391 ierr = PCBDDCSetUpSolvers(pc);CHKERRQ(ierr); 1392 /* SetUp Scaling operator */ 1393 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 1394 } 1395 1396 if (pcbddc->dbg_flag) { 1397 ierr = PetscViewerASCIISubtractTab(pcbddc->dbg_viewer,2*pcbddc->current_level);CHKERRQ(ierr); 1398 } 1399 PetscFunctionReturn(0); 1400 } 1401 1402 /* -------------------------------------------------------------------------- */ 1403 /* 1404 PCApply_BDDC - Applies the BDDC operator to a vector. 1405 1406 Input Parameters: 1407 + pc - the preconditioner context 1408 - r - input vector (global) 1409 1410 Output Parameter: 1411 . z - output vector (global) 1412 1413 Application Interface Routine: PCApply() 1414 */ 1415 #undef __FUNCT__ 1416 #define __FUNCT__ "PCApply_BDDC" 1417 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 1418 { 1419 PC_IS *pcis = (PC_IS*)(pc->data); 1420 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1421 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1422 PetscErrorCode ierr; 1423 const PetscScalar one = 1.0; 1424 const PetscScalar m_one = -1.0; 1425 const PetscScalar zero = 0.0; 1426 1427 /* This code is similar to that provided in nn.c for PCNN 1428 NN interface preconditioner changed to BDDC 1429 Added support for M_3 preconditioner in the reference article (code is active if pcbddc->switch_static == PETSC_TRUE) */ 1430 1431 PetscFunctionBegin; 1432 if (!pcbddc->use_exact_dirichlet_trick) { 1433 ierr = VecCopy(r,z);CHKERRQ(ierr); 1434 /* First Dirichlet solve */ 1435 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1436 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1437 /* 1438 Assembling right hand side for BDDC operator 1439 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1440 - pcis->vec1_B the interface part of the global vector z 1441 */ 1442 if (n_D) { 1443 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1444 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1445 if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 1446 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1447 } else { 1448 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1449 } 1450 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1451 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1452 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1453 } else { 1454 if (pcbddc->switch_static) { 1455 ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr); 1456 } 1457 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1458 } 1459 1460 /* Apply interface preconditioner 1461 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1462 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_FALSE);CHKERRQ(ierr); 1463 1464 /* Apply transpose of partition of unity operator */ 1465 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1466 1467 /* Second Dirichlet solve and assembling of output */ 1468 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1469 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1470 if (n_B) { 1471 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1472 if (pcbddc->switch_static) { ierr = MatMultAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 1473 } else if (pcbddc->switch_static) { 1474 ierr = MatMult(pcis->A_II,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 1475 } 1476 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1477 if (!pcbddc->use_exact_dirichlet_trick) { 1478 if (pcbddc->switch_static) { 1479 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 1480 } else { 1481 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 1482 } 1483 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1484 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1485 } else { 1486 if (pcbddc->switch_static) { 1487 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 1488 } else { 1489 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1490 } 1491 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1492 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1493 } 1494 PetscFunctionReturn(0); 1495 } 1496 1497 /* -------------------------------------------------------------------------- */ 1498 /* 1499 PCApplyTranspose_BDDC - Applies the transpose of the BDDC operator to a vector. 1500 1501 Input Parameters: 1502 + pc - the preconditioner context 1503 - r - input vector (global) 1504 1505 Output Parameter: 1506 . z - output vector (global) 1507 1508 Application Interface Routine: PCApplyTranspose() 1509 */ 1510 #undef __FUNCT__ 1511 #define __FUNCT__ "PCApplyTranspose_BDDC" 1512 PetscErrorCode PCApplyTranspose_BDDC(PC pc,Vec r,Vec z) 1513 { 1514 PC_IS *pcis = (PC_IS*)(pc->data); 1515 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1516 PetscInt n_B = pcis->n_B, n_D = pcis->n - n_B; 1517 PetscErrorCode ierr; 1518 const PetscScalar one = 1.0; 1519 const PetscScalar m_one = -1.0; 1520 const PetscScalar zero = 0.0; 1521 1522 PetscFunctionBegin; 1523 if (!pcbddc->use_exact_dirichlet_trick) { 1524 ierr = VecCopy(r,z);CHKERRQ(ierr); 1525 /* First Dirichlet solve */ 1526 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1527 ierr = VecScatterEnd(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1528 /* 1529 Assembling right hand side for BDDC operator 1530 - pcis->vec1_D for the Dirichlet part (if needed, i.e. pcbddc->switch_static == PETSC_TRUE) 1531 - pcis->vec1_B the interface part of the global vector z 1532 */ 1533 if (n_D) { 1534 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1535 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 1536 if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 1537 ierr = MatMultTranspose(pcis->A_IB,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 1538 } else { 1539 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1540 } 1541 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1542 ierr = VecScatterEnd(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1543 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 1544 } else { 1545 if (pcbddc->switch_static) { 1546 ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr); 1547 } 1548 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 1549 } 1550 1551 /* Apply interface preconditioner 1552 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 1553 ierr = PCBDDCApplyInterfacePreconditioner(pc,PETSC_TRUE);CHKERRQ(ierr); 1554 1555 /* Apply transpose of partition of unity operator */ 1556 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 1557 1558 /* Second Dirichlet solve and assembling of output */ 1559 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1560 ierr = VecScatterEnd(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1561 if (n_B) { 1562 ierr = MatMultTranspose(pcis->A_BI,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 1563 if (pcbddc->switch_static) { ierr = MatMultTransposeAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 1564 } else if (pcbddc->switch_static) { 1565 ierr = MatMultTranspose(pcis->A_II,pcis->vec1_D,pcis->vec3_D);CHKERRQ(ierr); 1566 } 1567 ierr = KSPSolveTranspose(pcbddc->ksp_D,pcis->vec3_D,pcis->vec4_D);CHKERRQ(ierr); 1568 if (!pcbddc->use_exact_dirichlet_trick) { 1569 if (pcbddc->switch_static) { 1570 ierr = VecAXPBYPCZ(pcis->vec2_D,m_one,one,m_one,pcis->vec4_D,pcis->vec1_D);CHKERRQ(ierr); 1571 } else { 1572 ierr = VecAXPBY(pcis->vec2_D,m_one,m_one,pcis->vec4_D);CHKERRQ(ierr); 1573 } 1574 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1575 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1576 } else { 1577 if (pcbddc->switch_static) { 1578 ierr = VecAXPBY(pcis->vec4_D,one,m_one,pcis->vec1_D);CHKERRQ(ierr); 1579 } else { 1580 ierr = VecScale(pcis->vec4_D,m_one);CHKERRQ(ierr); 1581 } 1582 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1583 ierr = VecScatterEnd(pcis->global_to_D,pcis->vec4_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1584 } 1585 PetscFunctionReturn(0); 1586 } 1587 /* -------------------------------------------------------------------------- */ 1588 1589 #undef __FUNCT__ 1590 #define __FUNCT__ "PCDestroy_BDDC" 1591 PetscErrorCode PCDestroy_BDDC(PC pc) 1592 { 1593 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1594 PetscErrorCode ierr; 1595 1596 PetscFunctionBegin; 1597 /* free data created by PCIS */ 1598 ierr = PCISDestroy(pc);CHKERRQ(ierr); 1599 /* free BDDC custom data */ 1600 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 1601 /* destroy objects related to topography */ 1602 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 1603 /* free allocated graph structure */ 1604 ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr); 1605 /* free allocated sub schurs structure */ 1606 ierr = PetscFree(pcbddc->sub_schurs);CHKERRQ(ierr); 1607 /* destroy objects for scaling operator */ 1608 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 1609 ierr = PetscFree(pcbddc->deluxe_ctx);CHKERRQ(ierr); 1610 /* free solvers stuff */ 1611 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 1612 /* free global vectors needed in presolve */ 1613 ierr = VecDestroy(&pcbddc->temp_solution);CHKERRQ(ierr); 1614 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 1615 /* free stuff for change of basis hooks */ 1616 if (pcbddc->new_global_mat) { 1617 PCBDDCChange_ctx change_ctx; 1618 ierr = MatShellGetContext(pcbddc->new_global_mat,&change_ctx);CHKERRQ(ierr); 1619 ierr = MatDestroy(&change_ctx->original_mat);CHKERRQ(ierr); 1620 ierr = MatDestroy(&change_ctx->global_change);CHKERRQ(ierr); 1621 ierr = VecDestroyVecs(2,&change_ctx->work);CHKERRQ(ierr); 1622 ierr = PetscFree(change_ctx);CHKERRQ(ierr); 1623 } 1624 ierr = MatDestroy(&pcbddc->new_global_mat);CHKERRQ(ierr); 1625 /* remove functions */ 1626 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",NULL);CHKERRQ(ierr); 1627 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 1628 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 1629 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",NULL);CHKERRQ(ierr); 1630 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",NULL);CHKERRQ(ierr); 1631 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",NULL);CHKERRQ(ierr); 1632 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",NULL);CHKERRQ(ierr); 1633 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 1634 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 1635 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 1636 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 1637 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 1638 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",NULL);CHKERRQ(ierr); 1639 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 1640 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",NULL);CHKERRQ(ierr); 1641 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 1642 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",NULL);CHKERRQ(ierr); 1643 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 1644 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 1645 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 1646 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 1647 /* Free the private data structure */ 1648 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1649 PetscFunctionReturn(0); 1650 } 1651 /* -------------------------------------------------------------------------- */ 1652 1653 #undef __FUNCT__ 1654 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC" 1655 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 1656 { 1657 FETIDPMat_ctx mat_ctx; 1658 Vec copy_standard_rhs; 1659 PC_IS* pcis; 1660 PC_BDDC* pcbddc; 1661 PetscErrorCode ierr; 1662 1663 PetscFunctionBegin; 1664 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1665 pcis = (PC_IS*)mat_ctx->pc->data; 1666 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 1667 1668 /* 1669 change of basis for physical rhs if needed 1670 It also changes the rhs in case of dirichlet boundaries 1671 TODO: better management when FETIDP will have its own class 1672 */ 1673 ierr = VecDuplicate(standard_rhs,©_standard_rhs);CHKERRQ(ierr); 1674 ierr = VecCopy(standard_rhs,copy_standard_rhs);CHKERRQ(ierr); 1675 ierr = PCPreSolve_BDDC(mat_ctx->pc,NULL,copy_standard_rhs,NULL);CHKERRQ(ierr); 1676 /* store vectors for computation of fetidp final solution */ 1677 ierr = VecScatterBegin(pcis->global_to_D,copy_standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1678 ierr = VecScatterEnd(pcis->global_to_D,copy_standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1679 /* scale rhs since it should be unassembled */ 1680 /* TODO use counter scaling? (also below) */ 1681 ierr = VecScatterBegin(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1682 ierr = VecScatterEnd(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1683 /* Apply partition of unity */ 1684 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1685 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,copy_standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 1686 if (!pcbddc->switch_static) { 1687 /* compute partially subassembled Schur complement right-hand side */ 1688 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1689 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 1690 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 1691 ierr = VecSet(copy_standard_rhs,0.0);CHKERRQ(ierr); 1692 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,copy_standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1693 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,copy_standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1694 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,copy_standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 1695 ierr = VecScatterBegin(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1696 ierr = VecScatterEnd(pcis->global_to_B,copy_standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1697 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1698 } 1699 ierr = VecDestroy(©_standard_rhs);CHKERRQ(ierr); 1700 /* BDDC rhs */ 1701 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 1702 if (pcbddc->switch_static) { 1703 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1704 } 1705 /* apply BDDC */ 1706 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 1707 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 1708 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 1709 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 1710 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1711 ierr = VecScatterEnd(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1712 PetscFunctionReturn(0); 1713 } 1714 1715 #undef __FUNCT__ 1716 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS" 1717 /*@ 1718 PCBDDCMatFETIDPGetRHS - Compute the right-hand side for FETI-DP linear system using the physical right-hand side 1719 1720 Collective 1721 1722 Input Parameters: 1723 + fetidp_mat - the FETI-DP matrix object obtained by a call to PCBDDCCreateFETIDPOperators 1724 - standard_rhs - the right-hand side of the original linear system 1725 1726 Output Parameters: 1727 . fetidp_flux_rhs - the right-hand side for the FETI-DP linear system 1728 1729 Level: developer 1730 1731 Notes: 1732 1733 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetSolution 1734 @*/ 1735 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 1736 { 1737 FETIDPMat_ctx mat_ctx; 1738 PetscErrorCode ierr; 1739 1740 PetscFunctionBegin; 1741 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1742 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 1743 PetscFunctionReturn(0); 1744 } 1745 /* -------------------------------------------------------------------------- */ 1746 1747 #undef __FUNCT__ 1748 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC" 1749 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 1750 { 1751 FETIDPMat_ctx mat_ctx; 1752 PC_IS* pcis; 1753 PC_BDDC* pcbddc; 1754 PetscErrorCode ierr; 1755 1756 PetscFunctionBegin; 1757 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1758 pcis = (PC_IS*)mat_ctx->pc->data; 1759 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 1760 1761 /* apply B_delta^T */ 1762 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1763 ierr = VecScatterEnd (mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1764 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 1765 /* compute rhs for BDDC application */ 1766 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1767 if (pcbddc->switch_static) { 1768 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1769 } 1770 /* apply BDDC */ 1771 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc,PETSC_FALSE);CHKERRQ(ierr); 1772 /* put values into standard global vector */ 1773 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1774 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1775 if (!pcbddc->switch_static) { 1776 /* compute values into the interior if solved for the partially subassembled Schur complement */ 1777 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 1778 ierr = VecAXPY(mat_ctx->temp_solution_D,-1.0,pcis->vec1_D);CHKERRQ(ierr); 1779 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1780 } 1781 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1782 ierr = VecScatterEnd (pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1783 /* final change of basis if needed 1784 Is also sums the dirichlet part removed during RHS assembling */ 1785 ierr = PCPostSolve_BDDC(mat_ctx->pc,NULL,NULL,standard_sol);CHKERRQ(ierr); 1786 PetscFunctionReturn(0); 1787 } 1788 1789 #undef __FUNCT__ 1790 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution" 1791 /*@ 1792 PCBDDCMatFETIDPGetSolution - Compute the physical solution using the solution of the FETI-DP linear system 1793 1794 Collective 1795 1796 Input Parameters: 1797 + fetidp_mat - the FETI-DP matrix obtained by a call to PCBDDCCreateFETIDPOperators 1798 - fetidp_flux_sol - the solution of the FETI-DP linear system 1799 1800 Output Parameters: 1801 . standard_sol - the solution defined on the physical domain 1802 1803 Level: developer 1804 1805 Notes: 1806 1807 .seealso: PCBDDC, PCBDDCCreateFETIDPOperators, PCBDDCMatFETIDPGetRHS 1808 @*/ 1809 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 1810 { 1811 FETIDPMat_ctx mat_ctx; 1812 PetscErrorCode ierr; 1813 1814 PetscFunctionBegin; 1815 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1816 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 1817 PetscFunctionReturn(0); 1818 } 1819 /* -------------------------------------------------------------------------- */ 1820 1821 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec); 1822 extern PetscErrorCode FETIDPMatMultTranspose(Mat,Vec,Vec); 1823 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat); 1824 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec); 1825 extern PetscErrorCode FETIDPPCApplyTranspose(PC,Vec,Vec); 1826 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC); 1827 1828 #undef __FUNCT__ 1829 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC" 1830 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 1831 { 1832 1833 FETIDPMat_ctx fetidpmat_ctx; 1834 Mat newmat; 1835 FETIDPPC_ctx fetidppc_ctx; 1836 PC newpc; 1837 MPI_Comm comm; 1838 PetscErrorCode ierr; 1839 1840 PetscFunctionBegin; 1841 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 1842 /* FETIDP linear matrix */ 1843 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 1844 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 1845 ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 1846 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 1847 ierr = MatShellSetOperation(newmat,MATOP_MULT_TRANSPOSE,(void (*)(void))FETIDPMatMultTranspose);CHKERRQ(ierr); 1848 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 1849 ierr = MatSetUp(newmat);CHKERRQ(ierr); 1850 /* FETIDP preconditioner */ 1851 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 1852 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 1853 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 1854 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 1855 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 1856 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 1857 ierr = PCShellSetApplyTranspose(newpc,FETIDPPCApplyTranspose);CHKERRQ(ierr); 1858 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 1859 ierr = PCSetOperators(newpc,newmat,newmat);CHKERRQ(ierr); 1860 ierr = PCSetUp(newpc);CHKERRQ(ierr); 1861 /* return pointers for objects created */ 1862 *fetidp_mat=newmat; 1863 *fetidp_pc=newpc; 1864 PetscFunctionReturn(0); 1865 } 1866 1867 #undef __FUNCT__ 1868 #define __FUNCT__ "PCBDDCCreateFETIDPOperators" 1869 /*@ 1870 PCBDDCCreateFETIDPOperators - Create FETI-DP operators 1871 1872 Collective 1873 1874 Input Parameters: 1875 . pc - the BDDC preconditioning context (setup should have been called before) 1876 1877 Output Parameters: 1878 + fetidp_mat - shell FETI-DP matrix object 1879 - fetidp_pc - shell Dirichlet preconditioner for FETI-DP matrix 1880 1881 Options Database Keys: 1882 . -fetidp_fullyredundant <false> - use or not a fully redundant set of Lagrange multipliers 1883 1884 Level: developer 1885 1886 Notes: 1887 Currently the only operations provided for FETI-DP matrix are MatMult and MatMultTranspose 1888 1889 .seealso: PCBDDC, PCBDDCMatFETIDPGetRHS, PCBDDCMatFETIDPGetSolution 1890 @*/ 1891 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 1892 { 1893 PetscErrorCode ierr; 1894 1895 PetscFunctionBegin; 1896 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1897 if (pc->setupcalled) { 1898 ierr = PetscUseMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,Mat*,PC*),(pc,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 1899 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n"); 1900 PetscFunctionReturn(0); 1901 } 1902 /* -------------------------------------------------------------------------- */ 1903 /*MC 1904 PCBDDC - Balancing Domain Decomposition by Constraints. 1905 1906 An implementation of the BDDC preconditioner based on 1907 1908 .vb 1909 [1] C. R. Dohrmann. "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 1910 [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 1911 [3] J. Mandel, B. Sousedik, C. R. Dohrmann. "Multispace and Multilevel BDDC", http://arxiv.org/abs/0712.3977 1912 [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 1913 .ve 1914 1915 The matrix to be preconditioned (Pmat) must be of type MATIS. 1916 1917 Currently works with MATIS matrices with local matrices of type MATSEQAIJ, MATSEQBAIJ or MATSEQSBAIJ, either with real or complex numbers. 1918 1919 It also works with unsymmetric and indefinite problems. 1920 1921 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. 1922 1923 Approximate local solvers are automatically adapted for singular linear problems (see [1]) if the user has provided the nullspace using PCBDDCSetNullSpace() 1924 1925 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() 1926 Additional information on dofs can be provided by using PCBDDCSetDofsSplitting(), PCBDDCSetDirichletBoundaries(), PCBDDCSetNeumannBoundaries(), and PCBDDCSetPrimalVerticesLocalIS() 1927 1928 Constraints can be customized by attaching a MatNullSpace object to the MATIS matrix via MatSetNearNullSpace(). Non-singular modes are retained via SVD. 1929 1930 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. 1931 User defined change of basis can be passed to PCBDDC by using PCBDDCSetChangeOfBasisMat() 1932 1933 The PETSc implementation also supports multilevel BDDC [3]. Coarse grids are partitioned using a MatPartitioning object. 1934 1935 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. 1936 1937 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. 1938 Deluxe scaling is not supported yet for FETI-DP. 1939 1940 Options Database Keys (some of them, run with -h for a complete list): 1941 1942 . -pc_bddc_use_vertices <true> - use or not vertices in primal space 1943 . -pc_bddc_use_edges <true> - use or not edges in primal space 1944 . -pc_bddc_use_faces <false> - use or not faces in primal space 1945 . -pc_bddc_symmetric <true> - symmetric computation of primal basis functions. Specify false for unsymmetric problems 1946 . -pc_bddc_use_change_of_basis <false> - use change of basis approach (on edges only) 1947 . -pc_bddc_use_change_on_faces <false> - use change of basis approach on faces if change of basis has been requested 1948 . -pc_bddc_switch_static <false> - switches from M_2 (default) to M_3 operator (see reference article [1]) 1949 . -pc_bddc_levels <0> - maximum number of levels for multilevel 1950 . -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) 1951 . -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) 1952 . -pc_bddc_use_deluxe_scaling <false> - use deluxe scaling 1953 . -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) 1954 . -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) 1955 - -pc_bddc_check_level <0> - set verbosity level of debugging output 1956 1957 Options for Dirichlet, Neumann or coarse solver can be set with 1958 .vb 1959 -pc_bddc_dirichlet_ 1960 -pc_bddc_neumann_ 1961 -pc_bddc_coarse_ 1962 .ve 1963 e.g -pc_bddc_dirichlet_ksp_type richardson -pc_bddc_dirichlet_pc_type gamg. PCBDDC uses by default KPSPREONLY and PCLU. 1964 1965 When using a multilevel approach, solvers' options at the N-th level (N > 1) can be specified as 1966 .vb 1967 -pc_bddc_dirichlet_lN_ 1968 -pc_bddc_neumann_lN_ 1969 -pc_bddc_coarse_lN_ 1970 .ve 1971 Note that level number ranges from the finest (0) to the coarsest (N). 1972 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. 1973 .vb 1974 -pc_bddc_coarse_pc_bddc_adaptive_threshold 5 -pc_bddc_coarse_l1_pc_bddc_redistribute 3 1975 .ve 1976 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 1977 1978 Level: intermediate 1979 1980 Developer notes: 1981 1982 Contributed by Stefano Zampini 1983 1984 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 1985 M*/ 1986 1987 #undef __FUNCT__ 1988 #define __FUNCT__ "PCCreate_BDDC" 1989 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 1990 { 1991 PetscErrorCode ierr; 1992 PC_BDDC *pcbddc; 1993 1994 PetscFunctionBegin; 1995 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 1996 ierr = PetscNewLog(pc,&pcbddc);CHKERRQ(ierr); 1997 pc->data = (void*)pcbddc; 1998 1999 /* create PCIS data structure */ 2000 ierr = PCISCreate(pc);CHKERRQ(ierr); 2001 2002 /* BDDC customization */ 2003 pcbddc->use_local_adj = PETSC_TRUE; 2004 pcbddc->use_vertices = PETSC_TRUE; 2005 pcbddc->use_edges = PETSC_TRUE; 2006 pcbddc->use_faces = PETSC_FALSE; 2007 pcbddc->use_change_of_basis = PETSC_FALSE; 2008 pcbddc->use_change_on_faces = PETSC_FALSE; 2009 pcbddc->switch_static = PETSC_FALSE; 2010 pcbddc->use_nnsp_true = PETSC_FALSE; 2011 pcbddc->use_qr_single = PETSC_FALSE; 2012 pcbddc->symmetric_primal = PETSC_TRUE; 2013 pcbddc->dbg_flag = 0; 2014 /* private */ 2015 pcbddc->local_primal_size = 0; 2016 pcbddc->local_primal_size_cc = 0; 2017 pcbddc->local_primal_ref_node = 0; 2018 pcbddc->local_primal_ref_mult = 0; 2019 pcbddc->n_vertices = 0; 2020 pcbddc->primal_indices_local_idxs = 0; 2021 pcbddc->recompute_topography = PETSC_FALSE; 2022 pcbddc->coarse_size = -1; 2023 pcbddc->new_primal_space = PETSC_FALSE; 2024 pcbddc->new_primal_space_local = PETSC_FALSE; 2025 pcbddc->global_primal_indices = 0; 2026 pcbddc->onearnullspace = 0; 2027 pcbddc->onearnullvecs_state = 0; 2028 pcbddc->user_primal_vertices = 0; 2029 pcbddc->NullSpace = 0; 2030 pcbddc->temp_solution = 0; 2031 pcbddc->original_rhs = 0; 2032 pcbddc->local_mat = 0; 2033 pcbddc->ChangeOfBasisMatrix = 0; 2034 pcbddc->user_ChangeOfBasisMatrix = 0; 2035 pcbddc->new_global_mat = 0; 2036 pcbddc->coarse_vec = 0; 2037 pcbddc->coarse_ksp = 0; 2038 pcbddc->coarse_phi_B = 0; 2039 pcbddc->coarse_phi_D = 0; 2040 pcbddc->coarse_psi_B = 0; 2041 pcbddc->coarse_psi_D = 0; 2042 pcbddc->vec1_P = 0; 2043 pcbddc->vec1_R = 0; 2044 pcbddc->vec2_R = 0; 2045 pcbddc->local_auxmat1 = 0; 2046 pcbddc->local_auxmat2 = 0; 2047 pcbddc->R_to_B = 0; 2048 pcbddc->R_to_D = 0; 2049 pcbddc->ksp_D = 0; 2050 pcbddc->ksp_R = 0; 2051 pcbddc->NeumannBoundaries = 0; 2052 pcbddc->NeumannBoundariesLocal = 0; 2053 pcbddc->DirichletBoundaries = 0; 2054 pcbddc->DirichletBoundariesLocal = 0; 2055 pcbddc->user_provided_isfordofs = PETSC_FALSE; 2056 pcbddc->n_ISForDofs = 0; 2057 pcbddc->n_ISForDofsLocal = 0; 2058 pcbddc->ISForDofs = 0; 2059 pcbddc->ISForDofsLocal = 0; 2060 pcbddc->ConstraintMatrix = 0; 2061 pcbddc->use_exact_dirichlet_trick = PETSC_TRUE; 2062 pcbddc->coarse_loc_to_glob = 0; 2063 pcbddc->coarsening_ratio = 8; 2064 pcbddc->coarse_adj_red = 0; 2065 pcbddc->current_level = 0; 2066 pcbddc->max_levels = 0; 2067 pcbddc->use_coarse_estimates = PETSC_FALSE; 2068 pcbddc->redistribute_coarse = 0; 2069 pcbddc->coarse_subassembling = 0; 2070 pcbddc->coarse_subassembling_init = 0; 2071 2072 /* create local graph structure */ 2073 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 2074 2075 /* scaling */ 2076 pcbddc->work_scaling = 0; 2077 pcbddc->use_deluxe_scaling = PETSC_FALSE; 2078 pcbddc->faster_deluxe = PETSC_FALSE; 2079 2080 /* create sub schurs structure */ 2081 ierr = PCBDDCSubSchursCreate(&pcbddc->sub_schurs);CHKERRQ(ierr); 2082 pcbddc->sub_schurs_rebuild = PETSC_FALSE; 2083 pcbddc->sub_schurs_layers = -1; 2084 pcbddc->sub_schurs_use_useradj = PETSC_FALSE; 2085 2086 pcbddc->computed_rowadj = PETSC_FALSE; 2087 2088 /* adaptivity */ 2089 pcbddc->adaptive_threshold = 0.0; 2090 pcbddc->adaptive_nmax = 0; 2091 pcbddc->adaptive_nmin = 0; 2092 2093 /* function pointers */ 2094 pc->ops->apply = PCApply_BDDC; 2095 pc->ops->applytranspose = PCApplyTranspose_BDDC; 2096 pc->ops->setup = PCSetUp_BDDC; 2097 pc->ops->destroy = PCDestroy_BDDC; 2098 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 2099 pc->ops->view = PCView_BDDC; 2100 pc->ops->applyrichardson = 0; 2101 pc->ops->applysymmetricleft = 0; 2102 pc->ops->applysymmetricright = 0; 2103 pc->ops->presolve = PCPreSolve_BDDC; 2104 pc->ops->postsolve = PCPostSolve_BDDC; 2105 2106 /* composing function */ 2107 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetChangeOfBasisMat_C",PCBDDCSetChangeOfBasisMat_BDDC);CHKERRQ(ierr); 2108 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 2109 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 2110 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevel_C",PCBDDCSetLevel_BDDC);CHKERRQ(ierr); 2111 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetUseExactDirichlet_C",PCBDDCSetUseExactDirichlet_BDDC);CHKERRQ(ierr); 2112 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLevels_C",PCBDDCSetLevels_BDDC);CHKERRQ(ierr); 2113 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",PCBDDCSetNullSpace_BDDC);CHKERRQ(ierr); 2114 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2115 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundariesLocal_C",PCBDDCSetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2116 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2117 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundariesLocal_C",PCBDDCSetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2118 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 2119 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundariesLocal_C",PCBDDCGetDirichletBoundariesLocal_BDDC);CHKERRQ(ierr); 2120 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 2121 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundariesLocal_C",PCBDDCGetNeumannBoundariesLocal_BDDC);CHKERRQ(ierr); 2122 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 2123 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplittingLocal_C",PCBDDCSetDofsSplittingLocal_BDDC);CHKERRQ(ierr); 2124 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 2125 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 2126 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 2127 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 2128 PetscFunctionReturn(0); 2129 } 2130 2131