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