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