1 /* TODOLIST 2 DofSplitting and DM attached to pc? 3 Change SetNeumannBoundaries to SetNeumannBoundariesLocal and provide new SetNeumannBoundaries (same Dirichlet) 4 change how to deal with the coarse problem (PCBDDCSetCoarseEnvironment): 5 - simplify coarse problem structure -> PCBDDC or PCREDUDANT, nothing else -> same comm for all levels? 6 - remove coarse enums and allow use of PCBDDCGetCoarseKSP 7 - remove metis dependency -> use MatPartitioning for multilevel -> Assemble serial adjacency in PCBDDCAnalyzeInterface? 8 code refactoring: 9 - pick up better names for static functions 10 change options structure: 11 - insert BDDC into MG framework? 12 provide other ops? Ask to developers 13 remove all unused printf 14 man pages 15 */ 16 17 /* ---------------------------------------------------------------------------------------------------------------------------------------------- 18 Implementation of BDDC preconditioner based on: 19 C. Dohrmann "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 20 ---------------------------------------------------------------------------------------------------------------------------------------------- */ 21 22 #include "bddc.h" /*I "petscpc.h" I*/ /* includes for fortran wrappers */ 23 #include "bddcprivate.h" 24 #include <petscblaslapack.h> 25 26 /* prototypes for static functions contained in bddc.c */ 27 static PetscErrorCode PCBDDCSetUseExactDirichlet(PC,PetscBool); 28 static PetscErrorCode PCBDDCSetLevel(PC,PetscInt); 29 static PetscErrorCode PCBDDCCoarseSetUp(PC); 30 static PetscErrorCode PCBDDCSetUpCoarseEnvironment(PC,PetscScalar*); 31 32 /* -------------------------------------------------------------------------- */ 33 #undef __FUNCT__ 34 #define __FUNCT__ "PCSetFromOptions_BDDC" 35 PetscErrorCode PCSetFromOptions_BDDC(PC pc) 36 { 37 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 ierr = PetscOptionsHead("BDDC options");CHKERRQ(ierr); 42 /* Verbose debugging of main data structures */ 43 ierr = PetscOptionsInt("-pc_bddc_check_level" ,"Verbose (debugging) output for PCBDDC" ,"none",pcbddc->dbg_flag ,&pcbddc->dbg_flag ,NULL);CHKERRQ(ierr); 44 /* Some customization for default primal space */ 45 ierr = PetscOptionsBool("-pc_bddc_vertices_only" ,"Use only vertices in coarse space (i.e. discard constraints)","none",pcbddc->vertices_flag ,&pcbddc->vertices_flag ,NULL);CHKERRQ(ierr); 46 ierr = PetscOptionsBool("-pc_bddc_constraints_only","Use only constraints in coarse space (i.e. discard vertices)","none",pcbddc->constraints_flag,&pcbddc->constraints_flag,NULL);CHKERRQ(ierr); 47 ierr = PetscOptionsBool("-pc_bddc_faces_only" ,"Use only faces among constraints of coarse space (i.e. discard edges)" ,"none",pcbddc->faces_flag ,&pcbddc->faces_flag ,NULL);CHKERRQ(ierr); 48 ierr = PetscOptionsBool("-pc_bddc_edges_only" ,"Use only edges among constraints of coarse space (i.e. discard faces)" ,"none",pcbddc->edges_flag ,&pcbddc->edges_flag ,NULL);CHKERRQ(ierr); 49 /* Coarse solver context */ 50 static const char * const avail_coarse_problems[] = {"sequential","replicated","parallel","multilevel","CoarseProblemType","PC_BDDC_",0}; /*order of choiches depends on ENUM defined in bddc.h */ 51 ierr = PetscOptionsEnum("-pc_bddc_coarse_problem_type","Set coarse problem type","none",avail_coarse_problems,(PetscEnum)pcbddc->coarse_problem_type,(PetscEnum*)&pcbddc->coarse_problem_type,NULL);CHKERRQ(ierr); 52 /* Two different application of BDDC to the whole set of dofs, internal and interface */ 53 ierr = PetscOptionsBool("-pc_bddc_switch_preconditioning_type","Switch between M_2 (default) and M_3 preconditioners (as defined by Dohrmann)","none",pcbddc->inexact_prec_type,&pcbddc->inexact_prec_type,NULL);CHKERRQ(ierr); 54 ierr = PetscOptionsBool("-pc_bddc_use_change_of_basis","Use change of basis approach for primal space","none",pcbddc->use_change_of_basis,&pcbddc->use_change_of_basis,NULL);CHKERRQ(ierr); 55 ierr = PetscOptionsBool("-pc_bddc_use_change_on_faces","Use change of basis approach for face constraints","none",pcbddc->use_change_on_faces,&pcbddc->use_change_on_faces,NULL);CHKERRQ(ierr); 56 if (!pcbddc->use_change_of_basis) { 57 pcbddc->use_change_on_faces = PETSC_FALSE; 58 } 59 ierr = PetscOptionsInt("-pc_bddc_coarsening_ratio","Set coarsening ratio used in multilevel coarsening","none",pcbddc->coarsening_ratio,&pcbddc->coarsening_ratio,NULL);CHKERRQ(ierr); 60 ierr = PetscOptionsInt("-pc_bddc_max_levels","Set maximum number of levels for multilevel","none",pcbddc->max_levels,&pcbddc->max_levels,NULL);CHKERRQ(ierr); 61 ierr = PetscOptionsBool("-pc_bddc_use_deluxe_scaling","Use deluxe scaling for BDDC","none",pcbddc->use_deluxe_scaling,&pcbddc->use_deluxe_scaling,NULL);CHKERRQ(ierr); 62 ierr = PetscOptionsTail();CHKERRQ(ierr); 63 PetscFunctionReturn(0); 64 } 65 /* -------------------------------------------------------------------------- */ 66 #undef __FUNCT__ 67 #define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS_BDDC" 68 static PetscErrorCode PCBDDCSetPrimalVerticesLocalIS_BDDC(PC pc, IS PrimalVertices) 69 { 70 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 71 PetscErrorCode ierr; 72 73 PetscFunctionBegin; 74 ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr); 75 ierr = PetscObjectReference((PetscObject)PrimalVertices);CHKERRQ(ierr); 76 pcbddc->user_primal_vertices = PrimalVertices; 77 PetscFunctionReturn(0); 78 } 79 #undef __FUNCT__ 80 #define __FUNCT__ "PCBDDCSetPrimalVerticesLocalIS" 81 /*@ 82 PCBDDCSetPrimalVerticesLocalIS - Set user defined primal vertices in PCBDDC. 83 84 Not collective 85 86 Input Parameters: 87 + pc - the preconditioning context 88 - PrimalVertices - index sets of primal vertices in local numbering 89 90 Level: intermediate 91 92 Notes: 93 94 .seealso: PCBDDC 95 @*/ 96 PetscErrorCode PCBDDCSetPrimalVerticesLocalIS(PC pc, IS PrimalVertices) 97 { 98 PetscErrorCode ierr; 99 100 PetscFunctionBegin; 101 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 102 PetscValidHeaderSpecific(PrimalVertices,IS_CLASSID,2); 103 ierr = PetscTryMethod(pc,"PCBDDCSetPrimalVerticesLocalIS_C",(PC,IS),(pc,PrimalVertices));CHKERRQ(ierr); 104 PetscFunctionReturn(0); 105 } 106 /* -------------------------------------------------------------------------- */ 107 #undef __FUNCT__ 108 #define __FUNCT__ "PCBDDCSetCoarseProblemType_BDDC" 109 static PetscErrorCode PCBDDCSetCoarseProblemType_BDDC(PC pc, CoarseProblemType CPT) 110 { 111 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 112 113 PetscFunctionBegin; 114 pcbddc->coarse_problem_type = CPT; 115 PetscFunctionReturn(0); 116 } 117 118 #undef __FUNCT__ 119 #define __FUNCT__ "PCBDDCSetCoarseProblemType" 120 /*@ 121 PCBDDCSetCoarseProblemType - Set coarse problem type in PCBDDC. 122 123 Not collective 124 125 Input Parameters: 126 + pc - the preconditioning context 127 - CoarseProblemType - pick a better name and explain what this is 128 129 Level: intermediate 130 131 Notes: 132 Not collective but all procs must call with same arguments. 133 134 .seealso: PCBDDC 135 @*/ 136 PetscErrorCode PCBDDCSetCoarseProblemType(PC pc, CoarseProblemType CPT) 137 { 138 PetscErrorCode ierr; 139 140 PetscFunctionBegin; 141 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 142 ierr = PetscTryMethod(pc,"PCBDDCSetCoarseProblemType_C",(PC,CoarseProblemType),(pc,CPT));CHKERRQ(ierr); 143 PetscFunctionReturn(0); 144 } 145 /* -------------------------------------------------------------------------- */ 146 #undef __FUNCT__ 147 #define __FUNCT__ "PCBDDCSetCoarseningRatio_BDDC" 148 static PetscErrorCode PCBDDCSetCoarseningRatio_BDDC(PC pc,PetscInt k) 149 { 150 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 151 152 PetscFunctionBegin; 153 pcbddc->coarsening_ratio=k; 154 PetscFunctionReturn(0); 155 } 156 157 #undef __FUNCT__ 158 #define __FUNCT__ "PCBDDCSetCoarseningRatio" 159 /*@ 160 PCBDDCSetCoarseningRatio - Set coarsening ratio used in multilevel coarsening 161 162 Logically collective on PC 163 164 Input Parameters: 165 + pc - the preconditioning context 166 - k - coarsening ratio 167 168 Approximatively k subdomains at the finer level will be aggregated into a single subdomain at the coarser level. 169 170 Level: intermediate 171 172 Notes: 173 174 .seealso: PCBDDC 175 @*/ 176 PetscErrorCode PCBDDCSetCoarseningRatio(PC pc,PetscInt k) 177 { 178 PetscErrorCode ierr; 179 180 PetscFunctionBegin; 181 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 182 ierr = PetscTryMethod(pc,"PCBDDCSetCoarseningRatio_C",(PC,PetscInt),(pc,k));CHKERRQ(ierr); 183 PetscFunctionReturn(0); 184 } 185 /* -------------------------------------------------------------------------- */ 186 187 #undef __FUNCT__ 188 #define __FUNCT__ "PCBDDCSetMaxLevels_BDDC" 189 static PetscErrorCode PCBDDCSetMaxLevels_BDDC(PC pc,PetscInt max_levels) 190 { 191 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 192 193 PetscFunctionBegin; 194 pcbddc->max_levels=max_levels; 195 PetscFunctionReturn(0); 196 } 197 198 #undef __FUNCT__ 199 #define __FUNCT__ "PCBDDCSetMaxLevels" 200 /*@ 201 PCBDDCSetMaxLevels - Sets the maximum number of levels within the multilevel approach. 202 203 Logically collective on PC 204 205 Input Parameters: 206 + pc - the preconditioning context 207 - max_levels - the maximum number of levels 208 209 Default value is 1, i.e. coarse problem will be solved inexactly with one application 210 of PCBDDC preconditioner if the multilevel approach is requested. 211 212 Level: intermediate 213 214 Notes: 215 216 .seealso: PCBDDC 217 @*/ 218 PetscErrorCode PCBDDCSetMaxLevels(PC pc,PetscInt max_levels) 219 { 220 PetscErrorCode ierr; 221 222 PetscFunctionBegin; 223 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 224 ierr = PetscTryMethod(pc,"PCBDDCSetMaxLevels_C",(PC,PetscInt),(pc,max_levels));CHKERRQ(ierr); 225 PetscFunctionReturn(0); 226 } 227 /* -------------------------------------------------------------------------- */ 228 229 #undef __FUNCT__ 230 #define __FUNCT__ "PCBDDCSetNullSpace_BDDC" 231 static PetscErrorCode PCBDDCSetNullSpace_BDDC(PC pc,MatNullSpace NullSpace) 232 { 233 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 234 PetscErrorCode ierr; 235 236 PetscFunctionBegin; 237 ierr = PetscObjectReference((PetscObject)NullSpace);CHKERRQ(ierr); 238 ierr = MatNullSpaceDestroy(&pcbddc->NullSpace);CHKERRQ(ierr); 239 pcbddc->NullSpace=NullSpace; 240 PetscFunctionReturn(0); 241 } 242 243 #undef __FUNCT__ 244 #define __FUNCT__ "PCBDDCSetNullSpace" 245 /*@ 246 PCBDDCSetNullSpace - Set NullSpace of global operator of BDDC preconditioned mat. 247 248 Logically collective on PC and MatNullSpace 249 250 Input Parameters: 251 + pc - the preconditioning context 252 - NullSpace - Null space of the linear operator to be preconditioned. 253 254 Level: intermediate 255 256 Notes: 257 258 .seealso: PCBDDC 259 @*/ 260 PetscErrorCode PCBDDCSetNullSpace(PC pc,MatNullSpace NullSpace) 261 { 262 PetscErrorCode ierr; 263 264 PetscFunctionBegin; 265 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 266 PetscValidHeaderSpecific(NullSpace,MAT_NULLSPACE_CLASSID,2); 267 ierr = PetscTryMethod(pc,"PCBDDCSetNullSpace_C",(PC,MatNullSpace),(pc,NullSpace));CHKERRQ(ierr); 268 PetscFunctionReturn(0); 269 } 270 /* -------------------------------------------------------------------------- */ 271 272 #undef __FUNCT__ 273 #define __FUNCT__ "PCBDDCSetDirichletBoundaries_BDDC" 274 static PetscErrorCode PCBDDCSetDirichletBoundaries_BDDC(PC pc,IS DirichletBoundaries) 275 { 276 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 277 PetscErrorCode ierr; 278 279 PetscFunctionBegin; 280 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 281 ierr = PetscObjectReference((PetscObject)DirichletBoundaries);CHKERRQ(ierr); 282 pcbddc->DirichletBoundaries=DirichletBoundaries; 283 PetscFunctionReturn(0); 284 } 285 286 #undef __FUNCT__ 287 #define __FUNCT__ "PCBDDCSetDirichletBoundaries" 288 /*@ 289 PCBDDCSetDirichletBoundaries - Set index set defining subdomain part (in local ordering) 290 of Dirichlet boundaries for the global problem. 291 292 Not collective 293 294 Input Parameters: 295 + pc - the preconditioning context 296 - DirichletBoundaries - sequential index set defining the subdomain part of Dirichlet boundaries (can be NULL) 297 298 Level: intermediate 299 300 Notes: 301 302 .seealso: PCBDDC 303 @*/ 304 PetscErrorCode PCBDDCSetDirichletBoundaries(PC pc,IS DirichletBoundaries) 305 { 306 PetscErrorCode ierr; 307 308 PetscFunctionBegin; 309 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 310 PetscValidHeaderSpecific(DirichletBoundaries,IS_CLASSID,2); 311 ierr = PetscTryMethod(pc,"PCBDDCSetDirichletBoundaries_C",(PC,IS),(pc,DirichletBoundaries));CHKERRQ(ierr); 312 PetscFunctionReturn(0); 313 } 314 /* -------------------------------------------------------------------------- */ 315 316 #undef __FUNCT__ 317 #define __FUNCT__ "PCBDDCSetNeumannBoundaries_BDDC" 318 static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc,IS NeumannBoundaries) 319 { 320 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 321 PetscErrorCode ierr; 322 323 PetscFunctionBegin; 324 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 325 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 326 pcbddc->NeumannBoundaries=NeumannBoundaries; 327 PetscFunctionReturn(0); 328 } 329 330 #undef __FUNCT__ 331 #define __FUNCT__ "PCBDDCSetNeumannBoundaries" 332 /*@ 333 PCBDDCSetNeumannBoundaries - Set index set defining subdomain part (in local ordering) 334 of Neumann boundaries for the global problem. 335 336 Not collective 337 338 Input Parameters: 339 + pc - the preconditioning context 340 - NeumannBoundaries - sequential index set defining the subdomain part of Neumann boundaries (can be NULL) 341 342 Level: intermediate 343 344 Notes: 345 346 .seealso: PCBDDC 347 @*/ 348 PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc,IS NeumannBoundaries) 349 { 350 PetscErrorCode ierr; 351 352 PetscFunctionBegin; 353 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 354 PetscValidHeaderSpecific(NeumannBoundaries,IS_CLASSID,2); 355 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundaries_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 356 PetscFunctionReturn(0); 357 } 358 /* -------------------------------------------------------------------------- */ 359 360 #undef __FUNCT__ 361 #define __FUNCT__ "PCBDDCGetDirichletBoundaries_BDDC" 362 static PetscErrorCode PCBDDCGetDirichletBoundaries_BDDC(PC pc,IS *DirichletBoundaries) 363 { 364 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 365 366 PetscFunctionBegin; 367 *DirichletBoundaries = pcbddc->DirichletBoundaries; 368 PetscFunctionReturn(0); 369 } 370 371 #undef __FUNCT__ 372 #define __FUNCT__ "PCBDDCGetDirichletBoundaries" 373 /*@ 374 PCBDDCGetDirichletBoundaries - Get index set defining subdomain part (in local ordering) 375 of Dirichlet boundaries for the global problem. 376 377 Not collective 378 379 Input Parameters: 380 + pc - the preconditioning context 381 382 Output Parameters: 383 + DirichletBoundaries - index set defining the subdomain part of Dirichlet boundaries 384 385 Level: intermediate 386 387 Notes: 388 389 .seealso: PCBDDC 390 @*/ 391 PetscErrorCode PCBDDCGetDirichletBoundaries(PC pc,IS *DirichletBoundaries) 392 { 393 PetscErrorCode ierr; 394 395 PetscFunctionBegin; 396 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 397 ierr = PetscUseMethod(pc,"PCBDDCGetDirichletBoundaries_C",(PC,IS*),(pc,DirichletBoundaries));CHKERRQ(ierr); 398 PetscFunctionReturn(0); 399 } 400 /* -------------------------------------------------------------------------- */ 401 402 #undef __FUNCT__ 403 #define __FUNCT__ "PCBDDCGetNeumannBoundaries_BDDC" 404 static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc,IS *NeumannBoundaries) 405 { 406 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 407 408 PetscFunctionBegin; 409 *NeumannBoundaries = pcbddc->NeumannBoundaries; 410 PetscFunctionReturn(0); 411 } 412 413 #undef __FUNCT__ 414 #define __FUNCT__ "PCBDDCGetNeumannBoundaries" 415 /*@ 416 PCBDDCGetNeumannBoundaries - Get index set defining subdomain part (in local ordering) 417 of Neumann boundaries for the global problem. 418 419 Not collective 420 421 Input Parameters: 422 + pc - the preconditioning context 423 424 Output Parameters: 425 + NeumannBoundaries - index set defining the subdomain part of Neumann boundaries 426 427 Level: intermediate 428 429 Notes: 430 431 .seealso: PCBDDC 432 @*/ 433 PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc,IS *NeumannBoundaries) 434 { 435 PetscErrorCode ierr; 436 437 PetscFunctionBegin; 438 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 439 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundaries_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 440 PetscFunctionReturn(0); 441 } 442 /* -------------------------------------------------------------------------- */ 443 444 #undef __FUNCT__ 445 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph_BDDC" 446 static PetscErrorCode PCBDDCSetLocalAdjacencyGraph_BDDC(PC pc, PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 447 { 448 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 449 PCBDDCGraph mat_graph = pcbddc->mat_graph; 450 PetscErrorCode ierr; 451 452 PetscFunctionBegin; 453 /* free old CSR */ 454 ierr = PCBDDCGraphResetCSR(mat_graph);CHKERRQ(ierr); 455 /* TODO: PCBDDCGraphSetAdjacency */ 456 /* get CSR into graph structure */ 457 if (copymode == PETSC_COPY_VALUES) { 458 ierr = PetscMalloc((nvtxs+1)*sizeof(PetscInt),&mat_graph->xadj);CHKERRQ(ierr); 459 ierr = PetscMalloc(xadj[nvtxs]*sizeof(PetscInt),&mat_graph->adjncy);CHKERRQ(ierr); 460 ierr = PetscMemcpy(mat_graph->xadj,xadj,(nvtxs+1)*sizeof(PetscInt));CHKERRQ(ierr); 461 ierr = PetscMemcpy(mat_graph->adjncy,adjncy,xadj[nvtxs]*sizeof(PetscInt));CHKERRQ(ierr); 462 } else if (copymode == PETSC_OWN_POINTER) { 463 mat_graph->xadj = (PetscInt*)xadj; 464 mat_graph->adjncy = (PetscInt*)adjncy; 465 } 466 mat_graph->nvtxs_csr = nvtxs; 467 PetscFunctionReturn(0); 468 } 469 470 #undef __FUNCT__ 471 #define __FUNCT__ "PCBDDCSetLocalAdjacencyGraph" 472 /*@ 473 PCBDDCSetLocalAdjacencyGraph - Set CSR graph of local matrix for use of PCBDDC. 474 475 Not collective 476 477 Input Parameters: 478 + pc - the preconditioning context 479 - nvtxs - number of local vertices of the graph 480 - xadj, adjncy - the CSR graph 481 - copymode - either PETSC_COPY_VALUES or PETSC_OWN_POINTER. In the former case the user must free the array passed in; 482 in the latter case, memory must be obtained with PetscMalloc. 483 484 Level: intermediate 485 486 Notes: 487 488 .seealso: PCBDDC 489 @*/ 490 PetscErrorCode PCBDDCSetLocalAdjacencyGraph(PC pc,PetscInt nvtxs,const PetscInt xadj[],const PetscInt adjncy[], PetscCopyMode copymode) 491 { 492 void (*f)(void) = 0; 493 PetscErrorCode ierr; 494 495 PetscFunctionBegin; 496 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 497 PetscValidIntPointer(xadj,3); 498 PetscValidIntPointer(xadj,4); 499 if (copymode != PETSC_COPY_VALUES && copymode != PETSC_OWN_POINTER) { 500 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unsupported copy mode %d in %s\n",copymode,__FUNCT__); 501 } 502 ierr = PetscTryMethod(pc,"PCBDDCSetLocalAdjacencyGraph_C",(PC,PetscInt,const PetscInt[],const PetscInt[],PetscCopyMode),(pc,nvtxs,xadj,adjncy,copymode));CHKERRQ(ierr); 503 /* free arrays if PCBDDC is not the PC type */ 504 ierr = PetscObjectQueryFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",&f);CHKERRQ(ierr); 505 if (!f && copymode == PETSC_OWN_POINTER) { 506 ierr = PetscFree(xadj);CHKERRQ(ierr); 507 ierr = PetscFree(adjncy);CHKERRQ(ierr); 508 } 509 PetscFunctionReturn(0); 510 } 511 /* -------------------------------------------------------------------------- */ 512 513 #undef __FUNCT__ 514 #define __FUNCT__ "PCBDDCSetDofsSplitting_BDDC" 515 static PetscErrorCode PCBDDCSetDofsSplitting_BDDC(PC pc,PetscInt n_is, IS ISForDofs[]) 516 { 517 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 518 PetscInt i; 519 PetscErrorCode ierr; 520 521 PetscFunctionBegin; 522 /* Destroy ISes if they were already set */ 523 for (i=0;i<pcbddc->n_ISForDofs;i++) { 524 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 525 } 526 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 527 /* allocate space then set */ 528 ierr = PetscMalloc(n_is*sizeof(IS),&pcbddc->ISForDofs);CHKERRQ(ierr); 529 for (i=0;i<n_is;i++) { 530 ierr = PetscObjectReference((PetscObject)ISForDofs[i]);CHKERRQ(ierr); 531 pcbddc->ISForDofs[i]=ISForDofs[i]; 532 } 533 pcbddc->n_ISForDofs=n_is; 534 PetscFunctionReturn(0); 535 } 536 537 #undef __FUNCT__ 538 #define __FUNCT__ "PCBDDCSetDofsSplitting" 539 /*@ 540 PCBDDCSetDofsSplitting - Set index sets defining fields of local mat. 541 542 Not collective 543 544 Input Parameters: 545 + pc - the preconditioning context 546 - n - number of index sets defining the fields 547 - IS[] - array of IS describing the fields 548 549 Level: intermediate 550 551 Notes: 552 553 .seealso: PCBDDC 554 @*/ 555 PetscErrorCode PCBDDCSetDofsSplitting(PC pc,PetscInt n_is, IS ISForDofs[]) 556 { 557 PetscErrorCode ierr; 558 559 PetscFunctionBegin; 560 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 561 ierr = PetscTryMethod(pc,"PCBDDCSetDofsSplitting_C",(PC,PetscInt,IS[]),(pc,n_is,ISForDofs));CHKERRQ(ierr); 562 PetscFunctionReturn(0); 563 } 564 /* -------------------------------------------------------------------------- */ 565 #undef __FUNCT__ 566 #define __FUNCT__ "PCPreSolve_BDDC" 567 /* -------------------------------------------------------------------------- */ 568 /* 569 PCPreSolve_BDDC - Changes the right hand side and (if necessary) the initial 570 guess if a transformation of basis approach has been selected. 571 572 Input Parameter: 573 + pc - the preconditioner contex 574 575 Application Interface Routine: PCPreSolve() 576 577 Notes: 578 The interface routine PCPreSolve() is not usually called directly by 579 the user, but instead is called by KSPSolve(). 580 */ 581 static PetscErrorCode PCPreSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 582 { 583 PetscErrorCode ierr; 584 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 585 PC_IS *pcis = (PC_IS*)(pc->data); 586 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 587 Mat temp_mat; 588 IS dirIS; 589 PetscInt dirsize,i,*is_indices; 590 PetscScalar *array_x,*array_diagonal; 591 Vec used_vec; 592 PetscBool guess_nonzero; 593 594 PetscFunctionBegin; 595 if (x) { 596 ierr = PetscObjectReference((PetscObject)x);CHKERRQ(ierr); 597 used_vec = x; 598 } else { 599 ierr = PetscObjectReference((PetscObject)pcbddc->temp_solution);CHKERRQ(ierr); 600 used_vec = pcbddc->temp_solution; 601 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 602 } 603 /* hack into ksp data structure PCPreSolve comes earlier in src/ksp/ksp/interface/itfunc.c */ 604 if (ksp) { 605 ierr = KSPGetInitialGuessNonzero(ksp,&guess_nonzero);CHKERRQ(ierr); 606 if ( !guess_nonzero ) { 607 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 608 } 609 } 610 /* store the original rhs */ 611 ierr = VecCopy(rhs,pcbddc->original_rhs);CHKERRQ(ierr); 612 613 /* Take into account zeroed rows -> change rhs and store solution removed */ 614 ierr = MatGetDiagonal(pc->pmat,pcis->vec1_global);CHKERRQ(ierr); 615 ierr = VecPointwiseDivide(pcis->vec1_global,rhs,pcis->vec1_global);CHKERRQ(ierr); 616 ierr = VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 617 ierr = VecScatterEnd(matis->ctx,pcis->vec1_global,pcis->vec2_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 618 ierr = VecScatterBegin(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 619 ierr = VecScatterEnd(matis->ctx,used_vec,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 620 ierr = PCBDDCGetDirichletBoundaries(pc,&dirIS);CHKERRQ(ierr); 621 if (dirIS) { 622 ierr = ISGetSize(dirIS,&dirsize);CHKERRQ(ierr); 623 ierr = VecGetArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 624 ierr = VecGetArray(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 625 ierr = ISGetIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 626 for (i=0; i<dirsize; i++) array_x[is_indices[i]] = array_diagonal[is_indices[i]]; 627 ierr = ISRestoreIndices(dirIS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 628 ierr = VecRestoreArray(pcis->vec2_N,&array_diagonal);CHKERRQ(ierr); 629 ierr = VecRestoreArray(pcis->vec1_N,&array_x);CHKERRQ(ierr); 630 } 631 ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 632 ierr = VecScatterEnd(matis->ctx,pcis->vec1_N,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 633 634 /* remove the computed solution from the rhs */ 635 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 636 ierr = MatMultAdd(pc->pmat,used_vec,rhs,rhs);CHKERRQ(ierr); 637 ierr = VecScale(used_vec,-1.0);CHKERRQ(ierr); 638 639 /* store partially computed solution and set initial guess */ 640 if (x) { 641 ierr = VecCopy(used_vec,pcbddc->temp_solution);CHKERRQ(ierr); 642 ierr = VecSet(used_vec,0.0);CHKERRQ(ierr); 643 if (pcbddc->use_exact_dirichlet && !pcbddc->coarse_psi_B) { 644 ierr = VecScatterBegin(pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 645 ierr = VecScatterEnd (pcis->global_to_D,rhs,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 646 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 647 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 648 ierr = VecScatterEnd (pcis->global_to_D,pcis->vec2_D,used_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 649 if (ksp) { 650 ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr); 651 } 652 } 653 } 654 655 /* rhs change of basis */ 656 if (pcbddc->use_change_of_basis) { 657 /* swap pointers for local matrices */ 658 temp_mat = matis->A; 659 matis->A = pcbddc->local_mat; 660 pcbddc->local_mat = temp_mat; 661 /* Get local rhs and apply transformation of basis */ 662 ierr = VecScatterBegin(pcis->global_to_B,rhs,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 663 ierr = VecScatterEnd (pcis->global_to_B,rhs,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 664 /* from original basis to modified basis */ 665 ierr = MatMultTranspose(pcbddc->ChangeOfBasisMatrix,pcis->vec1_B,pcis->vec2_B);CHKERRQ(ierr); 666 /* put back modified values into the global vec using INSERT_VALUES copy mode */ 667 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec2_B,rhs,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 668 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec2_B,rhs,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 669 } 670 if (ksp && pcbddc->NullSpace) { 671 ierr = MatNullSpaceRemove(pcbddc->NullSpace,used_vec);CHKERRQ(ierr); 672 ierr = MatNullSpaceRemove(pcbddc->NullSpace,rhs);CHKERRQ(ierr); 673 } 674 ierr = VecDestroy(&used_vec);CHKERRQ(ierr); 675 PetscFunctionReturn(0); 676 } 677 /* -------------------------------------------------------------------------- */ 678 #undef __FUNCT__ 679 #define __FUNCT__ "PCPostSolve_BDDC" 680 /* -------------------------------------------------------------------------- */ 681 /* 682 PCPostSolve_BDDC - Changes the computed solution if a transformation of basis 683 approach has been selected. Also, restores rhs to its original state. 684 685 Input Parameter: 686 + pc - the preconditioner contex 687 688 Application Interface Routine: PCPostSolve() 689 690 Notes: 691 The interface routine PCPostSolve() is not usually called directly by 692 the user, but instead is called by KSPSolve(). 693 */ 694 static PetscErrorCode PCPostSolve_BDDC(PC pc, KSP ksp, Vec rhs, Vec x) 695 { 696 PetscErrorCode ierr; 697 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 698 PC_IS *pcis = (PC_IS*)(pc->data); 699 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 700 Mat temp_mat; 701 702 PetscFunctionBegin; 703 if (pcbddc->use_change_of_basis && x) { 704 /* swap pointers for local matrices */ 705 temp_mat = matis->A; 706 matis->A = pcbddc->local_mat; 707 pcbddc->local_mat = temp_mat; 708 /* Get Local boundary and apply transformation of basis to solution vector */ 709 ierr = VecScatterBegin(pcis->global_to_B,x,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 710 ierr = VecScatterEnd (pcis->global_to_B,x,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 711 /* from modified basis to original basis */ 712 ierr = MatMult(pcbddc->ChangeOfBasisMatrix,pcis->vec1_B,pcis->vec2_B);CHKERRQ(ierr); 713 /* put back modified values into the global vec using INSERT_VALUES copy mode */ 714 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec2_B,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 715 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec2_B,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 716 } 717 /* add solution removed in presolve */ 718 if (x) { 719 ierr = VecAXPY(x,1.0,pcbddc->temp_solution);CHKERRQ(ierr); 720 } 721 /* restore rhs to its original state */ 722 if (rhs) { 723 ierr = VecCopy(pcbddc->original_rhs,rhs);CHKERRQ(ierr); 724 } 725 PetscFunctionReturn(0); 726 } 727 /* -------------------------------------------------------------------------- */ 728 #undef __FUNCT__ 729 #define __FUNCT__ "PCSetUp_BDDC" 730 /* -------------------------------------------------------------------------- */ 731 /* 732 PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner 733 by setting data structures and options. 734 735 Input Parameter: 736 + pc - the preconditioner context 737 738 Application Interface Routine: PCSetUp() 739 740 Notes: 741 The interface routine PCSetUp() is not usually called directly by 742 the user, but instead is called by PCApply() if necessary. 743 */ 744 PetscErrorCode PCSetUp_BDDC(PC pc) 745 { 746 PetscErrorCode ierr; 747 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 748 MatStructure flag; 749 PetscBool computeis,computetopography,computesolvers; 750 751 PetscFunctionBegin; 752 /* the following lines of code should be replaced by a better logic between PCIS, PCNN, PCBDDC and other nonoverlapping preconditioners */ 753 /* For BDDC we need to define a local "Neumann" problem different to that defined in PCISSetup 754 So, we set to pcnone the Neumann problem of pcis in order to avoid unneeded computation 755 Also, we decide to directly build the (same) Dirichlet problem */ 756 ierr = PetscOptionsSetValue("-is_localN_pc_type","none");CHKERRQ(ierr); 757 ierr = PetscOptionsSetValue("-is_localD_pc_type","none");CHKERRQ(ierr); 758 /* Get stdout for dbg */ 759 if (pcbddc->dbg_flag && !pcbddc->dbg_viewer) { 760 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&pcbddc->dbg_viewer);CHKERRQ(ierr); 761 ierr = PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr); 762 } 763 /* first attempt to split work */ 764 if (pc->setupcalled) { 765 computeis = PETSC_FALSE; 766 ierr = PCGetOperators(pc,NULL,NULL,&flag);CHKERRQ(ierr); 767 if (flag == SAME_PRECONDITIONER) { 768 computetopography = PETSC_FALSE; 769 computesolvers = PETSC_FALSE; 770 } else if (flag == SAME_NONZERO_PATTERN) { 771 computetopography = PETSC_FALSE; 772 computesolvers = PETSC_TRUE; 773 } else { /* DIFFERENT_NONZERO_PATTERN */ 774 computetopography = PETSC_TRUE; 775 computesolvers = PETSC_TRUE; 776 } 777 } else { 778 computeis = PETSC_TRUE; 779 computetopography = PETSC_TRUE; 780 computesolvers = PETSC_TRUE; 781 } 782 /* Set up all the "iterative substructuring" common block */ 783 if (computeis) { 784 ierr = PCISSetUp(pc);CHKERRQ(ierr); 785 } 786 /* Analyze interface and set up local constraint and change of basis matrices */ 787 if (computetopography) { 788 /* reset data */ 789 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 790 ierr = PCBDDCAnalyzeInterface(pc);CHKERRQ(ierr); 791 ierr = PCBDDCConstraintsSetUp(pc);CHKERRQ(ierr); 792 } 793 if (computesolvers) { 794 /* reset data */ 795 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 796 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 797 /* Create coarse and local stuffs used for evaluating action of preconditioner */ 798 ierr = PCBDDCCoarseSetUp(pc);CHKERRQ(ierr); 799 ierr = PCBDDCScalingSetUp(pc);CHKERRQ(ierr); 800 } 801 PetscFunctionReturn(0); 802 } 803 804 /* -------------------------------------------------------------------------- */ 805 /* 806 PCApply_BDDC - Applies the BDDC preconditioner to a vector. 807 808 Input Parameters: 809 . pc - the preconditioner context 810 . r - input vector (global) 811 812 Output Parameter: 813 . z - output vector (global) 814 815 Application Interface Routine: PCApply() 816 */ 817 #undef __FUNCT__ 818 #define __FUNCT__ "PCApply_BDDC" 819 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 820 { 821 PC_IS *pcis = (PC_IS*)(pc->data); 822 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 823 PetscErrorCode ierr; 824 const PetscScalar one = 1.0; 825 const PetscScalar m_one = -1.0; 826 const PetscScalar zero = 0.0; 827 828 /* This code is similar to that provided in nn.c for PCNN 829 NN interface preconditioner changed to BDDC 830 Added support for M_3 preconditioner in the reference article (code is active if pcbddc->inexact_prec_type = PETSC_TRUE) */ 831 832 PetscFunctionBegin; 833 if (!pcbddc->use_exact_dirichlet || pcbddc->coarse_psi_B) { 834 /* First Dirichlet solve */ 835 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 836 ierr = VecScatterEnd (pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 837 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 838 /* 839 Assembling right hand side for BDDC operator 840 - pcis->vec1_D for the Dirichlet part (if needed, i.e. prec_flag=PETSC_TRUE) 841 - pcis->vec1_B the interface part of the global vector z 842 */ 843 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 844 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 845 if (pcbddc->inexact_prec_type) { ierr = MatMultAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 846 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 847 ierr = VecCopy(r,z);CHKERRQ(ierr); 848 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 849 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 850 ierr = PCBDDCScalingRestriction(pc,z,pcis->vec1_B);CHKERRQ(ierr); 851 } else { 852 ierr = VecSet(pcis->vec1_D,zero);CHKERRQ(ierr); 853 ierr = VecSet(pcis->vec2_D,zero);CHKERRQ(ierr); 854 ierr = PCBDDCScalingRestriction(pc,r,pcis->vec1_B);CHKERRQ(ierr); 855 } 856 857 /* Apply interface preconditioner 858 input/output vecs: pcis->vec1_B and pcis->vec1_D */ 859 ierr = PCBDDCApplyInterfacePreconditioner(pc);CHKERRQ(ierr); 860 861 /* Apply transpose of partition of unity operator */ 862 ierr = PCBDDCScalingExtension(pc,pcis->vec1_B,z);CHKERRQ(ierr); 863 864 /* Second Dirichlet solve and assembling of output */ 865 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 866 ierr = VecScatterEnd (pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 867 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 868 if (pcbddc->inexact_prec_type) { ierr = MatMultAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 869 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcbddc->vec4_D);CHKERRQ(ierr); 870 ierr = VecScale(pcbddc->vec4_D,m_one);CHKERRQ(ierr); 871 if (pcbddc->inexact_prec_type) { ierr = VecAXPY (pcbddc->vec4_D,one,pcis->vec1_D);CHKERRQ(ierr); } 872 ierr = VecAXPY (pcis->vec2_D,one,pcbddc->vec4_D);CHKERRQ(ierr); 873 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 874 ierr = VecScatterEnd (pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 875 PetscFunctionReturn(0); 876 } 877 /* -------------------------------------------------------------------------- */ 878 879 #undef __FUNCT__ 880 #define __FUNCT__ "PCDestroy_BDDC" 881 PetscErrorCode PCDestroy_BDDC(PC pc) 882 { 883 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 884 PetscErrorCode ierr; 885 886 PetscFunctionBegin; 887 /* free data created by PCIS */ 888 ierr = PCISDestroy(pc);CHKERRQ(ierr); 889 /* free BDDC custom data */ 890 ierr = PCBDDCResetCustomization(pc);CHKERRQ(ierr); 891 /* destroy objects related to topography */ 892 ierr = PCBDDCResetTopography(pc);CHKERRQ(ierr); 893 /* free allocated graph structure */ 894 ierr = PetscFree(pcbddc->mat_graph);CHKERRQ(ierr); 895 /* free data for scaling operator */ 896 ierr = PCBDDCScalingDestroy(pc);CHKERRQ(ierr); 897 /* free solvers stuff */ 898 ierr = PCBDDCResetSolvers(pc);CHKERRQ(ierr); 899 /* remove functions */ 900 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",NULL);CHKERRQ(ierr); 901 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",NULL);CHKERRQ(ierr); 902 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetMaxLevels_C",NULL);CHKERRQ(ierr); 903 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",NULL);CHKERRQ(ierr); 904 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 905 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 906 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",NULL);CHKERRQ(ierr); 907 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",NULL);CHKERRQ(ierr); 908 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseProblemType_C",NULL);CHKERRQ(ierr); 909 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",NULL);CHKERRQ(ierr); 910 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",NULL);CHKERRQ(ierr); 911 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",NULL);CHKERRQ(ierr); 912 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",NULL);CHKERRQ(ierr); 913 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",NULL);CHKERRQ(ierr); 914 /* Free the private data structure */ 915 ierr = PetscFree(pc->data);CHKERRQ(ierr); 916 PetscFunctionReturn(0); 917 } 918 /* -------------------------------------------------------------------------- */ 919 920 #undef __FUNCT__ 921 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS_BDDC" 922 static PetscErrorCode PCBDDCMatFETIDPGetRHS_BDDC(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 923 { 924 FETIDPMat_ctx mat_ctx; 925 PC_IS* pcis; 926 PC_BDDC* pcbddc; 927 PetscErrorCode ierr; 928 929 PetscFunctionBegin; 930 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 931 pcis = (PC_IS*)mat_ctx->pc->data; 932 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 933 934 /* change of basis for physical rhs if needed 935 It also changes the rhs in case of dirichlet boundaries */ 936 (*mat_ctx->pc->ops->presolve)(mat_ctx->pc,NULL,standard_rhs,NULL); 937 /* store vectors for computation of fetidp final solution */ 938 ierr = VecScatterBegin(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 939 ierr = VecScatterEnd(pcis->global_to_D,standard_rhs,mat_ctx->temp_solution_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 940 /* scale rhs since it should be unassembled */ 941 /* TODO use counter scaling? (also below) */ 942 ierr = VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 943 ierr = VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 944 /* Apply partition of unity */ 945 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 946 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 947 if (!pcbddc->inexact_prec_type) { 948 /* compute partially subassembled Schur complement right-hand side */ 949 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 950 ierr = MatMult(pcis->A_BI,pcis->vec1_D,pcis->vec1_B);CHKERRQ(ierr); 951 ierr = VecAXPY(mat_ctx->temp_solution_B,-1.0,pcis->vec1_B);CHKERRQ(ierr); 952 ierr = VecSet(standard_rhs,0.0);CHKERRQ(ierr); 953 ierr = VecScatterBegin(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 954 ierr = VecScatterEnd(pcis->global_to_B,mat_ctx->temp_solution_B,standard_rhs,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 955 /* ierr = PCBDDCScalingRestriction(mat_ctx->pc,standard_rhs,mat_ctx->temp_solution_B);CHKERRQ(ierr); */ 956 ierr = VecScatterBegin(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 957 ierr = VecScatterEnd(pcis->global_to_B,standard_rhs,mat_ctx->temp_solution_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 958 ierr = VecPointwiseMult(mat_ctx->temp_solution_B,pcis->D,mat_ctx->temp_solution_B);CHKERRQ(ierr); 959 } 960 /* BDDC rhs */ 961 ierr = VecCopy(mat_ctx->temp_solution_B,pcis->vec1_B);CHKERRQ(ierr); 962 if (pcbddc->inexact_prec_type) { 963 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 964 } 965 /* apply BDDC */ 966 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc);CHKERRQ(ierr); 967 /* Application of B_delta and assembling of rhs for fetidp fluxes */ 968 ierr = VecSet(fetidp_flux_rhs,0.0);CHKERRQ(ierr); 969 ierr = MatMult(mat_ctx->B_delta,pcis->vec1_B,mat_ctx->lambda_local);CHKERRQ(ierr); 970 ierr = VecScatterBegin(mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 971 ierr = VecScatterEnd (mat_ctx->l2g_lambda,mat_ctx->lambda_local,fetidp_flux_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 972 /* restore original rhs */ 973 ierr = VecCopy(pcbddc->original_rhs,standard_rhs);CHKERRQ(ierr); 974 PetscFunctionReturn(0); 975 } 976 977 #undef __FUNCT__ 978 #define __FUNCT__ "PCBDDCMatFETIDPGetRHS" 979 /*@ 980 PCBDDCMatFETIDPGetRHS - Get rhs for FETIDP linear system. 981 982 Collective 983 984 Input Parameters: 985 + fetidp_mat - the FETIDP mat obtained by a call to PCBDDCCreateFETIDPOperators 986 + standard_rhs - the rhs of your linear system 987 988 Output Parameters: 989 + fetidp_flux_rhs - the rhs of the FETIDP linear system 990 991 Level: developer 992 993 Notes: 994 995 .seealso: PCBDDC 996 @*/ 997 PetscErrorCode PCBDDCMatFETIDPGetRHS(Mat fetidp_mat, Vec standard_rhs, Vec fetidp_flux_rhs) 998 { 999 FETIDPMat_ctx mat_ctx; 1000 PetscErrorCode ierr; 1001 1002 PetscFunctionBegin; 1003 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1004 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetRHS_C",(Mat,Vec,Vec),(fetidp_mat,standard_rhs,fetidp_flux_rhs));CHKERRQ(ierr); 1005 PetscFunctionReturn(0); 1006 } 1007 /* -------------------------------------------------------------------------- */ 1008 1009 #undef __FUNCT__ 1010 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution_BDDC" 1011 static PetscErrorCode PCBDDCMatFETIDPGetSolution_BDDC(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 1012 { 1013 FETIDPMat_ctx mat_ctx; 1014 PC_IS* pcis; 1015 PC_BDDC* pcbddc; 1016 PetscErrorCode ierr; 1017 1018 PetscFunctionBegin; 1019 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1020 pcis = (PC_IS*)mat_ctx->pc->data; 1021 pcbddc = (PC_BDDC*)mat_ctx->pc->data; 1022 1023 /* apply B_delta^T */ 1024 ierr = VecScatterBegin(mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1025 ierr = VecScatterEnd (mat_ctx->l2g_lambda,fetidp_flux_sol,mat_ctx->lambda_local,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1026 ierr = MatMultTranspose(mat_ctx->B_delta,mat_ctx->lambda_local,pcis->vec1_B);CHKERRQ(ierr); 1027 /* compute rhs for BDDC application */ 1028 ierr = VecAYPX(pcis->vec1_B,-1.0,mat_ctx->temp_solution_B);CHKERRQ(ierr); 1029 if (pcbddc->inexact_prec_type) { 1030 ierr = VecCopy(mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1031 } 1032 /* apply BDDC */ 1033 ierr = PCBDDCApplyInterfacePreconditioner(mat_ctx->pc);CHKERRQ(ierr); 1034 /* put values into standard global vector */ 1035 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1036 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1037 if (!pcbddc->inexact_prec_type) { 1038 /* compute values into the interior if solved for the partially subassembled Schur complement */ 1039 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec1_D);CHKERRQ(ierr); 1040 ierr = VecAXPY(mat_ctx->temp_solution_D,-1.0,pcis->vec1_D);CHKERRQ(ierr); 1041 ierr = KSPSolve(pcbddc->ksp_D,mat_ctx->temp_solution_D,pcis->vec1_D);CHKERRQ(ierr); 1042 } 1043 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1044 ierr = VecScatterEnd (pcis->global_to_D,pcis->vec1_D,standard_sol,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1045 /* final change of basis if needed 1046 Is also sums the dirichlet part removed during RHS assembling */ 1047 (*mat_ctx->pc->ops->postsolve)(mat_ctx->pc,NULL,NULL,standard_sol); 1048 PetscFunctionReturn(0); 1049 1050 } 1051 1052 #undef __FUNCT__ 1053 #define __FUNCT__ "PCBDDCMatFETIDPGetSolution" 1054 /*@ 1055 PCBDDCMatFETIDPGetSolution - Get Solution for FETIDP linear system. 1056 1057 Collective 1058 1059 Input Parameters: 1060 + fetidp_mat - the FETIDP mat obtained by a call to PCBDDCCreateFETIDPOperators 1061 + fetidp_flux_sol - the solution of the FETIDP linear system 1062 1063 Output Parameters: 1064 + standard_sol - the solution on the global domain 1065 1066 Level: developer 1067 1068 Notes: 1069 1070 .seealso: PCBDDC 1071 @*/ 1072 PetscErrorCode PCBDDCMatFETIDPGetSolution(Mat fetidp_mat, Vec fetidp_flux_sol, Vec standard_sol) 1073 { 1074 FETIDPMat_ctx mat_ctx; 1075 PetscErrorCode ierr; 1076 1077 PetscFunctionBegin; 1078 ierr = MatShellGetContext(fetidp_mat,&mat_ctx);CHKERRQ(ierr); 1079 ierr = PetscTryMethod(mat_ctx->pc,"PCBDDCMatFETIDPGetSolution_C",(Mat,Vec,Vec),(fetidp_mat,fetidp_flux_sol,standard_sol));CHKERRQ(ierr); 1080 PetscFunctionReturn(0); 1081 } 1082 /* -------------------------------------------------------------------------- */ 1083 1084 extern PetscErrorCode FETIDPMatMult(Mat,Vec,Vec); 1085 extern PetscErrorCode PCBDDCDestroyFETIDPMat(Mat); 1086 extern PetscErrorCode FETIDPPCApply(PC,Vec,Vec); 1087 extern PetscErrorCode PCBDDCDestroyFETIDPPC(PC); 1088 1089 #undef __FUNCT__ 1090 #define __FUNCT__ "PCBDDCCreateFETIDPOperators_BDDC" 1091 static PetscErrorCode PCBDDCCreateFETIDPOperators_BDDC(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 1092 { 1093 1094 FETIDPMat_ctx fetidpmat_ctx; 1095 Mat newmat; 1096 FETIDPPC_ctx fetidppc_ctx; 1097 PC newpc; 1098 MPI_Comm comm; 1099 PetscErrorCode ierr; 1100 1101 PetscFunctionBegin; 1102 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 1103 /* FETIDP linear matrix */ 1104 ierr = PCBDDCCreateFETIDPMatContext(pc,&fetidpmat_ctx);CHKERRQ(ierr); 1105 ierr = PCBDDCSetupFETIDPMatContext(fetidpmat_ctx);CHKERRQ(ierr); 1106 ierr = MatCreateShell(comm,PETSC_DECIDE,PETSC_DECIDE,fetidpmat_ctx->n_lambda,fetidpmat_ctx->n_lambda,fetidpmat_ctx,&newmat);CHKERRQ(ierr); 1107 ierr = MatShellSetOperation(newmat,MATOP_MULT,(void (*)(void))FETIDPMatMult);CHKERRQ(ierr); 1108 ierr = MatShellSetOperation(newmat,MATOP_DESTROY,(void (*)(void))PCBDDCDestroyFETIDPMat);CHKERRQ(ierr); 1109 ierr = MatSetUp(newmat);CHKERRQ(ierr); 1110 /* FETIDP preconditioner */ 1111 ierr = PCBDDCCreateFETIDPPCContext(pc,&fetidppc_ctx);CHKERRQ(ierr); 1112 ierr = PCBDDCSetupFETIDPPCContext(newmat,fetidppc_ctx);CHKERRQ(ierr); 1113 ierr = PCCreate(comm,&newpc);CHKERRQ(ierr); 1114 ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr); 1115 ierr = PCShellSetContext(newpc,fetidppc_ctx);CHKERRQ(ierr); 1116 ierr = PCShellSetApply(newpc,FETIDPPCApply);CHKERRQ(ierr); 1117 ierr = PCShellSetDestroy(newpc,PCBDDCDestroyFETIDPPC);CHKERRQ(ierr); 1118 ierr = PCSetOperators(newpc,newmat,newmat,SAME_PRECONDITIONER);CHKERRQ(ierr); 1119 ierr = PCSetUp(newpc);CHKERRQ(ierr); 1120 /* return pointers for objects created */ 1121 *fetidp_mat=newmat; 1122 *fetidp_pc=newpc; 1123 PetscFunctionReturn(0); 1124 } 1125 1126 #undef __FUNCT__ 1127 #define __FUNCT__ "PCBDDCCreateFETIDPOperators" 1128 /*@ 1129 PCBDDCCreateFETIDPOperators - Create operators for FETIDP. 1130 1131 Collective 1132 1133 Input Parameters: 1134 + pc - the BDDC preconditioning context (setup must be already called) 1135 1136 Level: developer 1137 1138 Notes: 1139 1140 .seealso: PCBDDC 1141 @*/ 1142 PetscErrorCode PCBDDCCreateFETIDPOperators(PC pc, Mat *fetidp_mat, PC *fetidp_pc) 1143 { 1144 PetscErrorCode ierr; 1145 1146 PetscFunctionBegin; 1147 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1148 if (pc->setupcalled) { 1149 ierr = PetscTryMethod(pc,"PCBDDCCreateFETIDPOperators_C",(PC,Mat*,PC*),(pc,fetidp_mat,fetidp_pc));CHKERRQ(ierr); 1150 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"You must call PCSetup_BDDC() first \n"); 1151 PetscFunctionReturn(0); 1152 } 1153 /* -------------------------------------------------------------------------- */ 1154 /*MC 1155 PCBDDC - Balancing Domain Decomposition by Constraints. 1156 1157 Options Database Keys: 1158 . -pcbddc ??? - 1159 1160 Level: intermediate 1161 1162 Notes: The matrix used with this preconditioner must be of type MATIS 1163 1164 Unlike more 'conventional' interface preconditioners, this iterates over ALL the 1165 degrees of freedom, NOT just those on the interface (this allows the use of approximate solvers 1166 on the subdomains). 1167 1168 Options for the coarse grid preconditioner can be set with - 1169 Options for the Dirichlet subproblem can be set with - 1170 Options for the Neumann subproblem can be set with - 1171 1172 Contributed by Stefano Zampini 1173 1174 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 1175 M*/ 1176 1177 #undef __FUNCT__ 1178 #define __FUNCT__ "PCCreate_BDDC" 1179 PETSC_EXTERN PetscErrorCode PCCreate_BDDC(PC pc) 1180 { 1181 PetscErrorCode ierr; 1182 PC_BDDC *pcbddc; 1183 1184 PetscFunctionBegin; 1185 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 1186 ierr = PetscNewLog(pc,PC_BDDC,&pcbddc);CHKERRQ(ierr); 1187 pc->data = (void*)pcbddc; 1188 1189 /* create PCIS data structure */ 1190 ierr = PCISCreate(pc);CHKERRQ(ierr); 1191 1192 /* BDDC specific */ 1193 pcbddc->user_primal_vertices = 0; 1194 pcbddc->NullSpace = 0; 1195 pcbddc->temp_solution = 0; 1196 pcbddc->original_rhs = 0; 1197 pcbddc->local_mat = 0; 1198 pcbddc->ChangeOfBasisMatrix = 0; 1199 pcbddc->use_change_of_basis = PETSC_TRUE; 1200 pcbddc->use_change_on_faces = PETSC_FALSE; 1201 pcbddc->coarse_vec = 0; 1202 pcbddc->coarse_rhs = 0; 1203 pcbddc->coarse_ksp = 0; 1204 pcbddc->coarse_phi_B = 0; 1205 pcbddc->coarse_phi_D = 0; 1206 pcbddc->coarse_psi_B = 0; 1207 pcbddc->coarse_psi_D = 0; 1208 pcbddc->vec1_P = 0; 1209 pcbddc->vec1_R = 0; 1210 pcbddc->vec2_R = 0; 1211 pcbddc->local_auxmat1 = 0; 1212 pcbddc->local_auxmat2 = 0; 1213 pcbddc->R_to_B = 0; 1214 pcbddc->R_to_D = 0; 1215 pcbddc->ksp_D = 0; 1216 pcbddc->ksp_R = 0; 1217 pcbddc->local_primal_indices = 0; 1218 pcbddc->inexact_prec_type = PETSC_FALSE; 1219 pcbddc->NeumannBoundaries = 0; 1220 pcbddc->ISForDofs = 0; 1221 pcbddc->ConstraintMatrix = 0; 1222 pcbddc->use_nnsp_true = PETSC_FALSE; 1223 pcbddc->local_primal_sizes = 0; 1224 pcbddc->local_primal_displacements = 0; 1225 pcbddc->coarse_loc_to_glob = 0; 1226 pcbddc->dbg_flag = 0; 1227 pcbddc->coarsening_ratio = 8; 1228 pcbddc->use_exact_dirichlet = PETSC_TRUE; 1229 pcbddc->current_level = 0; 1230 pcbddc->max_levels = 1; 1231 pcbddc->replicated_local_primal_indices = 0; 1232 pcbddc->replicated_local_primal_values = 0; 1233 1234 /* create local graph structure */ 1235 ierr = PCBDDCGraphCreate(&pcbddc->mat_graph);CHKERRQ(ierr); 1236 1237 /* scaling */ 1238 pcbddc->use_deluxe_scaling = PETSC_FALSE; 1239 pcbddc->work_scaling = 0; 1240 1241 /* function pointers */ 1242 pc->ops->apply = PCApply_BDDC; 1243 pc->ops->applytranspose = 0; 1244 pc->ops->setup = PCSetUp_BDDC; 1245 pc->ops->destroy = PCDestroy_BDDC; 1246 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 1247 pc->ops->view = 0; 1248 pc->ops->applyrichardson = 0; 1249 pc->ops->applysymmetricleft = 0; 1250 pc->ops->applysymmetricright = 0; 1251 pc->ops->presolve = PCPreSolve_BDDC; 1252 pc->ops->postsolve = PCPostSolve_BDDC; 1253 1254 /* composing function */ 1255 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetPrimalVerticesLocalIS_C",PCBDDCSetPrimalVerticesLocalIS_BDDC);CHKERRQ(ierr); 1256 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseningRatio_C",PCBDDCSetCoarseningRatio_BDDC);CHKERRQ(ierr); 1257 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetMaxLevels_C",PCBDDCSetMaxLevels_BDDC);CHKERRQ(ierr); 1258 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNullSpace_C",PCBDDCSetNullSpace_BDDC);CHKERRQ(ierr); 1259 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDirichletBoundaries_C",PCBDDCSetDirichletBoundaries_BDDC);CHKERRQ(ierr); 1260 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C",PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 1261 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetDirichletBoundaries_C",PCBDDCGetDirichletBoundaries_BDDC);CHKERRQ(ierr); 1262 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C",PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 1263 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetCoarseProblemType_C",PCBDDCSetCoarseProblemType_BDDC);CHKERRQ(ierr); 1264 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetDofsSplitting_C",PCBDDCSetDofsSplitting_BDDC);CHKERRQ(ierr); 1265 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCSetLocalAdjacencyGraph_C",PCBDDCSetLocalAdjacencyGraph_BDDC);CHKERRQ(ierr); 1266 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCCreateFETIDPOperators_C",PCBDDCCreateFETIDPOperators_BDDC);CHKERRQ(ierr); 1267 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetRHS_C",PCBDDCMatFETIDPGetRHS_BDDC);CHKERRQ(ierr); 1268 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCBDDCMatFETIDPGetSolution_C",PCBDDCMatFETIDPGetSolution_BDDC);CHKERRQ(ierr); 1269 PetscFunctionReturn(0); 1270 } 1271 1272 /* -------------------------------------------------------------------------- */ 1273 /* All static functions from now on */ 1274 /* -------------------------------------------------------------------------- */ 1275 1276 #undef __FUNCT__ 1277 #define __FUNCT__ "PCBDDCSetUseExactDirichlet" 1278 static PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool use) 1279 { 1280 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1281 1282 PetscFunctionBegin; 1283 pcbddc->use_exact_dirichlet=use; 1284 PetscFunctionReturn(0); 1285 } 1286 1287 #undef __FUNCT__ 1288 #define __FUNCT__ "PCBDDCSetLevel" 1289 static PetscErrorCode PCBDDCSetLevel(PC pc,PetscInt level) 1290 { 1291 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1292 1293 PetscFunctionBegin; 1294 pcbddc->current_level=level; 1295 PetscFunctionReturn(0); 1296 } 1297 1298 /* -------------------------------------------------------------------------- */ 1299 #undef __FUNCT__ 1300 #define __FUNCT__ "PCBDDCCoarseSetUp" 1301 static PetscErrorCode PCBDDCCoarseSetUp(PC pc) 1302 { 1303 PetscErrorCode ierr; 1304 1305 PC_IS* pcis = (PC_IS*)(pc->data); 1306 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 1307 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 1308 IS is_R_local; 1309 VecType impVecType; 1310 MatType impMatType; 1311 PetscInt n_R=0; 1312 PetscInt n_D=0; 1313 PetscInt n_B=0; 1314 PetscScalar zero=0.0; 1315 PetscScalar one=1.0; 1316 PetscScalar m_one=-1.0; 1317 PetscScalar* array; 1318 PetscScalar *coarse_submat_vals; 1319 PetscInt *idx_R_local; 1320 PetscReal *coarsefunctions_errors,*constraints_errors; 1321 /* auxiliary indices */ 1322 PetscInt i,j,k; 1323 /* for verbose output of bddc */ 1324 PetscViewer viewer=pcbddc->dbg_viewer; 1325 PetscInt dbg_flag=pcbddc->dbg_flag; 1326 /* for counting coarse dofs */ 1327 PetscInt n_vertices,n_constraints; 1328 PetscInt size_of_constraint; 1329 PetscInt *row_cmat_indices; 1330 PetscScalar *row_cmat_values; 1331 PetscInt *vertices; 1332 1333 PetscFunctionBegin; 1334 /* Set Non-overlapping dimensions */ 1335 n_B = pcis->n_B; n_D = pcis->n - n_B; 1336 1337 /* transform local matrices if needed */ 1338 if (pcbddc->use_change_of_basis) { 1339 Mat change_mat_all; 1340 PetscInt *nnz,*is_indices,*temp_indices; 1341 1342 ierr = PetscMalloc(pcis->n*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 1343 ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1344 for (i=0;i<n_D;i++) nnz[is_indices[i]] = 1; 1345 ierr = ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1346 ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1347 k=1; 1348 for (i=0;i<n_B;i++) { 1349 ierr = MatGetRow(pcbddc->ChangeOfBasisMatrix,i,&j,NULL,NULL);CHKERRQ(ierr); 1350 nnz[is_indices[i]]=j; 1351 if (k < j) k = j; 1352 ierr = MatRestoreRow(pcbddc->ChangeOfBasisMatrix,i,&j,NULL,NULL);CHKERRQ(ierr); 1353 } 1354 ierr = ISRestoreIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1355 /* assemble change of basis matrix on the whole set of local dofs */ 1356 ierr = PetscMalloc(k*sizeof(PetscInt),&temp_indices);CHKERRQ(ierr); 1357 ierr = MatCreate(PETSC_COMM_SELF,&change_mat_all);CHKERRQ(ierr); 1358 ierr = MatSetSizes(change_mat_all,pcis->n,pcis->n,pcis->n,pcis->n);CHKERRQ(ierr); 1359 ierr = MatSetType(change_mat_all,MATSEQAIJ);CHKERRQ(ierr); 1360 ierr = MatSeqAIJSetPreallocation(change_mat_all,0,nnz);CHKERRQ(ierr); 1361 ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1362 for (i=0;i<n_D;i++) { 1363 ierr = MatSetValue(change_mat_all,is_indices[i],is_indices[i],1.0,INSERT_VALUES);CHKERRQ(ierr); 1364 } 1365 ierr = ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1366 ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1367 for (i=0;i<n_B;i++) { 1368 ierr = MatGetRow(pcbddc->ChangeOfBasisMatrix,i,&j,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);CHKERRQ(ierr); 1369 for (k=0; k<j; k++) temp_indices[k]=is_indices[row_cmat_indices[k]]; 1370 ierr = MatSetValues(change_mat_all,1,&is_indices[i],j,temp_indices,row_cmat_values,INSERT_VALUES);CHKERRQ(ierr); 1371 ierr = MatRestoreRow(pcbddc->ChangeOfBasisMatrix,i,&j,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);CHKERRQ(ierr); 1372 } 1373 ierr = MatAssemblyBegin(change_mat_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1374 ierr = MatAssemblyEnd(change_mat_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1375 /* TODO: HOW TO WORK WITH BAIJ? PtAP not provided */ 1376 ierr = MatGetBlockSize(matis->A,&i);CHKERRQ(ierr); 1377 if (i==1) { 1378 ierr = MatPtAP(matis->A,change_mat_all,MAT_INITIAL_MATRIX,1.0,&pcbddc->local_mat);CHKERRQ(ierr); 1379 } else { 1380 Mat work_mat; 1381 ierr = MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&work_mat);CHKERRQ(ierr); 1382 ierr = MatPtAP(work_mat,change_mat_all,MAT_INITIAL_MATRIX,1.0,&pcbddc->local_mat);CHKERRQ(ierr); 1383 ierr = MatDestroy(&work_mat);CHKERRQ(ierr); 1384 } 1385 ierr = MatDestroy(&change_mat_all);CHKERRQ(ierr); 1386 ierr = PetscFree(nnz);CHKERRQ(ierr); 1387 ierr = PetscFree(temp_indices);CHKERRQ(ierr); 1388 } else { 1389 /* without change of basis, the local matrix is unchanged */ 1390 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 1391 pcbddc->local_mat = matis->A; 1392 } 1393 /* need to rebuild PCIS matrices during SNES or TS -> TODO move this to PCIS code */ 1394 ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr); 1395 ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr); 1396 ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr); 1397 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_IB);CHKERRQ(ierr); 1398 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_BI);CHKERRQ(ierr); 1399 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,MAT_INITIAL_MATRIX,&pcis->A_BB);CHKERRQ(ierr); 1400 /* Change global null space passed in by the user if change of basis has been requested */ 1401 if (pcbddc->NullSpace && pcbddc->use_change_of_basis) { 1402 ierr = PCBDDCNullSpaceAdaptGlobal(pc);CHKERRQ(ierr); 1403 } 1404 1405 /* Set types for local objects needed by BDDC precondtioner */ 1406 impMatType = MATSEQDENSE; 1407 impVecType = VECSEQ; 1408 /* get vertex indices from constraint matrix */ 1409 ierr = PCBDDCGetPrimalVerticesLocalIdx(pc,&n_vertices,&vertices);CHKERRQ(ierr); 1410 /* Set number of constraints */ 1411 n_constraints = pcbddc->local_primal_size-n_vertices; 1412 /* Dohrmann's notation: dofs splitted in R (Remaining: all dofs but the vertices) and V (Vertices) */ 1413 ierr = VecSet(pcis->vec1_N,one);CHKERRQ(ierr); 1414 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1415 for (i=0;i<n_vertices;i++) array[vertices[i]] = zero; 1416 ierr = PetscMalloc((pcis->n-n_vertices)*sizeof(PetscInt),&idx_R_local);CHKERRQ(ierr); 1417 for (i=0, n_R=0; i<pcis->n; i++) { 1418 if (array[i] == one) { 1419 idx_R_local[n_R] = i; 1420 n_R++; 1421 } 1422 } 1423 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1424 ierr = PetscFree(vertices);CHKERRQ(ierr); 1425 if (dbg_flag) { 1426 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 1427 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1428 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d local dimensions\n",PetscGlobalRank);CHKERRQ(ierr); 1429 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local_size = %d, dirichlet_size = %d, boundary_size = %d\n",pcis->n,n_D,n_B);CHKERRQ(ierr); 1430 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"r_size = %d, v_size = %d, constraints = %d, local_primal_size = %d\n",n_R,n_vertices,n_constraints,pcbddc->local_primal_size);CHKERRQ(ierr); 1431 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"pcbddc->n_vertices = %d, pcbddc->n_constraints = %d\n",pcbddc->n_vertices,pcbddc->n_constraints);CHKERRQ(ierr); 1432 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1433 } 1434 1435 /* Allocate needed vectors */ 1436 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 1437 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);CHKERRQ(ierr); 1438 ierr = VecDuplicate(pcis->vec1_D,&pcbddc->vec4_D);CHKERRQ(ierr); 1439 ierr = VecCreate(PETSC_COMM_SELF,&pcbddc->vec1_R);CHKERRQ(ierr); 1440 ierr = VecSetSizes(pcbddc->vec1_R,n_R,n_R);CHKERRQ(ierr); 1441 ierr = VecSetType(pcbddc->vec1_R,impVecType);CHKERRQ(ierr); 1442 ierr = VecDuplicate(pcbddc->vec1_R,&pcbddc->vec2_R);CHKERRQ(ierr); 1443 ierr = VecCreate(PETSC_COMM_SELF,&pcbddc->vec1_P);CHKERRQ(ierr); 1444 ierr = VecSetSizes(pcbddc->vec1_P,pcbddc->local_primal_size,pcbddc->local_primal_size);CHKERRQ(ierr); 1445 ierr = VecSetType(pcbddc->vec1_P,impVecType);CHKERRQ(ierr); 1446 1447 /* Creating some index sets needed */ 1448 /* For submatrices */ 1449 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_R,idx_R_local,PETSC_OWN_POINTER,&is_R_local);CHKERRQ(ierr); 1450 1451 /* For VecScatters pcbddc->R_to_B and (optionally) pcbddc->R_to_D */ 1452 { 1453 IS is_aux1,is_aux2; 1454 PetscInt *aux_array1; 1455 PetscInt *aux_array2; 1456 PetscInt *idx_I_local; 1457 1458 ierr = PetscMalloc((pcis->n_B-n_vertices)*sizeof(PetscInt),&aux_array1);CHKERRQ(ierr); 1459 ierr = PetscMalloc((pcis->n_B-n_vertices)*sizeof(PetscInt),&aux_array2);CHKERRQ(ierr); 1460 1461 ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&idx_I_local);CHKERRQ(ierr); 1462 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1463 for (i=0; i<n_D; i++) array[idx_I_local[i]] = 0; 1464 ierr = ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&idx_I_local);CHKERRQ(ierr); 1465 for (i=0, j=0; i<n_R; i++) { 1466 if (array[idx_R_local[i]] == one) { 1467 aux_array1[j] = i; 1468 j++; 1469 } 1470 } 1471 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1472 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,aux_array1,PETSC_COPY_VALUES,&is_aux1);CHKERRQ(ierr); 1473 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1474 ierr = VecScatterEnd (pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1475 ierr = VecGetArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1476 for (i=0, j=0; i<n_B; i++) { 1477 if (array[i] == one) { 1478 aux_array2[j] = i; j++; 1479 } 1480 } 1481 ierr = VecRestoreArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1482 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,aux_array2,PETSC_COPY_VALUES,&is_aux2);CHKERRQ(ierr); 1483 ierr = VecScatterCreate(pcbddc->vec1_R,is_aux1,pcis->vec1_B,is_aux2,&pcbddc->R_to_B);CHKERRQ(ierr); 1484 ierr = PetscFree(aux_array1);CHKERRQ(ierr); 1485 ierr = PetscFree(aux_array2);CHKERRQ(ierr); 1486 ierr = ISDestroy(&is_aux1);CHKERRQ(ierr); 1487 ierr = ISDestroy(&is_aux2);CHKERRQ(ierr); 1488 1489 if (pcbddc->inexact_prec_type || dbg_flag ) { 1490 ierr = PetscMalloc(n_D*sizeof(PetscInt),&aux_array1);CHKERRQ(ierr); 1491 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1492 for (i=0, j=0; i<n_R; i++) { 1493 if (array[idx_R_local[i]] == zero) { 1494 aux_array1[j] = i; 1495 j++; 1496 } 1497 } 1498 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1499 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,aux_array1,PETSC_COPY_VALUES,&is_aux1);CHKERRQ(ierr); 1500 ierr = VecScatterCreate(pcbddc->vec1_R,is_aux1,pcis->vec1_D,(IS)0,&pcbddc->R_to_D);CHKERRQ(ierr); 1501 ierr = PetscFree(aux_array1);CHKERRQ(ierr); 1502 ierr = ISDestroy(&is_aux1);CHKERRQ(ierr); 1503 } 1504 } 1505 1506 /* Creating PC contexts for local Dirichlet and Neumann problems */ 1507 { 1508 Mat A_RR; 1509 PC pc_temp; 1510 MatStructure matstruct; 1511 /* Matrix for Dirichlet problem is A_II */ 1512 /* HACK (TODO) A_II can be changed between nonlinear iterations */ 1513 if (pc->setupcalled) { /* we dont need to rebuild dirichlet problem the first time we build BDDC */ 1514 ierr = PCGetOperators(pc,NULL,NULL,&matstruct);CHKERRQ(ierr); 1515 if (matstruct == SAME_NONZERO_PATTERN) { 1516 ierr = MatGetSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,MAT_REUSE_MATRIX,&pcis->A_II);CHKERRQ(ierr); 1517 } else { 1518 ierr = MatDestroy(&pcis->A_II);CHKERRQ(ierr); 1519 ierr = MatGetSubMatrix(matis->A,pcis->is_I_local,pcis->is_I_local,MAT_INITIAL_MATRIX,&pcis->A_II);CHKERRQ(ierr); 1520 } 1521 } 1522 ierr = KSPCreate(PETSC_COMM_SELF,&pcbddc->ksp_D);CHKERRQ(ierr); 1523 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_D,(PetscObject)pc,1);CHKERRQ(ierr); 1524 ierr = KSPSetOperators(pcbddc->ksp_D,pcis->A_II,pcis->A_II,SAME_PRECONDITIONER);CHKERRQ(ierr); 1525 ierr = KSPSetType(pcbddc->ksp_D,KSPPREONLY);CHKERRQ(ierr); 1526 ierr = KSPSetOptionsPrefix(pcbddc->ksp_D,"dirichlet_");CHKERRQ(ierr); 1527 /* default */ 1528 ierr = KSPGetPC(pcbddc->ksp_D,&pc_temp);CHKERRQ(ierr); 1529 ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr); 1530 /* Allow user's customization */ 1531 ierr = KSPSetFromOptions(pcbddc->ksp_D);CHKERRQ(ierr); 1532 /* umfpack interface has a bug when matrix dimension is zero */ 1533 if (!n_D) { 1534 ierr = PCSetType(pc_temp,PCNONE);CHKERRQ(ierr); 1535 } 1536 /* Set Up KSP for Dirichlet problem of BDDC */ 1537 ierr = KSPSetUp(pcbddc->ksp_D);CHKERRQ(ierr); 1538 /* set ksp_D into pcis data */ 1539 ierr = KSPDestroy(&pcis->ksp_D);CHKERRQ(ierr); 1540 ierr = PetscObjectReference((PetscObject)pcbddc->ksp_D);CHKERRQ(ierr); 1541 pcis->ksp_D = pcbddc->ksp_D; 1542 /* Matrix for Neumann problem is A_RR -> we need to create it */ 1543 ierr = MatGetSubMatrix(pcbddc->local_mat,is_R_local,is_R_local,MAT_INITIAL_MATRIX,&A_RR);CHKERRQ(ierr); 1544 ierr = KSPCreate(PETSC_COMM_SELF,&pcbddc->ksp_R);CHKERRQ(ierr); 1545 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_R,(PetscObject)pc,1);CHKERRQ(ierr); 1546 ierr = KSPSetOperators(pcbddc->ksp_R,A_RR,A_RR,SAME_PRECONDITIONER);CHKERRQ(ierr); 1547 ierr = KSPSetType(pcbddc->ksp_R,KSPPREONLY);CHKERRQ(ierr); 1548 ierr = KSPSetOptionsPrefix(pcbddc->ksp_R,"neumann_");CHKERRQ(ierr); 1549 /* default */ 1550 ierr = KSPGetPC(pcbddc->ksp_R,&pc_temp);CHKERRQ(ierr); 1551 ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr); 1552 /* Allow user's customization */ 1553 ierr = KSPSetFromOptions(pcbddc->ksp_R);CHKERRQ(ierr); 1554 /* umfpack interface has a bug when matrix dimension is zero */ 1555 if (!n_R) { 1556 ierr = PCSetType(pc_temp,PCNONE);CHKERRQ(ierr); 1557 } 1558 /* Set Up KSP for Neumann problem of BDDC */ 1559 ierr = KSPSetUp(pcbddc->ksp_R);CHKERRQ(ierr); 1560 /* check Dirichlet and Neumann solvers and adapt them if a nullspace correction is needed */ 1561 { 1562 Vec temp_vec; 1563 PetscReal value; 1564 PetscInt use_exact,use_exact_reduced; 1565 1566 ierr = VecDuplicate(pcis->vec1_D,&temp_vec);CHKERRQ(ierr); 1567 ierr = VecSetRandom(pcis->vec1_D,NULL);CHKERRQ(ierr); 1568 ierr = MatMult(pcis->A_II,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 1569 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec2_D,temp_vec);CHKERRQ(ierr); 1570 ierr = VecAXPY(temp_vec,m_one,pcis->vec1_D);CHKERRQ(ierr); 1571 ierr = VecNorm(temp_vec,NORM_INFINITY,&value);CHKERRQ(ierr); 1572 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 1573 use_exact = 1; 1574 if (PetscAbsReal(value) > 1.e-4) use_exact = 0; 1575 ierr = MPI_Allreduce(&use_exact,&use_exact_reduced,1,MPIU_INT,MPI_LAND,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 1576 ierr = PCBDDCSetUseExactDirichlet(pc,(PetscBool)use_exact_reduced);CHKERRQ(ierr); 1577 if (dbg_flag) { 1578 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1579 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 1580 ierr = PetscViewerASCIIPrintf(viewer,"Checking solution of Dirichlet and Neumann problems\n");CHKERRQ(ierr); 1581 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d infinity error for Dirichlet solve = % 1.14e \n",PetscGlobalRank,value);CHKERRQ(ierr); 1582 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1583 } 1584 if (n_D && pcbddc->NullSpace && !use_exact_reduced && !pcbddc->inexact_prec_type) { 1585 ierr = PCBDDCNullSpaceAssembleCorrection(pc,pcis->is_I_local); 1586 } 1587 ierr = VecDuplicate(pcbddc->vec1_R,&temp_vec);CHKERRQ(ierr); 1588 ierr = VecSetRandom(pcbddc->vec1_R,NULL);CHKERRQ(ierr); 1589 ierr = MatMult(A_RR,pcbddc->vec1_R,pcbddc->vec2_R);CHKERRQ(ierr); 1590 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec2_R,temp_vec);CHKERRQ(ierr); 1591 ierr = VecAXPY(temp_vec,m_one,pcbddc->vec1_R);CHKERRQ(ierr); 1592 ierr = VecNorm(temp_vec,NORM_INFINITY,&value);CHKERRQ(ierr); 1593 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 1594 use_exact = 1; 1595 if (PetscAbsReal(value) > 1.e-4) use_exact = 0; 1596 ierr = MPI_Allreduce(&use_exact,&use_exact_reduced,1,MPIU_INT,MPI_LAND,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 1597 if (dbg_flag) { 1598 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d infinity error for Neumann solve = % 1.14e \n",PetscGlobalRank,value);CHKERRQ(ierr); 1599 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1600 } 1601 if (n_R && pcbddc->NullSpace && !use_exact_reduced) { /* is it the right logic? */ 1602 ierr = PCBDDCNullSpaceAssembleCorrection(pc,is_R_local); 1603 } 1604 } 1605 /* free Neumann problem's matrix */ 1606 ierr = MatDestroy(&A_RR);CHKERRQ(ierr); 1607 } 1608 1609 /* Assemble all remaining stuff needed to apply BDDC */ 1610 { 1611 Mat A_RV,A_VR,A_VV; 1612 Mat M1; 1613 Mat C_CR; 1614 Mat AUXMAT; 1615 Vec vec1_C; 1616 Vec vec2_C; 1617 Vec vec1_V; 1618 Vec vec2_V; 1619 IS is_C_local,is_V_local,is_aux1; 1620 ISLocalToGlobalMapping BtoNmap; 1621 PetscInt *nnz; 1622 PetscInt *idx_V_B; 1623 PetscInt *auxindices; 1624 PetscInt index; 1625 PetscScalar* array2; 1626 MatFactorInfo matinfo; 1627 PetscBool setsym=PETSC_FALSE,issym=PETSC_FALSE; 1628 1629 /* Allocating some extra storage just to be safe */ 1630 ierr = PetscMalloc (pcis->n*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 1631 ierr = PetscMalloc (pcis->n*sizeof(PetscInt),&auxindices);CHKERRQ(ierr); 1632 for (i=0;i<pcis->n;i++) auxindices[i]=i; 1633 1634 ierr = PCBDDCGetPrimalVerticesLocalIdx(pc,&n_vertices,&vertices);CHKERRQ(ierr); 1635 /* vertices in boundary numbering */ 1636 ierr = PetscMalloc(n_vertices*sizeof(PetscInt),&idx_V_B);CHKERRQ(ierr); 1637 ierr = ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&BtoNmap);CHKERRQ(ierr); 1638 ierr = ISGlobalToLocalMappingApply(BtoNmap,IS_GTOLM_DROP,n_vertices,vertices,&i,idx_V_B);CHKERRQ(ierr); 1639 ierr = ISLocalToGlobalMappingDestroy(&BtoNmap);CHKERRQ(ierr); 1640 if (i != n_vertices) { 1641 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Error in boundary numbering for BDDC vertices! %d != %d\n",n_vertices,i); 1642 } 1643 1644 /* some work vectors on vertices and/or constraints */ 1645 if (n_vertices) { 1646 ierr = VecCreate(PETSC_COMM_SELF,&vec1_V);CHKERRQ(ierr); 1647 ierr = VecSetSizes(vec1_V,n_vertices,n_vertices);CHKERRQ(ierr); 1648 ierr = VecSetType(vec1_V,impVecType);CHKERRQ(ierr); 1649 ierr = VecDuplicate(vec1_V,&vec2_V);CHKERRQ(ierr); 1650 } 1651 if (n_constraints) { 1652 ierr = VecCreate(PETSC_COMM_SELF,&vec1_C);CHKERRQ(ierr); 1653 ierr = VecSetSizes(vec1_C,n_constraints,n_constraints);CHKERRQ(ierr); 1654 ierr = VecSetType(vec1_C,impVecType);CHKERRQ(ierr); 1655 ierr = VecDuplicate(vec1_C,&vec2_C);CHKERRQ(ierr); 1656 ierr = VecDuplicate(vec1_C,&pcbddc->vec1_C);CHKERRQ(ierr); 1657 } 1658 /* Precompute stuffs needed for preprocessing and application of BDDC*/ 1659 if (n_constraints) { 1660 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->local_auxmat2);CHKERRQ(ierr); 1661 ierr = MatSetSizes(pcbddc->local_auxmat2,n_R,n_constraints,n_R,n_constraints);CHKERRQ(ierr); 1662 ierr = MatSetType(pcbddc->local_auxmat2,impMatType);CHKERRQ(ierr); 1663 ierr = MatSeqDenseSetPreallocation(pcbddc->local_auxmat2,NULL);CHKERRQ(ierr); 1664 1665 /* Create Constraint matrix on R nodes: C_{CR} */ 1666 ierr = ISCreateStride(PETSC_COMM_SELF,n_constraints,n_vertices,1,&is_C_local);CHKERRQ(ierr); 1667 ierr = MatGetSubMatrix(pcbddc->ConstraintMatrix,is_C_local,is_R_local,MAT_INITIAL_MATRIX,&C_CR);CHKERRQ(ierr); 1668 ierr = ISDestroy(&is_C_local);CHKERRQ(ierr); 1669 1670 /* Assemble local_auxmat2 = - A_{RR}^{-1} C^T_{CR} needed by BDDC application */ 1671 for (i=0;i<n_constraints;i++) { 1672 ierr = VecSet(pcbddc->vec1_R,zero);CHKERRQ(ierr); 1673 /* Get row of constraint matrix in R numbering */ 1674 ierr = VecGetArray(pcbddc->vec1_R,&array);CHKERRQ(ierr); 1675 ierr = MatGetRow(C_CR,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);CHKERRQ(ierr); 1676 for (j=0;j<size_of_constraint;j++) array[row_cmat_indices[j]] = -row_cmat_values[j]; 1677 ierr = MatRestoreRow(C_CR,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);CHKERRQ(ierr); 1678 ierr = VecRestoreArray(pcbddc->vec1_R,&array);CHKERRQ(ierr); 1679 1680 /* Solve for row of constraint matrix in R numbering */ 1681 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec2_R);CHKERRQ(ierr); 1682 1683 /* Set values */ 1684 ierr = VecGetArray(pcbddc->vec2_R,&array);CHKERRQ(ierr); 1685 ierr = MatSetValues(pcbddc->local_auxmat2,n_R,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1686 ierr = VecRestoreArray(pcbddc->vec2_R,&array);CHKERRQ(ierr); 1687 } 1688 ierr = MatAssemblyBegin(pcbddc->local_auxmat2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1689 ierr = MatAssemblyEnd(pcbddc->local_auxmat2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1690 1691 /* Assemble AUXMAT = ( LUFactor )( -C_{CR} A_{RR}^{-1} C^T_{CR} )^{-1} */ 1692 ierr = MatMatMult(C_CR,pcbddc->local_auxmat2,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AUXMAT);CHKERRQ(ierr); 1693 ierr = MatFactorInfoInitialize(&matinfo);CHKERRQ(ierr); 1694 ierr = ISCreateStride(PETSC_COMM_SELF,n_constraints,0,1,&is_aux1);CHKERRQ(ierr); 1695 ierr = MatLUFactor(AUXMAT,is_aux1,is_aux1,&matinfo);CHKERRQ(ierr); 1696 ierr = ISDestroy(&is_aux1);CHKERRQ(ierr); 1697 1698 /* Assemble explicitly M1 = ( C_{CR} A_{RR}^{-1} C^T_{CR} )^{-1} needed in preproc */ 1699 ierr = MatCreate(PETSC_COMM_SELF,&M1);CHKERRQ(ierr); 1700 ierr = MatSetSizes(M1,n_constraints,n_constraints,n_constraints,n_constraints);CHKERRQ(ierr); 1701 ierr = MatSetType(M1,impMatType);CHKERRQ(ierr); 1702 ierr = MatSeqDenseSetPreallocation(M1,NULL);CHKERRQ(ierr); 1703 for (i=0;i<n_constraints;i++) { 1704 ierr = VecSet(vec1_C,zero);CHKERRQ(ierr); 1705 ierr = VecSetValue(vec1_C,i,one,INSERT_VALUES);CHKERRQ(ierr); 1706 ierr = VecAssemblyBegin(vec1_C);CHKERRQ(ierr); 1707 ierr = VecAssemblyEnd(vec1_C);CHKERRQ(ierr); 1708 ierr = MatSolve(AUXMAT,vec1_C,vec2_C);CHKERRQ(ierr); 1709 ierr = VecScale(vec2_C,m_one);CHKERRQ(ierr); 1710 ierr = VecGetArray(vec2_C,&array);CHKERRQ(ierr); 1711 ierr = MatSetValues(M1,n_constraints,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1712 ierr = VecRestoreArray(vec2_C,&array);CHKERRQ(ierr); 1713 } 1714 ierr = MatAssemblyBegin(M1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1715 ierr = MatAssemblyEnd(M1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1716 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 1717 /* Assemble local_auxmat1 = M1*C_{CR} needed by BDDC application in KSP and in preproc */ 1718 ierr = MatMatMult(M1,C_CR,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pcbddc->local_auxmat1);CHKERRQ(ierr); 1719 1720 } 1721 1722 /* Get submatrices from subdomain matrix */ 1723 if (n_vertices) { 1724 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_vertices,vertices,PETSC_COPY_VALUES,&is_V_local);CHKERRQ(ierr); 1725 ierr = MatGetSubMatrix(pcbddc->local_mat,is_R_local,is_V_local,MAT_INITIAL_MATRIX,&A_RV);CHKERRQ(ierr); 1726 ierr = MatGetSubMatrix(pcbddc->local_mat,is_V_local,is_R_local,MAT_INITIAL_MATRIX,&A_VR);CHKERRQ(ierr); 1727 ierr = MatGetSubMatrix(pcbddc->local_mat,is_V_local,is_V_local,MAT_INITIAL_MATRIX,&A_VV);CHKERRQ(ierr); 1728 ierr = ISDestroy(&is_V_local);CHKERRQ(ierr); 1729 } 1730 1731 /* Matrix of coarse basis functions (local) */ 1732 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_phi_B);CHKERRQ(ierr); 1733 ierr = MatSetSizes(pcbddc->coarse_phi_B,n_B,pcbddc->local_primal_size,n_B,pcbddc->local_primal_size);CHKERRQ(ierr); 1734 ierr = MatSetType(pcbddc->coarse_phi_B,impMatType);CHKERRQ(ierr); 1735 ierr = MatSeqDenseSetPreallocation(pcbddc->coarse_phi_B,NULL);CHKERRQ(ierr); 1736 if (pcbddc->inexact_prec_type || dbg_flag ) { 1737 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_phi_D);CHKERRQ(ierr); 1738 ierr = MatSetSizes(pcbddc->coarse_phi_D,n_D,pcbddc->local_primal_size,n_D,pcbddc->local_primal_size);CHKERRQ(ierr); 1739 ierr = MatSetType(pcbddc->coarse_phi_D,impMatType);CHKERRQ(ierr); 1740 ierr = MatSeqDenseSetPreallocation(pcbddc->coarse_phi_D,NULL);CHKERRQ(ierr); 1741 } 1742 1743 if (dbg_flag) { 1744 ierr = PetscMalloc(2*pcbddc->local_primal_size*sizeof(*coarsefunctions_errors),&coarsefunctions_errors);CHKERRQ(ierr); 1745 ierr = PetscMalloc(2*pcbddc->local_primal_size*sizeof(*constraints_errors),&constraints_errors);CHKERRQ(ierr); 1746 } 1747 /* Subdomain contribution (Non-overlapping) to coarse matrix */ 1748 ierr = PetscMalloc ((pcbddc->local_primal_size)*(pcbddc->local_primal_size)*sizeof(PetscScalar),&coarse_submat_vals);CHKERRQ(ierr); 1749 1750 /* We are now ready to evaluate coarse basis functions and subdomain contribution to coarse problem */ 1751 for (i=0;i<n_vertices;i++){ 1752 ierr = VecSet(vec1_V,zero);CHKERRQ(ierr); 1753 ierr = VecSetValue(vec1_V,i,one,INSERT_VALUES);CHKERRQ(ierr); 1754 ierr = VecAssemblyBegin(vec1_V);CHKERRQ(ierr); 1755 ierr = VecAssemblyEnd(vec1_V);CHKERRQ(ierr); 1756 /* solution of saddle point problem */ 1757 ierr = MatMult(A_RV,vec1_V,pcbddc->vec1_R);CHKERRQ(ierr); 1758 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec1_R);CHKERRQ(ierr); 1759 ierr = VecScale(pcbddc->vec1_R,m_one);CHKERRQ(ierr); 1760 if (n_constraints) { 1761 ierr = MatMult(pcbddc->local_auxmat1,pcbddc->vec1_R,vec1_C);CHKERRQ(ierr); 1762 ierr = MatMultAdd(pcbddc->local_auxmat2,vec1_C,pcbddc->vec1_R,pcbddc->vec1_R);CHKERRQ(ierr); 1763 ierr = VecScale(vec1_C,m_one);CHKERRQ(ierr); 1764 } 1765 ierr = MatMult(A_VR,pcbddc->vec1_R,vec2_V);CHKERRQ(ierr); 1766 ierr = MatMultAdd(A_VV,vec1_V,vec2_V,vec2_V);CHKERRQ(ierr); 1767 1768 /* Set values in coarse basis function and subdomain part of coarse_mat */ 1769 /* coarse basis functions */ 1770 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1771 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1772 ierr = VecScatterEnd (pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1773 ierr = VecGetArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1774 ierr = MatSetValues(pcbddc->coarse_phi_B,n_B,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1775 ierr = VecRestoreArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1776 ierr = MatSetValue(pcbddc->coarse_phi_B,idx_V_B[i],i,one,INSERT_VALUES);CHKERRQ(ierr); 1777 if ( pcbddc->inexact_prec_type || dbg_flag ) { 1778 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1779 ierr = VecScatterEnd(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1780 ierr = VecGetArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1781 ierr = MatSetValues(pcbddc->coarse_phi_D,n_D,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1782 ierr = VecRestoreArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1783 } 1784 /* subdomain contribution to coarse matrix */ 1785 ierr = VecGetArray(vec2_V,&array);CHKERRQ(ierr); 1786 for (j=0; j<n_vertices; j++) coarse_submat_vals[i*pcbddc->local_primal_size+j] = array[j]; /* WARNING -> column major ordering */ 1787 ierr = VecRestoreArray(vec2_V,&array);CHKERRQ(ierr); 1788 if (n_constraints) { 1789 ierr = VecGetArray(vec1_C,&array);CHKERRQ(ierr); 1790 for (j=0; j<n_constraints; j++) coarse_submat_vals[i*pcbddc->local_primal_size+j+n_vertices] = array[j]; /* WARNING -> column major ordering */ 1791 ierr = VecRestoreArray(vec1_C,&array);CHKERRQ(ierr); 1792 } 1793 1794 if ( dbg_flag ) { 1795 /* assemble subdomain vector on nodes */ 1796 ierr = VecSet(pcis->vec1_N,zero);CHKERRQ(ierr); 1797 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1798 ierr = VecGetArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1799 for (j=0;j<n_R;j++) array[idx_R_local[j]] = array2[j]; 1800 array[ vertices[i] ] = one; 1801 ierr = VecRestoreArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1802 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1803 /* assemble subdomain vector of lagrange multipliers (i.e. primal nodes) */ 1804 ierr = VecSet(pcbddc->vec1_P,zero);CHKERRQ(ierr); 1805 ierr = VecGetArray(pcbddc->vec1_P,&array2);CHKERRQ(ierr); 1806 ierr = VecGetArray(vec2_V,&array);CHKERRQ(ierr); 1807 for (j=0;j<n_vertices;j++) array2[j]=array[j]; 1808 ierr = VecRestoreArray(vec2_V,&array);CHKERRQ(ierr); 1809 if (n_constraints) { 1810 ierr = VecGetArray(vec1_C,&array);CHKERRQ(ierr); 1811 for (j=0;j<n_constraints;j++) array2[j+n_vertices]=array[j]; 1812 ierr = VecRestoreArray(vec1_C,&array);CHKERRQ(ierr); 1813 } 1814 ierr = VecRestoreArray(pcbddc->vec1_P,&array2);CHKERRQ(ierr); 1815 ierr = VecScale(pcbddc->vec1_P,m_one);CHKERRQ(ierr); 1816 /* check saddle point solution */ 1817 ierr = MatMult(pcbddc->local_mat,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1818 ierr = MatMultTransposeAdd(pcbddc->ConstraintMatrix,pcbddc->vec1_P,pcis->vec2_N,pcis->vec2_N);CHKERRQ(ierr); 1819 ierr = VecNorm(pcis->vec2_N,NORM_INFINITY,&coarsefunctions_errors[i]);CHKERRQ(ierr); 1820 ierr = MatMult(pcbddc->ConstraintMatrix,pcis->vec1_N,pcbddc->vec1_P);CHKERRQ(ierr); 1821 ierr = VecGetArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1822 array[i]=array[i]+m_one; /* shift by the identity matrix */ 1823 ierr = VecRestoreArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1824 ierr = VecNorm(pcbddc->vec1_P,NORM_INFINITY,&constraints_errors[i]);CHKERRQ(ierr); 1825 } 1826 } 1827 1828 for (i=0;i<n_constraints;i++){ 1829 ierr = VecSet(vec2_C,zero);CHKERRQ(ierr); 1830 ierr = VecSetValue(vec2_C,i,m_one,INSERT_VALUES);CHKERRQ(ierr); 1831 ierr = VecAssemblyBegin(vec2_C);CHKERRQ(ierr); 1832 ierr = VecAssemblyEnd(vec2_C);CHKERRQ(ierr); 1833 /* solution of saddle point problem */ 1834 ierr = MatMult(M1,vec2_C,vec1_C);CHKERRQ(ierr); 1835 ierr = MatMult(pcbddc->local_auxmat2,vec1_C,pcbddc->vec1_R);CHKERRQ(ierr); 1836 ierr = VecScale(vec1_C,m_one);CHKERRQ(ierr); 1837 if (n_vertices) { ierr = MatMult(A_VR,pcbddc->vec1_R,vec2_V);CHKERRQ(ierr); } 1838 /* Set values in coarse basis function and subdomain part of coarse_mat */ 1839 /* coarse basis functions */ 1840 index=i+n_vertices; 1841 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1842 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1843 ierr = VecScatterEnd(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1844 ierr = VecGetArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1845 ierr = MatSetValues(pcbddc->coarse_phi_B,n_B,auxindices,1,&index,array,INSERT_VALUES);CHKERRQ(ierr); 1846 ierr = VecRestoreArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1847 if ( pcbddc->inexact_prec_type || dbg_flag ) { 1848 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1849 ierr = VecScatterEnd(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1850 ierr = VecGetArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1851 ierr = MatSetValues(pcbddc->coarse_phi_D,n_D,auxindices,1,&index,array,INSERT_VALUES);CHKERRQ(ierr); 1852 ierr = VecRestoreArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1853 } 1854 /* subdomain contribution to coarse matrix */ 1855 if (n_vertices) { 1856 ierr = VecGetArray(vec2_V,&array);CHKERRQ(ierr); 1857 for (j=0; j<n_vertices; j++) coarse_submat_vals[index*pcbddc->local_primal_size+j]=array[j]; /* WARNING -> column major ordering */ 1858 ierr = VecRestoreArray(vec2_V,&array);CHKERRQ(ierr); 1859 } 1860 ierr = VecGetArray(vec1_C,&array);CHKERRQ(ierr); 1861 for (j=0; j<n_constraints; j++) coarse_submat_vals[index*pcbddc->local_primal_size+j+n_vertices]=array[j]; /* WARNING -> column major ordering */ 1862 ierr = VecRestoreArray(vec1_C,&array);CHKERRQ(ierr); 1863 1864 if ( dbg_flag ) { 1865 /* assemble subdomain vector on nodes */ 1866 ierr = VecSet(pcis->vec1_N,zero);CHKERRQ(ierr); 1867 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1868 ierr = VecGetArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1869 for (j=0;j<n_R;j++) array[idx_R_local[j]] = array2[j]; 1870 ierr = VecRestoreArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1871 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1872 /* assemble subdomain vector of lagrange multipliers */ 1873 ierr = VecSet(pcbddc->vec1_P,zero);CHKERRQ(ierr); 1874 ierr = VecGetArray(pcbddc->vec1_P,&array2);CHKERRQ(ierr); 1875 if ( n_vertices) { 1876 ierr = VecGetArray(vec2_V,&array);CHKERRQ(ierr); 1877 for (j=0;j<n_vertices;j++) array2[j]=-array[j]; 1878 ierr = VecRestoreArray(vec2_V,&array);CHKERRQ(ierr); 1879 } 1880 ierr = VecGetArray(vec1_C,&array);CHKERRQ(ierr); 1881 for (j=0;j<n_constraints;j++) {array2[j+n_vertices]=-array[j];} 1882 ierr = VecRestoreArray(vec1_C,&array);CHKERRQ(ierr); 1883 ierr = VecRestoreArray(pcbddc->vec1_P,&array2);CHKERRQ(ierr); 1884 /* check saddle point solution */ 1885 ierr = MatMult(pcbddc->local_mat,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1886 ierr = MatMultTransposeAdd(pcbddc->ConstraintMatrix,pcbddc->vec1_P,pcis->vec2_N,pcis->vec2_N);CHKERRQ(ierr); 1887 ierr = VecNorm(pcis->vec2_N,NORM_INFINITY,&coarsefunctions_errors[index]);CHKERRQ(ierr); 1888 ierr = MatMult(pcbddc->ConstraintMatrix,pcis->vec1_N,pcbddc->vec1_P);CHKERRQ(ierr); 1889 ierr = VecGetArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1890 array[index]=array[index]+m_one; /* shift by the identity matrix */ 1891 ierr = VecRestoreArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1892 ierr = VecNorm(pcbddc->vec1_P,NORM_INFINITY,&constraints_errors[index]);CHKERRQ(ierr); 1893 } 1894 } 1895 ierr = MatAssemblyBegin(pcbddc->coarse_phi_B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1896 ierr = MatAssemblyEnd(pcbddc->coarse_phi_B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1897 if ( pcbddc->inexact_prec_type || dbg_flag ) { 1898 ierr = MatAssemblyBegin(pcbddc->coarse_phi_D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1899 ierr = MatAssemblyEnd(pcbddc->coarse_phi_D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1900 } 1901 /* compute other basis functions for non-symmetric problems */ 1902 ierr = MatIsSymmetricKnown(pc->pmat,&setsym,&issym);CHKERRQ(ierr); 1903 if ( !setsym || (setsym && !issym) ) { 1904 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_psi_B);CHKERRQ(ierr); 1905 ierr = MatSetSizes(pcbddc->coarse_psi_B,n_B,pcbddc->local_primal_size,n_B,pcbddc->local_primal_size);CHKERRQ(ierr); 1906 ierr = MatSetType(pcbddc->coarse_psi_B,impMatType);CHKERRQ(ierr); 1907 ierr = MatSeqDenseSetPreallocation(pcbddc->coarse_psi_B,NULL);CHKERRQ(ierr); 1908 if (pcbddc->inexact_prec_type || dbg_flag ) { 1909 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_psi_D);CHKERRQ(ierr); 1910 ierr = MatSetSizes(pcbddc->coarse_psi_D,n_D,pcbddc->local_primal_size,n_D,pcbddc->local_primal_size);CHKERRQ(ierr); 1911 ierr = MatSetType(pcbddc->coarse_psi_D,impMatType);CHKERRQ(ierr); 1912 ierr = MatSeqDenseSetPreallocation(pcbddc->coarse_psi_D,NULL);CHKERRQ(ierr); 1913 } 1914 for (i=0;i<pcbddc->local_primal_size;i++) { 1915 if (n_constraints) { 1916 ierr = VecSet(vec1_C,zero);CHKERRQ(ierr); 1917 ierr = VecGetArray(vec1_C,&array);CHKERRQ(ierr); 1918 for (j=0;j<n_constraints;j++) { 1919 array[j]=coarse_submat_vals[(j+n_vertices)*pcbddc->local_primal_size+i]; 1920 } 1921 ierr = VecRestoreArray(vec1_C,&array);CHKERRQ(ierr); 1922 } 1923 if (i<n_vertices) { 1924 ierr = VecSet(vec1_V,zero);CHKERRQ(ierr); 1925 ierr = VecSetValue(vec1_V,i,m_one,INSERT_VALUES);CHKERRQ(ierr); 1926 ierr = VecAssemblyBegin(vec1_V);CHKERRQ(ierr); 1927 ierr = VecAssemblyEnd(vec1_V);CHKERRQ(ierr); 1928 ierr = MatMultTranspose(A_VR,vec1_V,pcbddc->vec1_R);CHKERRQ(ierr); 1929 if (n_constraints) { 1930 ierr = MatMultTransposeAdd(C_CR,vec1_C,pcbddc->vec1_R,pcbddc->vec1_R);CHKERRQ(ierr); 1931 } 1932 } else { 1933 ierr = MatMultTranspose(C_CR,vec1_C,pcbddc->vec1_R);CHKERRQ(ierr); 1934 } 1935 ierr = KSPSolveTranspose(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec1_R);CHKERRQ(ierr); 1936 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1937 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1938 ierr = VecScatterEnd(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1939 ierr = VecGetArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1940 ierr = MatSetValues(pcbddc->coarse_psi_B,n_B,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1941 ierr = VecRestoreArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1942 if (i<n_vertices) { 1943 ierr = MatSetValue(pcbddc->coarse_psi_B,idx_V_B[i],i,one,INSERT_VALUES);CHKERRQ(ierr); 1944 } 1945 if ( pcbddc->inexact_prec_type || dbg_flag ) { 1946 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1947 ierr = VecScatterEnd(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1948 ierr = VecGetArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1949 ierr = MatSetValues(pcbddc->coarse_psi_D,n_D,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1950 ierr = VecRestoreArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1951 } 1952 1953 if ( dbg_flag ) { 1954 /* assemble subdomain vector on nodes */ 1955 ierr = VecSet(pcis->vec1_N,zero);CHKERRQ(ierr); 1956 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1957 ierr = VecGetArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1958 for (j=0;j<n_R;j++) array[idx_R_local[j]] = array2[j]; 1959 if (i<n_vertices) array[vertices[i]] = one; 1960 ierr = VecRestoreArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1961 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1962 /* assemble subdomain vector of lagrange multipliers */ 1963 ierr = VecGetArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1964 for (j=0;j<pcbddc->local_primal_size;j++) { 1965 array[j]=-coarse_submat_vals[j*pcbddc->local_primal_size+i]; 1966 } 1967 ierr = VecRestoreArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1968 /* check saddle point solution */ 1969 ierr = MatMultTranspose(pcbddc->local_mat,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1970 ierr = MatMultTransposeAdd(pcbddc->ConstraintMatrix,pcbddc->vec1_P,pcis->vec2_N,pcis->vec2_N);CHKERRQ(ierr); 1971 ierr = VecNorm(pcis->vec2_N,NORM_INFINITY,&coarsefunctions_errors[i+pcbddc->local_primal_size]);CHKERRQ(ierr); 1972 ierr = MatMult(pcbddc->ConstraintMatrix,pcis->vec1_N,pcbddc->vec1_P);CHKERRQ(ierr); 1973 ierr = VecGetArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1974 array[i]=array[i]+m_one; /* shift by the identity matrix */ 1975 ierr = VecRestoreArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1976 ierr = VecNorm(pcbddc->vec1_P,NORM_INFINITY,&constraints_errors[i+pcbddc->local_primal_size]);CHKERRQ(ierr); 1977 } 1978 } 1979 ierr = MatAssemblyBegin(pcbddc->coarse_psi_B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1980 ierr = MatAssemblyEnd(pcbddc->coarse_psi_B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1981 if ( pcbddc->inexact_prec_type || dbg_flag ) { 1982 ierr = MatAssemblyBegin(pcbddc->coarse_psi_D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1983 ierr = MatAssemblyEnd(pcbddc->coarse_psi_D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1984 } 1985 } 1986 ierr = PetscFree(idx_V_B);CHKERRQ(ierr); 1987 /* Checking coarse_sub_mat and coarse basis functios */ 1988 /* Symmetric case : It should be \Phi^{(j)^T} A^{(j)} \Phi^{(j)}=coarse_sub_mat */ 1989 /* Non-symmetric case : It should be \Psi^{(j)^T} A^{(j)} \Phi^{(j)}=coarse_sub_mat */ 1990 if (dbg_flag) { 1991 Mat coarse_sub_mat; 1992 Mat TM1,TM2,TM3,TM4; 1993 Mat coarse_phi_D,coarse_phi_B; 1994 Mat coarse_psi_D,coarse_psi_B; 1995 Mat A_II,A_BB,A_IB,A_BI; 1996 MatType checkmattype=MATSEQAIJ; 1997 PetscReal real_value; 1998 1999 ierr = MatConvert(pcis->A_II,checkmattype,MAT_INITIAL_MATRIX,&A_II);CHKERRQ(ierr); 2000 ierr = MatConvert(pcis->A_IB,checkmattype,MAT_INITIAL_MATRIX,&A_IB);CHKERRQ(ierr); 2001 ierr = MatConvert(pcis->A_BI,checkmattype,MAT_INITIAL_MATRIX,&A_BI);CHKERRQ(ierr); 2002 ierr = MatConvert(pcis->A_BB,checkmattype,MAT_INITIAL_MATRIX,&A_BB);CHKERRQ(ierr); 2003 ierr = MatConvert(pcbddc->coarse_phi_D,checkmattype,MAT_INITIAL_MATRIX,&coarse_phi_D);CHKERRQ(ierr); 2004 ierr = MatConvert(pcbddc->coarse_phi_B,checkmattype,MAT_INITIAL_MATRIX,&coarse_phi_B);CHKERRQ(ierr); 2005 if (pcbddc->coarse_psi_B) { 2006 ierr = MatConvert(pcbddc->coarse_psi_D,checkmattype,MAT_INITIAL_MATRIX,&coarse_psi_D);CHKERRQ(ierr); 2007 ierr = MatConvert(pcbddc->coarse_psi_B,checkmattype,MAT_INITIAL_MATRIX,&coarse_psi_B);CHKERRQ(ierr); 2008 } 2009 ierr = MatCreateSeqDense(PETSC_COMM_SELF,pcbddc->local_primal_size,pcbddc->local_primal_size,coarse_submat_vals,&coarse_sub_mat);CHKERRQ(ierr); 2010 2011 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 2012 ierr = PetscViewerASCIIPrintf(viewer,"Check coarse sub mat and local basis functions\n");CHKERRQ(ierr); 2013 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2014 if (pcbddc->coarse_psi_B) { 2015 ierr = MatMatMult(A_II,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 2016 ierr = MatTransposeMatMult(coarse_psi_D,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM1);CHKERRQ(ierr); 2017 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 2018 ierr = MatMatMult(A_BB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 2019 ierr = MatTransposeMatMult(coarse_psi_B,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM2);CHKERRQ(ierr); 2020 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 2021 ierr = MatMatMult(A_IB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 2022 ierr = MatTransposeMatMult(coarse_psi_D,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM3);CHKERRQ(ierr); 2023 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 2024 ierr = MatMatMult(A_BI,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 2025 ierr = MatTransposeMatMult(coarse_psi_B,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM4);CHKERRQ(ierr); 2026 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 2027 } else { 2028 ierr = MatPtAP(A_II,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&TM1);CHKERRQ(ierr); 2029 ierr = MatPtAP(A_BB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&TM2);CHKERRQ(ierr); 2030 ierr = MatMatMult(A_IB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 2031 ierr = MatTransposeMatMult(coarse_phi_D,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM3);CHKERRQ(ierr); 2032 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 2033 ierr = MatMatMult(A_BI,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 2034 ierr = MatTransposeMatMult(coarse_phi_B,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM4);CHKERRQ(ierr); 2035 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 2036 } 2037 ierr = MatAXPY(TM1,one,TM2,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 2038 ierr = MatAXPY(TM1,one,TM3,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 2039 ierr = MatAXPY(TM1,one,TM4,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 2040 ierr = MatConvert(TM1,MATSEQDENSE,MAT_REUSE_MATRIX,&TM1);CHKERRQ(ierr); 2041 ierr = MatAXPY(TM1,m_one,coarse_sub_mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2042 ierr = MatNorm(TM1,NORM_INFINITY,&real_value);CHKERRQ(ierr); 2043 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"----------------------------------\n");CHKERRQ(ierr); 2044 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d \n",PetscGlobalRank);CHKERRQ(ierr); 2045 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"matrix error = % 1.14e\n",real_value);CHKERRQ(ierr); 2046 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"coarse functions (phi) errors\n");CHKERRQ(ierr); 2047 for (i=0;i<pcbddc->local_primal_size;i++) { 2048 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local %02d-th function error = % 1.14e\n",i,coarsefunctions_errors[i]);CHKERRQ(ierr); 2049 } 2050 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"constraints (phi) errors\n");CHKERRQ(ierr); 2051 for (i=0;i<pcbddc->local_primal_size;i++) { 2052 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local %02d-th function error = % 1.14e\n",i,constraints_errors[i]);CHKERRQ(ierr); 2053 } 2054 if (pcbddc->coarse_psi_B) { 2055 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"coarse functions (psi) errors\n");CHKERRQ(ierr); 2056 for (i=pcbddc->local_primal_size;i<2*pcbddc->local_primal_size;i++) { 2057 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local %02d-th function error = % 1.14e\n",i-pcbddc->local_primal_size,coarsefunctions_errors[i]);CHKERRQ(ierr); 2058 } 2059 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"constraints (psi) errors\n");CHKERRQ(ierr); 2060 for (i=pcbddc->local_primal_size;i<2*pcbddc->local_primal_size;i++) { 2061 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local %02d-th function error = % 1.14e\n",i-pcbddc->local_primal_size,constraints_errors[i]);CHKERRQ(ierr); 2062 } 2063 } 2064 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2065 ierr = MatDestroy(&A_II);CHKERRQ(ierr); 2066 ierr = MatDestroy(&A_BB);CHKERRQ(ierr); 2067 ierr = MatDestroy(&A_IB);CHKERRQ(ierr); 2068 ierr = MatDestroy(&A_BI);CHKERRQ(ierr); 2069 ierr = MatDestroy(&TM1);CHKERRQ(ierr); 2070 ierr = MatDestroy(&TM2);CHKERRQ(ierr); 2071 ierr = MatDestroy(&TM3);CHKERRQ(ierr); 2072 ierr = MatDestroy(&TM4);CHKERRQ(ierr); 2073 ierr = MatDestroy(&coarse_phi_D);CHKERRQ(ierr); 2074 ierr = MatDestroy(&coarse_phi_B);CHKERRQ(ierr); 2075 if (pcbddc->coarse_psi_B) { 2076 ierr = MatDestroy(&coarse_psi_D);CHKERRQ(ierr); 2077 ierr = MatDestroy(&coarse_psi_B);CHKERRQ(ierr); 2078 } 2079 ierr = MatDestroy(&coarse_sub_mat);CHKERRQ(ierr); 2080 ierr = PetscFree(coarsefunctions_errors);CHKERRQ(ierr); 2081 ierr = PetscFree(constraints_errors);CHKERRQ(ierr); 2082 } 2083 /* free memory */ 2084 if (n_vertices) { 2085 ierr = PetscFree(vertices);CHKERRQ(ierr); 2086 ierr = VecDestroy(&vec1_V);CHKERRQ(ierr); 2087 ierr = VecDestroy(&vec2_V);CHKERRQ(ierr); 2088 ierr = MatDestroy(&A_RV);CHKERRQ(ierr); 2089 ierr = MatDestroy(&A_VR);CHKERRQ(ierr); 2090 ierr = MatDestroy(&A_VV);CHKERRQ(ierr); 2091 } 2092 if (n_constraints) { 2093 ierr = VecDestroy(&vec1_C);CHKERRQ(ierr); 2094 ierr = VecDestroy(&vec2_C);CHKERRQ(ierr); 2095 ierr = MatDestroy(&M1);CHKERRQ(ierr); 2096 ierr = MatDestroy(&C_CR);CHKERRQ(ierr); 2097 } 2098 ierr = PetscFree(auxindices);CHKERRQ(ierr); 2099 ierr = PetscFree(nnz);CHKERRQ(ierr); 2100 /* create coarse matrix and data structures for message passing associated actual choice of coarse problem type */ 2101 ierr = PCBDDCSetUpCoarseEnvironment(pc,coarse_submat_vals);CHKERRQ(ierr); 2102 ierr = PetscFree(coarse_submat_vals);CHKERRQ(ierr); 2103 } 2104 /* free memory */ 2105 ierr = ISDestroy(&is_R_local);CHKERRQ(ierr); 2106 2107 PetscFunctionReturn(0); 2108 } 2109 2110 /* -------------------------------------------------------------------------- */ 2111 2112 /* BDDC requires metis 5.0.1 for multilevel */ 2113 #if defined(PETSC_HAVE_METIS) 2114 #include "metis.h" 2115 #define MetisInt idx_t 2116 #define MetisScalar real_t 2117 #endif 2118 2119 #undef __FUNCT__ 2120 #define __FUNCT__ "PCBDDCSetUpCoarseEnvironment" 2121 static PetscErrorCode PCBDDCSetUpCoarseEnvironment(PC pc,PetscScalar* coarse_submat_vals) 2122 { 2123 2124 2125 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 2126 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 2127 PC_IS *pcis = (PC_IS*)pc->data; 2128 MPI_Comm prec_comm; 2129 MPI_Comm coarse_comm; 2130 2131 MatNullSpace CoarseNullSpace; 2132 2133 /* common to all choiches */ 2134 PetscScalar *temp_coarse_mat_vals; 2135 PetscScalar *ins_coarse_mat_vals; 2136 PetscInt *ins_local_primal_indices; 2137 PetscMPIInt *localsizes2,*localdispl2; 2138 PetscMPIInt size_prec_comm; 2139 PetscMPIInt rank_prec_comm; 2140 PetscMPIInt active_rank=MPI_PROC_NULL; 2141 PetscMPIInt master_proc=0; 2142 PetscInt ins_local_primal_size; 2143 /* specific to MULTILEVEL_BDDC */ 2144 PetscMPIInt *ranks_recv=0; 2145 PetscMPIInt count_recv=0; 2146 PetscMPIInt rank_coarse_proc_send_to=-1; 2147 PetscMPIInt coarse_color = MPI_UNDEFINED; 2148 ISLocalToGlobalMapping coarse_ISLG; 2149 /* some other variables */ 2150 PetscErrorCode ierr; 2151 MatType coarse_mat_type; 2152 PCType coarse_pc_type; 2153 KSPType coarse_ksp_type; 2154 PC pc_temp; 2155 PetscInt i,j,k; 2156 PetscInt max_it_coarse_ksp=1; /* don't increase this value */ 2157 /* verbose output viewer */ 2158 PetscViewer viewer=pcbddc->dbg_viewer; 2159 PetscInt dbg_flag=pcbddc->dbg_flag; 2160 2161 PetscInt offset,offset2; 2162 PetscMPIInt im_active,active_procs; 2163 PetscInt *dnz,*onz; 2164 2165 PetscBool setsym,issym=PETSC_FALSE; 2166 2167 PetscFunctionBegin; 2168 ierr = PetscObjectGetComm((PetscObject)pc,&prec_comm);CHKERRQ(ierr); 2169 ins_local_primal_indices = 0; 2170 ins_coarse_mat_vals = 0; 2171 localsizes2 = 0; 2172 localdispl2 = 0; 2173 temp_coarse_mat_vals = 0; 2174 coarse_ISLG = 0; 2175 2176 ierr = MPI_Comm_size(prec_comm,&size_prec_comm);CHKERRQ(ierr); 2177 ierr = MPI_Comm_rank(prec_comm,&rank_prec_comm);CHKERRQ(ierr); 2178 ierr = MatIsSymmetricKnown(pc->pmat,&setsym,&issym);CHKERRQ(ierr); 2179 2180 /* Assign global numbering to coarse dofs */ 2181 { 2182 PetscInt *auxlocal_primal,*aux_idx; 2183 PetscMPIInt mpi_local_primal_size; 2184 PetscScalar coarsesum,*array; 2185 2186 mpi_local_primal_size = (PetscMPIInt)pcbddc->local_primal_size; 2187 2188 /* Construct needed data structures for message passing */ 2189 j = 0; 2190 if (rank_prec_comm == 0 || pcbddc->coarse_problem_type == REPLICATED_BDDC || pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 2191 j = size_prec_comm; 2192 } 2193 ierr = PetscMalloc(j*sizeof(PetscMPIInt),&pcbddc->local_primal_sizes);CHKERRQ(ierr); 2194 ierr = PetscMalloc(j*sizeof(PetscMPIInt),&pcbddc->local_primal_displacements);CHKERRQ(ierr); 2195 /* Gather local_primal_size information for all processes */ 2196 if (pcbddc->coarse_problem_type == REPLICATED_BDDC || pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 2197 ierr = MPI_Allgather(&mpi_local_primal_size,1,MPIU_INT,&pcbddc->local_primal_sizes[0],1,MPIU_INT,prec_comm);CHKERRQ(ierr); 2198 } else { 2199 ierr = MPI_Gather(&mpi_local_primal_size,1,MPIU_INT,&pcbddc->local_primal_sizes[0],1,MPIU_INT,0,prec_comm);CHKERRQ(ierr); 2200 } 2201 pcbddc->replicated_primal_size = 0; 2202 for (i=0; i<j; i++) { 2203 pcbddc->local_primal_displacements[i] = pcbddc->replicated_primal_size ; 2204 pcbddc->replicated_primal_size += pcbddc->local_primal_sizes[i]; 2205 } 2206 2207 /* First let's count coarse dofs. 2208 This code fragment assumes that the number of local constraints per connected component 2209 is not greater than the number of nodes defined for the connected component 2210 (otherwise we will surely have linear dependence between constraints and thus a singular coarse problem) */ 2211 ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(PetscInt),&auxlocal_primal);CHKERRQ(ierr); 2212 ierr = PCBDDCGetPrimalVerticesLocalIdx(pc,&i,&aux_idx);CHKERRQ(ierr); 2213 ierr = PetscMemcpy(auxlocal_primal,aux_idx,i*sizeof(PetscInt));CHKERRQ(ierr); 2214 ierr = PetscFree(aux_idx);CHKERRQ(ierr); 2215 ierr = PCBDDCGetPrimalConstraintsLocalIdx(pc,&j,&aux_idx);CHKERRQ(ierr); 2216 ierr = PetscMemcpy(&auxlocal_primal[i],aux_idx,j*sizeof(PetscInt));CHKERRQ(ierr); 2217 ierr = PetscFree(aux_idx);CHKERRQ(ierr); 2218 /* Compute number of coarse dofs */ 2219 ierr = PCBDDCSubsetNumbering(prec_comm,matis->mapping,pcbddc->local_primal_size,auxlocal_primal,NULL,&pcbddc->coarse_size,&pcbddc->local_primal_indices);CHKERRQ(ierr); 2220 2221 if (dbg_flag) { 2222 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2223 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 2224 ierr = PetscViewerASCIIPrintf(viewer,"Check coarse indices\n");CHKERRQ(ierr); 2225 ierr = VecSet(pcis->vec1_N,0.0);CHKERRQ(ierr); 2226 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 2227 for (i=0;i<pcbddc->local_primal_size;i++) array[auxlocal_primal[i]]=1.0; 2228 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 2229 ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr); 2230 ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2231 ierr = VecScatterEnd (matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2232 ierr = VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2233 ierr = VecScatterEnd (matis->ctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2234 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 2235 for (i=0;i<pcis->n;i++) { 2236 if (array[i] == 1.0) { 2237 ierr = ISLocalToGlobalMappingApply(matis->mapping,1,&i,&j);CHKERRQ(ierr); 2238 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d: WRONG COARSE INDEX %d (local %d)\n",PetscGlobalRank,j,i);CHKERRQ(ierr); 2239 } 2240 } 2241 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2242 for (i=0;i<pcis->n;i++) { 2243 if (PetscRealPart(array[i]) > 0.0) array[i] = 1.0/PetscRealPart(array[i]); 2244 } 2245 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 2246 ierr = VecSet(pcis->vec1_global,0.0);CHKERRQ(ierr); 2247 ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2248 ierr = VecScatterEnd (matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2249 ierr = VecSum(pcis->vec1_global,&coarsesum);CHKERRQ(ierr); 2250 ierr = PetscViewerASCIIPrintf(viewer,"Size of coarse problem SHOULD be %lf\n",coarsesum);CHKERRQ(ierr); 2251 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2252 } 2253 ierr = PetscFree(auxlocal_primal);CHKERRQ(ierr); 2254 } 2255 2256 if (dbg_flag) { 2257 ierr = PetscViewerASCIIPrintf(viewer,"Size of coarse problem is %d\n",pcbddc->coarse_size);CHKERRQ(ierr); 2258 if (dbg_flag > 1) { 2259 ierr = PetscViewerASCIIPrintf(viewer,"Distribution of local primal indices\n");CHKERRQ(ierr); 2260 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2261 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d\n",PetscGlobalRank);CHKERRQ(ierr); 2262 for (i=0;i<pcbddc->local_primal_size;i++) { 2263 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local_primal_indices[%d]=%d \n",i,pcbddc->local_primal_indices[i]); 2264 } 2265 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2266 } 2267 } 2268 2269 im_active = 0; 2270 if (pcis->n) im_active = 1; 2271 ierr = MPI_Allreduce(&im_active,&active_procs,1,MPIU_INT,MPI_SUM,prec_comm);CHKERRQ(ierr); 2272 2273 /* adapt coarse problem type */ 2274 #if defined(PETSC_HAVE_METIS) 2275 if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 2276 if (pcbddc->current_level < pcbddc->max_levels) { 2277 if ( (active_procs/pcbddc->coarsening_ratio) < 2 ) { 2278 if (dbg_flag) { 2279 ierr = PetscViewerASCIIPrintf(viewer,"Not enough active processes on level %d (active %d,ratio %d). Parallel direct solve for coarse problem\n",pcbddc->current_level,active_procs,pcbddc->coarsening_ratio);CHKERRQ(ierr); 2280 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2281 } 2282 pcbddc->coarse_problem_type = PARALLEL_BDDC; 2283 } 2284 } else { 2285 if (dbg_flag) { 2286 ierr = PetscViewerASCIIPrintf(viewer,"Max number of levels reached. Using parallel direct solve for coarse problem\n",pcbddc->max_levels,active_procs,pcbddc->coarsening_ratio);CHKERRQ(ierr); 2287 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2288 } 2289 pcbddc->coarse_problem_type = PARALLEL_BDDC; 2290 } 2291 } 2292 #else 2293 pcbddc->coarse_problem_type = PARALLEL_BDDC; 2294 #endif 2295 2296 switch(pcbddc->coarse_problem_type){ 2297 2298 case(MULTILEVEL_BDDC): /* we define a coarse mesh where subdomains are elements */ 2299 { 2300 #if defined(PETSC_HAVE_METIS) 2301 /* we need additional variables */ 2302 MetisInt n_subdomains,n_parts,objval,ncon,faces_nvtxs; 2303 MetisInt *metis_coarse_subdivision; 2304 MetisInt options[METIS_NOPTIONS]; 2305 PetscMPIInt size_coarse_comm,rank_coarse_comm; 2306 PetscMPIInt procs_jumps_coarse_comm; 2307 PetscMPIInt *coarse_subdivision; 2308 PetscMPIInt *total_count_recv; 2309 PetscMPIInt *total_ranks_recv; 2310 PetscMPIInt *displacements_recv; 2311 PetscMPIInt *my_faces_connectivity; 2312 PetscMPIInt *petsc_faces_adjncy; 2313 MetisInt *faces_adjncy; 2314 MetisInt *faces_xadj; 2315 PetscMPIInt *number_of_faces; 2316 PetscMPIInt *faces_displacements; 2317 PetscInt *array_int; 2318 PetscMPIInt my_faces=0; 2319 PetscMPIInt total_faces=0; 2320 PetscInt ranks_stretching_ratio; 2321 2322 /* define some quantities */ 2323 pcbddc->coarse_communications_type = SCATTERS_BDDC; 2324 coarse_mat_type = MATIS; 2325 coarse_pc_type = PCBDDC; 2326 coarse_ksp_type = KSPRICHARDSON; 2327 2328 /* details of coarse decomposition */ 2329 n_subdomains = active_procs; 2330 n_parts = n_subdomains/pcbddc->coarsening_ratio; 2331 ranks_stretching_ratio = size_prec_comm/active_procs; 2332 procs_jumps_coarse_comm = pcbddc->coarsening_ratio*ranks_stretching_ratio; 2333 2334 #if 0 2335 PetscMPIInt *old_ranks; 2336 PetscInt *new_ranks,*jj,*ii; 2337 MatPartitioning mat_part; 2338 IS coarse_new_decomposition,is_numbering; 2339 PetscViewer viewer_test; 2340 MPI_Comm test_coarse_comm; 2341 PetscMPIInt test_coarse_color; 2342 Mat mat_adj; 2343 /* Create new communicator for coarse problem splitting the old one */ 2344 /* procs with coarse_color = MPI_UNDEFINED will have coarse_comm = MPI_COMM_NULL (from mpi standards) 2345 key = rank_prec_comm -> keep same ordering of ranks from the old to the new communicator */ 2346 test_coarse_color = ( im_active ? 0 : MPI_UNDEFINED ); 2347 test_coarse_comm = MPI_COMM_NULL; 2348 ierr = MPI_Comm_split(prec_comm,test_coarse_color,rank_prec_comm,&test_coarse_comm);CHKERRQ(ierr); 2349 if (im_active) { 2350 ierr = PetscMalloc(n_subdomains*sizeof(PetscMPIInt),&old_ranks); 2351 ierr = PetscMalloc(size_prec_comm*sizeof(PetscInt),&new_ranks); 2352 ierr = MPI_Comm_rank(test_coarse_comm,&rank_coarse_comm);CHKERRQ(ierr); 2353 ierr = MPI_Comm_size(test_coarse_comm,&j);CHKERRQ(ierr); 2354 ierr = MPI_Allgather(&rank_prec_comm,1,MPIU_INT,old_ranks,1,MPIU_INT,test_coarse_comm);CHKERRQ(ierr); 2355 for (i=0; i<size_prec_comm; i++) new_ranks[i] = -1; 2356 for (i=0; i<n_subdomains; i++) new_ranks[old_ranks[i]] = i; 2357 ierr = PetscViewerASCIIOpen(test_coarse_comm,"test_mat_part.out",&viewer_test);CHKERRQ(ierr); 2358 k = pcis->n_neigh-1; 2359 ierr = PetscMalloc(2*sizeof(PetscInt),&ii); 2360 ii[0]=0; 2361 ii[1]=k; 2362 ierr = PetscMalloc(k*sizeof(PetscInt),&jj); 2363 for (i=0; i<k; i++) jj[i]=new_ranks[pcis->neigh[i+1]]; 2364 ierr = PetscSortInt(k,jj);CHKERRQ(ierr); 2365 ierr = MatCreateMPIAdj(test_coarse_comm,1,n_subdomains,ii,jj,NULL,&mat_adj);CHKERRQ(ierr); 2366 ierr = MatView(mat_adj,viewer_test);CHKERRQ(ierr); 2367 ierr = MatPartitioningCreate(test_coarse_comm,&mat_part);CHKERRQ(ierr); 2368 ierr = MatPartitioningSetAdjacency(mat_part,mat_adj);CHKERRQ(ierr); 2369 ierr = MatPartitioningSetFromOptions(mat_part);CHKERRQ(ierr); 2370 printf("Setting Nparts %d\n",n_parts); 2371 ierr = MatPartitioningSetNParts(mat_part,n_parts);CHKERRQ(ierr); 2372 ierr = MatPartitioningView(mat_part,viewer_test);CHKERRQ(ierr); 2373 ierr = MatPartitioningApply(mat_part,&coarse_new_decomposition);CHKERRQ(ierr); 2374 ierr = ISView(coarse_new_decomposition,viewer_test);CHKERRQ(ierr); 2375 ierr = ISPartitioningToNumbering(coarse_new_decomposition,&is_numbering);CHKERRQ(ierr); 2376 ierr = ISView(is_numbering,viewer_test);CHKERRQ(ierr); 2377 ierr = PetscViewerDestroy(&viewer_test);CHKERRQ(ierr); 2378 ierr = ISDestroy(&coarse_new_decomposition);CHKERRQ(ierr); 2379 ierr = ISDestroy(&is_numbering);CHKERRQ(ierr); 2380 ierr = MatPartitioningDestroy(&mat_part);CHKERRQ(ierr); 2381 ierr = PetscFree(old_ranks);CHKERRQ(ierr); 2382 ierr = PetscFree(new_ranks);CHKERRQ(ierr); 2383 ierr = MPI_Comm_free(&test_coarse_comm);CHKERRQ(ierr); 2384 } 2385 #endif 2386 2387 /* build CSR graph of subdomains' connectivity */ 2388 ierr = PetscMalloc (pcis->n*sizeof(PetscInt),&array_int);CHKERRQ(ierr); 2389 ierr = PetscMemzero(array_int,pcis->n*sizeof(PetscInt));CHKERRQ(ierr); 2390 for (i=1;i<pcis->n_neigh;i++){/* i=1 so I don't count myself -> faces nodes counts to 1 */ 2391 for (j=0;j<pcis->n_shared[i];j++){ 2392 array_int[ pcis->shared[i][j] ]+=1; 2393 } 2394 } 2395 for (i=1;i<pcis->n_neigh;i++){ 2396 for (j=0;j<pcis->n_shared[i];j++){ 2397 if (array_int[ pcis->shared[i][j] ] > 0 ){ 2398 my_faces++; 2399 break; 2400 } 2401 } 2402 } 2403 2404 ierr = MPI_Reduce(&my_faces,&total_faces,1,MPIU_INT,MPI_SUM,master_proc,prec_comm);CHKERRQ(ierr); 2405 ierr = PetscMalloc (my_faces*sizeof(PetscInt),&my_faces_connectivity);CHKERRQ(ierr); 2406 my_faces=0; 2407 for (i=1;i<pcis->n_neigh;i++){ 2408 for (j=0;j<pcis->n_shared[i];j++){ 2409 if (array_int[ pcis->shared[i][j] ] > 0 ){ 2410 my_faces_connectivity[my_faces]=pcis->neigh[i]; 2411 my_faces++; 2412 break; 2413 } 2414 } 2415 } 2416 if (rank_prec_comm == master_proc) { 2417 ierr = PetscMalloc (total_faces*sizeof(PetscMPIInt),&petsc_faces_adjncy);CHKERRQ(ierr); 2418 ierr = PetscMalloc (size_prec_comm*sizeof(PetscMPIInt),&number_of_faces);CHKERRQ(ierr); 2419 ierr = PetscMalloc (total_faces*sizeof(MetisInt),&faces_adjncy);CHKERRQ(ierr); 2420 ierr = PetscMalloc ((n_subdomains+1)*sizeof(MetisInt),&faces_xadj);CHKERRQ(ierr); 2421 ierr = PetscMalloc ((size_prec_comm+1)*sizeof(PetscMPIInt),&faces_displacements);CHKERRQ(ierr); 2422 } 2423 ierr = MPI_Gather(&my_faces,1,MPIU_INT,&number_of_faces[0],1,MPIU_INT,master_proc,prec_comm);CHKERRQ(ierr); 2424 if (rank_prec_comm == master_proc) { 2425 faces_xadj[0]=0; 2426 faces_displacements[0]=0; 2427 j=0; 2428 for (i=1;i<size_prec_comm+1;i++) { 2429 faces_displacements[i]=faces_displacements[i-1]+number_of_faces[i-1]; 2430 if (number_of_faces[i-1]) { 2431 j++; 2432 faces_xadj[j]=faces_xadj[j-1]+number_of_faces[i-1]; 2433 } 2434 } 2435 } 2436 ierr = MPI_Gatherv(&my_faces_connectivity[0],my_faces,MPIU_INT,&petsc_faces_adjncy[0],number_of_faces,faces_displacements,MPIU_INT,master_proc,prec_comm);CHKERRQ(ierr); 2437 ierr = PetscFree(my_faces_connectivity);CHKERRQ(ierr); 2438 ierr = PetscFree(array_int);CHKERRQ(ierr); 2439 if (rank_prec_comm == master_proc) { 2440 for (i=0;i<total_faces;i++) faces_adjncy[i]=(MetisInt)(petsc_faces_adjncy[i]/ranks_stretching_ratio); /* cast to MetisInt */ 2441 ierr = PetscFree(faces_displacements);CHKERRQ(ierr); 2442 ierr = PetscFree(number_of_faces);CHKERRQ(ierr); 2443 ierr = PetscFree(petsc_faces_adjncy);CHKERRQ(ierr); 2444 } 2445 2446 if ( rank_prec_comm == master_proc ) { 2447 2448 PetscInt heuristic_for_metis=3; 2449 2450 ncon=1; 2451 faces_nvtxs=n_subdomains; 2452 /* partition graoh induced by face connectivity */ 2453 ierr = PetscMalloc (n_subdomains*sizeof(MetisInt),&metis_coarse_subdivision);CHKERRQ(ierr); 2454 ierr = METIS_SetDefaultOptions(options); 2455 /* we need a contiguous partition of the coarse mesh */ 2456 options[METIS_OPTION_CONTIG]=1; 2457 options[METIS_OPTION_NITER]=30; 2458 if (pcbddc->coarsening_ratio > 1) { 2459 if (n_subdomains>n_parts*heuristic_for_metis) { 2460 options[METIS_OPTION_IPTYPE]=METIS_IPTYPE_EDGE; 2461 options[METIS_OPTION_OBJTYPE]=METIS_OBJTYPE_CUT; 2462 ierr = METIS_PartGraphKway(&faces_nvtxs,&ncon,faces_xadj,faces_adjncy,NULL,NULL,NULL,&n_parts,NULL,NULL,options,&objval,metis_coarse_subdivision); 2463 if (ierr != METIS_OK) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in METIS_PartGraphKway (metis error code %D) called from PCBDDCSetUpCoarseEnvironment\n",ierr); 2464 } else { 2465 ierr = METIS_PartGraphRecursive(&faces_nvtxs,&ncon,faces_xadj,faces_adjncy,NULL,NULL,NULL,&n_parts,NULL,NULL,options,&objval,metis_coarse_subdivision); 2466 if (ierr != METIS_OK) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in METIS_PartGraphRecursive (metis error code %D) called from PCBDDCSetUpCoarseEnvironment\n",ierr); 2467 } 2468 } else { 2469 for (i=0;i<n_subdomains;i++) metis_coarse_subdivision[i]=i; 2470 } 2471 ierr = PetscFree(faces_xadj);CHKERRQ(ierr); 2472 ierr = PetscFree(faces_adjncy);CHKERRQ(ierr); 2473 ierr = PetscMalloc(size_prec_comm*sizeof(PetscMPIInt),&coarse_subdivision);CHKERRQ(ierr); 2474 2475 /* copy/cast values avoiding possible type conflicts between PETSc, MPI and METIS */ 2476 for (i=0;i<size_prec_comm;i++) coarse_subdivision[i]=MPI_PROC_NULL; 2477 for (i=0;i<n_subdomains;i++) coarse_subdivision[ranks_stretching_ratio*i]=(PetscInt)(metis_coarse_subdivision[i]); 2478 ierr = PetscFree(metis_coarse_subdivision);CHKERRQ(ierr); 2479 } 2480 2481 /* Create new communicator for coarse problem splitting the old one */ 2482 if ( !(rank_prec_comm%procs_jumps_coarse_comm) && rank_prec_comm < procs_jumps_coarse_comm*n_parts ){ 2483 coarse_color=0; /* for communicator splitting */ 2484 active_rank=rank_prec_comm; /* for insertion of matrix values */ 2485 } 2486 /* procs with coarse_color = MPI_UNDEFINED will have coarse_comm = MPI_COMM_NULL (from mpi standards) 2487 key = rank_prec_comm -> keep same ordering of ranks from the old to the new communicator */ 2488 ierr = MPI_Comm_split(prec_comm,coarse_color,rank_prec_comm,&coarse_comm);CHKERRQ(ierr); 2489 2490 if ( coarse_color == 0 ) { 2491 ierr = MPI_Comm_size(coarse_comm,&size_coarse_comm);CHKERRQ(ierr); 2492 ierr = MPI_Comm_rank(coarse_comm,&rank_coarse_comm);CHKERRQ(ierr); 2493 } else { 2494 rank_coarse_comm = MPI_PROC_NULL; 2495 } 2496 2497 /* master proc take care of arranging and distributing coarse information */ 2498 if (rank_coarse_comm == master_proc) { 2499 ierr = PetscMalloc (size_coarse_comm*sizeof(PetscMPIInt),&displacements_recv);CHKERRQ(ierr); 2500 ierr = PetscMalloc (size_coarse_comm*sizeof(PetscMPIInt),&total_count_recv);CHKERRQ(ierr); 2501 ierr = PetscMalloc (n_subdomains*sizeof(PetscMPIInt),&total_ranks_recv);CHKERRQ(ierr); 2502 /* some initializations */ 2503 displacements_recv[0]=0; 2504 ierr = PetscMemzero(total_count_recv,size_coarse_comm*sizeof(PetscMPIInt));CHKERRQ(ierr); 2505 /* count from how many processes the j-th process of the coarse decomposition will receive data */ 2506 for (j=0;j<size_coarse_comm;j++) { 2507 for (i=0;i<size_prec_comm;i++) { 2508 if (coarse_subdivision[i]==j) total_count_recv[j]++; 2509 } 2510 } 2511 /* displacements needed for scatterv of total_ranks_recv */ 2512 for (i=1; i<size_coarse_comm; i++) displacements_recv[i]=displacements_recv[i-1]+total_count_recv[i-1]; 2513 2514 /* Now fill properly total_ranks_recv -> each coarse process will receive the ranks (in prec_comm communicator) of its friend (sending) processes */ 2515 ierr = PetscMemzero(total_count_recv,size_coarse_comm*sizeof(PetscMPIInt));CHKERRQ(ierr); 2516 for (j=0;j<size_coarse_comm;j++) { 2517 for (i=0;i<size_prec_comm;i++) { 2518 if (coarse_subdivision[i]==j) { 2519 total_ranks_recv[displacements_recv[j]+total_count_recv[j]]=i; 2520 total_count_recv[j]+=1; 2521 } 2522 } 2523 } 2524 /*for (j=0;j<size_coarse_comm;j++) { 2525 printf("process %d in new rank will receive from %d processes (original ranks follows)\n",j,total_count_recv[j]); 2526 for (i=0;i<total_count_recv[j];i++) { 2527 printf("%d ",total_ranks_recv[displacements_recv[j]+i]); 2528 } 2529 printf("\n"); 2530 }*/ 2531 2532 /* identify new decomposition in terms of ranks in the old communicator */ 2533 for (i=0;i<n_subdomains;i++) { 2534 coarse_subdivision[ranks_stretching_ratio*i]=coarse_subdivision[ranks_stretching_ratio*i]*procs_jumps_coarse_comm; 2535 } 2536 /*printf("coarse_subdivision in old end new ranks\n"); 2537 for (i=0;i<size_prec_comm;i++) 2538 if (coarse_subdivision[i]!=MPI_PROC_NULL) { 2539 printf("%d=(%d %d), ",i,coarse_subdivision[i],coarse_subdivision[i]/procs_jumps_coarse_comm); 2540 } else { 2541 printf("%d=(%d %d), ",i,coarse_subdivision[i],coarse_subdivision[i]); 2542 } 2543 printf("\n");*/ 2544 } 2545 2546 /* Scatter new decomposition for send details */ 2547 ierr = MPI_Scatter(&coarse_subdivision[0],1,MPIU_INT,&rank_coarse_proc_send_to,1,MPIU_INT,master_proc,prec_comm);CHKERRQ(ierr); 2548 /* Scatter receiving details to members of coarse decomposition */ 2549 if ( coarse_color == 0) { 2550 ierr = MPI_Scatter(&total_count_recv[0],1,MPIU_INT,&count_recv,1,MPIU_INT,master_proc,coarse_comm);CHKERRQ(ierr); 2551 ierr = PetscMalloc (count_recv*sizeof(PetscMPIInt),&ranks_recv);CHKERRQ(ierr); 2552 ierr = MPI_Scatterv(&total_ranks_recv[0],total_count_recv,displacements_recv,MPIU_INT,&ranks_recv[0],count_recv,MPIU_INT,master_proc,coarse_comm);CHKERRQ(ierr); 2553 } 2554 2555 /*printf("I will send my matrix data to proc %d\n",rank_coarse_proc_send_to); 2556 if (coarse_color == 0) { 2557 printf("I will receive some matrix data from %d processes (ranks follows)\n",count_recv); 2558 for (i=0;i<count_recv;i++) 2559 printf("%d ",ranks_recv[i]); 2560 printf("\n"); 2561 }*/ 2562 2563 if (rank_prec_comm == master_proc) { 2564 ierr = PetscFree(coarse_subdivision);CHKERRQ(ierr); 2565 ierr = PetscFree(total_count_recv);CHKERRQ(ierr); 2566 ierr = PetscFree(total_ranks_recv);CHKERRQ(ierr); 2567 ierr = PetscFree(displacements_recv);CHKERRQ(ierr); 2568 } 2569 #endif 2570 break; 2571 } 2572 2573 case(REPLICATED_BDDC): 2574 2575 pcbddc->coarse_communications_type = GATHERS_BDDC; 2576 coarse_mat_type = MATSEQAIJ; 2577 coarse_pc_type = PCLU; 2578 coarse_ksp_type = KSPPREONLY; 2579 coarse_comm = PETSC_COMM_SELF; 2580 active_rank = rank_prec_comm; 2581 break; 2582 2583 case(PARALLEL_BDDC): 2584 2585 pcbddc->coarse_communications_type = SCATTERS_BDDC; 2586 coarse_mat_type = MATAIJ; 2587 coarse_pc_type = PCREDUNDANT; 2588 coarse_ksp_type = KSPPREONLY; 2589 coarse_comm = prec_comm; 2590 active_rank = rank_prec_comm; 2591 break; 2592 2593 case(SEQUENTIAL_BDDC): 2594 pcbddc->coarse_communications_type = GATHERS_BDDC; 2595 coarse_mat_type = MATAIJ; 2596 coarse_pc_type = PCLU; 2597 coarse_ksp_type = KSPPREONLY; 2598 coarse_comm = PETSC_COMM_SELF; 2599 active_rank = master_proc; 2600 break; 2601 } 2602 2603 switch(pcbddc->coarse_communications_type){ 2604 2605 case(SCATTERS_BDDC): 2606 { 2607 if (pcbddc->coarse_problem_type==MULTILEVEL_BDDC) { 2608 2609 IS coarse_IS; 2610 2611 if(pcbddc->coarsening_ratio == 1) { 2612 ins_local_primal_size = pcbddc->local_primal_size; 2613 ins_local_primal_indices = pcbddc->local_primal_indices; 2614 if (coarse_color == 0) { ierr = PetscFree(ranks_recv);CHKERRQ(ierr); } 2615 /* nonzeros */ 2616 ierr = PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&dnz);CHKERRQ(ierr); 2617 ierr = PetscMemzero(dnz,ins_local_primal_size*sizeof(PetscInt));CHKERRQ(ierr); 2618 for (i=0;i<ins_local_primal_size;i++) { 2619 dnz[i] = ins_local_primal_size; 2620 } 2621 } else { 2622 PetscMPIInt send_size; 2623 PetscMPIInt *send_buffer; 2624 PetscInt *aux_ins_indices; 2625 PetscInt ii,jj; 2626 MPI_Request *requests; 2627 2628 ierr = PetscMalloc(count_recv*sizeof(PetscMPIInt),&localdispl2);CHKERRQ(ierr); 2629 /* reusing pcbddc->local_primal_displacements and pcbddc->replicated_primal_size */ 2630 ierr = PetscFree(pcbddc->local_primal_displacements);CHKERRQ(ierr); 2631 ierr = PetscMalloc((count_recv+1)*sizeof(PetscMPIInt),&pcbddc->local_primal_displacements);CHKERRQ(ierr); 2632 pcbddc->replicated_primal_size = count_recv; 2633 j = 0; 2634 for (i=0;i<count_recv;i++) { 2635 pcbddc->local_primal_displacements[i] = j; 2636 j += pcbddc->local_primal_sizes[ranks_recv[i]]; 2637 } 2638 pcbddc->local_primal_displacements[count_recv] = j; 2639 ierr = PetscMalloc(j*sizeof(PetscMPIInt),&pcbddc->replicated_local_primal_indices);CHKERRQ(ierr); 2640 /* allocate auxiliary space */ 2641 ierr = PetscMalloc(count_recv*sizeof(PetscMPIInt),&localsizes2);CHKERRQ(ierr); 2642 ierr = PetscMalloc(pcbddc->coarse_size*sizeof(PetscInt),&aux_ins_indices);CHKERRQ(ierr); 2643 ierr = PetscMemzero(aux_ins_indices,pcbddc->coarse_size*sizeof(PetscInt));CHKERRQ(ierr); 2644 /* allocate stuffs for message massing */ 2645 ierr = PetscMalloc((count_recv+1)*sizeof(MPI_Request),&requests);CHKERRQ(ierr); 2646 for (i=0;i<count_recv+1;i++) { requests[i]=MPI_REQUEST_NULL; } 2647 /* send indices to be inserted */ 2648 for (i=0;i<count_recv;i++) { 2649 send_size = pcbddc->local_primal_sizes[ranks_recv[i]]; 2650 ierr = MPI_Irecv(&pcbddc->replicated_local_primal_indices[pcbddc->local_primal_displacements[i]],send_size,MPIU_INT,ranks_recv[i],999,prec_comm,&requests[i]);CHKERRQ(ierr); 2651 } 2652 if (rank_coarse_proc_send_to != MPI_PROC_NULL ) { 2653 send_size = pcbddc->local_primal_size; 2654 ierr = PetscMalloc(send_size*sizeof(PetscMPIInt),&send_buffer);CHKERRQ(ierr); 2655 for (i=0;i<send_size;i++) { 2656 send_buffer[i]=(PetscMPIInt)pcbddc->local_primal_indices[i]; 2657 } 2658 ierr = MPI_Isend(send_buffer,send_size,MPIU_INT,rank_coarse_proc_send_to,999,prec_comm,&requests[count_recv]);CHKERRQ(ierr); 2659 } 2660 ierr = MPI_Waitall(count_recv+1,requests,MPI_STATUSES_IGNORE);CHKERRQ(ierr); 2661 if (rank_coarse_proc_send_to != MPI_PROC_NULL ) { 2662 ierr = PetscFree(send_buffer);CHKERRQ(ierr); 2663 } 2664 j = 0; 2665 for (i=0;i<count_recv;i++) { 2666 ii = pcbddc->local_primal_displacements[i+1]-pcbddc->local_primal_displacements[i]; 2667 localsizes2[i] = ii*ii; 2668 localdispl2[i] = j; 2669 j += localsizes2[i]; 2670 jj = pcbddc->local_primal_displacements[i]; 2671 /* it counts the coarse subdomains sharing the coarse node */ 2672 for (k=0;k<ii;k++) { 2673 aux_ins_indices[pcbddc->replicated_local_primal_indices[jj+k]] += 1; 2674 } 2675 } 2676 /* temp_coarse_mat_vals used to store matrix values to be received */ 2677 ierr = PetscMalloc(j*sizeof(PetscScalar),&temp_coarse_mat_vals);CHKERRQ(ierr); 2678 /* evaluate how many values I will insert in coarse mat */ 2679 ins_local_primal_size = 0; 2680 for (i=0;i<pcbddc->coarse_size;i++) { 2681 if (aux_ins_indices[i]) { 2682 ins_local_primal_size++; 2683 } 2684 } 2685 /* evaluate indices I will insert in coarse mat */ 2686 ierr = PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);CHKERRQ(ierr); 2687 j = 0; 2688 for(i=0;i<pcbddc->coarse_size;i++) { 2689 if(aux_ins_indices[i]) { 2690 ins_local_primal_indices[j] = i; 2691 j++; 2692 } 2693 } 2694 /* processes partecipating in coarse problem receive matrix data from their friends */ 2695 for (i=0;i<count_recv;i++) { 2696 ierr = MPI_Irecv(&temp_coarse_mat_vals[localdispl2[i]],localsizes2[i],MPIU_SCALAR,ranks_recv[i],666,prec_comm,&requests[i]);CHKERRQ(ierr); 2697 } 2698 if (rank_coarse_proc_send_to != MPI_PROC_NULL ) { 2699 send_size = pcbddc->local_primal_size*pcbddc->local_primal_size; 2700 ierr = MPI_Isend(&coarse_submat_vals[0],send_size,MPIU_SCALAR,rank_coarse_proc_send_to,666,prec_comm,&requests[count_recv]);CHKERRQ(ierr); 2701 } 2702 ierr = MPI_Waitall(count_recv+1,requests,MPI_STATUSES_IGNORE);CHKERRQ(ierr); 2703 /* nonzeros */ 2704 ierr = PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&dnz);CHKERRQ(ierr); 2705 ierr = PetscMemzero(dnz,ins_local_primal_size*sizeof(PetscInt));CHKERRQ(ierr); 2706 /* use aux_ins_indices to realize a global to local mapping */ 2707 j=0; 2708 for(i=0;i<pcbddc->coarse_size;i++){ 2709 if(aux_ins_indices[i]==0){ 2710 aux_ins_indices[i]=-1; 2711 } else { 2712 aux_ins_indices[i]=j; 2713 j++; 2714 } 2715 } 2716 for (i=0;i<count_recv;i++) { 2717 j = pcbddc->local_primal_sizes[ranks_recv[i]]; 2718 for (k=0;k<j;k++) { 2719 dnz[aux_ins_indices[pcbddc->replicated_local_primal_indices[pcbddc->local_primal_displacements[i]+k]]] += j; 2720 } 2721 } 2722 /* check */ 2723 for (i=0;i<ins_local_primal_size;i++) { 2724 if (dnz[i] > ins_local_primal_size) { 2725 dnz[i] = ins_local_primal_size; 2726 } 2727 } 2728 ierr = PetscFree(requests);CHKERRQ(ierr); 2729 ierr = PetscFree(aux_ins_indices);CHKERRQ(ierr); 2730 if (coarse_color == 0) { ierr = PetscFree(ranks_recv);CHKERRQ(ierr); } 2731 } 2732 /* create local to global mapping needed by coarse MATIS */ 2733 if (coarse_comm != MPI_COMM_NULL ) {ierr = MPI_Comm_free(&coarse_comm);CHKERRQ(ierr);} 2734 coarse_comm = prec_comm; 2735 active_rank = rank_prec_comm; 2736 ierr = ISCreateGeneral(coarse_comm,ins_local_primal_size,ins_local_primal_indices,PETSC_COPY_VALUES,&coarse_IS);CHKERRQ(ierr); 2737 ierr = ISLocalToGlobalMappingCreateIS(coarse_IS,&coarse_ISLG);CHKERRQ(ierr); 2738 ierr = ISDestroy(&coarse_IS);CHKERRQ(ierr); 2739 } else if (pcbddc->coarse_problem_type==PARALLEL_BDDC) { 2740 /* arrays for values insertion */ 2741 ins_local_primal_size = pcbddc->local_primal_size; 2742 ierr = PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);CHKERRQ(ierr); 2743 ierr = PetscMalloc(ins_local_primal_size*ins_local_primal_size*sizeof(PetscScalar),&ins_coarse_mat_vals);CHKERRQ(ierr); 2744 for (j=0;j<ins_local_primal_size;j++){ 2745 ins_local_primal_indices[j]=pcbddc->local_primal_indices[j]; 2746 for (i=0;i<ins_local_primal_size;i++) { 2747 ins_coarse_mat_vals[j*ins_local_primal_size+i]=coarse_submat_vals[j*ins_local_primal_size+i]; 2748 } 2749 } 2750 } 2751 break; 2752 2753 } 2754 2755 case(GATHERS_BDDC): 2756 { 2757 2758 PetscMPIInt mysize,mysize2; 2759 PetscMPIInt *send_buffer; 2760 2761 if (rank_prec_comm==active_rank) { 2762 ierr = PetscMalloc ( pcbddc->replicated_primal_size*sizeof(PetscMPIInt),&pcbddc->replicated_local_primal_indices);CHKERRQ(ierr); 2763 ierr = PetscMalloc ( pcbddc->replicated_primal_size*sizeof(PetscScalar),&pcbddc->replicated_local_primal_values);CHKERRQ(ierr); 2764 ierr = PetscMalloc ( size_prec_comm*sizeof(PetscMPIInt),&localsizes2);CHKERRQ(ierr); 2765 ierr = PetscMalloc ( size_prec_comm*sizeof(PetscMPIInt),&localdispl2);CHKERRQ(ierr); 2766 /* arrays for values insertion */ 2767 for (i=0;i<size_prec_comm;i++) localsizes2[i]=pcbddc->local_primal_sizes[i]*pcbddc->local_primal_sizes[i]; 2768 localdispl2[0]=0; 2769 for (i=1;i<size_prec_comm;i++) localdispl2[i]=localsizes2[i-1]+localdispl2[i-1]; 2770 j=0; 2771 for (i=0;i<size_prec_comm;i++) j+=localsizes2[i]; 2772 ierr = PetscMalloc ( j*sizeof(PetscScalar),&temp_coarse_mat_vals);CHKERRQ(ierr); 2773 } 2774 2775 mysize=pcbddc->local_primal_size; 2776 mysize2=pcbddc->local_primal_size*pcbddc->local_primal_size; 2777 ierr = PetscMalloc(mysize*sizeof(PetscMPIInt),&send_buffer);CHKERRQ(ierr); 2778 for (i=0; i<mysize; i++) send_buffer[i]=(PetscMPIInt)pcbddc->local_primal_indices[i]; 2779 2780 if (pcbddc->coarse_problem_type == SEQUENTIAL_BDDC){ 2781 ierr = MPI_Gatherv(send_buffer,mysize,MPIU_INT,&pcbddc->replicated_local_primal_indices[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_INT,master_proc,prec_comm);CHKERRQ(ierr); 2782 ierr = MPI_Gatherv(&coarse_submat_vals[0],mysize2,MPIU_SCALAR,&temp_coarse_mat_vals[0],localsizes2,localdispl2,MPIU_SCALAR,master_proc,prec_comm);CHKERRQ(ierr); 2783 } else { 2784 ierr = MPI_Allgatherv(send_buffer,mysize,MPIU_INT,&pcbddc->replicated_local_primal_indices[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_INT,prec_comm);CHKERRQ(ierr); 2785 ierr = MPI_Allgatherv(&coarse_submat_vals[0],mysize2,MPIU_SCALAR,&temp_coarse_mat_vals[0],localsizes2,localdispl2,MPIU_SCALAR,prec_comm);CHKERRQ(ierr); 2786 } 2787 ierr = PetscFree(send_buffer);CHKERRQ(ierr); 2788 break; 2789 }/* switch on coarse problem and communications associated with finished */ 2790 } 2791 2792 /* Now create and fill up coarse matrix */ 2793 if ( rank_prec_comm == active_rank ) { 2794 2795 Mat matis_coarse_local_mat; 2796 2797 if (pcbddc->coarse_problem_type != MULTILEVEL_BDDC) { 2798 ierr = MatCreate(coarse_comm,&pcbddc->coarse_mat);CHKERRQ(ierr); 2799 ierr = MatSetSizes(pcbddc->coarse_mat,PETSC_DECIDE,PETSC_DECIDE,pcbddc->coarse_size,pcbddc->coarse_size);CHKERRQ(ierr); 2800 ierr = MatSetType(pcbddc->coarse_mat,coarse_mat_type);CHKERRQ(ierr); 2801 ierr = MatSetOptionsPrefix(pcbddc->coarse_mat,"coarse_");CHKERRQ(ierr); 2802 ierr = MatSetFromOptions(pcbddc->coarse_mat);CHKERRQ(ierr); 2803 ierr = MatSetUp(pcbddc->coarse_mat);CHKERRQ(ierr); 2804 ierr = MatSetOption(pcbddc->coarse_mat,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); /* local values stored in column major */ 2805 ierr = MatSetOption(pcbddc->coarse_mat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 2806 } else { 2807 ierr = MatCreateIS(coarse_comm,1,PETSC_DECIDE,PETSC_DECIDE,pcbddc->coarse_size,pcbddc->coarse_size,coarse_ISLG,&pcbddc->coarse_mat);CHKERRQ(ierr); 2808 ierr = MatSetUp(pcbddc->coarse_mat);CHKERRQ(ierr); 2809 ierr = MatISGetLocalMat(pcbddc->coarse_mat,&matis_coarse_local_mat);CHKERRQ(ierr); 2810 ierr = MatSetOptionsPrefix(pcbddc->coarse_mat,"coarse_");CHKERRQ(ierr); 2811 ierr = MatSetFromOptions(pcbddc->coarse_mat);CHKERRQ(ierr); 2812 ierr = MatSetUp(matis_coarse_local_mat);CHKERRQ(ierr); 2813 ierr = MatSetOption(matis_coarse_local_mat,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); /* local values stored in column major */ 2814 ierr = MatSetOption(matis_coarse_local_mat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 2815 } 2816 /* preallocation */ 2817 if (pcbddc->coarse_problem_type != MULTILEVEL_BDDC) { 2818 2819 PetscInt lrows,lcols,bs; 2820 2821 ierr = MatGetLocalSize(pcbddc->coarse_mat,&lrows,&lcols);CHKERRQ(ierr); 2822 ierr = MatPreallocateInitialize(coarse_comm,lrows,lcols,dnz,onz);CHKERRQ(ierr); 2823 ierr = MatGetBlockSize(pcbddc->coarse_mat,&bs);CHKERRQ(ierr); 2824 2825 if (pcbddc->coarse_problem_type == PARALLEL_BDDC) { 2826 2827 Vec vec_dnz,vec_onz; 2828 PetscScalar *my_dnz,*my_onz,*array; 2829 PetscInt *mat_ranges,*row_ownership; 2830 PetscInt coarse_index_row,coarse_index_col,owner; 2831 2832 ierr = VecCreate(prec_comm,&vec_dnz);CHKERRQ(ierr); 2833 ierr = VecSetBlockSize(vec_dnz,bs);CHKERRQ(ierr); 2834 ierr = VecSetSizes(vec_dnz,PETSC_DECIDE,pcbddc->coarse_size);CHKERRQ(ierr); 2835 ierr = VecSetType(vec_dnz,VECMPI);CHKERRQ(ierr); 2836 ierr = VecDuplicate(vec_dnz,&vec_onz);CHKERRQ(ierr); 2837 2838 ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(PetscScalar),&my_dnz);CHKERRQ(ierr); 2839 ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(PetscScalar),&my_onz);CHKERRQ(ierr); 2840 ierr = PetscMemzero(my_dnz,pcbddc->local_primal_size*sizeof(PetscScalar));CHKERRQ(ierr); 2841 ierr = PetscMemzero(my_onz,pcbddc->local_primal_size*sizeof(PetscScalar));CHKERRQ(ierr); 2842 2843 ierr = PetscMalloc(pcbddc->coarse_size*sizeof(PetscInt),&row_ownership);CHKERRQ(ierr); 2844 ierr = MatGetOwnershipRanges(pcbddc->coarse_mat,(const PetscInt**)&mat_ranges);CHKERRQ(ierr); 2845 for (i=0;i<size_prec_comm;i++) { 2846 for (j=mat_ranges[i];j<mat_ranges[i+1];j++) { 2847 row_ownership[j]=i; 2848 } 2849 } 2850 2851 for (i=0;i<pcbddc->local_primal_size;i++) { 2852 coarse_index_row = pcbddc->local_primal_indices[i]; 2853 owner = row_ownership[coarse_index_row]; 2854 for (j=i;j<pcbddc->local_primal_size;j++) { 2855 owner = row_ownership[coarse_index_row]; 2856 coarse_index_col = pcbddc->local_primal_indices[j]; 2857 if (coarse_index_col > mat_ranges[owner]-1 && coarse_index_col < mat_ranges[owner+1] ) { 2858 my_dnz[i] += 1.0; 2859 } else { 2860 my_onz[i] += 1.0; 2861 } 2862 if (i != j) { 2863 owner = row_ownership[coarse_index_col]; 2864 if (coarse_index_row > mat_ranges[owner]-1 && coarse_index_row < mat_ranges[owner+1] ) { 2865 my_dnz[j] += 1.0; 2866 } else { 2867 my_onz[j] += 1.0; 2868 } 2869 } 2870 } 2871 } 2872 ierr = VecSet(vec_dnz,0.0);CHKERRQ(ierr); 2873 ierr = VecSet(vec_onz,0.0);CHKERRQ(ierr); 2874 if (pcbddc->local_primal_size) { 2875 ierr = VecSetValues(vec_dnz,pcbddc->local_primal_size,pcbddc->local_primal_indices,my_dnz,ADD_VALUES);CHKERRQ(ierr); 2876 ierr = VecSetValues(vec_onz,pcbddc->local_primal_size,pcbddc->local_primal_indices,my_onz,ADD_VALUES);CHKERRQ(ierr); 2877 } 2878 ierr = VecAssemblyBegin(vec_dnz);CHKERRQ(ierr); 2879 ierr = VecAssemblyBegin(vec_onz);CHKERRQ(ierr); 2880 ierr = VecAssemblyEnd(vec_dnz);CHKERRQ(ierr); 2881 ierr = VecAssemblyEnd(vec_onz);CHKERRQ(ierr); 2882 j = mat_ranges[rank_prec_comm+1]-mat_ranges[rank_prec_comm]; 2883 ierr = VecGetArray(vec_dnz,&array);CHKERRQ(ierr); 2884 for (i=0; i<j; i++) dnz[i] = (PetscInt)PetscRealPart(array[i]); 2885 2886 ierr = VecRestoreArray(vec_dnz,&array);CHKERRQ(ierr); 2887 ierr = VecGetArray(vec_onz,&array);CHKERRQ(ierr); 2888 for (i=0;i<j;i++) onz[i] = (PetscInt)PetscRealPart(array[i]); 2889 2890 ierr = VecRestoreArray(vec_onz,&array);CHKERRQ(ierr); 2891 ierr = PetscFree(my_dnz);CHKERRQ(ierr); 2892 ierr = PetscFree(my_onz);CHKERRQ(ierr); 2893 ierr = PetscFree(row_ownership);CHKERRQ(ierr); 2894 ierr = VecDestroy(&vec_dnz);CHKERRQ(ierr); 2895 ierr = VecDestroy(&vec_onz);CHKERRQ(ierr); 2896 } else { 2897 for (k=0;k<size_prec_comm;k++){ 2898 offset=pcbddc->local_primal_displacements[k]; 2899 offset2=localdispl2[k]; 2900 ins_local_primal_size = pcbddc->local_primal_sizes[k]; 2901 ierr = PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);CHKERRQ(ierr); 2902 for (j=0;j<ins_local_primal_size;j++){ 2903 ins_local_primal_indices[j]=(PetscInt)pcbddc->replicated_local_primal_indices[offset+j]; 2904 } 2905 for (j=0;j<ins_local_primal_size;j++) { 2906 ierr = MatPreallocateSet(ins_local_primal_indices[j],ins_local_primal_size,ins_local_primal_indices,dnz,onz);CHKERRQ(ierr); 2907 } 2908 ierr = PetscFree(ins_local_primal_indices);CHKERRQ(ierr); 2909 } 2910 } 2911 2912 /* check */ 2913 for (i=0;i<lrows;i++) { 2914 if (dnz[i]>lcols) dnz[i]=lcols; 2915 if (onz[i]>pcbddc->coarse_size-lcols) onz[i]=pcbddc->coarse_size-lcols; 2916 } 2917 ierr = MatSeqAIJSetPreallocation(pcbddc->coarse_mat,0,dnz);CHKERRQ(ierr); 2918 ierr = MatMPIAIJSetPreallocation(pcbddc->coarse_mat,0,dnz,0,onz);CHKERRQ(ierr); 2919 ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 2920 } else { 2921 ierr = MatSeqAIJSetPreallocation(matis_coarse_local_mat,0,dnz);CHKERRQ(ierr); 2922 ierr = PetscFree(dnz);CHKERRQ(ierr); 2923 } 2924 /* insert values */ 2925 if (pcbddc->coarse_problem_type == PARALLEL_BDDC) { 2926 ierr = MatSetValues(pcbddc->coarse_mat,ins_local_primal_size,ins_local_primal_indices,ins_local_primal_size,ins_local_primal_indices,ins_coarse_mat_vals,ADD_VALUES);CHKERRQ(ierr); 2927 } else if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 2928 if (pcbddc->coarsening_ratio == 1) { 2929 ins_coarse_mat_vals = coarse_submat_vals; 2930 ierr = MatSetValues(pcbddc->coarse_mat,ins_local_primal_size,ins_local_primal_indices,ins_local_primal_size,ins_local_primal_indices,ins_coarse_mat_vals,INSERT_VALUES);CHKERRQ(ierr); 2931 } else { 2932 ierr = PetscFree(ins_local_primal_indices);CHKERRQ(ierr); 2933 for (k=0;k<pcbddc->replicated_primal_size;k++) { 2934 offset = pcbddc->local_primal_displacements[k]; 2935 offset2 = localdispl2[k]; 2936 ins_local_primal_size = pcbddc->local_primal_displacements[k+1]-pcbddc->local_primal_displacements[k]; 2937 ierr = PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);CHKERRQ(ierr); 2938 for (j=0;j<ins_local_primal_size;j++){ 2939 ins_local_primal_indices[j]=(PetscInt)pcbddc->replicated_local_primal_indices[offset+j]; 2940 } 2941 ins_coarse_mat_vals = &temp_coarse_mat_vals[offset2]; 2942 ierr = MatSetValues(pcbddc->coarse_mat,ins_local_primal_size,ins_local_primal_indices,ins_local_primal_size,ins_local_primal_indices,ins_coarse_mat_vals,ADD_VALUES);CHKERRQ(ierr); 2943 ierr = PetscFree(ins_local_primal_indices);CHKERRQ(ierr); 2944 } 2945 } 2946 ins_local_primal_indices = 0; 2947 ins_coarse_mat_vals = 0; 2948 } else { 2949 for (k=0;k<size_prec_comm;k++){ 2950 offset=pcbddc->local_primal_displacements[k]; 2951 offset2=localdispl2[k]; 2952 ins_local_primal_size = pcbddc->local_primal_sizes[k]; 2953 ierr = PetscMalloc(ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);CHKERRQ(ierr); 2954 for (j=0;j<ins_local_primal_size;j++){ 2955 ins_local_primal_indices[j]=(PetscInt)pcbddc->replicated_local_primal_indices[offset+j]; 2956 } 2957 ins_coarse_mat_vals = &temp_coarse_mat_vals[offset2]; 2958 ierr = MatSetValues(pcbddc->coarse_mat,ins_local_primal_size,ins_local_primal_indices,ins_local_primal_size,ins_local_primal_indices,ins_coarse_mat_vals,ADD_VALUES);CHKERRQ(ierr); 2959 ierr = PetscFree(ins_local_primal_indices);CHKERRQ(ierr); 2960 } 2961 ins_local_primal_indices = 0; 2962 ins_coarse_mat_vals = 0; 2963 } 2964 ierr = MatAssemblyBegin(pcbddc->coarse_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2965 ierr = MatAssemblyEnd(pcbddc->coarse_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2966 /* symmetry of coarse matrix */ 2967 if (issym) { 2968 ierr = MatSetOption(pcbddc->coarse_mat,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 2969 } 2970 ierr = MatGetVecs(pcbddc->coarse_mat,&pcbddc->coarse_vec,&pcbddc->coarse_rhs);CHKERRQ(ierr); 2971 } 2972 2973 /* create loc to glob scatters if needed */ 2974 if (pcbddc->coarse_communications_type == SCATTERS_BDDC) { 2975 IS local_IS,global_IS; 2976 ierr = ISCreateStride(PETSC_COMM_SELF,pcbddc->local_primal_size,0,1,&local_IS);CHKERRQ(ierr); 2977 ierr = ISCreateGeneral(PETSC_COMM_SELF,pcbddc->local_primal_size,pcbddc->local_primal_indices,PETSC_COPY_VALUES,&global_IS);CHKERRQ(ierr); 2978 ierr = VecScatterCreate(pcbddc->vec1_P,local_IS,pcbddc->coarse_vec,global_IS,&pcbddc->coarse_loc_to_glob);CHKERRQ(ierr); 2979 ierr = ISDestroy(&local_IS);CHKERRQ(ierr); 2980 ierr = ISDestroy(&global_IS);CHKERRQ(ierr); 2981 } 2982 2983 /* free memory no longer needed */ 2984 if (coarse_ISLG) { ierr = ISLocalToGlobalMappingDestroy(&coarse_ISLG);CHKERRQ(ierr); } 2985 if (ins_local_primal_indices) { ierr = PetscFree(ins_local_primal_indices);CHKERRQ(ierr); } 2986 if (ins_coarse_mat_vals) { ierr = PetscFree(ins_coarse_mat_vals);CHKERRQ(ierr); } 2987 if (localsizes2) { ierr = PetscFree(localsizes2);CHKERRQ(ierr); } 2988 if (localdispl2) { ierr = PetscFree(localdispl2);CHKERRQ(ierr); } 2989 if (temp_coarse_mat_vals) { ierr = PetscFree(temp_coarse_mat_vals);CHKERRQ(ierr); } 2990 2991 /* Compute coarse null space */ 2992 CoarseNullSpace = 0; 2993 if (pcbddc->NullSpace) { 2994 ierr = PCBDDCNullSpaceAssembleCoarse(pc,&CoarseNullSpace);CHKERRQ(ierr); 2995 } 2996 2997 /* KSP for coarse problem */ 2998 if (rank_prec_comm == active_rank) { 2999 PetscBool isbddc=PETSC_FALSE; 3000 3001 ierr = KSPCreate(coarse_comm,&pcbddc->coarse_ksp);CHKERRQ(ierr); 3002 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->coarse_ksp,(PetscObject)pc,1);CHKERRQ(ierr); 3003 ierr = KSPSetOperators(pcbddc->coarse_ksp,pcbddc->coarse_mat,pcbddc->coarse_mat,SAME_PRECONDITIONER);CHKERRQ(ierr); 3004 ierr = KSPSetTolerances(pcbddc->coarse_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,max_it_coarse_ksp);CHKERRQ(ierr); 3005 ierr = KSPSetType(pcbddc->coarse_ksp,coarse_ksp_type);CHKERRQ(ierr); 3006 ierr = KSPGetPC(pcbddc->coarse_ksp,&pc_temp);CHKERRQ(ierr); 3007 ierr = PCSetType(pc_temp,coarse_pc_type);CHKERRQ(ierr); 3008 /* Allow user's customization */ 3009 ierr = KSPSetOptionsPrefix(pcbddc->coarse_ksp,"coarse_");CHKERRQ(ierr); 3010 /* Set Up PC for coarse problem BDDC */ 3011 if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 3012 i = pcbddc->current_level+1; 3013 ierr = PCBDDCSetLevel(pc_temp,i);CHKERRQ(ierr); 3014 ierr = PCBDDCSetCoarseningRatio(pc_temp,pcbddc->coarsening_ratio);CHKERRQ(ierr); 3015 ierr = PCBDDCSetMaxLevels(pc_temp,pcbddc->max_levels);CHKERRQ(ierr); 3016 ierr = PCBDDCSetCoarseProblemType(pc_temp,MULTILEVEL_BDDC);CHKERRQ(ierr); 3017 if (CoarseNullSpace) { 3018 ierr = PCBDDCSetNullSpace(pc_temp,CoarseNullSpace);CHKERRQ(ierr); 3019 } 3020 if (dbg_flag) { 3021 ierr = PetscViewerASCIIPrintf(viewer,"----------------Level %d: Setting up level %d---------------\n",pcbddc->current_level,i);CHKERRQ(ierr); 3022 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 3023 } 3024 } else { 3025 if (CoarseNullSpace) { 3026 ierr = KSPSetNullSpace(pcbddc->coarse_ksp,CoarseNullSpace);CHKERRQ(ierr); 3027 } 3028 } 3029 ierr = KSPSetFromOptions(pcbddc->coarse_ksp);CHKERRQ(ierr); 3030 ierr = KSPSetUp(pcbddc->coarse_ksp);CHKERRQ(ierr); 3031 3032 ierr = KSPGetTolerances(pcbddc->coarse_ksp,NULL,NULL,NULL,&j);CHKERRQ(ierr); 3033 ierr = KSPGetPC(pcbddc->coarse_ksp,&pc_temp);CHKERRQ(ierr); 3034 ierr = PetscObjectTypeCompare((PetscObject)pc_temp,PCBDDC,&isbddc);CHKERRQ(ierr); 3035 if (j == 1) { 3036 ierr = KSPSetNormType(pcbddc->coarse_ksp,KSP_NORM_NONE);CHKERRQ(ierr); 3037 if (isbddc) { 3038 ierr = PCBDDCSetUseExactDirichlet(pc_temp,PETSC_FALSE);CHKERRQ(ierr); 3039 } 3040 } 3041 } 3042 /* Check coarse problem if requested */ 3043 if ( dbg_flag && rank_prec_comm == active_rank ) { 3044 KSP check_ksp; 3045 PC check_pc; 3046 Vec check_vec; 3047 PetscReal abs_infty_error,infty_error,lambda_min,lambda_max; 3048 KSPType check_ksp_type; 3049 3050 /* Create ksp object suitable for extreme eigenvalues' estimation */ 3051 ierr = KSPCreate(coarse_comm,&check_ksp);CHKERRQ(ierr); 3052 ierr = KSPSetOperators(check_ksp,pcbddc->coarse_mat,pcbddc->coarse_mat,SAME_PRECONDITIONER);CHKERRQ(ierr); 3053 ierr = KSPSetTolerances(check_ksp,1.e-12,1.e-12,PETSC_DEFAULT,pcbddc->coarse_size);CHKERRQ(ierr); 3054 if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 3055 if (issym) check_ksp_type = KSPCG; 3056 else check_ksp_type = KSPGMRES; 3057 ierr = KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);CHKERRQ(ierr); 3058 } else { 3059 check_ksp_type = KSPPREONLY; 3060 } 3061 ierr = KSPSetType(check_ksp,check_ksp_type);CHKERRQ(ierr); 3062 ierr = KSPGetPC(pcbddc->coarse_ksp,&check_pc);CHKERRQ(ierr); 3063 ierr = KSPSetPC(check_ksp,check_pc);CHKERRQ(ierr); 3064 ierr = KSPSetUp(check_ksp);CHKERRQ(ierr); 3065 /* create random vec */ 3066 ierr = VecDuplicate(pcbddc->coarse_vec,&check_vec);CHKERRQ(ierr); 3067 ierr = VecSetRandom(check_vec,NULL);CHKERRQ(ierr); 3068 if (CoarseNullSpace) { 3069 ierr = MatNullSpaceRemove(CoarseNullSpace,check_vec);CHKERRQ(ierr); 3070 } 3071 ierr = MatMult(pcbddc->coarse_mat,check_vec,pcbddc->coarse_rhs);CHKERRQ(ierr); 3072 /* solve coarse problem */ 3073 ierr = KSPSolve(check_ksp,pcbddc->coarse_rhs,pcbddc->coarse_vec);CHKERRQ(ierr); 3074 if (CoarseNullSpace) { 3075 ierr = MatNullSpaceRemove(CoarseNullSpace,pcbddc->coarse_vec);CHKERRQ(ierr); 3076 } 3077 /* check coarse problem residual error */ 3078 ierr = VecAXPY(check_vec,-1.0,pcbddc->coarse_vec);CHKERRQ(ierr); 3079 ierr = VecNorm(check_vec,NORM_INFINITY,&infty_error);CHKERRQ(ierr); 3080 ierr = MatMult(pcbddc->coarse_mat,check_vec,pcbddc->coarse_rhs);CHKERRQ(ierr); 3081 ierr = VecNorm(pcbddc->coarse_rhs,NORM_INFINITY,&abs_infty_error);CHKERRQ(ierr); 3082 ierr = VecDestroy(&check_vec);CHKERRQ(ierr); 3083 /* get eigenvalue estimation if inexact */ 3084 if (pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 3085 ierr = KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);CHKERRQ(ierr); 3086 ierr = KSPGetIterationNumber(check_ksp,&k);CHKERRQ(ierr); 3087 ierr = PetscViewerASCIIPrintf(viewer,"Coarse problem eigenvalues estimated with %d iterations of %s.\n",k,check_ksp_type);CHKERRQ(ierr); 3088 ierr = PetscViewerASCIIPrintf(viewer,"Coarse problem eigenvalues: % 1.14e %1.14e\n",lambda_min,lambda_max);CHKERRQ(ierr); 3089 } 3090 ierr = PetscViewerASCIIPrintf(viewer,"Coarse problem exact infty_error : %1.14e\n",infty_error);CHKERRQ(ierr); 3091 ierr = PetscViewerASCIIPrintf(viewer,"Coarse problem residual infty_error: %1.14e\n",abs_infty_error);CHKERRQ(ierr); 3092 ierr = KSPDestroy(&check_ksp);CHKERRQ(ierr); 3093 } 3094 if (dbg_flag) { 3095 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 3096 } 3097 ierr = MatNullSpaceDestroy(&CoarseNullSpace);CHKERRQ(ierr); 3098 3099 PetscFunctionReturn(0); 3100 } 3101 3102