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