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