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