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