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