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_N,&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_Gather(&mpi_local_primal_size,1,MPIU_INT,&pcbddc->local_primal_sizes[0],1,MPIU_INT,0,prec_comm);CHKERRQ(ierr); 1395 pcbddc->replicated_primal_size = 0; 1396 if(rank_prec_comm == 0) { 1397 for (i=0; i<size_prec_comm; i++) { 1398 pcbddc->local_primal_displacements[i] = pcbddc->replicated_primal_size ; 1399 pcbddc->replicated_primal_size += pcbddc->local_primal_sizes[i]; 1400 } 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\n",pcbddc->coarse_size);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 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 1528 /* this code has a bug (see below) for more then three levels -> I can solve it quickly */ 1529 1530 /* define some quantities */ 1531 pcbddc->coarse_communications_type = SCATTERS_BDDC; 1532 coarse_mat_type = MATIS; 1533 coarse_pc_type = PCBDDC; 1534 coarse_ksp_type = KSPRICHARDSON; 1535 1536 /* details of coarse decomposition */ 1537 n_subdomains = pcbddc->active_procs; 1538 n_parts = n_subdomains/pcbddc->coarsening_ratio; 1539 procs_jumps_coarse_comm = pcbddc->coarsening_ratio*(size_prec_comm/pcbddc->active_procs); 1540 1541 ierr = PetscMalloc (pcbddc->replicated_primal_size*sizeof(PetscMPIInt),&pcbddc->replicated_local_primal_indices);CHKERRQ(ierr); 1542 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); 1543 1544 /* build CSR graph of subdomains' connectivity through faces */ 1545 ierr = PetscMalloc (pcis->n*sizeof(PetscInt),&array_int);CHKERRQ(ierr); 1546 PetscMemzero(array_int,pcis->n*sizeof(PetscInt)); 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]); ///procs_jumps_coarse_comm); // cast to MetisInt 1602 printf("This is the face connectivity (%d)\n",procs_jumps_coarse_comm); 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 ncon=1; 1617 faces_nvtxs=n_subdomains; 1618 /* partition graoh induced by face connectivity */ 1619 ierr = PetscMalloc (n_subdomains*sizeof(MetisInt),&metis_coarse_subdivision);CHKERRQ(ierr); 1620 ierr = METIS_SetDefaultOptions(options); 1621 /* we need a contiguous partition of the coarse mesh */ 1622 options[METIS_OPTION_CONTIG]=1; 1623 options[METIS_OPTION_DBGLVL]=1; 1624 options[METIS_OPTION_OBJTYPE]=METIS_OBJTYPE_CUT; 1625 options[METIS_OPTION_IPTYPE]=METIS_IPTYPE_EDGE; 1626 options[METIS_OPTION_NITER]=30; 1627 //options[METIS_OPTION_NCUTS]=1; 1628 //printf("METIS PART GRAPH\n"); 1629 /* BUG: faces_xadj and faces_adjncy content must be adapted using the coarsening factor*/ 1630 ierr = METIS_PartGraphKway(&faces_nvtxs,&ncon,faces_xadj,faces_adjncy,NULL,NULL,NULL,&n_parts,NULL,NULL,options,&objval,metis_coarse_subdivision); 1631 //printf("OKOKOKOKOKOKOKOK\n"); 1632 if(ierr != METIS_OK) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in METIS_PartGraphKway (metis error code %D) called from PCBDDCSetupCoarseEnvironment\n",ierr); 1633 ierr = PetscFree(faces_xadj);CHKERRQ(ierr); 1634 ierr = PetscFree(faces_adjncy);CHKERRQ(ierr); 1635 coarse_subdivision = (PetscMPIInt*)calloc(size_prec_comm,sizeof(PetscMPIInt)); /* calloc for contiguous memory since we need to scatter these values later */ 1636 /* copy/cast values avoiding possible type conflicts between PETSc, MPI and METIS */ 1637 for(i=0;i<size_prec_comm;i++) coarse_subdivision[i]=MPI_PROC_NULL;; 1638 k=size_prec_comm/pcbddc->active_procs; 1639 for(i=0;i<n_subdomains;i++) coarse_subdivision[k*i]=(PetscInt)(metis_coarse_subdivision[i]); 1640 ierr = PetscFree(metis_coarse_subdivision);CHKERRQ(ierr); 1641 1642 } 1643 1644 /* Create new communicator for coarse problem splitting the old one */ 1645 if( !(rank_prec_comm%procs_jumps_coarse_comm) && rank_prec_comm < procs_jumps_coarse_comm*n_parts ){ 1646 coarse_color=0; //for communicator splitting 1647 active_rank=rank_prec_comm; //for insertion of matrix values 1648 } 1649 // procs with coarse_color = MPI_UNDEFINED will have coarse_comm = MPI_COMM_NULL (from mpi standards) 1650 // key = rank_prec_comm -> keep same ordering of ranks from the old to the new communicator 1651 ierr = MPI_Comm_split(prec_comm,coarse_color,rank_prec_comm,&coarse_comm);CHKERRQ(ierr); 1652 1653 if( coarse_color == 0 ) { 1654 ierr = MPI_Comm_size(coarse_comm,&size_coarse_comm);CHKERRQ(ierr); 1655 ierr = MPI_Comm_rank(coarse_comm,&rank_coarse_comm);CHKERRQ(ierr); 1656 //printf("Details of coarse comm\n"); 1657 //printf("size = %d, myrank = %d\n",size_coarse_comm,rank_coarse_comm); 1658 //printf("jumps = %d, coarse_color = %d, n_parts = %d\n",procs_jumps_coarse_comm,coarse_color,n_parts); 1659 } else { 1660 rank_coarse_comm = MPI_PROC_NULL; 1661 } 1662 1663 /* master proc take care of arranging and distributing coarse informations */ 1664 if(rank_coarse_comm == master_proc) { 1665 ierr = PetscMalloc (size_coarse_comm*sizeof(PetscMPIInt),&displacements_recv);CHKERRQ(ierr); 1666 //ierr = PetscMalloc (size_coarse_comm*sizeof(PetscMPIInt),&total_count_recv);CHKERRQ(ierr); 1667 //ierr = PetscMalloc (n_subdomains*sizeof(PetscMPIInt),&total_ranks_recv);CHKERRQ(ierr); 1668 total_count_recv = (PetscMPIInt*)calloc(size_prec_comm,sizeof(PetscMPIInt)); 1669 total_ranks_recv = (PetscMPIInt*)calloc(n_subdomains,sizeof(PetscMPIInt)); 1670 /* some initializations */ 1671 displacements_recv[0]=0; 1672 //PetscMemzero(total_count_recv,size_coarse_comm*sizeof(PetscMPIInt)); not needed -> calloc initializes to zero 1673 /* count from how many processes the j-th process of the coarse decomposition will receive data */ 1674 for(j=0;j<size_coarse_comm;j++) 1675 for(i=0;i<n_subdomains;i++) 1676 if(coarse_subdivision[i]==j) 1677 total_count_recv[j]++; 1678 /* displacements needed for scatterv of total_ranks_recv */ 1679 for(i=1;i<size_coarse_comm;i++) displacements_recv[i]=displacements_recv[i-1]+total_count_recv[i-1]; 1680 /* Now fill properly total_ranks_recv -> each coarse process will receive the ranks (in prec_comm communicator) of its friend (sending) processes */ 1681 ierr = PetscMemzero(total_count_recv,size_coarse_comm*sizeof(PetscMPIInt));CHKERRQ(ierr); 1682 for(j=0;j<size_coarse_comm;j++) { 1683 for(i=0;i<n_subdomains;i++) { 1684 if(coarse_subdivision[i]==j) { 1685 total_ranks_recv[displacements_recv[j]+total_count_recv[j]]=i; 1686 total_count_recv[j]=total_count_recv[j]+1; 1687 } 1688 } 1689 } 1690 //for(j=0;j<size_coarse_comm;j++) { 1691 // printf("process %d in new rank will receive from %d processes (ranks follows)\n",j,total_count_recv[j]); 1692 // for(i=0;i<total_count_recv[j];i++) { 1693 // printf("%d ",total_ranks_recv[displacements_recv[j]+i]); 1694 // } 1695 // printf("\n"); 1696 // } 1697 1698 /* identify new decomposition in terms of ranks in the old communicator */ 1699 k=size_prec_comm/pcbddc->active_procs; 1700 for(i=0;i<n_subdomains;i++) coarse_subdivision[k*i]=coarse_subdivision[k*i]*procs_jumps_coarse_comm; 1701 printf("coarse_subdivision in old end new ranks\n"); 1702 for(i=0;i<size_prec_comm;i++) 1703 printf("(%d %d) ",coarse_subdivision[i],coarse_subdivision[i]/procs_jumps_coarse_comm); 1704 printf("\n"); 1705 } 1706 1707 /* Scatter new decomposition for send details */ 1708 ierr = MPI_Scatter(&coarse_subdivision[0],1,MPIU_INT,&rank_coarse_proc_send_to,1,MPIU_INT,master_proc,prec_comm);CHKERRQ(ierr); 1709 /* Scatter receiving details to members of coarse decomposition */ 1710 if( coarse_color == 0) { 1711 ierr = MPI_Scatter(&total_count_recv[0],1,MPIU_INT,&count_recv,1,MPIU_INT,master_proc,coarse_comm);CHKERRQ(ierr); 1712 ierr = PetscMalloc (count_recv*sizeof(PetscMPIInt),&ranks_recv);CHKERRQ(ierr); 1713 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); 1714 } 1715 1716 //printf("I will send my matrix data to proc %d\n",rank_coarse_proc_send_to); 1717 //if(coarse_color == 0) { 1718 // printf("I will receive some matrix data from %d processes (ranks follows)\n",count_recv); 1719 // for(i=0;i<count_recv;i++) 1720 // printf("%d ",ranks_recv[i]); 1721 // printf("\n"); 1722 //} 1723 1724 if(rank_prec_comm == master_proc) { 1725 //ierr = PetscFree(coarse_subdivision);CHKERRQ(ierr); 1726 //ierr = PetscFree(total_count_recv);CHKERRQ(ierr); 1727 //ierr = PetscFree(total_ranks_recv);CHKERRQ(ierr); 1728 free(coarse_subdivision); 1729 free(total_count_recv); 1730 free(total_ranks_recv); 1731 ierr = PetscFree(displacements_recv);CHKERRQ(ierr); 1732 } 1733 break; 1734 } 1735 1736 case(REPLICATED_BDDC): 1737 1738 pcbddc->coarse_communications_type = GATHERS_BDDC; 1739 coarse_mat_type = MATSEQAIJ; 1740 coarse_pc_type = PCLU; 1741 coarse_ksp_type = KSPPREONLY; 1742 coarse_comm = PETSC_COMM_SELF; 1743 active_rank = rank_prec_comm; 1744 break; 1745 1746 case(PARALLEL_BDDC): 1747 1748 pcbddc->coarse_communications_type = SCATTERS_BDDC; 1749 coarse_mat_type = MATMPIAIJ; 1750 coarse_pc_type = PCREDUNDANT; 1751 coarse_ksp_type = KSPPREONLY; 1752 coarse_comm = prec_comm; 1753 active_rank = rank_prec_comm; 1754 break; 1755 1756 case(SEQUENTIAL_BDDC): 1757 pcbddc->coarse_communications_type = GATHERS_BDDC; 1758 coarse_mat_type = MATSEQAIJ; 1759 coarse_pc_type = PCLU; 1760 coarse_ksp_type = KSPPREONLY; 1761 coarse_comm = PETSC_COMM_SELF; 1762 active_rank = master_proc; 1763 break; 1764 } 1765 1766 switch(pcbddc->coarse_communications_type){ 1767 1768 case(SCATTERS_BDDC): 1769 { 1770 if(pcbddc->coarse_problem_type==MULTILEVEL_BDDC) { 1771 1772 PetscMPIInt send_size; 1773 PetscInt *aux_ins_indices; 1774 PetscInt ii,jj; 1775 MPI_Request *requests; 1776 1777 /* allocate auxiliary space */ 1778 ierr = PetscMalloc ( pcbddc->coarse_size*sizeof(PetscInt),&aux_ins_indices);CHKERRQ(ierr); 1779 ierr = PetscMemzero(aux_ins_indices,pcbddc->coarse_size*sizeof(PetscInt));CHKERRQ(ierr); 1780 /* allocate stuffs for message massing */ 1781 ierr = PetscMalloc ( (count_recv+1)*sizeof(MPI_Request),&requests);CHKERRQ(ierr); 1782 for(i=0;i<count_recv+1;i++) requests[i]=MPI_REQUEST_NULL; 1783 ierr = PetscMalloc ( count_recv*sizeof(PetscMPIInt),&localsizes2);CHKERRQ(ierr); 1784 ierr = PetscMalloc ( count_recv*sizeof(PetscMPIInt),&localdispl2);CHKERRQ(ierr); 1785 /* fill up quantities */ 1786 j=0; 1787 for(i=0;i<count_recv;i++){ 1788 ii = ranks_recv[i]; 1789 localsizes2[i]=pcbddc->local_primal_sizes[ii]*pcbddc->local_primal_sizes[ii]; 1790 localdispl2[i]=j; 1791 j+=localsizes2[i]; 1792 jj = pcbddc->local_primal_displacements[ii]; 1793 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 1794 } 1795 //printf("aux_ins_indices 1\n"); 1796 //for(i=0;i<pcbddc->coarse_size;i++) 1797 // printf("%d ",aux_ins_indices[i]); 1798 //printf("\n"); 1799 /* temp_coarse_mat_vals used to store temporarly received matrix values */ 1800 ierr = PetscMalloc ( j*sizeof(PetscScalar),&temp_coarse_mat_vals);CHKERRQ(ierr); 1801 /* evaluate how many values I will insert in coarse mat */ 1802 ins_local_primal_size=0; 1803 for(i=0;i<pcbddc->coarse_size;i++) 1804 if(aux_ins_indices[i]) 1805 ins_local_primal_size++; 1806 /* evaluate indices I will insert in coarse mat */ 1807 ierr = PetscMalloc ( ins_local_primal_size*sizeof(PetscInt),&ins_local_primal_indices);CHKERRQ(ierr); 1808 j=0; 1809 for(i=0;i<pcbddc->coarse_size;i++) 1810 if(aux_ins_indices[i]) 1811 ins_local_primal_indices[j++]=i; 1812 /* use aux_ins_indices to realize a global to local mapping */ 1813 j=0; 1814 for(i=0;i<pcbddc->coarse_size;i++){ 1815 if(aux_ins_indices[i]==0){ 1816 aux_ins_indices[i]=-1; 1817 } else { 1818 aux_ins_indices[i]=j; 1819 j++; 1820 } 1821 } 1822 1823 //printf("New details localsizes2 localdispl2\n"); 1824 //for(i=0;i<count_recv;i++) 1825 // printf("(%d %d) ",localsizes2[i],localdispl2[i]); 1826 //printf("\n"); 1827 //printf("aux_ins_indices 2\n"); 1828 //for(i=0;i<pcbddc->coarse_size;i++) 1829 // printf("%d ",aux_ins_indices[i]); 1830 //printf("\n"); 1831 //printf("ins_local_primal_indices\n"); 1832 //for(i=0;i<ins_local_primal_size;i++) 1833 // printf("%d ",ins_local_primal_indices[i]); 1834 //printf("\n"); 1835 //printf("coarse_submat_vals\n"); 1836 //for(i=0;i<pcbddc->local_primal_size;i++) 1837 // for(j=0;j<pcbddc->local_primal_size;j++) 1838 // printf("(%lf %d %d)\n",coarse_submat_vals[j*pcbddc->local_primal_size+i],pcbddc->local_primal_indices[i],pcbddc->local_primal_indices[j]); 1839 //printf("\n"); 1840 1841 /* processes partecipating in coarse problem receive matrix data from their friends */ 1842 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); 1843 if(rank_coarse_proc_send_to != MPI_PROC_NULL ) { 1844 send_size=pcbddc->local_primal_size*pcbddc->local_primal_size; 1845 ierr = MPI_Isend(&coarse_submat_vals[0],send_size,MPIU_SCALAR,rank_coarse_proc_send_to,666,prec_comm,&requests[count_recv]);CHKERRQ(ierr); 1846 } 1847 ierr = MPI_Waitall(count_recv+1,requests,MPI_STATUSES_IGNORE);CHKERRQ(ierr); 1848 1849 //if(coarse_color == 0) { 1850 // printf("temp_coarse_mat_vals\n"); 1851 // for(k=0;k<count_recv;k++){ 1852 // printf("---- %d ----\n",ranks_recv[k]); 1853 // for(i=0;i<pcbddc->local_primal_sizes[ranks_recv[k]];i++) 1854 // for(j=0;j<pcbddc->local_primal_sizes[ranks_recv[k]];j++) 1855 // 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]); 1856 // printf("\n"); 1857 // } 1858 //} 1859 /* calculate data to insert in coarse mat */ 1860 ierr = PetscMalloc ( ins_local_primal_size*ins_local_primal_size*sizeof(PetscScalar),&ins_coarse_mat_vals);CHKERRQ(ierr); 1861 PetscMemzero(ins_coarse_mat_vals,ins_local_primal_size*ins_local_primal_size*sizeof(PetscScalar)); 1862 1863 PetscMPIInt rr,kk,lps,lpd; 1864 PetscInt row_ind,col_ind; 1865 for(k=0;k<count_recv;k++){ 1866 rr = ranks_recv[k]; 1867 kk = localdispl2[k]; 1868 lps = pcbddc->local_primal_sizes[rr]; 1869 lpd = pcbddc->local_primal_displacements[rr]; 1870 //printf("Inserting the following indices (received from %d)\n",rr); 1871 for(j=0;j<lps;j++){ 1872 col_ind=aux_ins_indices[pcbddc->replicated_local_primal_indices[lpd+j]]; 1873 for(i=0;i<lps;i++){ 1874 row_ind=aux_ins_indices[pcbddc->replicated_local_primal_indices[lpd+i]]; 1875 //printf("%d %d\n",row_ind,col_ind); 1876 ins_coarse_mat_vals[col_ind*ins_local_primal_size+row_ind]+=temp_coarse_mat_vals[kk+j*lps+i]; 1877 } 1878 } 1879 } 1880 ierr = PetscFree(requests);CHKERRQ(ierr); 1881 ierr = PetscFree(aux_ins_indices);CHKERRQ(ierr); 1882 ierr = PetscFree(temp_coarse_mat_vals);CHKERRQ(ierr); 1883 if(coarse_color == 0) { ierr = PetscFree(ranks_recv);CHKERRQ(ierr); } 1884 1885 /* create local to global mapping needed by coarse MATIS */ 1886 { 1887 IS coarse_IS; 1888 if(coarse_comm != MPI_COMM_NULL ) ierr = MPI_Comm_free(&coarse_comm);CHKERRQ(ierr); 1889 coarse_comm = prec_comm; 1890 active_rank=rank_prec_comm; 1891 ierr = ISCreateGeneral(coarse_comm,ins_local_primal_size,ins_local_primal_indices,PETSC_COPY_VALUES,&coarse_IS);CHKERRQ(ierr); 1892 //ierr = ISSetBlockSize(coarse_IS,bs);CHKERRQ(ierr); 1893 ierr = ISLocalToGlobalMappingCreateIS(coarse_IS,&coarse_ISLG);CHKERRQ(ierr); 1894 ierr = ISDestroy(&coarse_IS);CHKERRQ(ierr); 1895 } 1896 } 1897 if(pcbddc->coarse_problem_type==PARALLEL_BDDC) { 1898 /* arrays for values insertion */ 1899 ins_local_primal_size = pcbddc->local_primal_size; 1900 ierr = PetscMalloc ( ins_local_primal_size*sizeof(PetscMPIInt),&ins_local_primal_indices);CHKERRQ(ierr); 1901 ierr = PetscMalloc ( ins_local_primal_size*ins_local_primal_size*sizeof(PetscScalar),&ins_coarse_mat_vals);CHKERRQ(ierr); 1902 for(j=0;j<ins_local_primal_size;j++){ 1903 ins_local_primal_indices[j]=pcbddc->local_primal_indices[j]; 1904 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]; 1905 } 1906 } 1907 break; 1908 1909 } 1910 1911 case(GATHERS_BDDC): 1912 { 1913 1914 PetscMPIInt mysize,mysize2; 1915 1916 if(rank_prec_comm==active_rank) { 1917 ierr = PetscMalloc ( pcbddc->replicated_primal_size*sizeof(PetscMPIInt),&pcbddc->replicated_local_primal_indices);CHKERRQ(ierr); 1918 pcbddc->replicated_local_primal_values = (PetscScalar*)calloc(pcbddc->replicated_primal_size,sizeof(PetscScalar)); 1919 ierr = PetscMalloc ( size_prec_comm*sizeof(PetscMPIInt),&localsizes2);CHKERRQ(ierr); 1920 ierr = PetscMalloc ( size_prec_comm*sizeof(PetscMPIInt),&localdispl2);CHKERRQ(ierr); 1921 /* arrays for values insertion */ 1922 ins_local_primal_size = pcbddc->coarse_size; 1923 ierr = PetscMalloc ( ins_local_primal_size*sizeof(PetscMPIInt),&ins_local_primal_indices);CHKERRQ(ierr); 1924 ierr = PetscMalloc ( ins_local_primal_size*ins_local_primal_size*sizeof(PetscScalar),&ins_coarse_mat_vals);CHKERRQ(ierr); 1925 for(i=0;i<size_prec_comm;i++) localsizes2[i]=pcbddc->local_primal_sizes[i]*pcbddc->local_primal_sizes[i]; 1926 localdispl2[0]=0; 1927 for(i=1;i<size_prec_comm;i++) localdispl2[i]=localsizes2[i-1]+localdispl2[i-1]; 1928 j=0; 1929 for(i=0;i<size_prec_comm;i++) j+=localsizes2[i]; 1930 ierr = PetscMalloc ( j*sizeof(PetscScalar),&temp_coarse_mat_vals);CHKERRQ(ierr); 1931 } 1932 1933 mysize=pcbddc->local_primal_size; 1934 mysize2=pcbddc->local_primal_size*pcbddc->local_primal_size; 1935 if(pcbddc->coarse_problem_type == SEQUENTIAL_BDDC){ 1936 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); 1937 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); 1938 } else { 1939 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); 1940 ierr = MPI_Allgatherv(&coarse_submat_vals[0],mysize2,MPIU_SCALAR,&temp_coarse_mat_vals[0],localsizes2,localdispl2,MPIU_SCALAR,prec_comm);CHKERRQ(ierr); 1941 } 1942 1943 /* free data structures no longer needed and allocate some space which will be needed in BDDC application */ 1944 if(rank_prec_comm==active_rank) { 1945 PetscInt offset,offset2,row_ind,col_ind; 1946 for(j=0;j<ins_local_primal_size;j++){ 1947 ins_local_primal_indices[j]=j; 1948 for(i=0;i<ins_local_primal_size;i++) ins_coarse_mat_vals[j*ins_local_primal_size+i]=0.0; 1949 } 1950 for(k=0;k<size_prec_comm;k++){ 1951 offset=pcbddc->local_primal_displacements[k]; 1952 offset2=localdispl2[k]; 1953 for(j=0;j<pcbddc->local_primal_sizes[k];j++){ 1954 col_ind=pcbddc->replicated_local_primal_indices[offset+j]; 1955 for(i=0;i<pcbddc->local_primal_sizes[k];i++){ 1956 row_ind=pcbddc->replicated_local_primal_indices[offset+i]; 1957 ins_coarse_mat_vals[col_ind*pcbddc->coarse_size+row_ind]+=temp_coarse_mat_vals[offset2+j*pcbddc->local_primal_sizes[k]+i]; 1958 } 1959 } 1960 } 1961 } 1962 break; 1963 }//switch on coarse problem and communications associated with finished 1964 } 1965 1966 /* Now create and fill up coarse matrix */ 1967 if( rank_prec_comm == active_rank ) { 1968 if(pcbddc->coarse_problem_type != MULTILEVEL_BDDC) { 1969 ierr = MatCreate(coarse_comm,&pcbddc->coarse_mat);CHKERRQ(ierr); 1970 ierr = MatSetSizes(pcbddc->coarse_mat,PETSC_DECIDE,PETSC_DECIDE,pcbddc->coarse_size,pcbddc->coarse_size);CHKERRQ(ierr); 1971 ierr = MatSetType(pcbddc->coarse_mat,coarse_mat_type);CHKERRQ(ierr); 1972 ierr = MatSetOption(pcbddc->coarse_mat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 1973 ierr = MatSetOption(pcbddc->coarse_mat,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); //local values stored in column major 1974 } else { 1975 Mat matis_coarse_local_mat; 1976 ierr = MatCreateIS(coarse_comm,PETSC_DECIDE,PETSC_DECIDE,pcbddc->coarse_size,pcbddc->coarse_size,coarse_ISLG,&pcbddc->coarse_mat);CHKERRQ(ierr); 1977 ierr = MatISGetLocalMat(pcbddc->coarse_mat,&matis_coarse_local_mat);CHKERRQ(ierr); 1978 ierr = MatSetOption(matis_coarse_local_mat,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); //local values stored in column major 1979 ierr = MatSetOption(matis_coarse_local_mat,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 1980 ierr = MatSetOption(matis_coarse_local_mat,MAT_USE_INODES,PETSC_FALSE);CHKERRQ(ierr); 1981 } 1982 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); 1983 ierr = MatAssemblyBegin(pcbddc->coarse_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1984 ierr = MatAssemblyEnd(pcbddc->coarse_mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1985 if(pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 1986 Mat matis_coarse_local_mat; 1987 ierr = MatISGetLocalMat(pcbddc->coarse_mat,&matis_coarse_local_mat);CHKERRQ(ierr); 1988 ierr = MatSetBlockSize(matis_coarse_local_mat,bs);CHKERRQ(ierr); 1989 } 1990 1991 ierr = MatGetVecs(pcbddc->coarse_mat,&pcbddc->coarse_vec,&pcbddc->coarse_rhs);CHKERRQ(ierr); 1992 /* Preconditioner for coarse problem */ 1993 ierr = KSPCreate(coarse_comm,&pcbddc->coarse_ksp);CHKERRQ(ierr); 1994 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->coarse_ksp,(PetscObject)pc,1);CHKERRQ(ierr); 1995 ierr = KSPSetOperators(pcbddc->coarse_ksp,pcbddc->coarse_mat,pcbddc->coarse_mat,SAME_PRECONDITIONER);CHKERRQ(ierr); 1996 ierr = KSPSetType(pcbddc->coarse_ksp,coarse_ksp_type);CHKERRQ(ierr); 1997 ierr = KSPGetPC(pcbddc->coarse_ksp,&pc_temp);CHKERRQ(ierr); 1998 ierr = PCSetType(pc_temp,coarse_pc_type);CHKERRQ(ierr); 1999 //ierr = PCSetOptionsPrefix(pcbddc->coarse_pc,"pc_bddc_coarse_");CHKERRQ(ierr); 2000 /* Allow user's customization */ 2001 ierr = KSPSetFromOptions(pcbddc->coarse_ksp);CHKERRQ(ierr); 2002 /* Set Up PC for coarse problem BDDC */ 2003 //if(pcbddc->dbg_flag && pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 2004 if(pcbddc->coarse_problem_type == MULTILEVEL_BDDC) { 2005 if(dbg_flag) { 2006 ierr = PetscViewerASCIIPrintf(viewer,"----------------Setting up a new level---------------\n");CHKERRQ(ierr); 2007 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2008 } 2009 ierr = PCBDDCSetCoarseProblemType(pc_temp,MULTILEVEL_BDDC);CHKERRQ(ierr); 2010 } 2011 ierr = KSPSetUp(pcbddc->coarse_ksp);CHKERRQ(ierr); 2012 } 2013 if(pcbddc->coarse_communications_type == SCATTERS_BDDC) { 2014 IS local_IS,global_IS; 2015 ierr = ISCreateStride(PETSC_COMM_SELF,pcbddc->local_primal_size,0,1,&local_IS);CHKERRQ(ierr); 2016 ierr = ISCreateGeneral(PETSC_COMM_SELF,pcbddc->local_primal_size,pcbddc->local_primal_indices,PETSC_COPY_VALUES,&global_IS);CHKERRQ(ierr); 2017 ierr = VecScatterCreate(pcbddc->vec1_P,local_IS,pcbddc->coarse_vec,global_IS,&pcbddc->coarse_loc_to_glob);CHKERRQ(ierr); 2018 ierr = ISDestroy(&local_IS);CHKERRQ(ierr); 2019 ierr = ISDestroy(&global_IS);CHKERRQ(ierr); 2020 } 2021 2022 2023 /* Check condition number of coarse problem */ 2024 if( pcbddc->coarse_problem_type == MULTILEVEL_BDDC && dbg_flag && rank_prec_comm == active_rank ) { 2025 PetscScalar m_one=-1.0; 2026 PetscReal infty_error,lambda_min,lambda_max; 2027 2028 ierr = KSPSetComputeSingularValues(pcbddc->coarse_ksp,PETSC_TRUE);CHKERRQ(ierr); 2029 ierr = VecSetRandom(pcbddc->coarse_rhs,PETSC_NULL);CHKERRQ(ierr); 2030 ierr = MatMult(pcbddc->coarse_mat,pcbddc->coarse_rhs,pcbddc->coarse_vec);CHKERRQ(ierr); 2031 ierr = MatMult(pcbddc->coarse_mat,pcbddc->coarse_vec,pcbddc->coarse_rhs);CHKERRQ(ierr); 2032 ierr = KSPSolve(pcbddc->coarse_ksp,pcbddc->coarse_rhs,pcbddc->coarse_rhs);CHKERRQ(ierr); 2033 ierr = KSPComputeExtremeSingularValues(pcbddc->coarse_ksp,&lambda_max,&lambda_min);CHKERRQ(ierr); 2034 ierr = VecAXPY(pcbddc->coarse_rhs,m_one,pcbddc->coarse_vec);CHKERRQ(ierr); 2035 ierr = VecNorm(pcbddc->coarse_rhs,NORM_INFINITY,&infty_error);CHKERRQ(ierr); 2036 ierr = PetscViewerASCIIPrintf(viewer,"Coarse problem eigenvalues: % 1.14e %1.14e\n",lambda_min,lambda_max);CHKERRQ(ierr); 2037 ierr = PetscViewerASCIIPrintf(viewer,"Coarse problem infty_error: %1.14e\n",infty_error);CHKERRQ(ierr); 2038 ierr = KSPSetComputeSingularValues(pcbddc->coarse_ksp,PETSC_FALSE);CHKERRQ(ierr); 2039 } 2040 2041 /* free data structures no longer needed */ 2042 if(coarse_ISLG) { ierr = ISLocalToGlobalMappingDestroy(&coarse_ISLG);CHKERRQ(ierr); } 2043 if(ins_local_primal_indices) { ierr = PetscFree(ins_local_primal_indices);CHKERRQ(ierr); } 2044 if(ins_coarse_mat_vals) { ierr = PetscFree(ins_coarse_mat_vals);CHKERRQ(ierr);} 2045 if(localsizes2) { ierr = PetscFree(localsizes2);CHKERRQ(ierr);} 2046 if(localdispl2) { ierr = PetscFree(localdispl2);CHKERRQ(ierr);} 2047 if(temp_coarse_mat_vals) { ierr = PetscFree(temp_coarse_mat_vals);CHKERRQ(ierr);} 2048 2049 PetscFunctionReturn(0); 2050 } 2051 2052 #undef __FUNCT__ 2053 #define __FUNCT__ "PCBDDCManageLocalBoundaries" 2054 static PetscErrorCode PCBDDCManageLocalBoundaries(PC pc) 2055 { 2056 2057 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 2058 PC_IS *pcis = (PC_IS*)pc->data; 2059 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 2060 PetscInt *distinct_values; 2061 PetscInt **neighbours_set; 2062 PetscInt bs,ierr,i,j,s,k,iindex; 2063 PetscInt total_counts; 2064 PetscBool flg_row; 2065 PCBDDCGraph mat_graph; 2066 PetscInt neumann_bsize; 2067 const PetscInt* neumann_nodes; 2068 Mat mat_adj; 2069 2070 PetscFunctionBegin; 2071 2072 ierr = MatGetBlockSize(matis->A,&bs);CHKERRQ(ierr); 2073 // allocate and initialize needed graph structure 2074 ierr = PetscMalloc(sizeof(*mat_graph),&mat_graph);CHKERRQ(ierr); 2075 ierr = MatConvert(matis->A,MATMPIADJ,MAT_INITIAL_MATRIX,&mat_adj);CHKERRQ(ierr); 2076 ierr = MatGetRowIJ(mat_adj,0,PETSC_FALSE,PETSC_FALSE,&mat_graph->nvtxs,&mat_graph->xadj,&mat_graph->adjncy,&flg_row);CHKERRQ(ierr); 2077 if(!flg_row) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatGetRowIJ called from PCBDDCManageLocalBoundaries.\n"); 2078 ierr = PetscMalloc(mat_graph->nvtxs*sizeof(PetscInt),&mat_graph->where);CHKERRQ(ierr); 2079 ierr = PetscMalloc(mat_graph->nvtxs*sizeof(PetscInt),&mat_graph->count);CHKERRQ(ierr); 2080 ierr = PetscMalloc(mat_graph->nvtxs*sizeof(PetscInt),&mat_graph->which_dof);CHKERRQ(ierr); 2081 ierr = PetscMalloc(mat_graph->nvtxs*sizeof(PetscInt),&mat_graph->queue);CHKERRQ(ierr); 2082 ierr = PetscMalloc((mat_graph->nvtxs+1)*sizeof(PetscInt),&mat_graph->cptr);CHKERRQ(ierr); 2083 ierr = PetscMalloc(mat_graph->nvtxs*sizeof(PetscBool),&mat_graph->touched);CHKERRQ(ierr); 2084 for(i=0;i<mat_graph->nvtxs;i++){ 2085 mat_graph->count[i]=0; 2086 mat_graph->touched[i]=PETSC_FALSE; 2087 } 2088 for(i=0;i<mat_graph->nvtxs/bs;i++) { 2089 for(s=0;s<bs;s++) { 2090 mat_graph->which_dof[i*bs+s]=s; 2091 } 2092 } 2093 //printf("nvtxs = %d\n",mat_graph->nvtxs); 2094 //printf("these are my IS data with n_neigh = %d\n",pcis->n_neigh); 2095 //for(i=0;i<pcis->n_neigh;i++){ 2096 // printf("number of shared nodes with rank %d is %d \n",pcis->neigh[i],pcis->n_shared[i]); 2097 // } 2098 2099 total_counts=0; 2100 for(i=0;i<pcis->n_neigh;i++){ 2101 s=pcis->n_shared[i]; 2102 total_counts+=s; 2103 //printf("computing neigh %d (rank = %d, n_shared = %d)\n",i,pcis->neigh[i],s); 2104 for(j=0;j<s;j++){ 2105 mat_graph->count[pcis->shared[i][j]] += 1; 2106 } 2107 } 2108 /* Take into account Neumann data incrementing number of sharing subdomains for all but faces nodes lying on the interface */ 2109 if(pcbddc->NeumannBoundaries) { 2110 ierr = ISGetLocalSize(pcbddc->NeumannBoundaries,&neumann_bsize);CHKERRQ(ierr); 2111 ierr = ISGetIndices(pcbddc->NeumannBoundaries,&neumann_nodes);CHKERRQ(ierr); 2112 for(i=0;i<neumann_bsize;i++){ 2113 iindex = neumann_nodes[i]; 2114 if(mat_graph->count[iindex] > 2){ 2115 mat_graph->count[iindex]+=1; 2116 total_counts++; 2117 } 2118 } 2119 } 2120 2121 ierr = PetscMalloc(mat_graph->nvtxs*sizeof(PetscInt*),&neighbours_set);CHKERRQ(ierr); 2122 if(mat_graph->nvtxs) ierr = PetscMalloc(total_counts*sizeof(PetscInt),&neighbours_set[0]);CHKERRQ(ierr); 2123 for(i=1;i<mat_graph->nvtxs;i++) neighbours_set[i]=neighbours_set[i-1]+mat_graph->count[i-1]; 2124 2125 for(i=0;i<mat_graph->nvtxs;i++) mat_graph->count[i]=0; 2126 for(i=0;i<pcis->n_neigh;i++){ 2127 s=pcis->n_shared[i]; 2128 for(j=0;j<s;j++) { 2129 k=pcis->shared[i][j]; 2130 neighbours_set[k][mat_graph->count[k]] = pcis->neigh[i]; 2131 mat_graph->count[k]+=1; 2132 } 2133 } 2134 /* set -1 fake neighbour */ 2135 if(pcbddc->NeumannBoundaries) { 2136 for(i=0;i<neumann_bsize;i++){ 2137 iindex = neumann_nodes[i]; 2138 if( mat_graph->count[iindex] > 2){ 2139 neighbours_set[iindex][mat_graph->count[iindex]] = -1; //An additional fake neighbour (with rank -1) 2140 mat_graph->count[iindex]+=1; 2141 } 2142 } 2143 ierr = ISRestoreIndices(pcbddc->NeumannBoundaries,&neumann_nodes);CHKERRQ(ierr); 2144 } 2145 2146 /* sort sharing subdomains */ 2147 for(i=0;i<mat_graph->nvtxs;i++) { ierr = PetscSortInt(mat_graph->count[i],neighbours_set[i]);CHKERRQ(ierr); } 2148 2149 /* Prepare for FindConnectedComponents 2150 Vertices will be eliminated later (if needed) */ 2151 PetscInt nodes_touched=0; 2152 for(i=0;i<mat_graph->nvtxs;i++){ 2153 if(!mat_graph->count[i]){ //internal nodes 2154 mat_graph->touched[i]=PETSC_TRUE; 2155 mat_graph->where[i]=0; 2156 nodes_touched++; 2157 } 2158 if(pcbddc->faces_flag){ 2159 if(mat_graph->count[i]>2){ //all but faces 2160 mat_graph->touched[i]=PETSC_TRUE; 2161 mat_graph->where[i]=0; 2162 nodes_touched++; 2163 } 2164 } 2165 if(pcbddc->edges_flag){ 2166 if(mat_graph->count[i]==2){ //faces 2167 mat_graph->touched[i]=PETSC_TRUE; 2168 mat_graph->where[i]=0; 2169 nodes_touched++; 2170 } 2171 } 2172 } 2173 2174 PetscInt rvalue=1; 2175 PetscBool same_set; 2176 mat_graph->ncmps = 0; 2177 while(nodes_touched<mat_graph->nvtxs) { 2178 // find first untouched node in local ordering 2179 i=0; 2180 while(mat_graph->touched[i]) i++; 2181 mat_graph->touched[i]=PETSC_TRUE; 2182 mat_graph->where[i]=rvalue; 2183 nodes_touched++; 2184 // now find other connected nodes shared by the same set of subdomains 2185 for(j=i+1;j<mat_graph->nvtxs;j++){ 2186 // check for same number of sharing subdomains and dof number 2187 if(mat_graph->count[i]==mat_graph->count[j] && mat_graph->which_dof[i] == mat_graph->which_dof[j] ){ 2188 // check for same set of sharing subdomains 2189 same_set=PETSC_TRUE; 2190 for(k=0;k<mat_graph->count[j];k++){ 2191 if(neighbours_set[i][k]!=neighbours_set[j][k]) { 2192 same_set=PETSC_FALSE; 2193 } 2194 } 2195 // OK, I found a friend of mine 2196 if(same_set) { 2197 mat_graph->where[j]=rvalue; 2198 mat_graph->touched[j]=PETSC_TRUE; 2199 nodes_touched++; 2200 } 2201 } 2202 } 2203 rvalue++; 2204 } 2205 // printf("where and count contains %d distinct values\n",rvalue); 2206 // for(j=0;j<mat_graph->nvtxs;j++) 2207 // printf("[%d %d %d]\n",j,mat_graph->where[j],mat_graph->count[j]); 2208 2209 if(mat_graph->nvtxs) { 2210 ierr = PetscFree(neighbours_set[0]);CHKERRQ(ierr); 2211 ierr = PetscFree(neighbours_set);CHKERRQ(ierr); 2212 } 2213 2214 rvalue--; 2215 ierr = PetscMalloc ( rvalue*sizeof(PetscInt),&distinct_values);CHKERRQ(ierr); 2216 for(i=0;i<rvalue;i++) distinct_values[i]=i+1; //initializiation 2217 if(rvalue) ierr = PCBDDCFindConnectedComponents(mat_graph, rvalue, distinct_values); 2218 //printf("total number of connected components %d \n",mat_graph->ncmps); 2219 //for (i=0; i<mat_graph->ncmps; i++) { 2220 // 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]]); 2221 //} 2222 PetscInt nfc=0; 2223 PetscInt nec=0; 2224 PetscInt nvc=0; 2225 for (i=0; i<mat_graph->ncmps; i++) { 2226 // sort each connected component (by local ordering) 2227 ierr = PetscSortInt(mat_graph->cptr[i+1]-mat_graph->cptr[i],&mat_graph->queue[mat_graph->cptr[i]]);CHKERRQ(ierr); 2228 // count edge and faces 2229 if( !pcbddc->vertices_flag ) { 2230 if( mat_graph->cptr[i+1]-mat_graph->cptr[i] > 1 ){ 2231 if(mat_graph->count[mat_graph->queue[mat_graph->cptr[i]]]==2){ 2232 nfc++; 2233 } else { 2234 nec++; 2235 } 2236 } 2237 } 2238 // count vertices 2239 if( !pcbddc->constraints_flag ){ 2240 if( mat_graph->cptr[i+1]-mat_graph->cptr[i] == 1 ){ 2241 nvc++; 2242 } 2243 } 2244 } 2245 2246 pcbddc->n_constraints = nec+nfc; 2247 pcbddc->n_vertices = nvc; 2248 2249 if(pcbddc->n_constraints){ 2250 /* allocate space for local constraints of BDDC */ 2251 ierr = PetscMalloc (pcbddc->n_constraints*sizeof(PetscInt*),&pcbddc->indices_to_constraint);CHKERRQ(ierr); 2252 ierr = PetscMalloc (pcbddc->n_constraints*sizeof(PetscScalar*),&pcbddc->quadrature_constraint);CHKERRQ(ierr); 2253 ierr = PetscMalloc (pcbddc->n_constraints*sizeof(PetscInt),&pcbddc->sizes_of_constraint);CHKERRQ(ierr); 2254 k=0; 2255 for (i=0; i<mat_graph->ncmps; i++) { 2256 if( mat_graph->cptr[i+1]-mat_graph->cptr[i] > 1 ){ 2257 pcbddc->sizes_of_constraint[k] = mat_graph->cptr[i+1]-mat_graph->cptr[i]; 2258 k++; 2259 } 2260 } 2261 // printf("check constraints %d (should be %d)\n",k,pcbddc->n_constraints); 2262 // for(i=0;i<k;i++) 2263 // printf("%d ",pcbddc->sizes_of_constraint[i]); 2264 // printf("\n"); 2265 k=0; 2266 for (i=0; i<pcbddc->n_constraints; i++) k+=pcbddc->sizes_of_constraint[i]; 2267 ierr = PetscMalloc (k*sizeof(PetscInt),&pcbddc->indices_to_constraint[0]);CHKERRQ(ierr); 2268 ierr = PetscMalloc (k*sizeof(PetscScalar),&pcbddc->quadrature_constraint[0]);CHKERRQ(ierr); 2269 for (i=1; i<pcbddc->n_constraints; i++) { 2270 pcbddc->indices_to_constraint[i] = pcbddc->indices_to_constraint[i-1] + pcbddc->sizes_of_constraint[i-1]; 2271 pcbddc->quadrature_constraint[i] = pcbddc->quadrature_constraint[i-1] + pcbddc->sizes_of_constraint[i-1]; 2272 } 2273 k=0; 2274 PetscScalar quad_value; 2275 for (i=0; i<mat_graph->ncmps; i++) { 2276 if( mat_graph->cptr[i+1]-mat_graph->cptr[i] > 1 ){ 2277 quad_value=1.0/( (PetscScalar) (mat_graph->cptr[i+1]-mat_graph->cptr[i]) ); 2278 for(j=0;j<mat_graph->cptr[i+1]-mat_graph->cptr[i];j++) { 2279 pcbddc->indices_to_constraint[k][j] = mat_graph->queue[mat_graph->cptr[i]+j]; 2280 pcbddc->quadrature_constraint[k][j] = quad_value; 2281 } 2282 k++; 2283 } 2284 } 2285 } 2286 if(pcbddc->n_vertices){ 2287 /* allocate space for local vertices of BDDC */ 2288 ierr = PetscMalloc (pcbddc->n_vertices*sizeof(PetscInt),&pcbddc->vertices);CHKERRQ(ierr); 2289 k=0; 2290 for (i=0; i<mat_graph->ncmps; i++) { 2291 if( mat_graph->cptr[i+1]-mat_graph->cptr[i] == 1 ){ 2292 pcbddc->vertices[k] = mat_graph->queue[mat_graph->cptr[i]]; 2293 k++; 2294 } 2295 } 2296 // sort vertex set (by local ordering) 2297 ierr = PetscSortInt(pcbddc->n_vertices,pcbddc->vertices);CHKERRQ(ierr); 2298 } 2299 2300 if(pcbddc->dbg_flag) { 2301 PetscViewer viewer=pcbddc->dbg_viewer; 2302 2303 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"--------------------------------------------------------------\n");CHKERRQ(ierr); 2304 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Details from PCBDDCManageLocalBoundaries for subdomain %04d\n",PetscGlobalRank);CHKERRQ(ierr); 2305 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"--------------------------------------------------------------\n");CHKERRQ(ierr); 2306 /* PetscViewerASCIISynchronizedPrintf(viewer,"Graph (adjacency structure) of local Neumann mat\n"); 2307 PetscViewerASCIISynchronizedPrintf(viewer,"--------------------------------------------------------------\n"); 2308 for(i=0;i<mat_graph->nvtxs;i++) { 2309 PetscViewerASCIISynchronizedPrintf(viewer,"Nodes connected to node number %d are %d\n",i,mat_graph->xadj[i+1]-mat_graph->xadj[i]); 2310 for(j=mat_graph->xadj[i];j<mat_graph->xadj[i+1];j++){ 2311 PetscViewerASCIISynchronizedPrintf(viewer,"%d ",mat_graph->adjncy[j]); 2312 } 2313 PetscViewerASCIISynchronizedPrintf(viewer,"\n--------------------------------------------------------------\n"); 2314 }*/ 2315 // TODO: APPLY Local to Global Mapping from IS object? 2316 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Matrix graph has %d connected components", mat_graph->ncmps);CHKERRQ(ierr); 2317 for(i=0;i<mat_graph->ncmps;i++) { 2318 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); 2319 for (j=mat_graph->cptr[i]; j<mat_graph->cptr[i+1]; j++){ 2320 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"%d ",mat_graph->queue[j]);CHKERRQ(ierr); 2321 } 2322 } 2323 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\n--------------------------------------------------------------\n");CHKERRQ(ierr); 2324 if( pcbddc->n_vertices ) { ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d detected %02d local vertices\n",PetscGlobalRank,pcbddc->n_vertices);CHKERRQ(ierr); } 2325 if( nfc ) { ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d detected %02d local faces\n",PetscGlobalRank,nfc);CHKERRQ(ierr); } 2326 if( nec ) { ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Subdomain %04d detected %02d local edges\n",PetscGlobalRank,nec);CHKERRQ(ierr); } 2327 if( pcbddc->n_vertices ) { ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Indices for subdomain vertices follows\n",PetscGlobalRank,pcbddc->n_vertices);CHKERRQ(ierr); } 2328 for(i=0;i<pcbddc->n_vertices;i++){ ierr = PetscViewerASCIISynchronizedPrintf(viewer,"%d ",pcbddc->vertices[i]);CHKERRQ(ierr); } 2329 if( pcbddc->n_vertices ) { ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\n");CHKERRQ(ierr); } 2330 /* if( pcbddc->n_constraints ) PetscViewerASCIISynchronizedPrintf(viewer,"Indices and quadrature constraints"); 2331 for(i=0;i<pcbddc->n_constraints;i++){ 2332 PetscViewerASCIISynchronizedPrintf(viewer,"\nConstraint number %d\n",i); 2333 for(j=0;j<pcbddc->sizes_of_constraint[i];j++) { 2334 PetscViewerASCIISynchronizedPrintf(viewer,"(%d, %f) ",pcbddc->indices_to_constraint[i][j],pcbddc->quadrature_constraint[i][j]); 2335 } 2336 } 2337 if( pcbddc->n_constraints ) PetscViewerASCIISynchronizedPrintf(viewer,"\n");*/ 2338 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 2339 } 2340 2341 // Restore CSR structure into sequantial matrix and free memory space no longer needed 2342 ierr = MatRestoreRowIJ(mat_adj,0,PETSC_FALSE,PETSC_TRUE,&mat_graph->nvtxs,&mat_graph->xadj,&mat_graph->adjncy,&flg_row);CHKERRQ(ierr); 2343 ierr = MatDestroy(&mat_adj);CHKERRQ(ierr); 2344 if(!flg_row) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called from PCBDDCManageLocalBoundaries.\n"); 2345 ierr = PetscFree(distinct_values);CHKERRQ(ierr); 2346 // Free graph structure 2347 if(mat_graph->nvtxs){ 2348 ierr = PetscFree(mat_graph->where);CHKERRQ(ierr); 2349 ierr = PetscFree(mat_graph->touched);CHKERRQ(ierr); 2350 ierr = PetscFree(mat_graph->which_dof);CHKERRQ(ierr); 2351 ierr = PetscFree(mat_graph->queue);CHKERRQ(ierr); 2352 ierr = PetscFree(mat_graph->cptr);CHKERRQ(ierr); 2353 ierr = PetscFree(mat_graph->count);CHKERRQ(ierr); 2354 } 2355 ierr = PetscFree(mat_graph);CHKERRQ(ierr); 2356 2357 PetscFunctionReturn(0); 2358 2359 } 2360 2361 /* -------------------------------------------------------------------------- */ 2362 2363 /* The following code has been adapted from function IsConnectedSubdomain contained 2364 in source file contig.c of METIS library (version 5.0.1) */ 2365 2366 #undef __FUNCT__ 2367 #define __FUNCT__ "PCBDDCFindConnectedComponents" 2368 static PetscErrorCode PCBDDCFindConnectedComponents(PCBDDCGraph graph, PetscInt n_dist, PetscInt *dist_vals) 2369 { 2370 PetscInt i, j, k, nvtxs, first, last, nleft, ncmps,pid,cum_queue,n,ncmps_pid; 2371 PetscInt *xadj, *adjncy, *where, *queue; 2372 PetscInt *cptr; 2373 PetscBool *touched; 2374 2375 PetscFunctionBegin; 2376 2377 nvtxs = graph->nvtxs; 2378 xadj = graph->xadj; 2379 adjncy = graph->adjncy; 2380 where = graph->where; 2381 touched = graph->touched; 2382 queue = graph->queue; 2383 cptr = graph->cptr; 2384 2385 for (i=0; i<nvtxs; i++) 2386 touched[i] = PETSC_FALSE; 2387 2388 cum_queue=0; 2389 ncmps=0; 2390 2391 for(n=0; n<n_dist; n++) { 2392 pid = dist_vals[n]; 2393 nleft = 0; 2394 for (i=0; i<nvtxs; i++) { 2395 if (where[i] == pid) 2396 nleft++; 2397 } 2398 for (i=0; i<nvtxs; i++) { 2399 if (where[i] == pid) 2400 break; 2401 } 2402 touched[i] = PETSC_TRUE; 2403 queue[cum_queue] = i; 2404 first = 0; last = 1; 2405 cptr[ncmps] = cum_queue; /* This actually points to queue */ 2406 ncmps_pid = 0; 2407 while (first != nleft) { 2408 if (first == last) { /* Find another starting vertex */ 2409 cptr[++ncmps] = first+cum_queue; 2410 ncmps_pid++; 2411 for (i=0; i<nvtxs; i++) { 2412 if (where[i] == pid && !touched[i]) 2413 break; 2414 } 2415 queue[cum_queue+last] = i; 2416 last++; 2417 touched[i] = PETSC_TRUE; 2418 } 2419 i = queue[cum_queue+first]; 2420 first++; 2421 for (j=xadj[i]; j<xadj[i+1]; j++) { 2422 k = adjncy[j]; 2423 if (where[k] == pid && !touched[k]) { 2424 queue[cum_queue+last] = k; 2425 last++; 2426 touched[k] = PETSC_TRUE; 2427 } 2428 } 2429 } 2430 cptr[++ncmps] = first+cum_queue; 2431 ncmps_pid++; 2432 cum_queue=cptr[ncmps]; 2433 } 2434 graph->ncmps = ncmps; 2435 2436 PetscFunctionReturn(0); 2437 } 2438 2439