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