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