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