1 /* TODOLIST 2 remove // commments 3 remove coarse enums -> allows use of PCBDDCGetCoarseKSP 4 code refactoring 5 remove metis dependency -> use MatPartitioning for multilevel 6 analyze MatSetNearNullSpace as suggested by Jed 7 options structure 8 */ 9 10 /* ---------------------------------------------------------------------------------------------------------------------------------------------- 11 Implementation of BDDC preconditioner based on: 12 C. Dohrmann "An approximate BDDC preconditioner", Numerical Linear Algebra with Applications Volume 14, Issue 2, pages 149-168, March 2007 13 ---------------------------------------------------------------------------------------------------------------------------------------------- */ 14 15 #include "bddc.h" /*I "petscpc.h" I*/ /* includes for fortran wrappers */ 16 17 /* -------------------------------------------------------------------------- */ 18 #undef __FUNCT__ 19 #define __FUNCT__ "PCSetFromOptions_BDDC" 20 PetscErrorCode PCSetFromOptions_BDDC(PC pc) 21 { 22 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 23 PetscErrorCode ierr; 24 25 PetscFunctionBegin; 26 ierr = PetscOptionsHead("BDDC options");CHKERRQ(ierr); 27 /* Verbose debugging of main data structures */ 28 ierr = PetscOptionsBool("-pc_bddc_check_all" ,"Verbose (debugging) output for PCBDDC" ,"none",pcbddc->dbg_flag ,&pcbddc->dbg_flag ,PETSC_NULL);CHKERRQ(ierr); 29 /* Some customization for default primal space */ 30 ierr = PetscOptionsBool("-pc_bddc_vertices_only" ,"Use vertices only in coarse space (i.e. discard constraints)","none",pcbddc->vertices_flag ,&pcbddc->vertices_flag ,PETSC_NULL);CHKERRQ(ierr); 31 ierr = PetscOptionsBool("-pc_bddc_constraints_only","Use constraints only in coarse space (i.e. discard vertices)","none",pcbddc->constraints_flag,&pcbddc->constraints_flag,PETSC_NULL);CHKERRQ(ierr); 32 ierr = PetscOptionsBool("-pc_bddc_faces_only" ,"Use faces only in coarse space (i.e. discard edges)" ,"none",pcbddc->faces_flag ,&pcbddc->faces_flag ,PETSC_NULL);CHKERRQ(ierr); 33 ierr = PetscOptionsBool("-pc_bddc_edges_only" ,"Use edges only in coarse space (i.e. discard faces)" ,"none",pcbddc->edges_flag ,&pcbddc->edges_flag ,PETSC_NULL);CHKERRQ(ierr); 34 /* Coarse solver context */ 35 static const char *avail_coarse_problems[] = {"sequential","replicated","parallel","multilevel",""}; //order of choiches depends on ENUM defined in bddc.h 36 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,PETSC_NULL);CHKERRQ(ierr); 37 /* Two different application of BDDC to the whole set of dofs, internal and interface */ 38 ierr = PetscOptionsBool("-pc_bddc_switch_preconditioning_type","Switch between M_2 (default) and M_3 preconditioners (as defined by Dohrmann)","none",pcbddc->prec_type,&pcbddc->prec_type,PETSC_NULL);CHKERRQ(ierr); 39 ierr = PetscOptionsInt("-pc_bddc_coarsening_ratio","Set coarsening ratio used in multilevel coarsening","none",pcbddc->coarsening_ratio,&pcbddc->coarsening_ratio,PETSC_NULL);CHKERRQ(ierr); 40 ierr = PetscOptionsTail();CHKERRQ(ierr); 41 PetscFunctionReturn(0); 42 } 43 44 /* -------------------------------------------------------------------------- */ 45 EXTERN_C_BEGIN 46 #undef __FUNCT__ 47 #define __FUNCT__ "PCBDDCSetCoarseProblemType_BDDC" 48 static PetscErrorCode PCBDDCSetCoarseProblemType_BDDC(PC pc, CoarseProblemType CPT) 49 { 50 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 51 52 PetscFunctionBegin; 53 pcbddc->coarse_problem_type = CPT; 54 PetscFunctionReturn(0); 55 } 56 EXTERN_C_END 57 58 /* -------------------------------------------------------------------------- */ 59 #undef __FUNCT__ 60 #define __FUNCT__ "PCBDDCSetCoarseProblemType" 61 /*@ 62 PCBDDCSetCoarseProblemType - brief explanation 63 64 Collective on PC 65 66 Input Parameters: 67 + pc - the preconditioning context 68 - CoarseProblemType - pick a better name and explain what this is 69 70 Level: intermediate 71 72 Notes: 73 usage notes, perhaps an example 74 75 .seealso: PCBDDC 76 @*/ 77 PetscErrorCode PCBDDCSetCoarseProblemType(PC pc, CoarseProblemType CPT) 78 { 79 PetscErrorCode ierr; 80 81 PetscFunctionBegin; 82 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 83 ierr = PetscTryMethod(pc,"PCBDDCSetCoarseProblemType_C",(PC,CoarseProblemType),(pc,CPT));CHKERRQ(ierr); 84 PetscFunctionReturn(0); 85 } 86 87 /* -------------------------------------------------------------------------- */ 88 EXTERN_C_BEGIN 89 #undef __FUNCT__ 90 #define __FUNCT__ "PCBDDCSetNeumannBoundaries_BDDC" 91 static PetscErrorCode PCBDDCSetNeumannBoundaries_BDDC(PC pc,IS NeumannBoundaries) 92 { 93 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 ierr = PetscObjectReference((PetscObject)NeumannBoundaries);CHKERRQ(ierr); 98 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 99 ierr = ISDuplicate(NeumannBoundaries,&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 100 ierr = ISCopy(NeumannBoundaries,pcbddc->NeumannBoundaries);CHKERRQ(ierr); 101 PetscFunctionReturn(0); 102 } 103 EXTERN_C_END 104 105 /* -------------------------------------------------------------------------- */ 106 #undef __FUNCT__ 107 #define __FUNCT__ "PCBDDCSetNeumannBoundaries" 108 /*@ 109 PCBDDCSetNeumannBoundaries - Set index set defining subdomain part of 110 Neumann boundaries for the global problem. 111 112 Logically Collective on PC 113 114 Input Parameters: 115 + pc - the preconditioning context 116 - NeumannBoundaries - index set defining the subdomain part of Neumann boundaries 117 118 Level: intermediate 119 120 Notes: 121 usage notes, perhaps an example 122 123 .seealso: PCBDDC 124 @*/ 125 PetscErrorCode PCBDDCSetNeumannBoundaries(PC pc,IS NeumannBoundaries) 126 { 127 PetscErrorCode ierr; 128 129 PetscFunctionBegin; 130 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 131 ierr = PetscTryMethod(pc,"PCBDDCSetNeumannBoundaries_C",(PC,IS),(pc,NeumannBoundaries));CHKERRQ(ierr); 132 PetscFunctionReturn(0); 133 } 134 /* -------------------------------------------------------------------------- */ 135 EXTERN_C_BEGIN 136 #undef __FUNCT__ 137 #define __FUNCT__ "PCBDDCGetNeumannBoundaries_BDDC" 138 static PetscErrorCode PCBDDCGetNeumannBoundaries_BDDC(PC pc,IS *NeumannBoundaries) 139 { 140 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 141 142 PetscFunctionBegin; 143 if(pcbddc->NeumannBoundaries) { 144 *NeumannBoundaries = pcbddc->NeumannBoundaries; 145 } else { 146 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Error in %s: Neumann boundaries not set!.\n",__FUNCT__); 147 } 148 PetscFunctionReturn(0); 149 } 150 EXTERN_C_END 151 152 /* -------------------------------------------------------------------------- */ 153 #undef __FUNCT__ 154 #define __FUNCT__ "PCBDDCGetNeumannBoundaries" 155 /*@ 156 PCBDDCGetNeumannBoundaries - Get index set defining subdomain part of 157 Neumann boundaries for the global problem. 158 159 Logically Collective on PC 160 161 Input Parameters: 162 + pc - the preconditioning context 163 164 Output Parameters: 165 + NeumannBoundaries - index set defining the subdomain part of Neumann boundaries 166 167 Level: intermediate 168 169 Notes: 170 usage notes, perhaps an example 171 172 .seealso: PCBDDC 173 @*/ 174 PetscErrorCode PCBDDCGetNeumannBoundaries(PC pc,IS *NeumannBoundaries) 175 { 176 PetscErrorCode ierr; 177 178 PetscFunctionBegin; 179 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 180 ierr = PetscUseMethod(pc,"PCBDDCGetNeumannBoundaries_C",(PC,IS*),(pc,NeumannBoundaries));CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183 184 #undef __FUNCT__ 185 #define __FUNCT__ "PCSetUp_BDDC" 186 /* -------------------------------------------------------------------------- */ 187 /* 188 PCSetUp_BDDC - Prepares for the use of the BDDC preconditioner 189 by setting data structures and options. 190 191 Input Parameter: 192 + pc - the preconditioner context 193 194 Application Interface Routine: PCSetUp() 195 196 Notes: 197 The interface routine PCSetUp() is not usually called directly by 198 the user, but instead is called by PCApply() if necessary. 199 */ 200 PetscErrorCode PCSetUp_BDDC(PC pc) 201 { 202 PetscErrorCode ierr; 203 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 204 PC_IS *pcis = (PC_IS*)(pc->data); 205 206 PetscFunctionBegin; 207 if (!pc->setupcalled) { 208 209 /* For BDDC we need to define a local "Neumann" to that defined in PCISSetup 210 So, we set to pcnone the Neumann problem of pcid in order to avoid unneeded computation 211 Also, we decide to directly build the (same) Dirichlet problem */ 212 ierr = PetscOptionsSetValue("-is_localN_pc_type","none");CHKERRQ(ierr); 213 ierr = PetscOptionsSetValue("-is_localD_pc_type","none");CHKERRQ(ierr); 214 /* Set up all the "iterative substructuring" common block */ 215 ierr = PCISSetUp(pc);CHKERRQ(ierr); 216 /* Destroy some PC_IS data which is not needed by BDDC */ 217 //if (pcis->ksp_D) {ierr = KSPDestroy(&pcis->ksp_D);CHKERRQ(ierr);} 218 //if (pcis->ksp_N) {ierr = KSPDestroy(&pcis->ksp_N);CHKERRQ(ierr);} 219 //if (pcis->vec2_B) {ierr = VecDestroy(&pcis->vec2_B);CHKERRQ(ierr);} 220 //if (pcis->vec3_B) {ierr = VecDestroy(&pcis->vec3_B);CHKERRQ(ierr);} 221 //pcis->ksp_D = 0; 222 //pcis->ksp_N = 0; 223 //pcis->vec2_B = 0; 224 //pcis->vec3_B = 0; 225 if(pcbddc->dbg_flag) { 226 ierr = PetscViewerASCIIGetStdout(((PetscObject)pc)->comm,&pcbddc->dbg_viewer);CHKERRQ(ierr); 227 ierr = PetscViewerASCIISynchronizedAllow(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr); 228 } 229 230 //TODO MOVE CODE FRAGMENT 231 PetscInt im_active=0; 232 if(pcis->n) im_active = 1; 233 ierr = MPI_Allreduce(&im_active,&pcbddc->active_procs,1,MPIU_INT,MPI_SUM,((PetscObject)pc)->comm);CHKERRQ(ierr); 234 //printf("Calling PCBDDC MANAGE with active procs %d (im_active = %d)\n",pcbddc->active_procs,im_active); 235 /* Set up grid quantities for BDDC */ 236 //TODO only active procs must call this -> remove synchronized print inside (the only point of sync) 237 ierr = PCBDDCManageLocalBoundaries(pc);CHKERRQ(ierr); 238 239 /* Create coarse and local stuffs used for evaluating action of preconditioner */ 240 ierr = PCBDDCCoarseSetUp(pc);CHKERRQ(ierr); 241 242 if ( !pcis->n_neigh ) pcis->ISLocalToGlobalMappingGetInfoWasCalled=PETSC_FALSE; //processes fakely involved in multilevel should not call ISLocalToGlobalMappingRestoreInfo 243 /* We no more need this matrix */ 244 //if (pcis->A_BB) {ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr);} 245 //pcis->A_BB = 0; 246 247 //printf("Using coarse problem type %d\n",pcbddc->coarse_problem_type); 248 //printf("Using coarse communications type %d\n",pcbddc->coarse_communications_type); 249 //printf("Using coarsening ratio %d\n",pcbddc->coarsening_ratio); 250 } 251 252 PetscFunctionReturn(0); 253 } 254 255 /* -------------------------------------------------------------------------- */ 256 /* 257 PCApply_BDDC - Applies the BDDC preconditioner to a vector. 258 259 Input Parameters: 260 . pc - the preconditioner context 261 . r - input vector (global) 262 263 Output Parameter: 264 . z - output vector (global) 265 266 Application Interface Routine: PCApply() 267 */ 268 #undef __FUNCT__ 269 #define __FUNCT__ "PCApply_BDDC" 270 PetscErrorCode PCApply_BDDC(PC pc,Vec r,Vec z) 271 { 272 PC_IS *pcis = (PC_IS*)(pc->data); 273 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 274 PetscErrorCode ierr; 275 PetscScalar one = 1.0; 276 PetscScalar m_one = -1.0; 277 278 /* This code is similar to that provided in nn.c for PCNN 279 NN interface preconditioner changed to BDDC 280 Added support for M_3 preconditioenr in the reference article (code is active if pcbddc->prec_type = PETSC_TRUE) */ 281 282 PetscFunctionBegin; 283 /* First Dirichlet solve */ 284 ierr = VecScatterBegin(pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 285 ierr = VecScatterEnd (pcis->global_to_D,r,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 286 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec1_D,pcis->vec2_D);CHKERRQ(ierr); 287 /* 288 Assembling right hand side for BDDC operator 289 - vec1_D for the Dirichlet part (if needed, i.e. prec_flag=PETSC_TRUE) 290 - the interface part of the global vector z 291 */ 292 //ierr = VecScatterBegin(pcis->global_to_B,r,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 293 //ierr = VecScatterEnd (pcis->global_to_B,r,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 294 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 295 ierr = MatMult(pcis->A_BI,pcis->vec2_D,pcis->vec1_B);CHKERRQ(ierr); 296 //ierr = MatMultAdd(pcis->A_BI,pcis->vec2_D,pcis->vec1_B,pcis->vec1_B);CHKERRQ(ierr); 297 if(pcbddc->prec_type) { ierr = MatMultAdd(pcis->A_II,pcis->vec2_D,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 298 ierr = VecScale(pcis->vec2_D,m_one);CHKERRQ(ierr); 299 ierr = VecCopy(r,z);CHKERRQ(ierr); 300 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 301 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 302 303 /* 304 Apply interface preconditioner 305 Results are stored in: 306 - vec1_D (if needed, i.e. with prec_type = PETSC_TRUE) 307 - the interface part of the global vector z 308 */ 309 ierr = PCBDDCApplyInterfacePreconditioner(pc,z);CHKERRQ(ierr); 310 311 /* Second Dirichlet solves and assembling of output */ 312 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 313 ierr = VecScatterEnd (pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 314 ierr = MatMult(pcis->A_IB,pcis->vec1_B,pcis->vec3_D);CHKERRQ(ierr); 315 if(pcbddc->prec_type) { ierr = MatMultAdd(pcis->A_II,pcis->vec1_D,pcis->vec3_D,pcis->vec3_D);CHKERRQ(ierr); } 316 ierr = KSPSolve(pcbddc->ksp_D,pcis->vec3_D,pcbddc->vec4_D);CHKERRQ(ierr); 317 ierr = VecScale(pcbddc->vec4_D,m_one);CHKERRQ(ierr); 318 if(pcbddc->prec_type) { ierr = VecAXPY (pcbddc->vec4_D,one,pcis->vec1_D);CHKERRQ(ierr); } 319 ierr = VecAXPY (pcis->vec2_D,one,pcbddc->vec4_D);CHKERRQ(ierr); 320 ierr = VecScatterBegin(pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 321 ierr = VecScatterEnd (pcis->global_to_D,pcis->vec2_D,z,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 322 323 PetscFunctionReturn(0); 324 325 } 326 327 /* -------------------------------------------------------------------------- */ 328 /* 329 PCBDDCApplyInterfacePreconditioner - Apply the BDDC preconditioner at the interface. 330 331 */ 332 #undef __FUNCT__ 333 #define __FUNCT__ "PCBDDCApplyInterfacePreconditioner" 334 static PetscErrorCode PCBDDCApplyInterfacePreconditioner(PC pc, Vec z) 335 { 336 PetscErrorCode ierr; 337 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 338 PC_IS* pcis = (PC_IS*) (pc->data); 339 PetscScalar zero = 0.0; 340 341 PetscFunctionBegin; 342 /* Get Local boundary and apply partition of unity */ 343 ierr = VecScatterBegin(pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 344 ierr = VecScatterEnd (pcis->global_to_B,z,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 345 ierr = VecPointwiseMult(pcis->vec1_B,pcis->D,pcis->vec1_B);CHKERRQ(ierr); 346 347 /* Application of PHI^T */ 348 ierr = MatMultTranspose(pcbddc->coarse_phi_B,pcis->vec1_B,pcbddc->vec1_P);CHKERRQ(ierr); 349 if(pcbddc->prec_type) { ierr = MatMultTransposeAdd(pcbddc->coarse_phi_D,pcis->vec1_D,pcbddc->vec1_P,pcbddc->vec1_P);CHKERRQ(ierr); } 350 351 /* Scatter data of coarse_rhs */ 352 if(pcbddc->coarse_rhs) ierr = VecSet(pcbddc->coarse_rhs,zero);CHKERRQ(ierr); 353 ierr = PCBDDCScatterCoarseDataBegin(pc,pcbddc->vec1_P,pcbddc->coarse_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 354 355 /* Local solution on R nodes */ 356 ierr = VecSet(pcbddc->vec1_R,zero);CHKERRQ(ierr); 357 ierr = VecScatterBegin(pcbddc->R_to_B,pcis->vec1_B,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 358 ierr = VecScatterEnd (pcbddc->R_to_B,pcis->vec1_B,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 359 if(pcbddc->prec_type) { 360 ierr = VecScatterBegin(pcbddc->R_to_D,pcis->vec1_D,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 361 ierr = VecScatterEnd (pcbddc->R_to_D,pcis->vec1_D,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 362 } 363 ierr = PCBDDCSolveSaddlePoint(pc);CHKERRQ(ierr); 364 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 365 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec2_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 366 ierr = VecScatterEnd (pcbddc->R_to_B,pcbddc->vec2_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 367 if(pcbddc->prec_type) { 368 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec2_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 369 ierr = VecScatterEnd (pcbddc->R_to_D,pcbddc->vec2_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 370 } 371 372 /* Coarse solution */ 373 ierr = PCBDDCScatterCoarseDataEnd(pc,pcbddc->vec1_P,pcbddc->coarse_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 374 if(pcbddc->coarse_rhs) ierr = KSPSolve(pcbddc->coarse_ksp,pcbddc->coarse_rhs,pcbddc->coarse_vec);CHKERRQ(ierr); 375 ierr = PCBDDCScatterCoarseDataBegin(pc,pcbddc->coarse_vec,pcbddc->vec1_P,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 376 ierr = PCBDDCScatterCoarseDataEnd (pc,pcbddc->coarse_vec,pcbddc->vec1_P,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 377 378 /* Sum contributions from two levels */ 379 /* Apply partition of unity and sum boundary values */ 380 ierr = MatMultAdd(pcbddc->coarse_phi_B,pcbddc->vec1_P,pcis->vec1_B,pcis->vec1_B);CHKERRQ(ierr); 381 ierr = VecPointwiseMult(pcis->vec1_B,pcis->D,pcis->vec1_B);CHKERRQ(ierr); 382 if(pcbddc->prec_type) { ierr = MatMultAdd(pcbddc->coarse_phi_D,pcbddc->vec1_P,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 383 ierr = VecSet(z,zero);CHKERRQ(ierr); 384 ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 385 ierr = VecScatterEnd (pcis->global_to_B,pcis->vec1_B,z,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 386 387 PetscFunctionReturn(0); 388 } 389 390 391 /* -------------------------------------------------------------------------- */ 392 /* 393 PCBDDCSolveSaddlePoint 394 395 */ 396 #undef __FUNCT__ 397 #define __FUNCT__ "PCBDDCSolveSaddlePoint" 398 static PetscErrorCode PCBDDCSolveSaddlePoint(PC pc) 399 { 400 PetscErrorCode ierr; 401 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 402 403 PetscFunctionBegin; 404 405 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec2_R);CHKERRQ(ierr); 406 if(pcbddc->n_constraints) { 407 ierr = MatMult(pcbddc->local_auxmat1,pcbddc->vec2_R,pcbddc->vec1_C);CHKERRQ(ierr); 408 ierr = MatMultAdd(pcbddc->local_auxmat2,pcbddc->vec1_C,pcbddc->vec2_R,pcbddc->vec2_R);CHKERRQ(ierr); 409 } 410 411 PetscFunctionReturn(0); 412 } 413 /* -------------------------------------------------------------------------- */ 414 /* 415 PCBDDCScatterCoarseDataBegin 416 417 */ 418 #undef __FUNCT__ 419 #define __FUNCT__ "PCBDDCScatterCoarseDataBegin" 420 static PetscErrorCode PCBDDCScatterCoarseDataBegin(PC pc,Vec vec_from, Vec vec_to, InsertMode imode, ScatterMode smode) 421 { 422 PetscErrorCode ierr; 423 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 424 425 PetscFunctionBegin; 426 427 switch(pcbddc->coarse_communications_type){ 428 case SCATTERS_BDDC: 429 ierr = VecScatterBegin(pcbddc->coarse_loc_to_glob,vec_from,vec_to,imode,smode);CHKERRQ(ierr); 430 break; 431 case GATHERS_BDDC: 432 break; 433 } 434 PetscFunctionReturn(0); 435 436 } 437 /* -------------------------------------------------------------------------- */ 438 /* 439 PCBDDCScatterCoarseDataEnd 440 441 */ 442 #undef __FUNCT__ 443 #define __FUNCT__ "PCBDDCScatterCoarseDataEnd" 444 static PetscErrorCode PCBDDCScatterCoarseDataEnd(PC pc,Vec vec_from, Vec vec_to, InsertMode imode, ScatterMode smode) 445 { 446 PetscErrorCode ierr; 447 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 448 PetscScalar* array_to; 449 PetscScalar* array_from; 450 MPI_Comm comm=((PetscObject)pc)->comm; 451 PetscInt i; 452 453 PetscFunctionBegin; 454 455 switch(pcbddc->coarse_communications_type){ 456 case SCATTERS_BDDC: 457 ierr = VecScatterEnd(pcbddc->coarse_loc_to_glob,vec_from,vec_to,imode,smode);CHKERRQ(ierr); 458 break; 459 case GATHERS_BDDC: 460 if(vec_from) VecGetArray(vec_from,&array_from); 461 if(vec_to) VecGetArray(vec_to,&array_to); 462 switch(pcbddc->coarse_problem_type){ 463 case SEQUENTIAL_BDDC: 464 if(smode == SCATTER_FORWARD) { 465 ierr = MPI_Gatherv(&array_from[0],pcbddc->local_primal_size,MPIU_SCALAR,&pcbddc->replicated_local_primal_values[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_SCALAR,0,comm);CHKERRQ(ierr); 466 if(vec_to) { 467 for(i=0;i<pcbddc->replicated_primal_size;i++) 468 array_to[pcbddc->replicated_local_primal_indices[i]]+=pcbddc->replicated_local_primal_values[i]; 469 } 470 } else { 471 if(vec_from) 472 for(i=0;i<pcbddc->replicated_primal_size;i++) 473 pcbddc->replicated_local_primal_values[i]=array_from[pcbddc->replicated_local_primal_indices[i]]; 474 ierr = MPI_Scatterv(&pcbddc->replicated_local_primal_values[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_SCALAR,&array_to[0],pcbddc->local_primal_size,MPIU_SCALAR,0,comm);CHKERRQ(ierr); 475 } 476 break; 477 case REPLICATED_BDDC: 478 if(smode == SCATTER_FORWARD) { 479 ierr = MPI_Allgatherv(&array_from[0],pcbddc->local_primal_size,MPIU_SCALAR,&pcbddc->replicated_local_primal_values[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_SCALAR,comm);CHKERRQ(ierr); 480 for(i=0;i<pcbddc->replicated_primal_size;i++) 481 array_to[pcbddc->replicated_local_primal_indices[i]]+=pcbddc->replicated_local_primal_values[i]; 482 } else { /* no communications needed for SCATTER_REVERSE since needed data is already present */ 483 for(i=0;i<pcbddc->local_primal_size;i++) 484 array_to[i]=array_from[pcbddc->local_primal_indices[i]]; 485 } 486 break; 487 case MULTILEVEL_BDDC: 488 break; 489 case PARALLEL_BDDC: 490 break; 491 } 492 if(vec_from) VecRestoreArray(vec_from,&array_from); 493 if(vec_to) VecRestoreArray(vec_to,&array_to); 494 break; 495 } 496 PetscFunctionReturn(0); 497 498 } 499 500 /* -------------------------------------------------------------------------- */ 501 /* 502 PCDestroy_BDDC - Destroys the private context for the NN preconditioner 503 that was created with PCCreate_BDDC(). 504 505 Input Parameter: 506 . pc - the preconditioner context 507 508 Application Interface Routine: PCDestroy() 509 */ 510 #undef __FUNCT__ 511 #define __FUNCT__ "PCDestroy_BDDC" 512 PetscErrorCode PCDestroy_BDDC(PC pc) 513 { 514 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 515 PetscErrorCode ierr; 516 517 PetscFunctionBegin; 518 /* free data created by PCIS */ 519 ierr = PCISDestroy(pc);CHKERRQ(ierr); 520 /* free BDDC data */ 521 ierr = VecDestroy(&pcbddc->coarse_vec);CHKERRQ(ierr); 522 ierr = VecDestroy(&pcbddc->coarse_rhs);CHKERRQ(ierr); 523 ierr = KSPDestroy(&pcbddc->coarse_ksp);CHKERRQ(ierr); 524 ierr = MatDestroy(&pcbddc->coarse_mat);CHKERRQ(ierr); 525 ierr = MatDestroy(&pcbddc->coarse_phi_B);CHKERRQ(ierr); 526 ierr = MatDestroy(&pcbddc->coarse_phi_D);CHKERRQ(ierr); 527 ierr = VecDestroy(&pcbddc->vec1_P);CHKERRQ(ierr); 528 ierr = VecDestroy(&pcbddc->vec1_C);CHKERRQ(ierr); 529 ierr = MatDestroy(&pcbddc->local_auxmat1);CHKERRQ(ierr); 530 ierr = MatDestroy(&pcbddc->local_auxmat2);CHKERRQ(ierr); 531 ierr = VecDestroy(&pcbddc->vec1_R);CHKERRQ(ierr); 532 ierr = VecDestroy(&pcbddc->vec2_R);CHKERRQ(ierr); 533 ierr = VecDestroy(&pcbddc->vec4_D);CHKERRQ(ierr); 534 ierr = VecScatterDestroy(&pcbddc->R_to_B);CHKERRQ(ierr); 535 ierr = VecScatterDestroy(&pcbddc->R_to_D);CHKERRQ(ierr); 536 ierr = VecScatterDestroy(&pcbddc->coarse_loc_to_glob);CHKERRQ(ierr); 537 ierr = KSPDestroy(&pcbddc->ksp_D);CHKERRQ(ierr); 538 ierr = KSPDestroy(&pcbddc->ksp_R);CHKERRQ(ierr); 539 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 540 ierr = PetscFree(pcbddc->vertices);CHKERRQ(ierr); 541 ierr = PetscFree(pcbddc->local_primal_indices);CHKERRQ(ierr); 542 ierr = PetscFree(pcbddc->replicated_local_primal_indices);CHKERRQ(ierr); 543 if (pcbddc->replicated_local_primal_values) { free(pcbddc->replicated_local_primal_values); } 544 ierr = PetscFree(pcbddc->local_primal_displacements);CHKERRQ(ierr); 545 ierr = PetscFree(pcbddc->local_primal_sizes);CHKERRQ(ierr); 546 if (pcbddc->n_constraints) { 547 ierr = PetscFree(pcbddc->indices_to_constraint[0]);CHKERRQ(ierr); 548 ierr = PetscFree(pcbddc->indices_to_constraint);CHKERRQ(ierr); 549 ierr = PetscFree(pcbddc->quadrature_constraint[0]);CHKERRQ(ierr); 550 ierr = PetscFree(pcbddc->quadrature_constraint);CHKERRQ(ierr); 551 ierr = PetscFree(pcbddc->sizes_of_constraint);CHKERRQ(ierr); 552 } 553 /* Free the private data structure that was hanging off the PC */ 554 ierr = PetscFree(pcbddc);CHKERRQ(ierr); 555 PetscFunctionReturn(0); 556 } 557 558 /* -------------------------------------------------------------------------- */ 559 /*MC 560 PCBDDC - Balancing Domain Decomposition by Constraints. 561 562 Options Database Keys: 563 . -pc_is_damp_fixed <fact> - 564 . -pc_is_remove_nullspace_fixed - 565 . -pc_is_set_damping_factor_floating <fact> - 566 . -pc_is_not_damp_floating - 567 568 Level: intermediate 569 570 Notes: The matrix used with this preconditioner must be of type MATIS 571 572 Unlike more 'conventional' interface preconditioners, this iterates over ALL the 573 degrees of freedom, NOT just those on the interface (this allows the use of approximate solvers 574 on the subdomains). 575 576 Options for the coarse grid preconditioner can be set with -pc_bddc_coarse_pc_xxx 577 Options for the Dirichlet subproblem can be set with -pc_bddc_localD_xxx 578 Options for the Neumann subproblem can be set with -pc_bddc_localN_xxx 579 580 Contributed by Stefano Zampini 581 582 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, MATIS 583 M*/ 584 EXTERN_C_BEGIN 585 #undef __FUNCT__ 586 #define __FUNCT__ "PCCreate_BDDC" 587 PetscErrorCode PCCreate_BDDC(PC pc) 588 { 589 PetscErrorCode ierr; 590 PC_BDDC *pcbddc; 591 592 PetscFunctionBegin; 593 /* Creates the private data structure for this preconditioner and attach it to the PC object. */ 594 ierr = PetscNewLog(pc,PC_BDDC,&pcbddc);CHKERRQ(ierr); 595 pc->data = (void*)pcbddc; 596 /* create PCIS data structure */ 597 ierr = PCISCreate(pc);CHKERRQ(ierr); 598 /* BDDC specific */ 599 pcbddc->coarse_vec = 0; 600 pcbddc->coarse_rhs = 0; 601 pcbddc->coarse_ksp = 0; 602 pcbddc->coarse_phi_B = 0; 603 pcbddc->coarse_phi_D = 0; 604 pcbddc->vec1_P = 0; 605 pcbddc->vec1_R = 0; 606 pcbddc->vec2_R = 0; 607 pcbddc->local_auxmat1 = 0; 608 pcbddc->local_auxmat2 = 0; 609 pcbddc->R_to_B = 0; 610 pcbddc->R_to_D = 0; 611 pcbddc->ksp_D = 0; 612 pcbddc->ksp_R = 0; 613 pcbddc->n_constraints = 0; 614 pcbddc->n_vertices = 0; 615 pcbddc->vertices = 0; 616 pcbddc->indices_to_constraint = 0; 617 pcbddc->quadrature_constraint = 0; 618 pcbddc->sizes_of_constraint = 0; 619 pcbddc->local_primal_indices = 0; 620 pcbddc->prec_type = PETSC_FALSE; 621 pcbddc->NeumannBoundaries = 0; 622 pcbddc->local_primal_sizes = 0; 623 pcbddc->local_primal_displacements = 0; 624 pcbddc->replicated_local_primal_indices = 0; 625 pcbddc->replicated_local_primal_values = 0; 626 pcbddc->coarse_loc_to_glob = 0; 627 pcbddc->dbg_flag = PETSC_FALSE; 628 pcbddc->vertices_flag = PETSC_FALSE; 629 pcbddc->constraints_flag = PETSC_FALSE; 630 pcbddc->faces_flag = PETSC_FALSE; 631 pcbddc->edges_flag = PETSC_FALSE; 632 pcbddc->coarsening_ratio = 8; 633 /* function pointers */ 634 pc->ops->apply = PCApply_BDDC; 635 pc->ops->applytranspose = 0; 636 pc->ops->setup = PCSetUp_BDDC; 637 pc->ops->destroy = PCDestroy_BDDC; 638 pc->ops->setfromoptions = PCSetFromOptions_BDDC; 639 pc->ops->view = 0; 640 pc->ops->applyrichardson = 0; 641 pc->ops->applysymmetricleft = 0; 642 pc->ops->applysymmetricright = 0; 643 /* composing function */ 644 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBDDCSetNeumannBoundaries_C","PCBDDCSetNeumannBoundaries_BDDC", 645 PCBDDCSetNeumannBoundaries_BDDC);CHKERRQ(ierr); 646 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBDDCGetNeumannBoundaries_C","PCBDDCGetNeumannBoundaries_BDDC", 647 PCBDDCGetNeumannBoundaries_BDDC);CHKERRQ(ierr); 648 ierr = PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBDDCSetCoarseProblemType_C","PCBDDCSetCoarseProblemType_BDDC", 649 PCBDDCSetCoarseProblemType_BDDC);CHKERRQ(ierr); 650 PetscFunctionReturn(0); 651 } 652 EXTERN_C_END 653 654 /* -------------------------------------------------------------------------- */ 655 /* 656 PCBDDCCoarseSetUp - 657 */ 658 #undef __FUNCT__ 659 #define __FUNCT__ "PCBDDCCoarseSetUp" 660 static PetscErrorCode PCBDDCCoarseSetUp(PC pc) 661 { 662 PetscErrorCode ierr; 663 664 PC_IS* pcis = (PC_IS*)(pc->data); 665 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 666 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 667 IS is_R_local; 668 IS is_V_local; 669 IS is_C_local; 670 IS is_aux1; 671 IS is_aux2; 672 const VecType impVecType; 673 const MatType impMatType; 674 PetscInt n_R=0; 675 PetscInt n_D=0; 676 PetscInt n_B=0; 677 PetscMPIInt totprocs; 678 PetscScalar zero=0.0; 679 PetscScalar one=1.0; 680 PetscScalar m_one=-1.0; 681 PetscScalar* array; 682 PetscScalar *coarse_submat_vals; 683 PetscInt *idx_R_local; 684 PetscInt *idx_V_B; 685 PetscScalar *coarsefunctions_errors; 686 PetscScalar *constraints_errors; 687 /* auxiliary indices */ 688 PetscInt s,i,j,k; 689 /* for verbose output of bddc */ 690 PetscViewer viewer=pcbddc->dbg_viewer; 691 PetscBool dbg_flag=pcbddc->dbg_flag; 692 693 PetscFunctionBegin; 694 695 /* Set Non-overlapping dimensions */ 696 n_B = pcis->n_B; n_D = pcis->n - n_B; 697 /* Set local primal size */ 698 pcbddc->local_primal_size = pcbddc->n_vertices + pcbddc->n_constraints; 699 /* Dohrmann's notation: dofs splitted in R (Remaining: all dofs but the vertices) and V (Vertices) */ 700 ierr = VecSet(pcis->vec1_N,one);CHKERRQ(ierr); 701 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 702 for (i=0;i<pcbddc->n_vertices;i++) { array[ pcbddc->vertices[i] ] = zero; } 703 ierr = PetscMalloc(( pcis->n - pcbddc->n_vertices )*sizeof(PetscInt),&idx_R_local);CHKERRQ(ierr); 704 for (i=0, n_R=0; i<pcis->n; i++) { if (array[i] == one) { idx_R_local[n_R] = i; n_R++; } } 705 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 706 if(dbg_flag) { 707 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 708 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 709 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d local dimensions\n",PetscGlobalRank);CHKERRQ(ierr); 710 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local_size = %d, dirichlet_size = %d, boundary_size = %d\n",pcis->n,n_D,n_B);CHKERRQ(ierr); 711 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"r_size = %d, v_size = %d, constraints = %d, local_primal_size = %d\n",n_R,pcbddc->n_vertices,pcbddc->n_constraints,pcbddc->local_primal_size);CHKERRQ(ierr); 712 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 713 } 714 /* Allocate needed vectors */ 715 /* Set Mat type for local matrices needed by BDDC precondtioner */ 716 // ierr = MatGetType(matis->A,&impMatType);CHKERRQ(ierr); 717 //ierr = VecGetType(pcis->vec1_N,&impVecType);CHKERRQ(ierr); 718 impMatType = MATSEQDENSE; 719 impVecType = VECSEQ; 720 ierr = VecDuplicate(pcis->vec1_D,&pcbddc->vec4_D);CHKERRQ(ierr); 721 ierr = VecDuplicate(pcis->vec1_N,&pcis->vec2_N);CHKERRQ(ierr); 722 ierr = VecCreate(PETSC_COMM_SELF,&pcbddc->vec1_R);CHKERRQ(ierr); 723 ierr = VecSetSizes(pcbddc->vec1_R,n_R,n_R);CHKERRQ(ierr); 724 ierr = VecSetType(pcbddc->vec1_R,impVecType);CHKERRQ(ierr); 725 ierr = VecDuplicate(pcbddc->vec1_R,&pcbddc->vec2_R);CHKERRQ(ierr); 726 ierr = VecCreate(PETSC_COMM_SELF,&pcbddc->vec1_P);CHKERRQ(ierr); 727 ierr = VecSetSizes(pcbddc->vec1_P,pcbddc->local_primal_size,pcbddc->local_primal_size);CHKERRQ(ierr); 728 ierr = VecSetType(pcbddc->vec1_P,impVecType);CHKERRQ(ierr); 729 730 /* Creating some index sets needed */ 731 /* For submatrices */ 732 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_R,idx_R_local,PETSC_COPY_VALUES,&is_R_local);CHKERRQ(ierr); 733 if(pcbddc->n_vertices) { ierr = ISCreateGeneral(PETSC_COMM_SELF,pcbddc->n_vertices,pcbddc->vertices,PETSC_COPY_VALUES,&is_V_local);CHKERRQ(ierr); } 734 if(pcbddc->n_constraints) { ierr = ISCreateStride (PETSC_COMM_SELF,pcbddc->n_constraints,pcbddc->n_vertices,1,&is_C_local);CHKERRQ(ierr); } 735 /* For VecScatters pcbddc->R_to_B and (optionally) pcbddc->R_to_D */ 736 { 737 PetscInt *aux_array1; 738 PetscInt *aux_array2; 739 PetscScalar value; 740 741 ierr = PetscMalloc( (pcis->n_B-pcbddc->n_vertices)*sizeof(PetscInt),&aux_array1);CHKERRQ(ierr); 742 ierr = PetscMalloc( (pcis->n_B-pcbddc->n_vertices)*sizeof(PetscInt),&aux_array2);CHKERRQ(ierr); 743 744 ierr = VecSet(pcis->vec1_global,zero);CHKERRQ(ierr); 745 ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 746 ierr = VecScatterEnd (matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 747 ierr = VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 748 ierr = VecScatterEnd (matis->ctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 749 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 750 ierr = VecScatterEnd (pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 751 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 752 for (i=0, s=0; i<n_R; i++) { if (array[idx_R_local[i]] > one) { aux_array1[s] = i; s++; } } 753 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 754 ierr = ISCreateGeneral(PETSC_COMM_SELF,s,aux_array1,PETSC_COPY_VALUES,&is_aux1);CHKERRQ(ierr); 755 ierr = VecGetArray(pcis->vec1_B,&array);CHKERRQ(ierr); 756 for (i=0, s=0; i<n_B; i++) { if (array[i] > one) { aux_array2[s] = i; s++; } } 757 ierr = VecRestoreArray(pcis->vec1_B,&array);CHKERRQ(ierr); 758 ierr = ISCreateGeneral(PETSC_COMM_SELF,s,aux_array2,PETSC_COPY_VALUES,&is_aux2);CHKERRQ(ierr); 759 ierr = VecScatterCreate(pcbddc->vec1_R,is_aux1,pcis->vec1_B,is_aux2,&pcbddc->R_to_B);CHKERRQ(ierr); 760 ierr = PetscFree(aux_array1);CHKERRQ(ierr); 761 ierr = PetscFree(aux_array2);CHKERRQ(ierr); 762 ierr = ISDestroy(&is_aux1);CHKERRQ(ierr); 763 ierr = ISDestroy(&is_aux2);CHKERRQ(ierr); 764 765 if(pcbddc->prec_type || dbg_flag ) { 766 ierr = PetscMalloc(n_D*sizeof(PetscInt),&aux_array1);CHKERRQ(ierr); 767 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 768 for (i=0, s=0; i<n_R; i++) { if (array[idx_R_local[i]] == one) { aux_array1[s] = i; s++; } } 769 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 770 ierr = ISCreateGeneral(PETSC_COMM_SELF,s,aux_array1,PETSC_COPY_VALUES,&is_aux1);CHKERRQ(ierr); 771 ierr = VecScatterCreate(pcbddc->vec1_R,is_aux1,pcis->vec1_D,(IS)0,&pcbddc->R_to_D);CHKERRQ(ierr); 772 ierr = PetscFree(aux_array1);CHKERRQ(ierr); 773 ierr = ISDestroy(&is_aux1);CHKERRQ(ierr); 774 } 775 776 /* Check scatters */ 777 if(dbg_flag) { 778 779 Vec vec_aux; 780 781 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 782 ierr = PetscViewerASCIIPrintf(viewer,"Checking pcbddc->R_to_B scatter\n");CHKERRQ(ierr); 783 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 784 ierr = VecSetRandom(pcbddc->vec1_R,PETSC_NULL);CHKERRQ(ierr); 785 ierr = VecSetRandom(pcis->vec1_B,PETSC_NULL);CHKERRQ(ierr); 786 ierr = VecDuplicate(pcbddc->vec1_R,&vec_aux);CHKERRQ(ierr); 787 ierr = VecCopy(pcbddc->vec1_R,vec_aux);CHKERRQ(ierr); 788 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 789 ierr = VecScatterEnd (pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 790 ierr = VecScatterBegin(pcbddc->R_to_B,pcis->vec1_B,vec_aux,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 791 ierr = VecScatterEnd (pcbddc->R_to_B,pcis->vec1_B,vec_aux,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 792 ierr = VecAXPY(vec_aux,m_one,pcbddc->vec1_R);CHKERRQ(ierr); 793 ierr = VecNorm(vec_aux,NORM_INFINITY,&value);CHKERRQ(ierr); 794 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d R_to_B FORWARD error = % 1.14e\n",PetscGlobalRank,value);CHKERRQ(ierr); 795 ierr = VecDestroy(&vec_aux);CHKERRQ(ierr); 796 797 ierr = VecSetRandom(pcbddc->vec1_R,PETSC_NULL);CHKERRQ(ierr); 798 ierr = VecSetRandom(pcis->vec1_B,PETSC_NULL);CHKERRQ(ierr); 799 ierr = VecDuplicate(pcis->vec1_B,&vec_aux);CHKERRQ(ierr); 800 ierr = VecCopy(pcis->vec1_B,vec_aux);CHKERRQ(ierr); 801 ierr = VecScatterBegin(pcbddc->R_to_B,pcis->vec1_B,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 802 ierr = VecScatterEnd (pcbddc->R_to_B,pcis->vec1_B,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 803 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,vec_aux,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 804 ierr = VecScatterEnd (pcbddc->R_to_B,pcbddc->vec1_R,vec_aux,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 805 ierr = VecAXPY(vec_aux,m_one,pcis->vec1_B);CHKERRQ(ierr); 806 ierr = VecNorm(vec_aux,NORM_INFINITY,&value);CHKERRQ(ierr); 807 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d R_to_B REVERSE error = % 1.14e\n",PetscGlobalRank,value);CHKERRQ(ierr); 808 ierr = VecDestroy(&vec_aux);CHKERRQ(ierr); 809 810 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 811 ierr = PetscViewerASCIIPrintf(viewer,"Checking pcbddc->R_to_D scatter\n");CHKERRQ(ierr); 812 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 813 814 ierr = VecSetRandom(pcbddc->vec1_R,PETSC_NULL);CHKERRQ(ierr); 815 ierr = VecSetRandom(pcis->vec1_D,PETSC_NULL);CHKERRQ(ierr); 816 ierr = VecDuplicate(pcbddc->vec1_R,&vec_aux);CHKERRQ(ierr); 817 ierr = VecCopy(pcbddc->vec1_R,vec_aux);CHKERRQ(ierr); 818 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 819 ierr = VecScatterEnd (pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 820 ierr = VecScatterBegin(pcbddc->R_to_D,pcis->vec1_D,vec_aux,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 821 ierr = VecScatterEnd (pcbddc->R_to_D,pcis->vec1_D,vec_aux,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 822 ierr = VecAXPY(vec_aux,m_one,pcbddc->vec1_R);CHKERRQ(ierr); 823 ierr = VecNorm(vec_aux,NORM_INFINITY,&value);CHKERRQ(ierr); 824 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d R_to_D FORWARD error = % 1.14e\n",PetscGlobalRank,value);CHKERRQ(ierr); 825 ierr = VecDestroy(&vec_aux);CHKERRQ(ierr); 826 827 ierr = VecSetRandom(pcbddc->vec1_R,PETSC_NULL);CHKERRQ(ierr); 828 ierr = VecSetRandom(pcis->vec1_D,PETSC_NULL);CHKERRQ(ierr); 829 ierr = VecDuplicate(pcis->vec1_D,&vec_aux);CHKERRQ(ierr); 830 ierr = VecCopy(pcis->vec1_D,vec_aux);CHKERRQ(ierr); 831 ierr = VecScatterBegin(pcbddc->R_to_D,pcis->vec1_D,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 832 ierr = VecScatterEnd (pcbddc->R_to_D,pcis->vec1_D,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 833 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,vec_aux,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 834 ierr = VecScatterEnd (pcbddc->R_to_D,pcbddc->vec1_R,vec_aux,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 835 ierr = VecAXPY(vec_aux,m_one,pcis->vec1_D);CHKERRQ(ierr); 836 ierr = VecNorm(vec_aux,NORM_INFINITY,&value);CHKERRQ(ierr); 837 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d R_to_D REVERSE error = % 1.14e\n",PetscGlobalRank,value);CHKERRQ(ierr); 838 ierr = VecDestroy(&vec_aux);CHKERRQ(ierr); 839 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 840 841 } 842 } 843 844 /* vertices in boundary numbering */ 845 if(pcbddc->n_vertices) { 846 ierr = VecSet(pcis->vec1_N,m_one);CHKERRQ(ierr); 847 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 848 for (i=0; i<pcbddc->n_vertices; i++) { array[ pcbddc->vertices[i] ] = i; } 849 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 850 ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 851 ierr = VecScatterEnd (pcis->N_to_B,pcis->vec1_N,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 852 ierr = PetscMalloc(pcbddc->n_vertices*sizeof(PetscInt),&idx_V_B);CHKERRQ(ierr); 853 ierr = VecGetArray(pcis->vec1_B,&array);CHKERRQ(ierr); 854 //printf("vertices in boundary numbering\n"); 855 for (i=0; i<pcbddc->n_vertices; i++) { 856 s=0; 857 while (array[s] != i ) {s++;} 858 idx_V_B[i]=s; 859 //printf("idx_V_B[%d]=%d\n",i,s); 860 } 861 ierr = VecRestoreArray(pcis->vec1_B,&array);CHKERRQ(ierr); 862 } 863 864 865 /* Creating PC contexts for local Dirichlet and Neumann problems */ 866 { 867 Mat A_RR; 868 PC pc_temp; 869 /* Matrix for Dirichlet problem is A_II -> we already have it from pcis.c code */ 870 ierr = KSPCreate(PETSC_COMM_SELF,&pcbddc->ksp_D);CHKERRQ(ierr); 871 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_D,(PetscObject)pc,1);CHKERRQ(ierr); 872 ierr = KSPSetOperators(pcbddc->ksp_D,pcis->A_II,pcis->A_II,SAME_PRECONDITIONER);CHKERRQ(ierr); 873 ierr = KSPSetType(pcbddc->ksp_D,KSPPREONLY);CHKERRQ(ierr); 874 //ierr = KSPSetOptionsPrefix(pcbddc->pc_D,"pc_bddc_localD_");CHKERRQ(ierr); 875 /* default */ 876 ierr = KSPGetPC(pcbddc->ksp_D,&pc_temp);CHKERRQ(ierr); 877 ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr); 878 /* Allow user's customization */ 879 ierr = KSPSetFromOptions(pcbddc->ksp_D);CHKERRQ(ierr); 880 /* Set Up KSP for Dirichlet problem of BDDC */ 881 ierr = KSPSetUp(pcbddc->ksp_D);CHKERRQ(ierr); 882 /* Matrix for Neumann problem is A_RR -> we need to create it */ 883 ierr = MatGetSubMatrix(matis->A,is_R_local,is_R_local,MAT_INITIAL_MATRIX,&A_RR);CHKERRQ(ierr); 884 ierr = KSPCreate(PETSC_COMM_SELF,&pcbddc->ksp_R);CHKERRQ(ierr); 885 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_R,(PetscObject)pc,1);CHKERRQ(ierr); 886 ierr = KSPSetOperators(pcbddc->ksp_R,A_RR,A_RR,SAME_PRECONDITIONER);CHKERRQ(ierr); 887 ierr = KSPSetType(pcbddc->ksp_R,KSPPREONLY);CHKERRQ(ierr); 888 //ierr = PCSetOptionsPrefix(pcbddc->pc_R,"pc_bddc_localN_");CHKERRQ(ierr); 889 /* default */ 890 ierr = KSPGetPC(pcbddc->ksp_R,&pc_temp);CHKERRQ(ierr); 891 ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr); 892 /* Allow user's customization */ 893 ierr = KSPSetFromOptions(pcbddc->ksp_R);CHKERRQ(ierr); 894 /* Set Up KSP for Neumann problem of BDDC */ 895 ierr = KSPSetUp(pcbddc->ksp_R);CHKERRQ(ierr); 896 /* check Neumann solve */ 897 if(pcbddc->dbg_flag) { 898 Vec temp_vec; 899 PetscScalar value; 900 901 ierr = VecDuplicate(pcbddc->vec1_R,&temp_vec);CHKERRQ(ierr); 902 ierr = VecSetRandom(pcbddc->vec1_R,PETSC_NULL);CHKERRQ(ierr); 903 ierr = MatMult(A_RR,pcbddc->vec1_R,pcbddc->vec2_R);CHKERRQ(ierr); 904 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec2_R,temp_vec);CHKERRQ(ierr); 905 ierr = VecAXPY(temp_vec,m_one,pcbddc->vec1_R);CHKERRQ(ierr); 906 ierr = VecNorm(temp_vec,NORM_INFINITY,&value);CHKERRQ(ierr); 907 ierr = VecDestroy(&temp_vec);CHKERRQ(ierr); 908 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 909 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 910 ierr = PetscViewerASCIIPrintf(viewer,"Checking solution of Neumann problem\n");CHKERRQ(ierr); 911 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d infinity error for Neumann solve = % 1.14e \n",PetscGlobalRank,value);CHKERRQ(ierr); 912 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 913 } 914 /* free Neumann problem's matrix */ 915 ierr = MatDestroy(&A_RR);CHKERRQ(ierr); 916 } 917 918 /* Assemble all remaining stuff needed to apply BDDC */ 919 { 920 Mat A_RV,A_VR,A_VV; 921 Mat M1,M2; 922 Mat C_CR; 923 Mat CMAT,AUXMAT; 924 Vec vec1_C; 925 Vec vec2_C; 926 Vec vec1_V; 927 Vec vec2_V; 928 PetscInt *nnz; 929 PetscInt *auxindices; 930 PetscInt index; 931 PetscScalar* array2; 932 MatFactorInfo matinfo; 933 934 /* Allocating some extra storage just to be safe */ 935 ierr = PetscMalloc (pcis->n*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 936 ierr = PetscMalloc (pcis->n*sizeof(PetscInt),&auxindices);CHKERRQ(ierr); 937 for(i=0;i<pcis->n;i++) {auxindices[i]=i;} 938 939 /* some work vectors on vertices and/or constraints */ 940 if(pcbddc->n_vertices) { 941 ierr = VecCreate(PETSC_COMM_SELF,&vec1_V);CHKERRQ(ierr); 942 ierr = VecSetSizes(vec1_V,pcbddc->n_vertices,pcbddc->n_vertices);CHKERRQ(ierr); 943 ierr = VecSetType(vec1_V,impVecType);CHKERRQ(ierr); 944 ierr = VecDuplicate(vec1_V,&vec2_V);CHKERRQ(ierr); 945 } 946 if(pcbddc->n_constraints) { 947 ierr = VecCreate(PETSC_COMM_SELF,&vec1_C);CHKERRQ(ierr); 948 ierr = VecSetSizes(vec1_C,pcbddc->n_constraints,pcbddc->n_constraints);CHKERRQ(ierr); 949 ierr = VecSetType(vec1_C,impVecType);CHKERRQ(ierr); 950 ierr = VecDuplicate(vec1_C,&vec2_C);CHKERRQ(ierr); 951 ierr = VecDuplicate(vec1_C,&pcbddc->vec1_C);CHKERRQ(ierr); 952 } 953 /* Create C matrix [I 0; 0 const] */ 954 ierr = MatCreate(PETSC_COMM_SELF,&CMAT);CHKERRQ(ierr); 955 ierr = MatSetType(CMAT,MATSEQAIJ);CHKERRQ(ierr); 956 ierr = MatSetSizes(CMAT,pcbddc->local_primal_size,pcis->n,pcbddc->local_primal_size,pcis->n);CHKERRQ(ierr); 957 /* nonzeros */ 958 for(i=0;i<pcbddc->n_vertices;i++) { nnz[i]= 1; } 959 for(i=0;i<pcbddc->n_constraints;i++) { nnz[i+pcbddc->n_vertices]=pcbddc->sizes_of_constraint[i];} 960 ierr = MatSeqAIJSetPreallocation(CMAT,0,nnz);CHKERRQ(ierr); 961 for(i=0;i<pcbddc->n_vertices;i++) { 962 ierr = MatSetValue(CMAT,i,pcbddc->vertices[i],1.0,INSERT_VALUES);CHKERRQ(ierr); 963 } 964 for(i=0;i<pcbddc->n_constraints;i++) { 965 index=i+pcbddc->n_vertices; 966 ierr = MatSetValues(CMAT,1,&index,pcbddc->sizes_of_constraint[i],pcbddc->indices_to_constraint[i],pcbddc->quadrature_constraint[i],INSERT_VALUES);CHKERRQ(ierr); 967 } 968 //if(dbg_flag) printf("CMAT assembling\n"); 969 ierr = MatAssemblyBegin(CMAT,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 970 ierr = MatAssemblyEnd(CMAT,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 971 //ierr = MatView(CMAT,PETSC_VIEWER_STDOUT_SELF); 972 973 /* Precompute stuffs needed for preprocessing and application of BDDC*/ 974 975 if(pcbddc->n_constraints) { 976 /* some work vectors */ 977 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->local_auxmat2);CHKERRQ(ierr); 978 ierr = MatSetSizes(pcbddc->local_auxmat2,n_R,pcbddc->n_constraints,n_R,pcbddc->n_constraints);CHKERRQ(ierr); 979 ierr = MatSetType(pcbddc->local_auxmat2,impMatType);CHKERRQ(ierr); 980 981 /* Assemble local_auxmat2 = - A_{RR}^{-1} C^T_{CR} needed by BDDC application */ 982 for(i=0;i<pcbddc->n_constraints;i++) { 983 ierr = VecSet(pcis->vec1_N,zero);CHKERRQ(ierr); 984 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 985 for(j=0;j<pcbddc->sizes_of_constraint[i];j++) { array[ pcbddc->indices_to_constraint[i][j] ] = - pcbddc->quadrature_constraint[i][j]; } 986 ierr = VecGetArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 987 for(j=0;j<n_R;j++) { array2[j] = array[ idx_R_local[j] ]; } 988 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 989 ierr = VecRestoreArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 990 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec2_R);CHKERRQ(ierr); 991 ierr = VecGetArray(pcbddc->vec2_R,&array);CHKERRQ(ierr); 992 index=i; 993 ierr = MatSetValues(pcbddc->local_auxmat2,n_R,auxindices,1,&index,array,INSERT_VALUES);CHKERRQ(ierr); 994 ierr = VecRestoreArray(pcbddc->vec2_R,&array);CHKERRQ(ierr); 995 } 996 //if(dbg_flag) printf("pcbddc->local_auxmat2 assembling\n"); 997 ierr = MatAssemblyBegin(pcbddc->local_auxmat2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 998 ierr = MatAssemblyEnd (pcbddc->local_auxmat2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 999 1000 /* Create Constraint matrix on R nodes: C_{CR} */ 1001 ierr = MatGetSubMatrix(CMAT,is_C_local,is_R_local,MAT_INITIAL_MATRIX,&C_CR);CHKERRQ(ierr); 1002 ierr = ISDestroy(&is_C_local);CHKERRQ(ierr); 1003 1004 /* Assemble AUXMAT = ( LUFactor )( -C_{CR} A_{RR}^{-1} C^T_{CR} )^{-1} */ 1005 ierr = MatMatMult(C_CR,pcbddc->local_auxmat2,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AUXMAT);CHKERRQ(ierr); 1006 ierr = MatFactorInfoInitialize(&matinfo);CHKERRQ(ierr); 1007 ierr = ISCreateStride (PETSC_COMM_SELF,pcbddc->n_constraints,0,1,&is_aux1);CHKERRQ(ierr); 1008 ierr = MatLUFactor(AUXMAT,is_aux1,is_aux1,&matinfo);CHKERRQ(ierr); 1009 ierr = ISDestroy(&is_aux1);CHKERRQ(ierr); 1010 1011 /* Assemble explicitly M1 = ( C_{CR} A_{RR}^{-1} C^T_{CR} )^{-1} needed in preproc (should be dense) */ 1012 ierr = MatCreate(PETSC_COMM_SELF,&M1);CHKERRQ(ierr); 1013 ierr = MatSetSizes(M1,pcbddc->n_constraints,pcbddc->n_constraints,pcbddc->n_constraints,pcbddc->n_constraints);CHKERRQ(ierr); 1014 ierr = MatSetType(M1,impMatType);CHKERRQ(ierr); 1015 for(i=0;i<pcbddc->n_constraints;i++) { 1016 ierr = VecSet(vec1_C,zero);CHKERRQ(ierr); 1017 ierr = VecSetValue(vec1_C,i,one,INSERT_VALUES);CHKERRQ(ierr); 1018 ierr = VecAssemblyBegin(vec1_C);CHKERRQ(ierr); 1019 ierr = VecAssemblyEnd(vec1_C);CHKERRQ(ierr); 1020 ierr = MatSolve(AUXMAT,vec1_C,vec2_C);CHKERRQ(ierr); 1021 ierr = VecScale(vec2_C,m_one);CHKERRQ(ierr); 1022 ierr = VecGetArray(vec2_C,&array);CHKERRQ(ierr); 1023 index=i; 1024 ierr = MatSetValues(M1,pcbddc->n_constraints,auxindices,1,&index,array,INSERT_VALUES);CHKERRQ(ierr); 1025 ierr = VecRestoreArray(vec2_C,&array);CHKERRQ(ierr); 1026 } 1027 //if(dbg_flag) printf("M1 assembling\n"); 1028 ierr = MatAssemblyBegin(M1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1029 ierr = MatAssemblyEnd (M1,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1030 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 1031 /* Assemble local_auxmat1 = M1*C_{CR} needed by BDDC application in KSP and in preproc */ 1032 //if(dbg_flag) printf("pcbddc->local_auxmat1 computing and assembling\n"); 1033 ierr = MatMatMult(M1,C_CR,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pcbddc->local_auxmat1);CHKERRQ(ierr); 1034 1035 } 1036 1037 /* Get submatrices from subdomain matrix */ 1038 if(pcbddc->n_vertices){ 1039 ierr = MatGetSubMatrix(matis->A,is_R_local,is_V_local,MAT_INITIAL_MATRIX,&A_RV);CHKERRQ(ierr); 1040 ierr = MatGetSubMatrix(matis->A,is_V_local,is_R_local,MAT_INITIAL_MATRIX,&A_VR);CHKERRQ(ierr); 1041 ierr = MatGetSubMatrix(matis->A,is_V_local,is_V_local,MAT_INITIAL_MATRIX,&A_VV);CHKERRQ(ierr); 1042 /* Assemble M2 = A_RR^{-1}A_RV */ 1043 ierr = MatCreate(PETSC_COMM_SELF,&M2);CHKERRQ(ierr); 1044 ierr = MatSetSizes(M2,n_R,pcbddc->n_vertices,n_R,pcbddc->n_vertices);CHKERRQ(ierr); 1045 ierr = MatSetType(M2,impMatType);CHKERRQ(ierr); 1046 for(i=0;i<pcbddc->n_vertices;i++) { 1047 ierr = VecSet(vec1_V,zero);CHKERRQ(ierr); 1048 ierr = VecSetValue(vec1_V,i,one,INSERT_VALUES);CHKERRQ(ierr); 1049 ierr = VecAssemblyBegin(vec1_V);CHKERRQ(ierr); 1050 ierr = VecAssemblyEnd(vec1_V);CHKERRQ(ierr); 1051 ierr = MatMult(A_RV,vec1_V,pcbddc->vec1_R);CHKERRQ(ierr); 1052 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec2_R);CHKERRQ(ierr); 1053 index=i; 1054 ierr = VecGetArray(pcbddc->vec2_R,&array);CHKERRQ(ierr); 1055 ierr = MatSetValues(M2,n_R,auxindices,1,&index,array,INSERT_VALUES);CHKERRQ(ierr); 1056 ierr = VecRestoreArray(pcbddc->vec2_R,&array);CHKERRQ(ierr); 1057 } 1058 //if(dbg_flag) printf("M2 assembling\n"); 1059 ierr = MatAssemblyBegin(M2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1060 ierr = MatAssemblyEnd (M2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1061 } 1062 1063 /* Matrix of coarse basis functions (local) */ 1064 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_phi_B);CHKERRQ(ierr); 1065 ierr = MatSetSizes(pcbddc->coarse_phi_B,n_B,pcbddc->local_primal_size,n_B,pcbddc->local_primal_size);CHKERRQ(ierr); 1066 ierr = MatSetType(pcbddc->coarse_phi_B,impMatType);CHKERRQ(ierr); 1067 if(pcbddc->prec_type || dbg_flag ) { 1068 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_phi_D);CHKERRQ(ierr); 1069 ierr = MatSetSizes(pcbddc->coarse_phi_D,n_D,pcbddc->local_primal_size,n_D,pcbddc->local_primal_size);CHKERRQ(ierr); 1070 ierr = MatSetType(pcbddc->coarse_phi_D,impMatType);CHKERRQ(ierr); 1071 } 1072 1073 /* Subdomain contribution (Non-overlapping) to coarse matrix */ 1074 if(dbg_flag) { 1075 ierr = PetscMalloc( pcbddc->local_primal_size*sizeof(PetscScalar),&coarsefunctions_errors);CHKERRQ(ierr); 1076 ierr = PetscMalloc( pcbddc->local_primal_size*sizeof(PetscScalar),&constraints_errors);CHKERRQ(ierr); 1077 } 1078 ierr = PetscMalloc ((pcbddc->local_primal_size)*(pcbddc->local_primal_size)*sizeof(PetscScalar),&coarse_submat_vals);CHKERRQ(ierr); 1079 1080 /* We are now ready to evaluate coarse basis functions and subdomain contribution to coarse problem */ 1081 for(i=0;i<pcbddc->n_vertices;i++){ 1082 ierr = VecSet(vec1_V,zero);CHKERRQ(ierr); 1083 ierr = VecSetValue(vec1_V,i,one,INSERT_VALUES);CHKERRQ(ierr); 1084 ierr = VecAssemblyBegin(vec1_V);CHKERRQ(ierr); 1085 ierr = VecAssemblyEnd(vec1_V);CHKERRQ(ierr); 1086 /* solution of saddle point problem */ 1087 ierr = MatMult(M2,vec1_V,pcbddc->vec1_R);CHKERRQ(ierr); 1088 ierr = VecScale(pcbddc->vec1_R,m_one);CHKERRQ(ierr); 1089 if(pcbddc->n_constraints) { 1090 ierr = MatMult(pcbddc->local_auxmat1,pcbddc->vec1_R,vec1_C);CHKERRQ(ierr); 1091 ierr = MatMultAdd(pcbddc->local_auxmat2,vec1_C,pcbddc->vec1_R,pcbddc->vec1_R);CHKERRQ(ierr); 1092 ierr = VecScale(vec1_C,m_one);CHKERRQ(ierr); 1093 } 1094 ierr = MatMult(A_VR,pcbddc->vec1_R,vec2_V);CHKERRQ(ierr); 1095 ierr = MatMultAdd(A_VV,vec1_V,vec2_V,vec2_V);CHKERRQ(ierr); 1096 1097 /* Set values in coarse basis function and subdomain part of coarse_mat */ 1098 /* coarse basis functions */ 1099 index=i; 1100 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1101 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1102 ierr = VecScatterEnd (pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1103 ierr = VecGetArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1104 ierr = MatSetValues(pcbddc->coarse_phi_B,n_B,auxindices,1,&index,array,INSERT_VALUES);CHKERRQ(ierr); 1105 ierr = VecRestoreArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1106 ierr = MatSetValue(pcbddc->coarse_phi_B,idx_V_B[i],i,one,INSERT_VALUES);CHKERRQ(ierr); 1107 if( pcbddc->prec_type || dbg_flag ) { 1108 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1109 ierr = VecScatterEnd(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1110 ierr = VecGetArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1111 ierr = MatSetValues(pcbddc->coarse_phi_D,n_D,auxindices,1,&index,array,INSERT_VALUES);CHKERRQ(ierr); 1112 ierr = VecRestoreArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1113 } 1114 /* subdomain contribution to coarse matrix */ 1115 ierr = VecGetArray(vec2_V,&array);CHKERRQ(ierr); 1116 for(j=0;j<pcbddc->n_vertices;j++) {coarse_submat_vals[i*pcbddc->local_primal_size+j]=array[j];} //WARNING -> column major ordering 1117 ierr = VecRestoreArray(vec2_V,&array);CHKERRQ(ierr); 1118 if( pcbddc->n_constraints) { 1119 ierr = VecGetArray(vec1_C,&array);CHKERRQ(ierr); 1120 for(j=0;j<pcbddc->n_constraints;j++) {coarse_submat_vals[i*pcbddc->local_primal_size+j+pcbddc->n_vertices]=array[j];} //WARNING -> column major ordering 1121 ierr = VecRestoreArray(vec1_C,&array);CHKERRQ(ierr); 1122 } 1123 1124 if( dbg_flag ) { 1125 /* assemble subdomain vector on nodes */ 1126 ierr = VecSet(pcis->vec1_N,zero);CHKERRQ(ierr); 1127 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1128 ierr = VecGetArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1129 for(j=0;j<n_R;j++) { array[idx_R_local[j]] = array2[j]; } 1130 array[ pcbddc->vertices[i] ] = one; 1131 ierr = VecRestoreArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1132 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1133 /* assemble subdomain vector of lagrange multipliers (i.e. primal nodes) */ 1134 ierr = VecSet(pcbddc->vec1_P,zero);CHKERRQ(ierr); 1135 ierr = VecGetArray(pcbddc->vec1_P,&array2);CHKERRQ(ierr); 1136 ierr = VecGetArray(vec2_V,&array);CHKERRQ(ierr); 1137 for(j=0;j<pcbddc->n_vertices;j++) { array2[j]=array[j]; } 1138 ierr = VecRestoreArray(vec2_V,&array);CHKERRQ(ierr); 1139 if(pcbddc->n_constraints) { 1140 ierr = VecGetArray(vec1_C,&array);CHKERRQ(ierr); 1141 for(j=0;j<pcbddc->n_constraints;j++) { array2[j+pcbddc->n_vertices]=array[j]; } 1142 ierr = VecRestoreArray(vec1_C,&array);CHKERRQ(ierr); 1143 } 1144 ierr = VecRestoreArray(pcbddc->vec1_P,&array2);CHKERRQ(ierr); 1145 ierr = VecScale(pcbddc->vec1_P,m_one);CHKERRQ(ierr); 1146 /* check saddle point solution */ 1147 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1148 ierr = MatMultTransposeAdd(CMAT,pcbddc->vec1_P,pcis->vec2_N,pcis->vec2_N);CHKERRQ(ierr); 1149 ierr = VecNorm(pcis->vec2_N,NORM_INFINITY,&coarsefunctions_errors[index]);CHKERRQ(ierr); 1150 ierr = MatMult(CMAT,pcis->vec1_N,pcbddc->vec1_P);CHKERRQ(ierr); 1151 ierr = VecGetArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1152 array[index]=array[index]+m_one; /* shift by the identity matrix */ 1153 ierr = VecRestoreArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1154 ierr = VecNorm(pcbddc->vec1_P,NORM_INFINITY,&constraints_errors[index]);CHKERRQ(ierr); 1155 } 1156 } 1157 1158 for(i=0;i<pcbddc->n_constraints;i++){ 1159 ierr = VecSet(vec2_C,zero);CHKERRQ(ierr); 1160 ierr = VecSetValue(vec2_C,i,m_one,INSERT_VALUES);CHKERRQ(ierr); 1161 ierr = VecAssemblyBegin(vec2_C);CHKERRQ(ierr); 1162 ierr = VecAssemblyEnd(vec2_C);CHKERRQ(ierr); 1163 /* solution of saddle point problem */ 1164 ierr = MatMult(M1,vec2_C,vec1_C);CHKERRQ(ierr); 1165 ierr = MatMult(pcbddc->local_auxmat2,vec1_C,pcbddc->vec1_R);CHKERRQ(ierr); 1166 ierr = VecScale(vec1_C,m_one);CHKERRQ(ierr); 1167 if(pcbddc->n_vertices) { ierr = MatMult(A_VR,pcbddc->vec1_R,vec2_V);CHKERRQ(ierr); } 1168 /* Set values in coarse basis function and subdomain part of coarse_mat */ 1169 /* coarse basis functions */ 1170 index=i+pcbddc->n_vertices; 1171 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1172 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1173 ierr = VecScatterEnd (pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1174 ierr = VecGetArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1175 ierr = MatSetValues(pcbddc->coarse_phi_B,n_B,auxindices,1,&index,array,INSERT_VALUES);CHKERRQ(ierr); 1176 ierr = VecRestoreArray(pcis->vec1_B,&array);CHKERRQ(ierr); 1177 if( pcbddc->prec_type || dbg_flag ) { 1178 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1179 ierr = VecScatterEnd(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1180 ierr = VecGetArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1181 ierr = MatSetValues(pcbddc->coarse_phi_D,n_D,auxindices,1,&index,array,INSERT_VALUES);CHKERRQ(ierr); 1182 ierr = VecRestoreArray(pcis->vec1_D,&array);CHKERRQ(ierr); 1183 } 1184 /* subdomain contribution to coarse matrix */ 1185 if(pcbddc->n_vertices) { 1186 ierr = VecGetArray(vec2_V,&array);CHKERRQ(ierr); 1187 for(j=0;j<pcbddc->n_vertices;j++) {coarse_submat_vals[index*pcbddc->local_primal_size+j]=array[j];} //WARNING -> column major ordering 1188 ierr = VecRestoreArray(vec2_V,&array);CHKERRQ(ierr); 1189 } 1190 ierr = VecGetArray(vec1_C,&array);CHKERRQ(ierr); 1191 for(j=0;j<pcbddc->n_constraints;j++) {coarse_submat_vals[index*pcbddc->local_primal_size+j+pcbddc->n_vertices]=array[j];} //WARNING -> column major ordering 1192 ierr = VecRestoreArray(vec1_C,&array);CHKERRQ(ierr); 1193 1194 if( dbg_flag ) { 1195 /* assemble subdomain vector on nodes */ 1196 ierr = VecSet(pcis->vec1_N,zero);CHKERRQ(ierr); 1197 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1198 ierr = VecGetArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1199 for(j=0;j<n_R;j++){ array[ idx_R_local[j] ] = array2[j]; } 1200 ierr = VecRestoreArray(pcbddc->vec1_R,&array2);CHKERRQ(ierr); 1201 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1202 /* assemble subdomain vector of lagrange multipliers */ 1203 ierr = VecSet(pcbddc->vec1_P,zero);CHKERRQ(ierr); 1204 ierr = VecGetArray(pcbddc->vec1_P,&array2);CHKERRQ(ierr); 1205 if( pcbddc->n_vertices) { 1206 ierr = VecGetArray(vec2_V,&array);CHKERRQ(ierr); 1207 for(j=0;j<pcbddc->n_vertices;j++) {array2[j]=-array[j];} 1208 ierr = VecRestoreArray(vec2_V,&array);CHKERRQ(ierr); 1209 } 1210 ierr = VecGetArray(vec1_C,&array);CHKERRQ(ierr); 1211 for(j=0;j<pcbddc->n_constraints;j++) {array2[j+pcbddc->n_vertices]=-array[j];} 1212 ierr = VecRestoreArray(vec1_C,&array);CHKERRQ(ierr); 1213 ierr = VecRestoreArray(pcbddc->vec1_P,&array2);CHKERRQ(ierr); 1214 /* check saddle point solution */ 1215 ierr = MatMult(matis->A,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 1216 ierr = MatMultTransposeAdd(CMAT,pcbddc->vec1_P,pcis->vec2_N,pcis->vec2_N);CHKERRQ(ierr); 1217 ierr = VecNorm(pcis->vec2_N,NORM_INFINITY,&coarsefunctions_errors[index]);CHKERRQ(ierr); 1218 ierr = MatMult(CMAT,pcis->vec1_N,pcbddc->vec1_P);CHKERRQ(ierr); 1219 ierr = VecGetArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1220 array[index]=array[index]+m_one; /* shift by the identity matrix */ 1221 ierr = VecRestoreArray(pcbddc->vec1_P,&array);CHKERRQ(ierr); 1222 ierr = VecNorm(pcbddc->vec1_P,NORM_INFINITY,&constraints_errors[index]);CHKERRQ(ierr); 1223 } 1224 } 1225 ierr = MatAssemblyBegin(pcbddc->coarse_phi_B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1226 ierr = MatAssemblyEnd (pcbddc->coarse_phi_B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1227 if( pcbddc->prec_type || dbg_flag ) { 1228 ierr = MatAssemblyBegin(pcbddc->coarse_phi_D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1229 ierr = MatAssemblyEnd (pcbddc->coarse_phi_D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1230 } 1231 /* Checking coarse_sub_mat and coarse basis functios */ 1232 /* It shuld be \Phi^{(j)^T} A^{(j)} \Phi^{(j)}=coarse_sub_mat */ 1233 if(dbg_flag) { 1234 1235 Mat coarse_sub_mat; 1236 Mat TM1,TM2,TM3,TM4; 1237 Mat coarse_phi_D,coarse_phi_B,A_II,A_BB,A_IB,A_BI; 1238 const MatType checkmattype; 1239 PetscScalar value; 1240 1241 ierr = MatGetType(matis->A,&checkmattype);CHKERRQ(ierr); 1242 ierr = MatGetType(pcis->A_II,&checkmattype);CHKERRQ(ierr); 1243 checkmattype = MATSEQAIJ; 1244 ierr = MatConvert(pcis->A_II,checkmattype,MAT_INITIAL_MATRIX,&A_II);CHKERRQ(ierr); 1245 ierr = MatConvert(pcis->A_IB,checkmattype,MAT_INITIAL_MATRIX,&A_IB);CHKERRQ(ierr); 1246 ierr = MatConvert(pcis->A_BI,checkmattype,MAT_INITIAL_MATRIX,&A_BI);CHKERRQ(ierr); 1247 ierr = MatConvert(pcis->A_BB,checkmattype,MAT_INITIAL_MATRIX,&A_BB);CHKERRQ(ierr); 1248 ierr = MatConvert(pcbddc->coarse_phi_D,checkmattype,MAT_INITIAL_MATRIX,&coarse_phi_D);CHKERRQ(ierr); 1249 ierr = MatConvert(pcbddc->coarse_phi_B,checkmattype,MAT_INITIAL_MATRIX,&coarse_phi_B);CHKERRQ(ierr); 1250 ierr = MatCreateSeqDense(PETSC_COMM_SELF,pcbddc->local_primal_size,pcbddc->local_primal_size,coarse_submat_vals,&coarse_sub_mat);CHKERRQ(ierr); 1251 ierr = MatConvert(coarse_sub_mat,checkmattype,MAT_REUSE_MATRIX,&coarse_sub_mat);CHKERRQ(ierr); 1252 1253 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 1254 ierr = PetscViewerASCIIPrintf(viewer,"Check coarse sub mat and local basis functions\n");CHKERRQ(ierr); 1255 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1256 ierr = MatPtAP(A_II,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&TM1);CHKERRQ(ierr); 1257 ierr = MatPtAP(A_BB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&TM2);CHKERRQ(ierr); 1258 ierr = MatMatMult(A_IB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 1259 ierr = MatTransposeMatMult(coarse_phi_D,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM3);CHKERRQ(ierr); 1260 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 1261 ierr = MatMatMult(A_BI,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 1262 ierr = MatTransposeMatMult(coarse_phi_B,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM4);CHKERRQ(ierr); 1263 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 1264 ierr = MatAXPY(TM1,one,TM2,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 1265 ierr = MatAXPY(TM1,one,TM3,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 1266 ierr = MatAXPY(TM1,one,TM4,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 1267 ierr = MatAXPY(TM1,m_one,coarse_sub_mat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 1268 ierr = MatNorm(TM1,NORM_INFINITY,&value);CHKERRQ(ierr); 1269 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"----------------------------------\n");CHKERRQ(ierr); 1270 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d \n",PetscGlobalRank);CHKERRQ(ierr); 1271 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"matrix error = % 1.14e\n",value);CHKERRQ(ierr); 1272 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"coarse functions errors\n");CHKERRQ(ierr); 1273 for(i=0;i<pcbddc->local_primal_size;i++) { ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local %02d-th function error = % 1.14e\n",i,coarsefunctions_errors[i]);CHKERRQ(ierr); } 1274 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"constraints errors\n");CHKERRQ(ierr); 1275 for(i=0;i<pcbddc->local_primal_size;i++) { ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local %02d-th function error = % 1.14e\n",i,constraints_errors[i]);CHKERRQ(ierr); } 1276 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1277 ierr = MatDestroy(&A_II);CHKERRQ(ierr); 1278 ierr = MatDestroy(&A_BB);CHKERRQ(ierr); 1279 ierr = MatDestroy(&A_IB);CHKERRQ(ierr); 1280 ierr = MatDestroy(&A_BI);CHKERRQ(ierr); 1281 ierr = MatDestroy(&TM1);CHKERRQ(ierr); 1282 ierr = MatDestroy(&TM2);CHKERRQ(ierr); 1283 ierr = MatDestroy(&TM3);CHKERRQ(ierr); 1284 ierr = MatDestroy(&TM4);CHKERRQ(ierr); 1285 ierr = MatDestroy(&coarse_phi_D);CHKERRQ(ierr); 1286 ierr = MatDestroy(&coarse_sub_mat);CHKERRQ(ierr); 1287 ierr = MatDestroy(&coarse_phi_B);CHKERRQ(ierr); 1288 ierr = PetscFree(coarsefunctions_errors);CHKERRQ(ierr); 1289 ierr = PetscFree(constraints_errors);CHKERRQ(ierr); 1290 } 1291 1292 /* create coarse matrix and data structures for message passing associated actual choice of coarse problem type */ 1293 ierr = PCBDDCSetupCoarseEnvironment(pc,coarse_submat_vals);CHKERRQ(ierr); 1294 /* free memory */ 1295 ierr = PetscFree(coarse_submat_vals);CHKERRQ(ierr); 1296 ierr = PetscFree(auxindices);CHKERRQ(ierr); 1297 ierr = PetscFree(nnz);CHKERRQ(ierr); 1298 if(pcbddc->n_vertices) { 1299 ierr = VecDestroy(&vec1_V);CHKERRQ(ierr); 1300 ierr = VecDestroy(&vec2_V);CHKERRQ(ierr); 1301 ierr = MatDestroy(&M2);CHKERRQ(ierr); 1302 ierr = MatDestroy(&A_RV);CHKERRQ(ierr); 1303 ierr = MatDestroy(&A_VR);CHKERRQ(ierr); 1304 ierr = MatDestroy(&A_VV);CHKERRQ(ierr); 1305 } 1306 if(pcbddc->n_constraints) { 1307 ierr = VecDestroy(&vec1_C);CHKERRQ(ierr); 1308 ierr = VecDestroy(&vec2_C);CHKERRQ(ierr); 1309 ierr = MatDestroy(&M1);CHKERRQ(ierr); 1310 ierr = MatDestroy(&C_CR);CHKERRQ(ierr); 1311 } 1312 ierr = MatDestroy(&CMAT);CHKERRQ(ierr); 1313 } 1314 /* free memory */ 1315 if(pcbddc->n_vertices) { 1316 ierr = PetscFree(idx_V_B);CHKERRQ(ierr); 1317 ierr = ISDestroy(&is_V_local);CHKERRQ(ierr); 1318 } 1319 ierr = PetscFree(idx_R_local);CHKERRQ(ierr); 1320 ierr = ISDestroy(&is_R_local);CHKERRQ(ierr); 1321 1322 PetscFunctionReturn(0); 1323 } 1324 1325 /* -------------------------------------------------------------------------- */ 1326 1327 #undef __FUNCT__ 1328 #define __FUNCT__ "PCBDDCSetupCoarseEnvironment" 1329 static PetscErrorCode PCBDDCSetupCoarseEnvironment(PC pc,PetscScalar* coarse_submat_vals) 1330 { 1331 1332 1333 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 1334 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1335 PC_IS *pcis = (PC_IS*)pc->data; 1336 MPI_Comm prec_comm = ((PetscObject)pc)->comm; 1337 MPI_Comm coarse_comm; 1338 1339 /* common to all choiches */ 1340 PetscScalar *temp_coarse_mat_vals; 1341 PetscScalar *ins_coarse_mat_vals; 1342 PetscInt *ins_local_primal_indices; 1343 PetscMPIInt *localsizes2,*localdispl2; 1344 PetscMPIInt size_prec_comm; 1345 PetscMPIInt rank_prec_comm; 1346 PetscMPIInt active_rank=MPI_PROC_NULL; 1347 PetscMPIInt master_proc=0; 1348 PetscInt ins_local_primal_size; 1349 /* specific to MULTILEVEL_BDDC */ 1350 PetscMPIInt *ranks_recv; 1351 PetscMPIInt count_recv=0; 1352 PetscMPIInt rank_coarse_proc_send_to; 1353 PetscMPIInt coarse_color = MPI_UNDEFINED; 1354 ISLocalToGlobalMapping coarse_ISLG; 1355 /* some other variables */ 1356 PetscErrorCode ierr; 1357 const MatType coarse_mat_type; 1358 const PCType coarse_pc_type; 1359 const KSPType coarse_ksp_type; 1360 PC pc_temp; 1361 PetscInt i,j,k,bs; 1362 /* verbose output viewer */ 1363 PetscViewer viewer=pcbddc->dbg_viewer; 1364 PetscBool dbg_flag=pcbddc->dbg_flag; 1365 1366 PetscFunctionBegin; 1367 1368 ins_local_primal_indices = 0; 1369 ins_coarse_mat_vals = 0; 1370 localsizes2 = 0; 1371 localdispl2 = 0; 1372 temp_coarse_mat_vals = 0; 1373 coarse_ISLG = 0; 1374 1375 ierr = MPI_Comm_size(prec_comm,&size_prec_comm);CHKERRQ(ierr); 1376 ierr = MPI_Comm_rank(prec_comm,&rank_prec_comm);CHKERRQ(ierr); 1377 ierr = MatGetBlockSize(matis->A,&bs);CHKERRQ(ierr); 1378 1379 /* Assign global numbering to coarse dofs */ 1380 { 1381 PetscScalar coarsesum,one=1.,zero=0.; 1382 PetscScalar *array; 1383 PetscMPIInt *auxlocal_primal; 1384 PetscMPIInt *auxglobal_primal; 1385 PetscMPIInt *all_auxglobal_primal; 1386 PetscMPIInt *all_auxglobal_primal_dummy; 1387 PetscMPIInt mpi_local_primal_size = (PetscMPIInt)pcbddc->local_primal_size; 1388 1389 /* Construct needed data structures for message passing */ 1390 ierr = PetscMalloc(mpi_local_primal_size*sizeof(PetscMPIInt),&pcbddc->local_primal_indices);CHKERRQ(ierr); 1391 ierr = PetscMalloc(size_prec_comm*sizeof(PetscMPIInt),&pcbddc->local_primal_sizes);CHKERRQ(ierr); 1392 ierr = PetscMalloc(size_prec_comm*sizeof(PetscMPIInt),&pcbddc->local_primal_displacements);CHKERRQ(ierr); 1393 /* Gather local_primal_size information for all processes */ 1394 ierr = MPI_Allgather(&mpi_local_primal_size,1,MPIU_INT,&pcbddc->local_primal_sizes[0],1,MPIU_INT,prec_comm);CHKERRQ(ierr); 1395 pcbddc->replicated_primal_size = 0; 1396 for (i=0; i<size_prec_comm; i++) { 1397 pcbddc->local_primal_displacements[i] = pcbddc->replicated_primal_size ; 1398 pcbddc->replicated_primal_size += pcbddc->local_primal_sizes[i]; 1399 } 1400 if(rank_prec_comm == 0) { 1401 /* allocate some auxiliary space */ 1402 ierr = PetscMalloc(pcbddc->replicated_primal_size*sizeof(*all_auxglobal_primal),&all_auxglobal_primal);CHKERRQ(ierr); 1403 ierr = PetscMalloc(pcbddc->replicated_primal_size*sizeof(*all_auxglobal_primal_dummy),&all_auxglobal_primal_dummy);CHKERRQ(ierr); 1404 } 1405 ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(PetscMPIInt),&auxlocal_primal);CHKERRQ(ierr); 1406 ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(PetscMPIInt),&auxglobal_primal);CHKERRQ(ierr); 1407 1408 /* First let's count coarse dofs: note that we allow to have a constraint on a subdomain and not its counterpart on the neighbour subdomain (if user wants) 1409 This code fragment assumes that the number of local constraints per connected component 1410 is not greater than the number of nodes defined for the connected component 1411 (otherwise we will surely have linear dependence between constraints and thus a singular coarse problem) */ 1412 /* auxlocal_primal : primal indices in local nodes numbering (internal and interface) */ 1413 ierr = VecSet(pcis->vec1_N,zero);CHKERRQ(ierr); 1414 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1415 for(i=0;i<pcbddc->n_vertices;i++) { 1416 array[ pcbddc->vertices[i] ] = one; 1417 auxlocal_primal[i] = pcbddc->vertices[i]; 1418 } 1419 for(i=0;i<pcbddc->n_constraints;i++) { 1420 for (j=0; j<pcbddc->sizes_of_constraint[i]; j++) { 1421 k = pcbddc->indices_to_constraint[i][j]; 1422 if( array[k] == zero ) { 1423 array[k] = one; 1424 auxlocal_primal[i+pcbddc->n_vertices] = k; 1425 break; 1426 } 1427 } 1428 } 1429 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1430 ierr = VecSet(pcis->vec1_global,zero);CHKERRQ(ierr); 1431 ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1432 ierr = VecScatterEnd (matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1433 ierr = VecScatterBegin(matis->ctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1434 ierr = VecScatterEnd (matis->ctx,pcis->vec1_global,pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1435 ierr = VecGetArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1436 for(i=0;i<pcis->n;i++) { if(array[i] > 0.0) array[i] = one/array[i]; } 1437 ierr = VecRestoreArray(pcis->vec1_N,&array);CHKERRQ(ierr); 1438 ierr = VecSet(pcis->vec1_global,zero);CHKERRQ(ierr); 1439 ierr = VecScatterBegin(matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1440 ierr = VecScatterEnd (matis->ctx,pcis->vec1_N,pcis->vec1_global,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1441 1442 ierr = VecSum(pcis->vec1_global,&coarsesum);CHKERRQ(ierr); 1443 pcbddc->coarse_size = (PetscInt) coarsesum; 1444 if(dbg_flag) { 1445 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 1446 ierr = PetscViewerASCIIPrintf(viewer,"Size of coarse problem = %d (%f)\n",pcbddc->coarse_size,coarsesum);CHKERRQ(ierr); 1447 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1448 } 1449 /* Now assign them a global numbering */ 1450 /* auxglobal_primal contains indices in global nodes numbering (internal and interface) */ 1451 ierr = ISLocalToGlobalMappingApply(matis->mapping,pcbddc->local_primal_size,auxlocal_primal,auxglobal_primal);CHKERRQ(ierr); 1452 /* all_auxglobal_primal contains all primal nodes indices in global nodes numbering (internal and interface) */ 1453 ierr = MPI_Gatherv(&auxglobal_primal[0],pcbddc->local_primal_size,MPIU_INT,&all_auxglobal_primal[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_INT,0,prec_comm);CHKERRQ(ierr); 1454 1455 /* After this block all_auxglobal_primal should contains one copy of each primal node's indices in global nodes numbering */ 1456 /* It implements a function similar to PetscSortRemoveDupsInt */ 1457 if(rank_prec_comm==0) { 1458 /* dummy argument since PetscSortMPIInt doesn't exist! */ 1459 ierr = PetscSortMPIIntWithArray(pcbddc->replicated_primal_size,all_auxglobal_primal,all_auxglobal_primal_dummy);CHKERRQ(ierr); 1460 k=1; 1461 j=all_auxglobal_primal[0]; /* first dof in global numbering */ 1462 for(i=1;i< pcbddc->replicated_primal_size ;i++) { 1463 if(j != all_auxglobal_primal[i] ) { 1464 all_auxglobal_primal[k]=all_auxglobal_primal[i]; 1465 k++; 1466 j=all_auxglobal_primal[i]; 1467 } 1468 } 1469 } else { 1470 ierr = PetscMalloc(pcbddc->coarse_size*sizeof(PetscMPIInt),&all_auxglobal_primal);CHKERRQ(ierr); 1471 } 1472 /* We only need to broadcast the indices from 0 to pcbddc->coarse_size. Remaning elements of array all_aux_global_primal are garbage. */ 1473 ierr = MPI_Bcast(all_auxglobal_primal,pcbddc->coarse_size,MPIU_INT,0,prec_comm);CHKERRQ(ierr); 1474 1475 /* Now get global coarse numbering of local primal nodes */ 1476 for(i=0;i<pcbddc->local_primal_size;i++) { 1477 k=0; 1478 while( all_auxglobal_primal[k] != auxglobal_primal[i] ) { k++;} 1479 pcbddc->local_primal_indices[i]=k; 1480 } 1481 if(dbg_flag) { 1482 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 1483 ierr = PetscViewerASCIIPrintf(viewer,"Distribution of local primal indices\n");CHKERRQ(ierr); 1484 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1485 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d\n",PetscGlobalRank);CHKERRQ(ierr); 1486 for(i=0;i<pcbddc->local_primal_size;i++) { 1487 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"local_primal_indices[%d]=%d \n",i,pcbddc->local_primal_indices[i]); 1488 } 1489 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 1490 } 1491 /* free allocated memory */ 1492 ierr = PetscFree(auxlocal_primal);CHKERRQ(ierr); 1493 ierr = PetscFree(auxglobal_primal);CHKERRQ(ierr); 1494 ierr = PetscFree(all_auxglobal_primal);CHKERRQ(ierr); 1495 if(rank_prec_comm == 0) { 1496 ierr = PetscFree(all_auxglobal_primal_dummy);CHKERRQ(ierr); 1497 } 1498 } 1499 1500 /* adapt coarse problem type */ 1501 if(pcbddc->coarse_problem_type == MULTILEVEL_BDDC && pcbddc->active_procs < MIN_PROCS_FOR_BDDC ) 1502 pcbddc->coarse_problem_type = PARALLEL_BDDC; 1503 1504 switch(pcbddc->coarse_problem_type){ 1505 1506 case(MULTILEVEL_BDDC): //we define a coarse mesh where subdomains are elements 1507 { 1508 /* we need additional variables */ 1509 MetisInt n_subdomains,n_parts,objval,ncon,faces_nvtxs; 1510 MetisInt *metis_coarse_subdivision; 1511 MetisInt options[METIS_NOPTIONS]; 1512 PetscMPIInt size_coarse_comm,rank_coarse_comm; 1513 PetscMPIInt procs_jumps_coarse_comm; 1514 PetscMPIInt *coarse_subdivision; 1515 PetscMPIInt *total_count_recv; 1516 PetscMPIInt *total_ranks_recv; 1517 PetscMPIInt *displacements_recv; 1518 PetscMPIInt *my_faces_connectivity; 1519 PetscMPIInt *petsc_faces_adjncy; 1520 MetisInt *faces_adjncy; 1521 MetisInt *faces_xadj; 1522 PetscMPIInt *number_of_faces; 1523 PetscMPIInt *faces_displacements; 1524 PetscInt *array_int; 1525 PetscMPIInt my_faces=0; 1526 PetscMPIInt total_faces=0; 1527 PetscInt ranks_stretching_ratio; 1528 1529 /* define some quantities */ 1530 pcbddc->coarse_communications_type = SCATTERS_BDDC; 1531 coarse_mat_type = MATIS; 1532 coarse_pc_type = PCBDDC; 1533 coarse_ksp_type = KSPRICHARDSON; 1534 1535 /* details of coarse decomposition */ 1536 n_subdomains = pcbddc->active_procs; 1537 n_parts = n_subdomains/pcbddc->coarsening_ratio; 1538 ranks_stretching_ratio = size_prec_comm/pcbddc->active_procs; 1539 procs_jumps_coarse_comm = pcbddc->coarsening_ratio*ranks_stretching_ratio; 1540 1541 printf("Coarse algorithm details: \n"); 1542 printf("n_subdomains %d, n_parts %d\nstretch %d,jumps %d,coarse_ratio %d\nlevel should be %d\n",n_subdomains,n_parts,ranks_stretching_ratio,procs_jumps_coarse_comm,pcbddc->coarsening_ratio,ranks_stretching_ratio/pcbddc->coarsening_ratio); 1543 1544 /* build CSR graph of subdomains' connectivity through faces */ 1545 ierr = PetscMalloc (pcis->n*sizeof(PetscInt),&array_int);CHKERRQ(ierr); 1546 ierr = PetscMemzero(array_int,pcis->n*sizeof(PetscInt));CHKERRQ(ierr); 1547 for(i=1;i<pcis->n_neigh;i++){/* i=1 so I don't count myself -> faces nodes counts to 1 */ 1548 for(j=0;j<pcis->n_shared[i];j++){ 1549 array_int[ pcis->shared[i][j] ]+=1; 1550 } 1551 } 1552 for(i=1;i<pcis->n_neigh;i++){ 1553 for(j=0;j<pcis->n_shared[i];j++){ 1554 if(array_int[ pcis->shared[i][j] ] == 1 ){ 1555 my_faces++; 1556 break; 1557 } 1558 } 1559 } 1560 //printf("I found %d faces.\n",my_faces); 1561 1562 ierr = MPI_Reduce(&my_faces,&total_faces,1,MPIU_INT,MPI_SUM,master_proc,prec_comm);CHKERRQ(ierr); 1563 ierr = PetscMalloc (my_faces*sizeof(PetscInt),&my_faces_connectivity);CHKERRQ(ierr); 1564 my_faces=0; 1565 for(i=1;i<pcis->n_neigh;i++){ 1566 for(j=0;j<pcis->n_shared[i];j++){ 1567 if(array_int[ pcis->shared[i][j] ] == 1 ){ 1568 my_faces_connectivity[my_faces]=pcis->neigh[i]; 1569 my_faces++; 1570 break; 1571 } 1572 } 1573 } 1574 if(rank_prec_comm == master_proc) { 1575 //printf("I found %d total faces.\n",total_faces); 1576 ierr = PetscMalloc (total_faces*sizeof(PetscMPIInt),&petsc_faces_adjncy);CHKERRQ(ierr); 1577 ierr = PetscMalloc (size_prec_comm*sizeof(PetscMPIInt),&number_of_faces);CHKERRQ(ierr); 1578 ierr = PetscMalloc (total_faces*sizeof(MetisInt),&faces_adjncy);CHKERRQ(ierr); 1579 ierr = PetscMalloc ((n_subdomains+1)*sizeof(MetisInt),&faces_xadj);CHKERRQ(ierr); 1580 ierr = PetscMalloc ((size_prec_comm+1)*sizeof(PetscMPIInt),&faces_displacements);CHKERRQ(ierr); 1581 } 1582 ierr = MPI_Gather(&my_faces,1,MPIU_INT,&number_of_faces[0],1,MPIU_INT,master_proc,prec_comm);CHKERRQ(ierr); 1583 if(rank_prec_comm == master_proc) { 1584 faces_xadj[0]=0; 1585 faces_displacements[0]=0; 1586 j=0; 1587 for(i=1;i<size_prec_comm+1;i++) { 1588 faces_displacements[i]=faces_displacements[i-1]+number_of_faces[i-1]; 1589 if(number_of_faces[i-1]) { 1590 j++; 1591 faces_xadj[j]=faces_xadj[j-1]+number_of_faces[i-1]; 1592 } 1593 } 1594 printf("The J I count is %d and should be %d\n",j,n_subdomains); 1595 printf("Total faces seem %d and should be %d\n",faces_xadj[j],total_faces); 1596 } 1597 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); 1598 ierr = PetscFree(my_faces_connectivity);CHKERRQ(ierr); 1599 ierr = PetscFree(array_int);CHKERRQ(ierr); 1600 if(rank_prec_comm == master_proc) { 1601 for(i=0;i<total_faces;i++) faces_adjncy[i]=(MetisInt)(petsc_faces_adjncy[i]/ranks_stretching_ratio); /* cast to MetisInt */ 1602 printf("This is the face connectivity (actual ranks)\n"); 1603 for(i=0;i<n_subdomains;i++){ 1604 printf("proc %d is connected with \n",i); 1605 for(j=faces_xadj[i];j<faces_xadj[i+1];j++) 1606 printf("%d ",faces_adjncy[j]); 1607 printf("\n"); 1608 } 1609 ierr = PetscFree(faces_displacements);CHKERRQ(ierr); 1610 ierr = PetscFree(number_of_faces);CHKERRQ(ierr); 1611 ierr = PetscFree(petsc_faces_adjncy);CHKERRQ(ierr); 1612 } 1613 1614 if( rank_prec_comm == master_proc ) { 1615 1616 PetscInt heuristic_for_metis=3; 1617 1618 ncon=1; 1619 faces_nvtxs=n_subdomains; 1620 /* partition graoh induced by face connectivity */ 1621 ierr = PetscMalloc (n_subdomains*sizeof(MetisInt),&metis_coarse_subdivision);CHKERRQ(ierr); 1622 ierr = METIS_SetDefaultOptions(options); 1623 /* we need a contiguous partition of the coarse mesh */ 1624 options[METIS_OPTION_CONTIG]=1; 1625 options[METIS_OPTION_DBGLVL]=1; 1626 options[METIS_OPTION_NITER]=30; 1627 //options[METIS_OPTION_NCUTS]=1; 1628 printf("METIS PART GRAPH\n"); 1629 if(n_subdomains>n_parts*heuristic_for_metis) { 1630 printf("Using Kway\n"); 1631 options[METIS_OPTION_IPTYPE]=METIS_IPTYPE_EDGE; 1632 options[METIS_OPTION_OBJTYPE]=METIS_OBJTYPE_CUT; 1633 ierr = METIS_PartGraphKway(&faces_nvtxs,&ncon,faces_xadj,faces_adjncy,NULL,NULL,NULL,&n_parts,NULL,NULL,options,&objval,metis_coarse_subdivision); 1634 } else { 1635 printf("Using Recursive\n"); 1636 ierr = METIS_PartGraphRecursive(&faces_nvtxs,&ncon,faces_xadj,faces_adjncy,NULL,NULL,NULL,&n_parts,NULL,NULL,options,&objval,metis_coarse_subdivision); 1637 } 1638 if(ierr != METIS_OK) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in METIS_PartGraphKway (metis error code %D) called from PCBDDCSetupCoarseEnvironment\n",ierr); 1639 printf("Partition done!\n"); 1640 ierr = PetscFree(faces_xadj);CHKERRQ(ierr); 1641 ierr = PetscFree(faces_adjncy);CHKERRQ(ierr); 1642 coarse_subdivision = (PetscMPIInt*)calloc(size_prec_comm,sizeof(PetscMPIInt)); /* calloc for contiguous memory since we need to scatter these values later */ 1643 /* copy/cast values avoiding possible type conflicts between PETSc, MPI and METIS */ 1644 for(i=0;i<size_prec_comm;i++) coarse_subdivision[i]=MPI_PROC_NULL; 1645 for(i=0;i<n_subdomains;i++) coarse_subdivision[ranks_stretching_ratio*i]=(PetscInt)(metis_coarse_subdivision[i]); 1646 ierr = PetscFree(metis_coarse_subdivision);CHKERRQ(ierr); 1647 } 1648 1649 /* Create new communicator for coarse problem splitting the old one */ 1650 if( !(rank_prec_comm%procs_jumps_coarse_comm) && rank_prec_comm < procs_jumps_coarse_comm*n_parts ){ 1651 coarse_color=0; //for communicator splitting 1652 active_rank=rank_prec_comm; //for insertion of matrix values 1653 } 1654 // procs with coarse_color = MPI_UNDEFINED will have coarse_comm = MPI_COMM_NULL (from mpi standards) 1655 // key = rank_prec_comm -> keep same ordering of ranks from the old to the new communicator 1656 ierr = MPI_Comm_split(prec_comm,coarse_color,rank_prec_comm,&coarse_comm);CHKERRQ(ierr); 1657 1658 if( coarse_color == 0 ) { 1659 ierr = MPI_Comm_size(coarse_comm,&size_coarse_comm);CHKERRQ(ierr); 1660 ierr = MPI_Comm_rank(coarse_comm,&rank_coarse_comm);CHKERRQ(ierr); 1661 printf("Details of coarse comm\n"); 1662 printf("size = %d, myrank = %d\n",size_coarse_comm,rank_coarse_comm); 1663 printf("jumps = %d, coarse_color = %d, n_parts = %d\n",procs_jumps_coarse_comm,coarse_color,n_parts); 1664 } else { 1665 rank_coarse_comm = MPI_PROC_NULL; 1666 } 1667 1668 /* master proc take care of arranging and distributing coarse informations */ 1669 if(rank_coarse_comm == master_proc) { 1670 ierr = PetscMalloc (size_coarse_comm*sizeof(PetscMPIInt),&displacements_recv);CHKERRQ(ierr); 1671 //ierr = PetscMalloc (size_coarse_comm*sizeof(PetscMPIInt),&total_count_recv);CHKERRQ(ierr); 1672 //ierr = PetscMalloc (n_subdomains*sizeof(PetscMPIInt),&total_ranks_recv);CHKERRQ(ierr); 1673 total_count_recv = (PetscMPIInt*)calloc(size_prec_comm,sizeof(PetscMPIInt)); 1674 total_ranks_recv = (PetscMPIInt*)calloc(n_subdomains,sizeof(PetscMPIInt)); 1675 /* some initializations */ 1676 displacements_recv[0]=0; 1677 //PetscMemzero(total_count_recv,size_coarse_comm*sizeof(PetscMPIInt)); not needed -> calloc initializes to zero 1678 /* count from how many processes the j-th process of the coarse decomposition will receive data */ 1679 for(j=0;j<size_coarse_comm;j++) 1680 for(i=0;i<size_prec_comm;i++) 1681 if(coarse_subdivision[i]==j) 1682 total_count_recv[j]++; 1683 /* displacements needed for scatterv of total_ranks_recv */ 1684 for(i=1;i<size_coarse_comm;i++) displacements_recv[i]=displacements_recv[i-1]+total_count_recv[i-1]; 1685 /* Now fill properly total_ranks_recv -> each coarse process will receive the ranks (in prec_comm communicator) of its friend (sending) processes */ 1686 ierr = PetscMemzero(total_count_recv,size_coarse_comm*sizeof(PetscMPIInt));CHKERRQ(ierr); 1687 for(j=0;j<size_coarse_comm;j++) { 1688 for(i=0;i<size_prec_comm;i++) { 1689 if(coarse_subdivision[i]==j) { 1690 total_ranks_recv[displacements_recv[j]+total_count_recv[j]]=i; 1691 total_count_recv[j]+=1; 1692 } 1693 } 1694 } 1695 for(j=0;j<size_coarse_comm;j++) { 1696 printf("process %d in new rank will receive from %d processes (original ranks follows)\n",j,total_count_recv[j]); 1697 for(i=0;i<total_count_recv[j];i++) { 1698 printf("%d ",total_ranks_recv[displacements_recv[j]+i]); 1699 } 1700 printf("\n"); 1701 } 1702 1703 /* identify new decomposition in terms of ranks in the old communicator */ 1704 for(i=0;i<n_subdomains;i++) coarse_subdivision[ranks_stretching_ratio*i]=coarse_subdivision[ranks_stretching_ratio*i]*procs_jumps_coarse_comm; 1705 printf("coarse_subdivision in old end new ranks\n"); 1706 for(i=0;i<size_prec_comm;i++) 1707 if(coarse_subdivision[i]!=MPI_PROC_NULL) { 1708 printf("%d=(%d %d), ",i,coarse_subdivision[i],coarse_subdivision[i]/procs_jumps_coarse_comm); 1709 } else { 1710 printf("%d=(%d %d), ",i,coarse_subdivision[i],coarse_subdivision[i]); 1711 } 1712 printf("\n"); 1713 } 1714 1715 /* Scatter new decomposition for send details */ 1716 ierr = MPI_Scatter(&coarse_subdivision[0],1,MPIU_INT,&rank_coarse_proc_send_to,1,MPIU_INT,master_proc,prec_comm);CHKERRQ(ierr); 1717 /* Scatter receiving details to members of coarse decomposition */ 1718 if( coarse_color == 0) { 1719 ierr = MPI_Scatter(&total_count_recv[0],1,MPIU_INT,&count_recv,1,MPIU_INT,master_proc,coarse_comm);CHKERRQ(ierr); 1720 ierr = PetscMalloc (count_recv*sizeof(PetscMPIInt),&ranks_recv);CHKERRQ(ierr); 1721 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); 1722 } 1723 1724 //printf("I will send my matrix data to proc %d\n",rank_coarse_proc_send_to); 1725 //if(coarse_color == 0) { 1726 // printf("I will receive some matrix data from %d processes (ranks follows)\n",count_recv); 1727 // for(i=0;i<count_recv;i++) 1728 // printf("%d ",ranks_recv[i]); 1729 // printf("\n"); 1730 //} 1731 1732 if(rank_prec_comm == master_proc) { 1733 //ierr = PetscFree(coarse_subdivision);CHKERRQ(ierr); 1734 //ierr = PetscFree(total_count_recv);CHKERRQ(ierr); 1735 //ierr = PetscFree(total_ranks_recv);CHKERRQ(ierr); 1736 free(coarse_subdivision); 1737 free(total_count_recv); 1738 free(total_ranks_recv); 1739 ierr = PetscFree(displacements_recv);CHKERRQ(ierr); 1740 } 1741 break; 1742 } 1743 1744 case(REPLICATED_BDDC): 1745 1746 pcbddc->coarse_communications_type = GATHERS_BDDC; 1747 coarse_mat_type = MATSEQAIJ; 1748 coarse_pc_type = PCLU; 1749 coarse_ksp_type = KSPPREONLY; 1750 coarse_comm = PETSC_COMM_SELF; 1751 active_rank = rank_prec_comm; 1752 break; 1753 1754 case(PARALLEL_BDDC): 1755 1756 pcbddc->coarse_communications_type = SCATTERS_BDDC; 1757 coarse_mat_type = MATMPIAIJ; 1758 coarse_pc_type = PCREDUNDANT; 1759 coarse_ksp_type = KSPPREONLY; 1760 coarse_comm = prec_comm; 1761 active_rank = rank_prec_comm; 1762 break; 1763 1764 case(SEQUENTIAL_BDDC): 1765 pcbddc->coarse_communications_type = GATHERS_BDDC; 1766 coarse_mat_type = MATSEQAIJ; 1767 coarse_pc_type = PCLU; 1768 coarse_ksp_type = KSPPREONLY; 1769 coarse_comm = PETSC_COMM_SELF; 1770 active_rank = master_proc; 1771 break; 1772 } 1773 1774 switch(pcbddc->coarse_communications_type){ 1775 1776 case(SCATTERS_BDDC): 1777 { 1778 if(pcbddc->coarse_problem_type==MULTILEVEL_BDDC) { 1779 1780 PetscMPIInt send_size; 1781 PetscInt *aux_ins_indices; 1782 PetscInt ii,jj; 1783 MPI_Request *requests; 1784 1785 /* allocate auxiliary space */ 1786 ierr = PetscMalloc (pcbddc->replicated_primal_size*sizeof(PetscMPIInt),&pcbddc->replicated_local_primal_indices);CHKERRQ(ierr); 1787 ierr = MPI_Allgatherv(&pcbddc->local_primal_indices[0],pcbddc->local_primal_size,MPIU_INT,&pcbddc->replicated_local_primal_indices[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_INT,prec_comm);CHKERRQ(ierr); 1788 ierr = PetscMalloc ( pcbddc->coarse_size*sizeof(PetscInt),&aux_ins_indices);CHKERRQ(ierr); 1789 ierr = PetscMemzero(aux_ins_indices,pcbddc->coarse_size*sizeof(PetscInt));CHKERRQ(ierr); 1790 /* allocate stuffs for message massing */ 1791 ierr = PetscMalloc ( (count_recv+1)*sizeof(MPI_Request),&requests);CHKERRQ(ierr); 1792 for(i=0;i<count_recv+1;i++) requests[i]=MPI_REQUEST_NULL; 1793 ierr = PetscMalloc ( count_recv*sizeof(PetscMPIInt),&localsizes2);CHKERRQ(ierr); 1794 ierr = PetscMalloc ( count_recv*sizeof(PetscMPIInt),&localdispl2);CHKERRQ(ierr); 1795 /* fill up quantities */ 1796 j=0; 1797 for(i=0;i<count_recv;i++){ 1798 ii = ranks_recv[i]; 1799 localsizes2[i]=pcbddc->local_primal_sizes[ii]*pcbddc->local_primal_sizes[ii]; 1800 localdispl2[i]=j; 1801 j+=localsizes2[i]; 1802 jj = pcbddc->local_primal_displacements[ii]; 1803 for(k=0;k<pcbddc->local_primal_sizes[ii];k++) aux_ins_indices[pcbddc->replicated_local_primal_indices[jj+k]]+=1; // it counts the coarse subdomains sharing the coarse node 1804 } 1805 //printf("aux_ins_indices 1\n"); 1806 //for(i=0;i<pcbddc->coarse_size;i++) 1807 // printf("%d ",aux_ins_indices[i]); 1808 //printf("\n"); 1809 /* temp_coarse_mat_vals used to store temporarly received matrix values */ 1810 ierr = PetscMalloc ( j*sizeof(PetscScalar),&temp_coarse_mat_vals);CHKERRQ(ierr); 1811 /* evaluate how many values I will insert in coarse mat */ 1812 ins_local_primal_size=0; 1813 for(i=0;i<pcbddc->coarse_size;i++) 1814 if(aux_ins_indices[i]) 1815 ins_local_primal_size++; 1816 /* evaluate indices I will insert in coarse mat */ 1817 ierr = PetscMalloc ( ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);CHKERRQ(ierr); 1818 j=0; 1819 for(i=0;i<pcbddc->coarse_size;i++) 1820 if(aux_ins_indices[i]) 1821 ins_local_primal_indices[j++]=i; 1822 /* use aux_ins_indices to realize a global to local mapping */ 1823 j=0; 1824 for(i=0;i<pcbddc->coarse_size;i++){ 1825 if(aux_ins_indices[i]==0){ 1826 aux_ins_indices[i]=-1; 1827 } else { 1828 aux_ins_indices[i]=j; 1829 j++; 1830 } 1831 } 1832 1833 //printf("New details localsizes2 localdispl2\n"); 1834 //for(i=0;i<count_recv;i++) 1835 // printf("(%d %d) ",localsizes2[i],localdispl2[i]); 1836 //printf("\n"); 1837 //printf("aux_ins_indices 2\n"); 1838 //for(i=0;i<pcbddc->coarse_size;i++) 1839 // printf("%d ",aux_ins_indices[i]); 1840 //printf("\n"); 1841 //printf("ins_local_primal_indices\n"); 1842 //for(i=0;i<ins_local_primal_size;i++) 1843 // printf("%d ",ins_local_primal_indices[i]); 1844 //printf("\n"); 1845 //printf("coarse_submat_vals\n"); 1846 //for(i=0;i<pcbddc->local_primal_size;i++) 1847 // for(j=0;j<pcbddc->local_primal_size;j++) 1848 // printf("(%lf %d %d)\n",coarse_submat_vals[j*pcbddc->local_primal_size+i],pcbddc->local_primal_indices[i],pcbddc->local_primal_indices[j]); 1849 //printf("\n"); 1850 1851 /* processes partecipating in coarse problem receive matrix data from their friends */ 1852 for(i=0;i<count_recv;i++) ierr = MPI_Irecv(&temp_coarse_mat_vals[localdispl2[i]],localsizes2[i],MPIU_SCALAR,ranks_recv[i],666,prec_comm,&requests[i]);CHKERRQ(ierr); 1853 if(rank_coarse_proc_send_to != MPI_PROC_NULL ) { 1854 send_size=pcbddc->local_primal_size*pcbddc->local_primal_size; 1855 ierr = MPI_Isend(&coarse_submat_vals[0],send_size,MPIU_SCALAR,rank_coarse_proc_send_to,666,prec_comm,&requests[count_recv]);CHKERRQ(ierr); 1856 } 1857 ierr = MPI_Waitall(count_recv+1,requests,MPI_STATUSES_IGNORE);CHKERRQ(ierr); 1858 1859 //if(coarse_color == 0) { 1860 // printf("temp_coarse_mat_vals\n"); 1861 // for(k=0;k<count_recv;k++){ 1862 // printf("---- %d ----\n",ranks_recv[k]); 1863 // for(i=0;i<pcbddc->local_primal_sizes[ranks_recv[k]];i++) 1864 // for(j=0;j<pcbddc->local_primal_sizes[ranks_recv[k]];j++) 1865 // printf("(%lf %d %d)\n",temp_coarse_mat_vals[localdispl2[k]+j*pcbddc->local_primal_sizes[ranks_recv[k]]+i],pcbddc->replicated_local_primal_indices[pcbddc->local_primal_displacements[ranks_recv[k]]+i],pcbddc->replicated_local_primal_indices[pcbddc->local_primal_displacements[ranks_recv[k]]+j]); 1866 // printf("\n"); 1867 // } 1868 //} 1869 /* calculate data to insert in coarse mat */ 1870 ierr = PetscMalloc ( ins_local_primal_size*ins_local_primal_size*sizeof(PetscScalar),&ins_coarse_mat_vals);CHKERRQ(ierr); 1871 PetscMemzero(ins_coarse_mat_vals,ins_local_primal_size*ins_local_primal_size*sizeof(PetscScalar)); 1872 1873 PetscMPIInt rr,kk,lps,lpd; 1874 PetscInt row_ind,col_ind; 1875 for(k=0;k<count_recv;k++){ 1876 rr = ranks_recv[k]; 1877 kk = localdispl2[k]; 1878 lps = pcbddc->local_primal_sizes[rr]; 1879 lpd = pcbddc->local_primal_displacements[rr]; 1880 //printf("Inserting the following indices (received from %d)\n",rr); 1881 for(j=0;j<lps;j++){ 1882 col_ind=aux_ins_indices[pcbddc->replicated_local_primal_indices[lpd+j]]; 1883 for(i=0;i<lps;i++){ 1884 row_ind=aux_ins_indices[pcbddc->replicated_local_primal_indices[lpd+i]]; 1885 //printf("%d %d\n",row_ind,col_ind); 1886 ins_coarse_mat_vals[col_ind*ins_local_primal_size+row_ind]+=temp_coarse_mat_vals[kk+j*lps+i]; 1887 } 1888 } 1889 } 1890 ierr = PetscFree(requests);CHKERRQ(ierr); 1891 ierr = PetscFree(aux_ins_indices);CHKERRQ(ierr); 1892 ierr = PetscFree(temp_coarse_mat_vals);CHKERRQ(ierr); 1893 if(coarse_color == 0) { ierr = PetscFree(ranks_recv);CHKERRQ(ierr); } 1894 1895 /* create local to global mapping needed by coarse MATIS */ 1896 { 1897 IS coarse_IS; 1898 if(coarse_comm != MPI_COMM_NULL ) ierr = MPI_Comm_free(&coarse_comm);CHKERRQ(ierr); 1899 coarse_comm = prec_comm; 1900 active_rank=rank_prec_comm; 1901 ierr = ISCreateGeneral(coarse_comm,ins_local_primal_size,ins_local_primal_indices,PETSC_COPY_VALUES,&coarse_IS);CHKERRQ(ierr); 1902 //ierr = ISSetBlockSize(coarse_IS,bs);CHKERRQ(ierr); 1903 ierr = ISLocalToGlobalMappingCreateIS(coarse_IS,&coarse_ISLG);CHKERRQ(ierr); 1904 ierr = ISDestroy(&coarse_IS);CHKERRQ(ierr); 1905 } 1906 } 1907 if(pcbddc->coarse_problem_type==PARALLEL_BDDC) { 1908 /* arrays for values insertion */ 1909 ins_local_primal_size = pcbddc->local_primal_size; 1910 ierr = PetscMalloc ( ins_local_primal_size*sizeof(PetscMPIInt),&ins_local_primal_indices);CHKERRQ(ierr); 1911 ierr = PetscMalloc ( ins_local_primal_size*ins_local_primal_size*sizeof(PetscScalar),&ins_coarse_mat_vals);CHKERRQ(ierr); 1912 for(j=0;j<ins_local_primal_size;j++){ 1913 ins_local_primal_indices[j]=pcbddc->local_primal_indices[j]; 1914 for(i=0;i<ins_local_primal_size;i++) ins_coarse_mat_vals[j*ins_local_primal_size+i]=coarse_submat_vals[j*ins_local_primal_size+i]; 1915 } 1916 } 1917 break; 1918 1919 } 1920 1921 case(GATHERS_BDDC): 1922 { 1923 1924 PetscMPIInt mysize,mysize2; 1925 1926 if(rank_prec_comm==active_rank) { 1927 ierr = PetscMalloc ( pcbddc->replicated_primal_size*sizeof(PetscMPIInt),&pcbddc->replicated_local_primal_indices);CHKERRQ(ierr); 1928 pcbddc->replicated_local_primal_values = (PetscScalar*)calloc(pcbddc->replicated_primal_size,sizeof(PetscScalar)); 1929 ierr = PetscMalloc ( size_prec_comm*sizeof(PetscMPIInt),&localsizes2);CHKERRQ(ierr); 1930 ierr = PetscMalloc ( size_prec_comm*sizeof(PetscMPIInt),&localdispl2);CHKERRQ(ierr); 1931 /* arrays for values insertion */ 1932 ins_local_primal_size = pcbddc->coarse_size; 1933 ierr = PetscMalloc ( ins_local_primal_size*sizeof(PetscMPIInt),&ins_local_primal_indices);CHKERRQ(ierr); 1934 ierr = PetscMalloc ( ins_local_primal_size*ins_local_primal_size*sizeof(PetscScalar),&ins_coarse_mat_vals);CHKERRQ(ierr); 1935 for(i=0;i<size_prec_comm;i++) localsizes2[i]=pcbddc->local_primal_sizes[i]*pcbddc->local_primal_sizes[i]; 1936 localdispl2[0]=0; 1937 for(i=1;i<size_prec_comm;i++) localdispl2[i]=localsizes2[i-1]+localdispl2[i-1]; 1938 j=0; 1939 for(i=0;i<size_prec_comm;i++) j+=localsizes2[i]; 1940 ierr = PetscMalloc ( j*sizeof(PetscScalar),&temp_coarse_mat_vals);CHKERRQ(ierr); 1941 } 1942 1943 mysize=pcbddc->local_primal_size; 1944 mysize2=pcbddc->local_primal_size*pcbddc->local_primal_size; 1945 if(pcbddc->coarse_problem_type == SEQUENTIAL_BDDC){ 1946 ierr = MPI_Gatherv(&pcbddc->local_primal_indices[0],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); 1947 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); 1948 } else { 1949 ierr = MPI_Allgatherv(&pcbddc->local_primal_indices[0],mysize,MPIU_INT,&pcbddc->replicated_local_primal_indices[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_INT,prec_comm);CHKERRQ(ierr); 1950 ierr = MPI_Allgatherv(&coarse_submat_vals[0],mysize2,MPIU_SCALAR,&temp_coarse_mat_vals[0],localsizes2,localdispl2,MPIU_SCALAR,prec_comm);CHKERRQ(ierr); 1951 } 1952 1953 /* free data structures no longer needed and allocate some space which will be needed in BDDC application */ 1954 if(rank_prec_comm==active_rank) { 1955 PetscInt offset,offset2,row_ind,col_ind; 1956 for(j=0;j<ins_local_primal_size;j++){ 1957 ins_local_primal_indices[j]=j; 1958 for(i=0;i<ins_local_primal_size;i++) ins_coarse_mat_vals[j*ins_local_primal_size+i]=0.0; 1959 } 1960 for(k=0;k<size_prec_comm;k++){ 1961 offset=pcbddc->local_primal_displacements[k]; 1962 offset2=localdispl2[k]; 1963 for(j=0;j<pcbddc->local_primal_sizes[k];j++){ 1964 col_ind=pcbddc->replicated_local_primal_indices[offset+j]; 1965 for(i=0;i<pcbddc->local_primal_sizes[k];i++){ 1966 row_ind=pcbddc->replicated_local_primal_indices[offset+i]; 1967 ins_coarse_mat_vals[col_ind*pcbddc->coarse_size+row_ind]+=temp_coarse_mat_vals[offset2+j*pcbddc->local_primal_sizes[k]+i]; 1968 } 1969 } 1970 } 1971 } 1972 break; 1973 }//switch on coarse problem and communications associated with finished 1974 } 1975 1976 /* Now create and fill up coarse matrix */ 1977 if( rank_prec_comm == active_rank ) { 1978 if(pcbddc->coarse_problem_type != MULTILEVEL_BDDC) { 1979 ierr = MatCreate(coarse_comm,&pcbddc->coarse_mat);CHKERRQ(ierr); 1980 ierr = MatSetSizes(pcbddc->coarse_mat,PETSC_DECIDE,PETSC_DECIDE,pcbddc->coarse_size,pcbddc->coarse_size);CHKERRQ(ierr); 1981 ierr = MatSetType(pcbddc->coarse_mat,coarse_mat_type);CHKERRQ(ierr); 1982 ierr = MatSetOption(pcbddc->coarse_mat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 1983 ierr = MatSetOption(pcbddc->coarse_mat,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); //local values stored in column major 1984 } else { 1985 Mat matis_coarse_local_mat; 1986 ierr = MatCreateIS(coarse_comm,PETSC_DECIDE,PETSC_DECIDE,pcbddc->coarse_size,pcbddc->coarse_size,coarse_ISLG,&pcbddc->coarse_mat);CHKERRQ(ierr); 1987 ierr = MatISGetLocalMat(pcbddc->coarse_mat,&matis_coarse_local_mat);CHKERRQ(ierr); 1988 ierr = MatSetOption(matis_coarse_local_mat,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); //local values stored in column major 1989 //ierr = MatSetOption(matis_coarse_local_mat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 1990 } 1991 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); 1992 ierr = MatAssemblyBegin(pcbddc->coarse_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1993 ierr = MatAssemblyEnd(pcbddc->coarse_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1994 if(pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 1995 Mat matis_coarse_local_mat; 1996 ierr = MatISGetLocalMat(pcbddc->coarse_mat,&matis_coarse_local_mat);CHKERRQ(ierr); 1997 ierr = MatSetBlockSize(matis_coarse_local_mat,bs);CHKERRQ(ierr); 1998 } 1999 2000 ierr = MatGetVecs(pcbddc->coarse_mat,&pcbddc->coarse_vec,&pcbddc->coarse_rhs);CHKERRQ(ierr); 2001 /* Preconditioner for coarse problem */ 2002 ierr = KSPCreate(coarse_comm,&pcbddc->coarse_ksp);CHKERRQ(ierr); 2003 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->coarse_ksp,(PetscObject)pc,1);CHKERRQ(ierr); 2004 ierr = KSPSetOperators(pcbddc->coarse_ksp,pcbddc->coarse_mat,pcbddc->coarse_mat,SAME_PRECONDITIONER);CHKERRQ(ierr); 2005 ierr = KSPSetTolerances(pcbddc->coarse_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,1);CHKERRQ(ierr); 2006 ierr = KSPSetType(pcbddc->coarse_ksp,coarse_ksp_type);CHKERRQ(ierr); 2007 ierr = KSPGetPC(pcbddc->coarse_ksp,&pc_temp);CHKERRQ(ierr); 2008 ierr = PCSetType(pc_temp,coarse_pc_type);CHKERRQ(ierr); 2009 //ierr = PCSetOptionsPrefix(pcbddc->coarse_pc,"pc_bddc_coarse_");CHKERRQ(ierr); 2010 /* Allow user's customization */ 2011 ierr = KSPSetFromOptions(pcbddc->coarse_ksp);CHKERRQ(ierr); 2012 /* Set Up PC for coarse problem BDDC */ 2013 //if(pcbddc->dbg_flag && pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 2014 if(pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 2015 if(dbg_flag) { 2016 ierr = PetscViewerASCIIPrintf(viewer,"----------------Setting up a new level---------------\n");CHKERRQ(ierr); 2017 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2018 } 2019 ierr = PCBDDCSetCoarseProblemType(pc_temp,MULTILEVEL_BDDC);CHKERRQ(ierr); 2020 } 2021 ierr = KSPSetUp(pcbddc->coarse_ksp);CHKERRQ(ierr); 2022 if(pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 2023 if(dbg_flag) { 2024 ierr = PetscViewerASCIIPrintf(viewer,"----------------New level set------------------------\n");CHKERRQ(ierr); 2025 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2026 } 2027 } 2028 } 2029 if(pcbddc->coarse_communications_type == SCATTERS_BDDC) { 2030 IS local_IS,global_IS; 2031 ierr = ISCreateStride(PETSC_COMM_SELF,pcbddc->local_primal_size,0,1,&local_IS);CHKERRQ(ierr); 2032 ierr = ISCreateGeneral(PETSC_COMM_SELF,pcbddc->local_primal_size,pcbddc->local_primal_indices,PETSC_COPY_VALUES,&global_IS);CHKERRQ(ierr); 2033 ierr = VecScatterCreate(pcbddc->vec1_P,local_IS,pcbddc->coarse_vec,global_IS,&pcbddc->coarse_loc_to_glob);CHKERRQ(ierr); 2034 ierr = ISDestroy(&local_IS);CHKERRQ(ierr); 2035 ierr = ISDestroy(&global_IS);CHKERRQ(ierr); 2036 } 2037 2038 2039 /* Check condition number of coarse problem */ 2040 if( pcbddc->coarse_problem_type == MULTILEVEL_BDDC && dbg_flag && rank_prec_comm == active_rank ) { 2041 PetscScalar m_one=-1.0; 2042 PetscReal infty_error,lambda_min,lambda_max,kappa_2; 2043 2044 /* change coarse ksp object to an iterative method suitable for extreme eigenvalues' estimation */ 2045 ierr = KSPSetComputeSingularValues(pcbddc->coarse_ksp,PETSC_TRUE);CHKERRQ(ierr); 2046 ierr = KSPSetType(pcbddc->coarse_ksp,KSPGMRES);CHKERRQ(ierr); 2047 ierr = KSPSetTolerances(pcbddc->coarse_ksp,1.e-8,1.e-8,PETSC_DEFAULT,pcbddc->coarse_size);CHKERRQ(ierr); 2048 ierr = KSPSetUp(pcbddc->coarse_ksp);CHKERRQ(ierr); 2049 ierr = VecSetRandom(pcbddc->coarse_rhs,PETSC_NULL);CHKERRQ(ierr); 2050 ierr = MatMult(pcbddc->coarse_mat,pcbddc->coarse_rhs,pcbddc->coarse_vec);CHKERRQ(ierr); 2051 ierr = MatMult(pcbddc->coarse_mat,pcbddc->coarse_vec,pcbddc->coarse_rhs);CHKERRQ(ierr); 2052 ierr = KSPSolve(pcbddc->coarse_ksp,pcbddc->coarse_rhs,pcbddc->coarse_rhs);CHKERRQ(ierr); 2053 ierr = KSPComputeExtremeSingularValues(pcbddc->coarse_ksp,&lambda_max,&lambda_min);CHKERRQ(ierr); 2054 kappa_2=lambda_max/lambda_min; 2055 ierr = KSPGetIterationNumber(pcbddc->coarse_ksp,&k);CHKERRQ(ierr); 2056 ierr = VecAXPY(pcbddc->coarse_rhs,m_one,pcbddc->coarse_vec);CHKERRQ(ierr); 2057 ierr = VecNorm(pcbddc->coarse_rhs,NORM_INFINITY,&infty_error);CHKERRQ(ierr); 2058 ierr = PetscViewerASCIIPrintf(viewer,"Coarse problem condition number estimated with %d iterations of gmres is: % 1.14e\n",k,kappa_2);CHKERRQ(ierr); 2059 ierr = PetscViewerASCIIPrintf(viewer,"Coarse problem eigenvalues: % 1.14e %1.14e\n",lambda_min,lambda_max);CHKERRQ(ierr); 2060 ierr = PetscViewerASCIIPrintf(viewer,"Coarse problem infty_error: %1.14e\n",infty_error);CHKERRQ(ierr); 2061 /* restore coarse ksp to default values */ 2062 ierr = KSPSetComputeSingularValues(pcbddc->coarse_ksp,PETSC_FALSE);CHKERRQ(ierr); 2063 ierr = KSPSetType(pcbddc->coarse_ksp,coarse_ksp_type);CHKERRQ(ierr); 2064 ierr = KSPSetTolerances(pcbddc->coarse_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,1);CHKERRQ(ierr); 2065 ierr = KSPSetFromOptions(pcbddc->coarse_ksp);CHKERRQ(ierr); 2066 ierr = KSPSetUp(pcbddc->coarse_ksp);CHKERRQ(ierr); 2067 } 2068 2069 /* free data structures no longer needed */ 2070 if(coarse_ISLG) { ierr = ISLocalToGlobalMappingDestroy(&coarse_ISLG);CHKERRQ(ierr); } 2071 if(ins_local_primal_indices) { ierr = PetscFree(ins_local_primal_indices);CHKERRQ(ierr); } 2072 if(ins_coarse_mat_vals) { ierr = PetscFree(ins_coarse_mat_vals);CHKERRQ(ierr);} 2073 if(localsizes2) { ierr = PetscFree(localsizes2);CHKERRQ(ierr);} 2074 if(localdispl2) { ierr = PetscFree(localdispl2);CHKERRQ(ierr);} 2075 if(temp_coarse_mat_vals) { ierr = PetscFree(temp_coarse_mat_vals);CHKERRQ(ierr);} 2076 2077 PetscFunctionReturn(0); 2078 } 2079 2080 #undef __FUNCT__ 2081 #define __FUNCT__ "PCBDDCManageLocalBoundaries" 2082 static PetscErrorCode PCBDDCManageLocalBoundaries(PC pc) 2083 { 2084 2085 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 2086 PC_IS *pcis = (PC_IS*)pc->data; 2087 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 2088 PetscInt *distinct_values; 2089 PetscInt **neighbours_set; 2090 PetscInt bs,ierr,i,j,s,k,iindex; 2091 PetscInt total_counts; 2092 PetscBool flg_row; 2093 PCBDDCGraph mat_graph; 2094 PetscInt neumann_bsize; 2095 const PetscInt* neumann_nodes; 2096 Mat mat_adj; 2097 2098 PetscFunctionBegin; 2099 2100 ierr = MatGetBlockSize(matis->A,&bs);CHKERRQ(ierr); 2101 // allocate and initialize needed graph structure 2102 ierr = PetscMalloc(sizeof(*mat_graph),&mat_graph);CHKERRQ(ierr); 2103 ierr = MatConvert(matis->A,MATMPIADJ,MAT_INITIAL_MATRIX,&mat_adj);CHKERRQ(ierr); 2104 ierr = MatGetRowIJ(mat_adj,0,PETSC_FALSE,PETSC_FALSE,&mat_graph->nvtxs,&mat_graph->xadj,&mat_graph->adjncy,&flg_row);CHKERRQ(ierr); 2105 if(!flg_row) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatGetRowIJ called from PCBDDCManageLocalBoundaries.\n"); 2106 ierr = PetscMalloc(mat_graph->nvtxs*sizeof(PetscInt),&mat_graph->where);CHKERRQ(ierr); 2107 ierr = PetscMemzero(mat_graph->where,mat_graph->nvtxs*sizeof(PetscInt));CHKERRQ(ierr); 2108 ierr = PetscMalloc(mat_graph->nvtxs*sizeof(PetscInt),&mat_graph->count);CHKERRQ(ierr); 2109 ierr = PetscMemzero(mat_graph->count,mat_graph->nvtxs*sizeof(PetscInt));CHKERRQ(ierr); 2110 ierr = PetscMalloc(mat_graph->nvtxs*sizeof(PetscInt),&mat_graph->which_dof);CHKERRQ(ierr); 2111 ierr = PetscMemzero(mat_graph->which_dof,mat_graph->nvtxs*sizeof(PetscInt));CHKERRQ(ierr); 2112 ierr = PetscMalloc(mat_graph->nvtxs*sizeof(PetscInt),&mat_graph->queue);CHKERRQ(ierr); 2113 ierr = PetscMemzero(mat_graph->queue,mat_graph->nvtxs*sizeof(PetscInt));CHKERRQ(ierr); 2114 ierr = PetscMalloc((mat_graph->nvtxs+1)*sizeof(PetscInt),&mat_graph->cptr);CHKERRQ(ierr); 2115 ierr = PetscMemzero(mat_graph->cptr,(mat_graph->nvtxs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2116 ierr = PetscMalloc(mat_graph->nvtxs*sizeof(PetscBool),&mat_graph->touched);CHKERRQ(ierr); 2117 for(i=0;i<mat_graph->nvtxs;i++){mat_graph->touched[i]=PETSC_FALSE;} 2118 for(i=0;i<mat_graph->nvtxs/bs;i++) { 2119 for(s=0;s<bs;s++) { 2120 mat_graph->which_dof[i*bs+s]=s; 2121 } 2122 } 2123 //printf("nvtxs = %d\n",mat_graph->nvtxs); 2124 //printf("these are my IS data with n_neigh = %d\n",pcis->n_neigh); 2125 //for(i=0;i<pcis->n_neigh;i++){ 2126 // printf("number of shared nodes with rank %d is %d \n",pcis->neigh[i],pcis->n_shared[i]); 2127 // } 2128 2129 total_counts=0; 2130 for(i=0;i<pcis->n_neigh;i++){ 2131 s=pcis->n_shared[i]; 2132 total_counts+=s; 2133 //printf("computing neigh %d (rank = %d, n_shared = %d)\n",i,pcis->neigh[i],s); 2134 for(j=0;j<s;j++){ 2135 mat_graph->count[pcis->shared[i][j]] += 1; 2136 } 2137 } 2138 /* Take into account Neumann data incrementing number of sharing subdomains for all but faces nodes lying on the interface */ 2139 if(pcbddc->NeumannBoundaries) { 2140 ierr = ISGetLocalSize(pcbddc->NeumannBoundaries,&neumann_bsize);CHKERRQ(ierr); 2141 ierr = ISGetIndices(pcbddc->NeumannBoundaries,&neumann_nodes);CHKERRQ(ierr); 2142 for(i=0;i<neumann_bsize;i++){ 2143 iindex = neumann_nodes[i]; 2144 if(mat_graph->count[iindex] > 2){ 2145 mat_graph->count[iindex]+=1; 2146 total_counts++; 2147 } 2148 } 2149 } 2150 2151 ierr = PetscMalloc(mat_graph->nvtxs*sizeof(PetscInt*),&neighbours_set);CHKERRQ(ierr); 2152 if(mat_graph->nvtxs) ierr = PetscMalloc(total_counts*sizeof(PetscInt),&neighbours_set[0]);CHKERRQ(ierr); 2153 for(i=1;i<mat_graph->nvtxs;i++) neighbours_set[i]=neighbours_set[i-1]+mat_graph->count[i-1]; 2154 2155 for(i=0;i<mat_graph->nvtxs;i++) mat_graph->count[i]=0; 2156 for(i=0;i<pcis->n_neigh;i++){ 2157 s=pcis->n_shared[i]; 2158 for(j=0;j<s;j++) { 2159 k=pcis->shared[i][j]; 2160 neighbours_set[k][mat_graph->count[k]] = pcis->neigh[i]; 2161 mat_graph->count[k]+=1; 2162 } 2163 } 2164 /* set -1 fake neighbour */ 2165 if(pcbddc->NeumannBoundaries) { 2166 for(i=0;i<neumann_bsize;i++){ 2167 iindex = neumann_nodes[i]; 2168 if( mat_graph->count[iindex] > 2){ 2169 neighbours_set[iindex][mat_graph->count[iindex]] = -1; //An additional fake neighbour (with rank -1) 2170 mat_graph->count[iindex]+=1; 2171 } 2172 } 2173 ierr = ISRestoreIndices(pcbddc->NeumannBoundaries,&neumann_nodes);CHKERRQ(ierr); 2174 } 2175 2176 /* sort sharing subdomains */ 2177 for(i=0;i<mat_graph->nvtxs;i++) { ierr = PetscSortInt(mat_graph->count[i],neighbours_set[i]);CHKERRQ(ierr); } 2178 2179 /* Prepare for FindConnectedComponents 2180 Vertices will be eliminated later (if needed) */ 2181 PetscInt nodes_touched=0; 2182 for(i=0;i<mat_graph->nvtxs;i++){ 2183 if(!mat_graph->count[i]){ //internal nodes 2184 mat_graph->touched[i]=PETSC_TRUE; 2185 mat_graph->where[i]=0; 2186 nodes_touched++; 2187 } 2188 if(pcbddc->faces_flag){ 2189 if(mat_graph->count[i]>2){ //all but faces 2190 mat_graph->touched[i]=PETSC_TRUE; 2191 mat_graph->where[i]=0; 2192 nodes_touched++; 2193 } 2194 } 2195 if(pcbddc->edges_flag){ 2196 if(mat_graph->count[i]==2){ //faces 2197 mat_graph->touched[i]=PETSC_TRUE; 2198 mat_graph->where[i]=0; 2199 nodes_touched++; 2200 } 2201 } 2202 } 2203 2204 PetscInt rvalue=1; 2205 PetscBool same_set; 2206 mat_graph->ncmps = 0; 2207 while(nodes_touched<mat_graph->nvtxs) { 2208 // find first untouched node in local ordering 2209 i=0; 2210 while(mat_graph->touched[i]) i++; 2211 mat_graph->touched[i]=PETSC_TRUE; 2212 mat_graph->where[i]=rvalue; 2213 nodes_touched++; 2214 // now find other connected nodes shared by the same set of subdomains 2215 for(j=i+1;j<mat_graph->nvtxs;j++){ 2216 // check for same number of sharing subdomains and dof number 2217 if(mat_graph->count[i]==mat_graph->count[j] && mat_graph->which_dof[i] == mat_graph->which_dof[j] ){ 2218 // check for same set of sharing subdomains 2219 same_set=PETSC_TRUE; 2220 for(k=0;k<mat_graph->count[j];k++){ 2221 if(neighbours_set[i][k]!=neighbours_set[j][k]) { 2222 same_set=PETSC_FALSE; 2223 } 2224 } 2225 // OK, I found a friend of mine 2226 if(same_set) { 2227 mat_graph->where[j]=rvalue; 2228 mat_graph->touched[j]=PETSC_TRUE; 2229 nodes_touched++; 2230 } 2231 } 2232 } 2233 rvalue++; 2234 } 2235 // printf("where and count contains %d distinct values\n",rvalue); 2236 // for(j=0;j<mat_graph->nvtxs;j++) 2237 // printf("[%d %d %d]\n",j,mat_graph->where[j],mat_graph->count[j]); 2238 2239 if(mat_graph->nvtxs) { 2240 ierr = PetscFree(neighbours_set[0]);CHKERRQ(ierr); 2241 ierr = PetscFree(neighbours_set);CHKERRQ(ierr); 2242 } 2243 2244 rvalue--; 2245 ierr = PetscMalloc ( rvalue*sizeof(PetscInt),&distinct_values);CHKERRQ(ierr); 2246 for(i=0;i<rvalue;i++) distinct_values[i]=i+1; //initializiation 2247 if(rvalue) ierr = PCBDDCFindConnectedComponents(mat_graph, rvalue, distinct_values); 2248 //printf("total number of connected components %d \n",mat_graph->ncmps); 2249 //for (i=0; i<mat_graph->ncmps; i++) { 2250 // printf("[queue num %d] ptr %d, length %d, start index %d\n",i,mat_graph->cptr[i],mat_graph->cptr[i+1]-mat_graph->cptr[i],mat_graph->queue[mat_graph->cptr[i]]); 2251 //} 2252 PetscInt nfc=0; 2253 PetscInt nec=0; 2254 PetscInt nvc=0; 2255 for (i=0; i<mat_graph->ncmps; i++) { 2256 // sort each connected component (by local ordering) 2257 ierr = PetscSortInt(mat_graph->cptr[i+1]-mat_graph->cptr[i],&mat_graph->queue[mat_graph->cptr[i]]);CHKERRQ(ierr); 2258 // count edge and faces 2259 if( !pcbddc->vertices_flag ) { 2260 if( mat_graph->cptr[i+1]-mat_graph->cptr[i] > 1 ){ 2261 if(mat_graph->count[mat_graph->queue[mat_graph->cptr[i]]]==2){ 2262 nfc++; 2263 } else { 2264 nec++; 2265 } 2266 } 2267 } 2268 // count vertices 2269 if( !pcbddc->constraints_flag ){ 2270 if( mat_graph->cptr[i+1]-mat_graph->cptr[i] == 1 ){ 2271 nvc++; 2272 } 2273 } 2274 } 2275 2276 pcbddc->n_constraints = nec+nfc; 2277 pcbddc->n_vertices = nvc; 2278 2279 if(pcbddc->n_constraints){ 2280 /* allocate space for local constraints of BDDC */ 2281 ierr = PetscMalloc (pcbddc->n_constraints*sizeof(PetscInt*),&pcbddc->indices_to_constraint);CHKERRQ(ierr); 2282 ierr = PetscMalloc (pcbddc->n_constraints*sizeof(PetscScalar*),&pcbddc->quadrature_constraint);CHKERRQ(ierr); 2283 ierr = PetscMalloc (pcbddc->n_constraints*sizeof(PetscInt),&pcbddc->sizes_of_constraint);CHKERRQ(ierr); 2284 k=0; 2285 for (i=0; i<mat_graph->ncmps; i++) { 2286 if( mat_graph->cptr[i+1]-mat_graph->cptr[i] > 1 ){ 2287 pcbddc->sizes_of_constraint[k] = mat_graph->cptr[i+1]-mat_graph->cptr[i]; 2288 k++; 2289 } 2290 } 2291 // printf("check constraints %d (should be %d)\n",k,pcbddc->n_constraints); 2292 // for(i=0;i<k;i++) 2293 // printf("%d ",pcbddc->sizes_of_constraint[i]); 2294 // printf("\n"); 2295 k=0; 2296 for (i=0; i<pcbddc->n_constraints; i++) k+=pcbddc->sizes_of_constraint[i]; 2297 ierr = PetscMalloc (k*sizeof(PetscInt),&pcbddc->indices_to_constraint[0]);CHKERRQ(ierr); 2298 ierr = PetscMalloc (k*sizeof(PetscScalar),&pcbddc->quadrature_constraint[0]);CHKERRQ(ierr); 2299 for (i=1; i<pcbddc->n_constraints; i++) { 2300 pcbddc->indices_to_constraint[i] = pcbddc->indices_to_constraint[i-1] + pcbddc->sizes_of_constraint[i-1]; 2301 pcbddc->quadrature_constraint[i] = pcbddc->quadrature_constraint[i-1] + pcbddc->sizes_of_constraint[i-1]; 2302 } 2303 k=0; 2304 PetscScalar quad_value; 2305 for (i=0; i<mat_graph->ncmps; i++) { 2306 if( mat_graph->cptr[i+1]-mat_graph->cptr[i] > 1 ){ 2307 quad_value=1.0/( (PetscScalar) (mat_graph->cptr[i+1]-mat_graph->cptr[i]) ); 2308 for(j=0;j<mat_graph->cptr[i+1]-mat_graph->cptr[i];j++) { 2309 pcbddc->indices_to_constraint[k][j] = mat_graph->queue[mat_graph->cptr[i]+j]; 2310 pcbddc->quadrature_constraint[k][j] = quad_value; 2311 } 2312 k++; 2313 } 2314 } 2315 } 2316 if(pcbddc->n_vertices){ 2317 /* allocate space for local vertices of BDDC */ 2318 ierr = PetscMalloc (pcbddc->n_vertices*sizeof(PetscInt),&pcbddc->vertices);CHKERRQ(ierr); 2319 k=0; 2320 for (i=0; i<mat_graph->ncmps; i++) { 2321 if( mat_graph->cptr[i+1]-mat_graph->cptr[i] == 1 ){ 2322 pcbddc->vertices[k] = mat_graph->queue[mat_graph->cptr[i]]; 2323 k++; 2324 } 2325 } 2326 // sort vertex set (by local ordering) 2327 ierr = PetscSortInt(pcbddc->n_vertices,pcbddc->vertices);CHKERRQ(ierr); 2328 } 2329 2330 if(pcbddc->dbg_flag) { 2331 PetscViewer viewer=pcbddc->dbg_viewer; 2332 PetscInt* queue_in_global_numbering; 2333 2334 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"--------------------------------------------------------------\n");CHKERRQ(ierr); 2335 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Details from PCBDDCManageLocalBoundaries for subdomain %04d\n",PetscGlobalRank);CHKERRQ(ierr); 2336 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"--------------------------------------------------------------\n");CHKERRQ(ierr); 2337 PetscViewerASCIISynchronizedPrintf(viewer,"Graph (adjacency structure) of local Neumann mat\n"); 2338 PetscViewerASCIISynchronizedPrintf(viewer,"--------------------------------------------------------------\n"); 2339 for(i=0;i<mat_graph->nvtxs;i++) { 2340 PetscViewerASCIISynchronizedPrintf(viewer,"Nodes connected to node number %d are %d\n",i,mat_graph->xadj[i+1]-mat_graph->xadj[i]); 2341 for(j=mat_graph->xadj[i];j<mat_graph->xadj[i+1];j++){ 2342 PetscViewerASCIISynchronizedPrintf(viewer,"%d ",mat_graph->adjncy[j]); 2343 } 2344 PetscViewerASCIISynchronizedPrintf(viewer,"\n--------------------------------------------------------------\n"); 2345 } 2346 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Matrix graph has %d connected components", mat_graph->ncmps);CHKERRQ(ierr); 2347 ierr = PetscMalloc(mat_graph->nvtxs*sizeof(PetscInt),&queue_in_global_numbering);CHKERRQ(ierr); 2348 ierr = ISLocalToGlobalMappingApply(matis->mapping,mat_graph->nvtxs,mat_graph->queue,queue_in_global_numbering);CHKERRQ(ierr); 2349 for(i=0;i<mat_graph->ncmps;i++) { 2350 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\nSize and count for connected component %02d : %04d %01d\n", i,mat_graph->cptr[i+1]-mat_graph->cptr[i],mat_graph->count[mat_graph->queue[mat_graph->cptr[i]]]);CHKERRQ(ierr); 2351 for (j=mat_graph->cptr[i]; j<mat_graph->cptr[i+1]; j++){ 2352 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"%d (%d), ",queue_in_global_numbering[j],mat_graph->queue[j]);CHKERRQ(ierr); 2353 } 2354 } 2355 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\n--------------------------------------------------------------\n");CHKERRQ(ierr); 2356 if( pcbddc->n_vertices ) { ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d detected %02d local vertices\n",PetscGlobalRank,pcbddc->n_vertices);CHKERRQ(ierr); } 2357 if( nfc ) { ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d detected %02d local faces\n",PetscGlobalRank,nfc);CHKERRQ(ierr); } 2358 if( nec ) { ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d detected %02d local edges\n",PetscGlobalRank,nec);CHKERRQ(ierr); } 2359 if( pcbddc->n_vertices ) { ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Indices for subdomain vertices follows\n",PetscGlobalRank,pcbddc->n_vertices);CHKERRQ(ierr); } 2360 ierr = ISLocalToGlobalMappingApply(matis->mapping,pcbddc->n_vertices,pcbddc->vertices,queue_in_global_numbering);CHKERRQ(ierr); 2361 for(i=0;i<pcbddc->n_vertices;i++){ ierr = PetscViewerASCIISynchronizedPrintf(viewer,"%d ",queue_in_global_numbering[i]);CHKERRQ(ierr); } 2362 if( pcbddc->n_vertices ) { ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\n");CHKERRQ(ierr); } 2363 /* if( pcbddc->n_constraints ) PetscViewerASCIISynchronizedPrintf(viewer,"Indices and quadrature constraints"); 2364 for(i=0;i<pcbddc->n_constraints;i++){ 2365 PetscViewerASCIISynchronizedPrintf(viewer,"\nConstraint number %d\n",i); 2366 for(j=0;j<pcbddc->sizes_of_constraint[i];j++) { 2367 PetscViewerASCIISynchronizedPrintf(viewer,"(%d, %f) ",pcbddc->indices_to_constraint[i][j],pcbddc->quadrature_constraint[i][j]); 2368 } 2369 } 2370 if( pcbddc->n_constraints ) PetscViewerASCIISynchronizedPrintf(viewer,"\n");*/ 2371 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2372 ierr = PetscFree(queue_in_global_numbering);CHKERRQ(ierr); 2373 } 2374 2375 // Restore CSR structure into sequantial matrix and free memory space no longer needed 2376 ierr = MatRestoreRowIJ(mat_adj,0,PETSC_FALSE,PETSC_TRUE,&mat_graph->nvtxs,&mat_graph->xadj,&mat_graph->adjncy,&flg_row);CHKERRQ(ierr); 2377 ierr = MatDestroy(&mat_adj);CHKERRQ(ierr); 2378 if(!flg_row) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called from PCBDDCManageLocalBoundaries.\n"); 2379 ierr = PetscFree(distinct_values);CHKERRQ(ierr); 2380 // Free graph structure 2381 if(mat_graph->nvtxs){ 2382 ierr = PetscFree(mat_graph->where);CHKERRQ(ierr); 2383 ierr = PetscFree(mat_graph->touched);CHKERRQ(ierr); 2384 ierr = PetscFree(mat_graph->which_dof);CHKERRQ(ierr); 2385 ierr = PetscFree(mat_graph->queue);CHKERRQ(ierr); 2386 ierr = PetscFree(mat_graph->cptr);CHKERRQ(ierr); 2387 ierr = PetscFree(mat_graph->count);CHKERRQ(ierr); 2388 } 2389 ierr = PetscFree(mat_graph);CHKERRQ(ierr); 2390 2391 PetscFunctionReturn(0); 2392 2393 } 2394 2395 /* -------------------------------------------------------------------------- */ 2396 2397 /* The following code has been adapted from function IsConnectedSubdomain contained 2398 in source file contig.c of METIS library (version 5.0.1) */ 2399 2400 #undef __FUNCT__ 2401 #define __FUNCT__ "PCBDDCFindConnectedComponents" 2402 static PetscErrorCode PCBDDCFindConnectedComponents(PCBDDCGraph graph, PetscInt n_dist, PetscInt *dist_vals) 2403 { 2404 PetscInt i, j, k, nvtxs, first, last, nleft, ncmps,pid,cum_queue,n,ncmps_pid; 2405 PetscInt *xadj, *adjncy, *where, *queue; 2406 PetscInt *cptr; 2407 PetscBool *touched; 2408 2409 PetscFunctionBegin; 2410 2411 nvtxs = graph->nvtxs; 2412 xadj = graph->xadj; 2413 adjncy = graph->adjncy; 2414 where = graph->where; 2415 touched = graph->touched; 2416 queue = graph->queue; 2417 cptr = graph->cptr; 2418 2419 for (i=0; i<nvtxs; i++) 2420 touched[i] = PETSC_FALSE; 2421 2422 cum_queue=0; 2423 ncmps=0; 2424 2425 for(n=0; n<n_dist; n++) { 2426 pid = dist_vals[n]; 2427 nleft = 0; 2428 for (i=0; i<nvtxs; i++) { 2429 if (where[i] == pid) 2430 nleft++; 2431 } 2432 for (i=0; i<nvtxs; i++) { 2433 if (where[i] == pid) 2434 break; 2435 } 2436 touched[i] = PETSC_TRUE; 2437 queue[cum_queue] = i; 2438 first = 0; last = 1; 2439 cptr[ncmps] = cum_queue; /* This actually points to queue */ 2440 ncmps_pid = 0; 2441 while (first != nleft) { 2442 if (first == last) { /* Find another starting vertex */ 2443 cptr[++ncmps] = first+cum_queue; 2444 ncmps_pid++; 2445 for (i=0; i<nvtxs; i++) { 2446 if (where[i] == pid && !touched[i]) 2447 break; 2448 } 2449 queue[cum_queue+last] = i; 2450 last++; 2451 touched[i] = PETSC_TRUE; 2452 } 2453 i = queue[cum_queue+first]; 2454 first++; 2455 for (j=xadj[i]; j<xadj[i+1]; j++) { 2456 k = adjncy[j]; 2457 if (where[k] == pid && !touched[k]) { 2458 queue[cum_queue+last] = k; 2459 last++; 2460 touched[k] = PETSC_TRUE; 2461 } 2462 } 2463 } 2464 cptr[++ncmps] = first+cum_queue; 2465 ncmps_pid++; 2466 cum_queue=cptr[ncmps]; 2467 } 2468 graph->ncmps = ncmps; 2469 2470 PetscFunctionReturn(0); 2471 } 2472 2473