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