1 #include <petsc/private/petscimpl.h> 2 #include <../src/ksp/pc/impls/bddc/bddcprivate.h> 3 #include <../src/ksp/pc/impls/bddc/bddcstructs.h> 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "PCBDDCGraphGetDirichletDofsB" 7 PetscErrorCode PCBDDCGraphGetDirichletDofsB(PCBDDCGraph graph, IS* dirdofs) 8 { 9 PetscErrorCode ierr; 10 11 PetscFunctionBegin; 12 if (graph->dirdofsB) { 13 ierr = PetscObjectReference((PetscObject)graph->dirdofsB);CHKERRQ(ierr); 14 } else if (graph->has_dirichlet) { 15 PetscInt i,size; 16 PetscInt *dirdofs_idxs; 17 18 size = 0; 19 for (i=0;i<graph->nvtxs;i++) { 20 if (graph->count[i] && graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) size++; 21 } 22 23 ierr = PetscMalloc1(size,&dirdofs_idxs);CHKERRQ(ierr); 24 size = 0; 25 for (i=0;i<graph->nvtxs;i++) { 26 if (graph->count[i] && graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) dirdofs_idxs[size++] = i; 27 } 28 ierr = ISCreateGeneral(PETSC_COMM_SELF,size,dirdofs_idxs,PETSC_OWN_POINTER,&graph->dirdofsB);CHKERRQ(ierr); 29 ierr = PetscObjectReference((PetscObject)graph->dirdofsB);CHKERRQ(ierr); 30 } 31 *dirdofs = graph->dirdofsB; 32 PetscFunctionReturn(0); 33 } 34 35 #undef __FUNCT__ 36 #define __FUNCT__ "PCBDDCGraphGetDirichletDofs" 37 PetscErrorCode PCBDDCGraphGetDirichletDofs(PCBDDCGraph graph, IS* dirdofs) 38 { 39 PetscErrorCode ierr; 40 41 PetscFunctionBegin; 42 if (graph->dirdofs) { 43 ierr = PetscObjectReference((PetscObject)graph->dirdofs);CHKERRQ(ierr); 44 } else if (graph->has_dirichlet) { 45 PetscInt i,size; 46 PetscInt *dirdofs_idxs; 47 48 size = 0; 49 for (i=0;i<graph->nvtxs;i++) { 50 if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) size++; 51 } 52 53 ierr = PetscMalloc1(size,&dirdofs_idxs);CHKERRQ(ierr); 54 size = 0; 55 for (i=0;i<graph->nvtxs;i++) { 56 if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK) dirdofs_idxs[size++] = i; 57 } 58 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)graph->l2gmap),size,dirdofs_idxs,PETSC_OWN_POINTER,&graph->dirdofs);CHKERRQ(ierr); 59 ierr = PetscObjectReference((PetscObject)graph->dirdofs);CHKERRQ(ierr); 60 } 61 *dirdofs = graph->dirdofs; 62 PetscFunctionReturn(0); 63 } 64 65 #undef __FUNCT__ 66 #define __FUNCT__ "PCBDDCGraphASCIIView" 67 PetscErrorCode PCBDDCGraphASCIIView(PCBDDCGraph graph, PetscInt verbosity_level, PetscViewer viewer) 68 { 69 PetscInt i,j,tabs; 70 PetscInt* queue_in_global_numbering; 71 PetscErrorCode ierr; 72 73 PetscFunctionBegin; 74 ierr = PetscViewerASCIIPushSynchronized(viewer);CHKERRQ(ierr); 75 ierr = PetscViewerASCIIGetTab(viewer,&tabs);CHKERRQ(ierr); 76 ierr = PetscViewerASCIIPrintf(viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 77 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 78 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Local BDDC graph for subdomain %04d\n",PetscGlobalRank);CHKERRQ(ierr); 79 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Number of vertices %d\n",graph->nvtxs);CHKERRQ(ierr); 80 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Custom minimal size %d\n",graph->custom_minimal_size);CHKERRQ(ierr); 81 if (verbosity_level > 2) { 82 for (i=0;i<graph->nvtxs;i++) { 83 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"%d:\n",i);CHKERRQ(ierr); 84 ierr = PetscViewerASCIISynchronizedPrintf(viewer," which_dof: %d\n",graph->which_dof[i]);CHKERRQ(ierr); 85 ierr = PetscViewerASCIISynchronizedPrintf(viewer," special_dof: %d\n",graph->special_dof[i]);CHKERRQ(ierr); 86 ierr = PetscViewerASCIISynchronizedPrintf(viewer," neighbours: %d\n",graph->count[i]);CHKERRQ(ierr); 87 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 88 if (graph->count[i]) { 89 ierr = PetscViewerASCIISynchronizedPrintf(viewer," set of neighbours:");CHKERRQ(ierr); 90 for (j=0;j<graph->count[i];j++) { 91 ierr = PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->neighbours_set[i][j]);CHKERRQ(ierr); 92 } 93 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\n");CHKERRQ(ierr); 94 } 95 ierr = PetscViewerASCIISetTab(viewer,tabs);CHKERRQ(ierr); 96 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 97 if (graph->mirrors) { 98 ierr = PetscViewerASCIISynchronizedPrintf(viewer," mirrors: %d\n",graph->mirrors[i]);CHKERRQ(ierr); 99 if (graph->mirrors[i]) { 100 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 101 ierr = PetscViewerASCIISynchronizedPrintf(viewer," set of mirrors:");CHKERRQ(ierr); 102 for (j=0;j<graph->mirrors[i];j++) { 103 ierr = PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->mirrors_set[i][j]);CHKERRQ(ierr); 104 } 105 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\n");CHKERRQ(ierr); 106 ierr = PetscViewerASCIISetTab(viewer,tabs);CHKERRQ(ierr); 107 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 108 } 109 } 110 if (verbosity_level > 3) { 111 if (graph->xadj && graph->adjncy) { 112 ierr = PetscViewerASCIISynchronizedPrintf(viewer," local adj list:");CHKERRQ(ierr); 113 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 114 for (j=graph->xadj[i];j<graph->xadj[i+1];j++) { 115 ierr = PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->adjncy[j]);CHKERRQ(ierr); 116 } 117 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\n");CHKERRQ(ierr); 118 ierr = PetscViewerASCIISetTab(viewer,tabs);CHKERRQ(ierr); 119 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 120 } else { 121 ierr = PetscViewerASCIISynchronizedPrintf(viewer," no adj info\n");CHKERRQ(ierr); 122 } 123 } 124 if (graph->n_local_subs) { 125 ierr = PetscViewerASCIISynchronizedPrintf(viewer," local sub id: %d\n",graph->local_subs[i]);CHKERRQ(ierr); 126 } 127 ierr = PetscViewerASCIISynchronizedPrintf(viewer," interface subset id: %d\n",graph->subset[i]);CHKERRQ(ierr); 128 if (graph->subset[i] && graph->subset_ncc) { 129 ierr = PetscViewerASCIISynchronizedPrintf(viewer," ncc for subset: %d\n",graph->subset_ncc[graph->subset[i]-1]);CHKERRQ(ierr); 130 } 131 } 132 } 133 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"Total number of connected components %d\n",graph->ncc);CHKERRQ(ierr); 134 ierr = PetscMalloc1(graph->cptr[graph->ncc],&queue_in_global_numbering);CHKERRQ(ierr); 135 ierr = ISLocalToGlobalMappingApply(graph->l2gmap,graph->cptr[graph->ncc],graph->queue,queue_in_global_numbering);CHKERRQ(ierr); 136 for (i=0;i<graph->ncc;i++) { 137 PetscInt node_num=graph->queue[graph->cptr[i]]; 138 PetscBool printcc = PETSC_FALSE; 139 ierr = PetscViewerASCIISynchronizedPrintf(viewer," %d (neighs:",i);CHKERRQ(ierr); 140 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);CHKERRQ(ierr); 141 for (j=0;j<graph->count[node_num];j++) { 142 ierr = PetscViewerASCIISynchronizedPrintf(viewer," %d",graph->neighbours_set[node_num][j]);CHKERRQ(ierr); 143 } 144 if (verbosity_level > 1) { 145 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"):");CHKERRQ(ierr); 146 if (graph->twodim || graph->count[node_num] > 1 || (graph->count[node_num] == 1 && graph->special_dof[node_num] == PCBDDCGRAPH_NEUMANN_MARK)) { 147 printcc = PETSC_TRUE; 148 } 149 if (printcc) { 150 for (j=graph->cptr[i];j<graph->cptr[i+1];j++) { 151 ierr = PetscViewerASCIISynchronizedPrintf(viewer," %d (%d)",graph->queue[j],queue_in_global_numbering[j]);CHKERRQ(ierr); 152 } 153 } 154 } else { 155 ierr = PetscViewerASCIISynchronizedPrintf(viewer,")");CHKERRQ(ierr); 156 } 157 ierr = PetscViewerASCIISynchronizedPrintf(viewer,"\n");CHKERRQ(ierr); 158 ierr = PetscViewerASCIISetTab(viewer,tabs);CHKERRQ(ierr); 159 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 160 } 161 ierr = PetscFree(queue_in_global_numbering);CHKERRQ(ierr); 162 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 163 PetscFunctionReturn(0); 164 } 165 166 #undef __FUNCT__ 167 #define __FUNCT__ "PCBDDCGraphRestoreCandidatesIS" 168 PetscErrorCode PCBDDCGraphRestoreCandidatesIS(PCBDDCGraph graph, PetscInt *n_faces, IS *FacesIS[], PetscInt *n_edges, IS *EdgesIS[], IS *VerticesIS) 169 { 170 PetscInt i; 171 PetscErrorCode ierr; 172 173 PetscFunctionBegin; 174 if (n_faces) { 175 if (FacesIS) { 176 for (i=0;i<*n_faces;i++) { 177 ierr = ISDestroy(&((*FacesIS)[i]));CHKERRQ(ierr); 178 } 179 ierr = PetscFree(*FacesIS);CHKERRQ(ierr); 180 } 181 *n_faces = 0; 182 } 183 if (n_edges) { 184 if (EdgesIS) { 185 for (i=0;i<*n_edges;i++) { 186 ierr = ISDestroy(&((*EdgesIS)[i]));CHKERRQ(ierr); 187 } 188 ierr = PetscFree(*EdgesIS);CHKERRQ(ierr); 189 } 190 *n_edges = 0; 191 } 192 if (VerticesIS) { 193 ierr = ISDestroy(VerticesIS);CHKERRQ(ierr); 194 } 195 PetscFunctionReturn(0); 196 } 197 198 #undef __FUNCT__ 199 #define __FUNCT__ "PCBDDCGraphGetCandidatesIS" 200 PetscErrorCode PCBDDCGraphGetCandidatesIS(PCBDDCGraph graph, PetscInt *n_faces, IS *FacesIS[], PetscInt *n_edges, IS *EdgesIS[], IS *VerticesIS) 201 { 202 IS *ISForFaces,*ISForEdges,ISForVertices; 203 PetscInt i,nfc,nec,nvc,*idx; 204 PetscErrorCode ierr; 205 206 PetscFunctionBegin; 207 /* loop on ccs to evalute number of faces, edges and vertices */ 208 nfc = 0; 209 nec = 0; 210 nvc = 0; 211 for (i=0;i<graph->ncc;i++) { 212 PetscInt repdof = graph->queue[graph->cptr[i]]; 213 if (graph->cptr[i+1]-graph->cptr[i] > graph->custom_minimal_size) { 214 if (graph->count[repdof] == 1 && graph->special_dof[repdof] != PCBDDCGRAPH_NEUMANN_MARK) { 215 nfc++; 216 } else { /* note that nec will be zero in 2d */ 217 nec++; 218 } 219 } else { 220 nvc += graph->cptr[i+1]-graph->cptr[i]; 221 } 222 } 223 224 /* check if we are in 2D or 3D */ 225 if (graph->twodim) { /* we are in a 2D case -> edges are shared by 2 subregions and faces don't exist */ 226 nec = nfc; 227 nfc = 0; 228 } 229 230 /* allocate IS arrays for faces, edges. Vertices need a single index set. */ 231 if (FacesIS) { 232 ierr = PetscMalloc1(nfc,&ISForFaces);CHKERRQ(ierr); 233 } 234 if (EdgesIS) { 235 ierr = PetscMalloc1(nec,&ISForEdges);CHKERRQ(ierr); 236 } 237 if (VerticesIS) { 238 ierr = PetscMalloc1(nvc,&idx);CHKERRQ(ierr); 239 } 240 241 /* loop on ccs to compute index sets for faces and edges */ 242 if (!graph->queue_sorted) { 243 PetscInt *queue_global; 244 245 ierr = PetscMalloc1(graph->cptr[graph->ncc],&queue_global);CHKERRQ(ierr); 246 ierr = ISLocalToGlobalMappingApply(graph->l2gmap,graph->cptr[graph->ncc],graph->queue,queue_global);CHKERRQ(ierr); 247 for (i=0;i<graph->ncc;i++) { 248 ierr = PetscSortIntWithArray(graph->cptr[i+1]-graph->cptr[i],&queue_global[graph->cptr[i]],&graph->queue[graph->cptr[i]]);CHKERRQ(ierr); 249 } 250 ierr = PetscFree(queue_global);CHKERRQ(ierr); 251 graph->queue_sorted = PETSC_TRUE; 252 } 253 nfc = 0; 254 nec = 0; 255 for (i=0;i<graph->ncc;i++) { 256 PetscInt repdof = graph->queue[graph->cptr[i]]; 257 if (graph->cptr[i+1]-graph->cptr[i] > graph->custom_minimal_size) { 258 if (graph->count[repdof] == 1 && graph->special_dof[repdof] != PCBDDCGRAPH_NEUMANN_MARK) { 259 if (graph->twodim) { 260 if (EdgesIS) { 261 ierr = ISCreateGeneral(PETSC_COMM_SELF,graph->cptr[i+1]-graph->cptr[i],&graph->queue[graph->cptr[i]],PETSC_USE_POINTER,&ISForEdges[nec]);CHKERRQ(ierr); 262 } 263 nec++; 264 } else { 265 if (FacesIS) { 266 ierr = ISCreateGeneral(PETSC_COMM_SELF,graph->cptr[i+1]-graph->cptr[i],&graph->queue[graph->cptr[i]],PETSC_USE_POINTER,&ISForFaces[nfc]);CHKERRQ(ierr); 267 } 268 nfc++; 269 } 270 } else { 271 if (EdgesIS) { 272 ierr = ISCreateGeneral(PETSC_COMM_SELF,graph->cptr[i+1]-graph->cptr[i],&graph->queue[graph->cptr[i]],PETSC_USE_POINTER,&ISForEdges[nec]);CHKERRQ(ierr); 273 } 274 nec++; 275 } 276 } 277 } 278 /* index set for vertices */ 279 if (VerticesIS) { 280 nvc = 0; 281 for (i=0;i<graph->ncc;i++) { 282 if (graph->cptr[i+1]-graph->cptr[i] <= graph->custom_minimal_size) { 283 PetscInt j; 284 285 for (j=graph->cptr[i];j<graph->cptr[i+1];j++) { 286 idx[nvc]=graph->queue[j]; 287 nvc++; 288 } 289 } 290 } 291 /* sort vertex set (by local ordering) */ 292 ierr = PetscSortInt(nvc,idx);CHKERRQ(ierr); 293 ierr = ISCreateGeneral(PETSC_COMM_SELF,nvc,idx,PETSC_OWN_POINTER,&ISForVertices);CHKERRQ(ierr); 294 } 295 /* get back info */ 296 if (n_faces) *n_faces = nfc; 297 if (FacesIS) *FacesIS = ISForFaces; 298 if (n_edges) *n_edges = nec; 299 if (EdgesIS) *EdgesIS = ISForEdges; 300 if (VerticesIS) *VerticesIS = ISForVertices; 301 PetscFunctionReturn(0); 302 } 303 304 #undef __FUNCT__ 305 #define __FUNCT__ "PCBDDCGraphComputeConnectedComponents" 306 PetscErrorCode PCBDDCGraphComputeConnectedComponents(PCBDDCGraph graph) 307 { 308 PetscBool adapt_interface_reduced; 309 MPI_Comm interface_comm; 310 PetscMPIInt size; 311 PetscInt i; 312 PetscBool twodim; 313 PetscErrorCode ierr; 314 315 PetscFunctionBegin; 316 /* compute connected components locally */ 317 ierr = PetscObjectGetComm((PetscObject)(graph->l2gmap),&interface_comm);CHKERRQ(ierr); 318 ierr = PCBDDCGraphComputeConnectedComponentsLocal(graph);CHKERRQ(ierr); 319 /* check consistency of connected components among neighbouring subdomains -> it adapt them in case it is needed */ 320 ierr = MPI_Comm_size(interface_comm,&size);CHKERRQ(ierr); 321 adapt_interface_reduced = PETSC_FALSE; 322 if (size > 1) { 323 PetscInt i; 324 PetscBool adapt_interface = PETSC_FALSE; 325 for (i=0;i<graph->n_subsets;i++) { 326 /* We are not sure that on a given subset of the local interface, 327 with two connected components, the latters be the same among sharing subdomains */ 328 if (graph->subset_ncc[i] > 1) { 329 adapt_interface = PETSC_TRUE; 330 break; 331 } 332 } 333 ierr = MPIU_Allreduce(&adapt_interface,&adapt_interface_reduced,1,MPIU_BOOL,MPI_LOR,interface_comm);CHKERRQ(ierr); 334 } 335 336 if (graph->n_subsets && adapt_interface_reduced) { 337 PetscBT subset_cc_adapt; 338 MPI_Request *send_requests,*recv_requests; 339 PetscInt *send_buffer,*recv_buffer; 340 PetscInt sum_requests,start_of_recv,start_of_send; 341 PetscInt *cum_recv_counts; 342 PetscInt *labels; 343 PetscInt ncc,cum_queue,mss,mns,j,k,s; 344 PetscInt **refine_buffer=NULL,*private_labels = NULL; 345 346 ierr = PetscMalloc1(graph->nvtxs,&labels);CHKERRQ(ierr); 347 ierr = PetscMemzero(labels,graph->nvtxs*sizeof(*labels));CHKERRQ(ierr); 348 for (i=0;i<graph->ncc;i++) 349 for (j=graph->cptr[i];j<graph->cptr[i+1];j++) 350 labels[graph->queue[j]] = i; 351 352 /* allocate some space */ 353 ierr = PetscMalloc1(graph->n_subsets+1,&cum_recv_counts);CHKERRQ(ierr); 354 ierr = PetscMemzero(cum_recv_counts,(graph->n_subsets+1)*sizeof(*cum_recv_counts));CHKERRQ(ierr); 355 356 /* first count how many neighbours per connected component I will receive from */ 357 cum_recv_counts[0] = 0; 358 for (i=0;i<graph->n_subsets;i++) cum_recv_counts[i+1] = cum_recv_counts[i]+graph->count[graph->subset_idxs[i][0]]; 359 ierr = PetscMalloc1(cum_recv_counts[graph->n_subsets],&recv_buffer);CHKERRQ(ierr); 360 ierr = PetscMalloc2(cum_recv_counts[graph->n_subsets],&send_requests,cum_recv_counts[graph->n_subsets],&recv_requests);CHKERRQ(ierr); 361 for (i=0;i<cum_recv_counts[graph->n_subsets];i++) { 362 send_requests[i] = MPI_REQUEST_NULL; 363 recv_requests[i] = MPI_REQUEST_NULL; 364 } 365 366 /* exchange with my neighbours the number of my connected components on the subset of interface */ 367 sum_requests = 0; 368 for (i=0;i<graph->n_subsets;i++) { 369 PetscMPIInt neigh,tag; 370 PetscInt count,*neighs; 371 372 count = graph->count[graph->subset_idxs[i][0]]; 373 neighs = graph->neighbours_set[graph->subset_idxs[i][0]]; 374 ierr = PetscMPIIntCast(2*graph->subset_ref_node[i],&tag);CHKERRQ(ierr); 375 for (k=0;k<count;k++) { 376 ierr = PetscMPIIntCast(neighs[k],&neigh);CHKERRQ(ierr); 377 ierr = MPI_Isend(&graph->subset_ncc[i],1,MPIU_INT,neigh,tag,interface_comm,&send_requests[sum_requests]);CHKERRQ(ierr); 378 ierr = MPI_Irecv(&recv_buffer[sum_requests],1,MPIU_INT,neigh,tag,interface_comm,&recv_requests[sum_requests]);CHKERRQ(ierr); 379 sum_requests++; 380 } 381 } 382 ierr = MPI_Waitall(sum_requests,recv_requests,MPI_STATUSES_IGNORE);CHKERRQ(ierr); 383 ierr = MPI_Waitall(sum_requests,send_requests,MPI_STATUSES_IGNORE);CHKERRQ(ierr); 384 385 /* determine the subsets I have to adapt (those having more than 1 cc) */ 386 ierr = PetscBTCreate(graph->n_subsets,&subset_cc_adapt);CHKERRQ(ierr); 387 ierr = PetscBTMemzero(graph->n_subsets,subset_cc_adapt);CHKERRQ(ierr); 388 for (i=0;i<graph->n_subsets;i++) { 389 if (graph->subset_ncc[i] > 1) { 390 ierr = PetscBTSet(subset_cc_adapt,i);CHKERRQ(ierr); 391 continue; 392 } 393 for (j=cum_recv_counts[i];j<cum_recv_counts[i+1];j++){ 394 if (recv_buffer[j] > 1) { 395 ierr = PetscBTSet(subset_cc_adapt,i); 396 break; 397 } 398 } 399 } 400 ierr = PetscFree(recv_buffer);CHKERRQ(ierr); 401 402 /* determine send/recv buffers sizes */ 403 j = 0; 404 mss = 0; 405 for (i=0;i<graph->n_subsets;i++) { 406 if (PetscBTLookup(subset_cc_adapt,i)) { 407 j += graph->subset_size[i]; 408 mss = PetscMax(graph->subset_size[i],mss); 409 } 410 } 411 k = 0; 412 mns = 0; 413 for (i=0;i<graph->n_subsets;i++) { 414 if (PetscBTLookup(subset_cc_adapt,i)) { 415 k += (cum_recv_counts[i+1]-cum_recv_counts[i])*graph->subset_size[i]; 416 mns = PetscMax(cum_recv_counts[i+1]-cum_recv_counts[i],mns); 417 } 418 } 419 ierr = PetscMalloc2(j,&send_buffer,k,&recv_buffer);CHKERRQ(ierr); 420 421 /* fill send buffer (order matters: subset_idxs ordered by global ordering) */ 422 j = 0; 423 for (i=0;i<graph->n_subsets;i++) 424 if (PetscBTLookup(subset_cc_adapt,i)) 425 for (k=0;k<graph->subset_size[i];k++) 426 send_buffer[j++] = labels[graph->subset_idxs[i][k]]; 427 428 /* now exchange the data */ 429 start_of_recv = 0; 430 start_of_send = 0; 431 sum_requests = 0; 432 for (i=0;i<graph->n_subsets;i++) { 433 if (PetscBTLookup(subset_cc_adapt,i)) { 434 PetscMPIInt neigh,tag; 435 PetscInt size_of_send = graph->subset_size[i]; 436 437 j = graph->subset_idxs[i][0]; 438 ierr = PetscMPIIntCast(2*graph->subset_ref_node[i]+1,&tag);CHKERRQ(ierr); 439 for (k=0;k<graph->count[j];k++) { 440 ierr = PetscMPIIntCast(graph->neighbours_set[j][k],&neigh);CHKERRQ(ierr); 441 ierr = MPI_Isend(&send_buffer[start_of_send],size_of_send,MPIU_INT,neigh,tag,interface_comm,&send_requests[sum_requests]);CHKERRQ(ierr); 442 ierr = MPI_Irecv(&recv_buffer[start_of_recv],size_of_send,MPIU_INT,neigh,tag,interface_comm,&recv_requests[sum_requests]);CHKERRQ(ierr); 443 start_of_recv += size_of_send; 444 sum_requests++; 445 } 446 start_of_send += size_of_send; 447 } 448 } 449 ierr = MPI_Waitall(sum_requests,recv_requests,MPI_STATUSES_IGNORE);CHKERRQ(ierr); 450 451 /* refine connected components */ 452 start_of_recv = 0; 453 /* allocate some temporary space */ 454 if (mss) { 455 ierr = PetscMalloc1(mss,&refine_buffer);CHKERRQ(ierr); 456 ierr = PetscMalloc2(mss*(mns+1),&refine_buffer[0],mss,&private_labels);CHKERRQ(ierr); 457 } 458 ncc = 0; 459 cum_queue = 0; 460 graph->cptr[0] = 0; 461 for (i=0;i<graph->n_subsets;i++) { 462 if (PetscBTLookup(subset_cc_adapt,i)) { 463 PetscInt subset_counter = 0; 464 PetscInt sharingprocs = cum_recv_counts[i+1]-cum_recv_counts[i]+1; /* count myself */ 465 PetscInt buffer_size = graph->subset_size[i]; 466 467 /* compute pointers */ 468 for (j=1;j<buffer_size;j++) refine_buffer[j] = refine_buffer[j-1] + sharingprocs; 469 /* analyze contributions from subdomains that share the i-th subset 470 The stricture of refine_buffer is suitable to find intersections of ccs among sharingprocs. 471 supposing the current subset is shared by 3 processes and has dimension 5 with global dofs 0,1,2,3,4 (local 0,4,3,1,2) 472 sharing procs connected components: 473 neigh 0: [0 1 4], [2 3], labels [4,7] (2 connected components) 474 neigh 1: [0 1], [2 3 4], labels [3 2] (2 connected components) 475 neigh 2: [0 4], [1], [2 3], labels [1 5 6] (3 connected components) 476 refine_buffer will be filled as: 477 [ 4, 3, 1; 478 4, 2, 1; 479 7, 2, 6; 480 4, 3, 5; 481 7, 2, 6; ]; 482 The connected components in local ordering are [0], [1], [2 3], [4] */ 483 /* fill temp_buffer */ 484 for (k=0;k<buffer_size;k++) refine_buffer[k][0] = labels[graph->subset_idxs[i][k]]; 485 for (j=0;j<sharingprocs-1;j++) { 486 for (k=0;k<buffer_size;k++) refine_buffer[k][j+1] = recv_buffer[start_of_recv+k]; 487 start_of_recv += buffer_size; 488 } 489 ierr = PetscMemzero(private_labels,buffer_size*sizeof(PetscInt));CHKERRQ(ierr); 490 for (j=0;j<buffer_size;j++) { 491 if (!private_labels[j]) { /* found a new cc */ 492 PetscBool same_set; 493 494 graph->cptr[ncc] = cum_queue; 495 ncc++; 496 subset_counter++; 497 private_labels[j] = subset_counter; 498 graph->queue[cum_queue++] = graph->subset_idxs[i][j]; 499 for (k=j+1;k<buffer_size;k++) { /* check for other nodes in new cc */ 500 same_set = PETSC_TRUE; 501 for (s=0;s<sharingprocs;s++) { 502 if (refine_buffer[j][s] != refine_buffer[k][s]) { 503 same_set = PETSC_FALSE; 504 break; 505 } 506 } 507 if (same_set) { 508 private_labels[k] = subset_counter; 509 graph->queue[cum_queue++] = graph->subset_idxs[i][k]; 510 } 511 } 512 } 513 } 514 graph->cptr[ncc] = cum_queue; 515 graph->subset_ncc[i] = subset_counter; 516 graph->queue_sorted = PETSC_FALSE; 517 } else { /* this subset does not need to be adapted */ 518 ierr = PetscMemcpy(graph->queue+cum_queue,graph->subset_idxs[i],graph->subset_size[i]*sizeof(PetscInt));CHKERRQ(ierr); 519 ncc++; 520 cum_queue += graph->subset_size[i]; 521 graph->cptr[ncc] = cum_queue; 522 } 523 } 524 graph->cptr[ncc] = cum_queue; 525 graph->ncc = ncc; 526 if (mss) { 527 ierr = PetscFree2(refine_buffer[0],private_labels);CHKERRQ(ierr); 528 ierr = PetscFree(refine_buffer);CHKERRQ(ierr); 529 } 530 ierr = PetscFree(labels);CHKERRQ(ierr); 531 ierr = MPI_Waitall(sum_requests,send_requests,MPI_STATUSES_IGNORE);CHKERRQ(ierr); 532 ierr = PetscFree2(send_requests,recv_requests);CHKERRQ(ierr); 533 ierr = PetscFree2(send_buffer,recv_buffer);CHKERRQ(ierr); 534 ierr = PetscFree(cum_recv_counts);CHKERRQ(ierr); 535 ierr = PetscBTDestroy(&subset_cc_adapt);CHKERRQ(ierr); 536 } 537 538 /* Determine if we are in 2D or 3D */ 539 twodim = PETSC_TRUE; 540 for (i=0;i<graph->ncc;i++) { 541 PetscInt repdof = graph->queue[graph->cptr[i]]; 542 if (graph->cptr[i+1]-graph->cptr[i] > graph->custom_minimal_size) { 543 if (graph->count[repdof] > 1 || graph->special_dof[repdof] == PCBDDCGRAPH_NEUMANN_MARK) { 544 twodim = PETSC_FALSE; 545 break; 546 } 547 } 548 } 549 ierr = MPIU_Allreduce(&twodim,&graph->twodim,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)graph->l2gmap));CHKERRQ(ierr); 550 PetscFunctionReturn(0); 551 } 552 553 554 #undef __FUNCT__ 555 #define __FUNCT__ "PCBDDCGraphComputeCC_Private" 556 PETSC_STATIC_INLINE PetscErrorCode PCBDDCGraphComputeCC_Private(PCBDDCGraph graph,PetscInt pid,PetscInt* queue_tip,PetscInt n_prev,PetscInt* n_added) 557 { 558 PetscInt i,j,n; 559 PetscInt *xadj = graph->xadj,*adjncy = graph->adjncy; 560 PetscBT touched = graph->touched; 561 PetscBool havecsr = (PetscBool)(xadj && adjncy); 562 PetscBool havesubs = (PetscBool)(!!graph->n_local_subs); 563 PetscErrorCode ierr; 564 565 PetscFunctionBegin; 566 n = 0; 567 if (havecsr && !havesubs) { 568 for (i=-n_prev;i<0;i++) { 569 PetscInt start_dof = queue_tip[i]; 570 for (j=xadj[start_dof];j<xadj[start_dof+1];j++) { 571 PetscInt dof = adjncy[j]; 572 if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid) { 573 ierr = PetscBTSet(touched,dof);CHKERRQ(ierr); 574 queue_tip[n] = dof; 575 n++; 576 } 577 } 578 } 579 } else if (havecsr && havesubs) { 580 PetscInt sid = graph->local_subs[queue_tip[-n_prev]]; 581 for (i=-n_prev;i<0;i++) { 582 PetscInt start_dof = queue_tip[i]; 583 for (j=xadj[start_dof];j<xadj[start_dof+1];j++) { 584 PetscInt dof = adjncy[j]; 585 if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid && graph->local_subs[dof] == sid) { 586 ierr = PetscBTSet(touched,dof);CHKERRQ(ierr); 587 queue_tip[n] = dof; 588 n++; 589 } 590 } 591 } 592 } else { /* sub info only */ 593 PetscInt sid = graph->local_subs[queue_tip[-n_prev]]; 594 for (j=0;j<graph->subset_size[pid-1];j++) { /* pid \in [1,graph->n_subsets] */ 595 PetscInt dof = graph->subset_idxs[pid-1][j]; 596 if (!PetscBTLookup(touched,dof) && graph->subset[dof] == pid && graph->local_subs[dof] == sid) { 597 ierr = PetscBTSet(touched,dof);CHKERRQ(ierr); 598 queue_tip[n] = dof; 599 n++; 600 } 601 } 602 } 603 *n_added = n; 604 PetscFunctionReturn(0); 605 } 606 607 #undef __FUNCT__ 608 #define __FUNCT__ "PCBDDCGraphComputeConnectedComponentsLocal" 609 PetscErrorCode PCBDDCGraphComputeConnectedComponentsLocal(PCBDDCGraph graph) 610 { 611 PetscInt ncc,cum_queue,n; 612 PetscMPIInt commsize; 613 PetscErrorCode ierr; 614 615 PetscFunctionBegin; 616 if (!graph->setupcalled) SETERRQ(PetscObjectComm((PetscObject)graph->l2gmap),PETSC_ERR_ORDER,"PCBDDCGraphSetUp should be called first"); 617 /* quiet return if there isn't any local info */ 618 if ((!graph->xadj || !graph->adjncy) && !graph->n_local_subs) { 619 PetscFunctionReturn(0); 620 } 621 622 /* reset any previous search of connected components */ 623 ierr = PetscBTMemzero(graph->nvtxs,graph->touched);CHKERRQ(ierr); 624 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)graph->l2gmap),&commsize);CHKERRQ(ierr); 625 if (commsize > graph->commsizelimit) { 626 PetscInt i; 627 for (i=0;i<graph->nvtxs;i++) { 628 if (graph->special_dof[i] == PCBDDCGRAPH_DIRICHLET_MARK || !graph->count[i]) { 629 ierr = PetscBTSet(graph->touched,i);CHKERRQ(ierr); 630 } 631 } 632 } 633 634 /* begin search for connected components */ 635 cum_queue = 0; 636 ncc = 0; 637 for (n=0;n<graph->n_subsets;n++) { 638 PetscInt pid = n+1; /* partition labeled by 0 is discarded */ 639 PetscInt found = 0,prev = 0,first = 0,ncc_pid = 0; 640 while (found != graph->subset_size[n]) { 641 PetscInt added = 0; 642 if (!prev) { /* search for new starting dof */ 643 while (PetscBTLookup(graph->touched,graph->subset_idxs[n][first])) first++; 644 ierr = PetscBTSet(graph->touched,graph->subset_idxs[n][first]);CHKERRQ(ierr); 645 graph->queue[cum_queue] = graph->subset_idxs[n][first]; 646 graph->cptr[ncc] = cum_queue; 647 prev = 1; 648 cum_queue++; 649 found++; 650 ncc_pid++; 651 ncc++; 652 } 653 ierr = PCBDDCGraphComputeCC_Private(graph,pid,graph->queue + cum_queue,prev,&added);CHKERRQ(ierr); 654 if (!added) { 655 graph->subset_ncc[n] = ncc_pid; 656 graph->cptr[ncc] = cum_queue; 657 } 658 prev = added; 659 found += added; 660 cum_queue += added; 661 if (added && found == graph->subset_size[n]) { 662 graph->subset_ncc[n] = ncc_pid; 663 graph->cptr[ncc] = cum_queue; 664 } 665 } 666 } 667 graph->ncc = ncc; 668 graph->queue_sorted = PETSC_FALSE; 669 PetscFunctionReturn(0); 670 } 671 672 #undef __FUNCT__ 673 #define __FUNCT__ "PCBDDCGraphSetUp" 674 PetscErrorCode PCBDDCGraphSetUp(PCBDDCGraph graph, PetscInt custom_minimal_size, IS neumann_is, IS dirichlet_is, PetscInt n_ISForDofs, IS ISForDofs[], IS custom_primal_vertices) 675 { 676 VecScatter scatter_ctx; 677 Vec local_vec,local_vec2,global_vec; 678 IS to,from,subset,subset_n; 679 MPI_Comm comm; 680 PetscScalar *array,*array2; 681 const PetscInt *is_indices; 682 PetscInt n_neigh,*neigh,*n_shared,**shared,*queue_global; 683 PetscInt i,j,k,s,total_counts,nodes_touched,is_size; 684 PetscMPIInt commsize; 685 PetscBool same_set,mirrors_found; 686 PetscErrorCode ierr; 687 688 PetscFunctionBegin; 689 graph->has_dirichlet = PETSC_FALSE; 690 if (dirichlet_is) { 691 PetscCheckSameComm(graph->l2gmap,1,dirichlet_is,4); 692 graph->has_dirichlet = PETSC_TRUE; 693 } 694 ierr = PetscObjectGetComm((PetscObject)(graph->l2gmap),&comm);CHKERRQ(ierr); 695 ierr = MPI_Comm_size(comm,&commsize);CHKERRQ(ierr); 696 697 /* custom_minimal_size */ 698 graph->custom_minimal_size = custom_minimal_size; 699 /* get info l2gmap and allocate work vectors */ 700 ierr = ISLocalToGlobalMappingGetInfo(graph->l2gmap,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr); 701 /* check if we have any local periodic nodes (periodic BCs) */ 702 mirrors_found = PETSC_FALSE; 703 if (graph->nvtxs && n_neigh) { 704 for (i=0; i<n_shared[0]; i++) graph->count[shared[0][i]] += 1; 705 for (i=0; i<n_shared[0]; i++) { 706 if (graph->count[shared[0][i]] > 1) { 707 mirrors_found = PETSC_TRUE; 708 break; 709 } 710 } 711 } 712 /* create some workspace objects */ 713 local_vec = NULL; 714 local_vec2 = NULL; 715 global_vec = NULL; 716 to = NULL; 717 from = NULL; 718 scatter_ctx = NULL; 719 if (n_ISForDofs || dirichlet_is || neumann_is || custom_primal_vertices) { 720 ierr = VecCreate(PETSC_COMM_SELF,&local_vec);CHKERRQ(ierr); 721 ierr = VecSetSizes(local_vec,PETSC_DECIDE,graph->nvtxs);CHKERRQ(ierr); 722 ierr = VecSetType(local_vec,VECSTANDARD);CHKERRQ(ierr); 723 ierr = VecDuplicate(local_vec,&local_vec2);CHKERRQ(ierr); 724 ierr = VecCreate(comm,&global_vec);CHKERRQ(ierr); 725 ierr = VecSetSizes(global_vec,PETSC_DECIDE,graph->nvtxs_global);CHKERRQ(ierr); 726 ierr = VecSetType(global_vec,VECSTANDARD);CHKERRQ(ierr); 727 ierr = ISCreateStride(PETSC_COMM_SELF,graph->nvtxs,0,1,&to);CHKERRQ(ierr); 728 ierr = ISLocalToGlobalMappingApplyIS(graph->l2gmap,to,&from);CHKERRQ(ierr); 729 ierr = VecScatterCreate(global_vec,from,local_vec,to,&scatter_ctx);CHKERRQ(ierr); 730 } else if (mirrors_found) { 731 ierr = ISCreateStride(PETSC_COMM_SELF,graph->nvtxs,0,1,&to);CHKERRQ(ierr); 732 ierr = ISLocalToGlobalMappingApplyIS(graph->l2gmap,to,&from);CHKERRQ(ierr); 733 } 734 /* compute local mirrors (if any) */ 735 if (mirrors_found) { 736 PetscInt *local_indices,*global_indices; 737 /* get arrays of local and global indices */ 738 ierr = PetscMalloc1(graph->nvtxs,&local_indices);CHKERRQ(ierr); 739 ierr = ISGetIndices(to,(const PetscInt**)&is_indices);CHKERRQ(ierr); 740 ierr = PetscMemcpy(local_indices,is_indices,graph->nvtxs*sizeof(PetscInt));CHKERRQ(ierr); 741 ierr = ISRestoreIndices(to,(const PetscInt**)&is_indices);CHKERRQ(ierr); 742 ierr = PetscMalloc1(graph->nvtxs,&global_indices);CHKERRQ(ierr); 743 ierr = ISGetIndices(from,(const PetscInt**)&is_indices);CHKERRQ(ierr); 744 ierr = PetscMemcpy(global_indices,is_indices,graph->nvtxs*sizeof(PetscInt));CHKERRQ(ierr); 745 ierr = ISRestoreIndices(from,(const PetscInt**)&is_indices);CHKERRQ(ierr); 746 /* allocate space for mirrors */ 747 ierr = PetscMalloc2(graph->nvtxs,&graph->mirrors,graph->nvtxs,&graph->mirrors_set);CHKERRQ(ierr); 748 ierr = PetscMemzero(graph->mirrors,graph->nvtxs*sizeof(PetscInt));CHKERRQ(ierr); 749 graph->mirrors_set[0] = 0; 750 751 k=0; 752 for (i=0;i<n_shared[0];i++) { 753 j=shared[0][i]; 754 if (graph->count[j] > 1) { 755 graph->mirrors[j]++; 756 k++; 757 } 758 } 759 /* allocate space for set of mirrors */ 760 ierr = PetscMalloc1(k,&graph->mirrors_set[0]);CHKERRQ(ierr); 761 for (i=1;i<graph->nvtxs;i++) 762 graph->mirrors_set[i]=graph->mirrors_set[i-1]+graph->mirrors[i-1]; 763 764 /* fill arrays */ 765 ierr = PetscMemzero(graph->mirrors,graph->nvtxs*sizeof(PetscInt));CHKERRQ(ierr); 766 for (j=0;j<n_shared[0];j++) { 767 i=shared[0][j]; 768 if (graph->count[i] > 1) 769 graph->mirrors_set[i][graph->mirrors[i]++]=global_indices[i]; 770 } 771 ierr = PetscSortIntWithArray(graph->nvtxs,global_indices,local_indices);CHKERRQ(ierr); 772 for (i=0;i<graph->nvtxs;i++) { 773 if (graph->mirrors[i] > 0) { 774 ierr = PetscFindInt(graph->mirrors_set[i][0],graph->nvtxs,global_indices,&k);CHKERRQ(ierr); 775 j = global_indices[k]; 776 while ( k > 0 && global_indices[k-1] == j) k--; 777 for (j=0;j<graph->mirrors[i];j++) { 778 graph->mirrors_set[i][j]=local_indices[k+j]; 779 } 780 ierr = PetscSortInt(graph->mirrors[i],graph->mirrors_set[i]);CHKERRQ(ierr); 781 } 782 } 783 ierr = PetscFree(local_indices);CHKERRQ(ierr); 784 ierr = PetscFree(global_indices);CHKERRQ(ierr); 785 } 786 ierr = PetscMemzero(graph->count,graph->nvtxs*sizeof(*graph->count));CHKERRQ(ierr); 787 ierr = ISDestroy(&to);CHKERRQ(ierr); 788 ierr = ISDestroy(&from);CHKERRQ(ierr); 789 790 /* Count total number of neigh per node */ 791 k = 0; 792 for (i=1;i<n_neigh;i++) { 793 k += n_shared[i]; 794 for (j=0;j<n_shared[i];j++) { 795 graph->count[shared[i][j]] += 1; 796 } 797 } 798 /* Allocate space for storing the set of neighbours for each node */ 799 if (graph->nvtxs) { 800 ierr = PetscMalloc1(k,&graph->neighbours_set[0]);CHKERRQ(ierr); 801 } 802 for (i=1;i<graph->nvtxs;i++) { /* dont count myself */ 803 graph->neighbours_set[i]=graph->neighbours_set[i-1]+graph->count[i-1]; 804 } 805 /* Get information for sharing subdomains */ 806 ierr = PetscMemzero(graph->count,graph->nvtxs*sizeof(*graph->count));CHKERRQ(ierr); 807 for (i=1;i<n_neigh;i++) { /* dont count myself */ 808 s = n_shared[i]; 809 for (j=0;j<s;j++) { 810 k = shared[i][j]; 811 graph->neighbours_set[k][graph->count[k]] = neigh[i]; 812 graph->count[k] += 1; 813 } 814 } 815 /* sort set of sharing subdomains */ 816 for (i=0;i<graph->nvtxs;i++) { 817 ierr = PetscSortRemoveDupsInt(&graph->count[i],graph->neighbours_set[i]);CHKERRQ(ierr); 818 } 819 /* free memory allocated by ISLocalToGlobalMappingGetInfo */ 820 ierr = ISLocalToGlobalMappingRestoreInfo(graph->l2gmap,&n_neigh,&neigh,&n_shared,&shared);CHKERRQ(ierr); 821 822 /* 823 Get info for dofs splitting 824 User can specify just a subset; an additional field is considered as a complementary field 825 */ 826 for (i=0;i<graph->nvtxs;i++) graph->which_dof[i] = n_ISForDofs; /* by default a dof belongs to the complement set */ 827 if (n_ISForDofs) { 828 ierr = VecSet(local_vec,-1.0);CHKERRQ(ierr); 829 } 830 for (i=0;i<n_ISForDofs;i++) { 831 ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr); 832 ierr = ISGetLocalSize(ISForDofs[i],&is_size);CHKERRQ(ierr); 833 ierr = ISGetIndices(ISForDofs[i],(const PetscInt**)&is_indices);CHKERRQ(ierr); 834 for (j=0;j<is_size;j++) { 835 if (is_indices[j] > -1 && is_indices[j] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */ 836 graph->which_dof[is_indices[j]] = i; 837 array[is_indices[j]] = 1.*i; 838 } 839 } 840 ierr = ISRestoreIndices(ISForDofs[i],(const PetscInt**)&is_indices);CHKERRQ(ierr); 841 ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr); 842 } 843 /* Check consistency among neighbours */ 844 if (n_ISForDofs) { 845 ierr = VecScatterBegin(scatter_ctx,local_vec,global_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 846 ierr = VecScatterEnd(scatter_ctx,local_vec,global_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 847 ierr = VecScatterBegin(scatter_ctx,global_vec,local_vec2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 848 ierr = VecScatterEnd(scatter_ctx,global_vec,local_vec2,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 849 ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr); 850 ierr = VecGetArray(local_vec2,&array2);CHKERRQ(ierr); 851 for (i=0;i<graph->nvtxs;i++){ 852 PetscInt field1,field2; 853 854 field1 = (PetscInt)PetscRealPart(array[i]); 855 field2 = (PetscInt)PetscRealPart(array2[i]); 856 if (field1 != field2) SETERRQ3(comm,PETSC_ERR_USER,"Local node %D have been assigned two different field ids %D and %D at the same time\n",i,field1,field2); 857 } 858 ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr); 859 ierr = VecRestoreArray(local_vec2,&array2);CHKERRQ(ierr); 860 } 861 /* Take into account Neumann nodes */ 862 if (neumann_is) { 863 ierr = VecSet(local_vec,0.0);CHKERRQ(ierr); 864 ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr); 865 ierr = ISGetLocalSize(neumann_is,&is_size);CHKERRQ(ierr); 866 ierr = ISGetIndices(neumann_is,(const PetscInt**)&is_indices);CHKERRQ(ierr); 867 for (i=0;i<is_size;i++) { 868 if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */ 869 array[is_indices[i]] = 1.0; 870 } 871 } 872 ierr = ISRestoreIndices(neumann_is,(const PetscInt**)&is_indices);CHKERRQ(ierr); 873 ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr); 874 /* Neumann nodes: impose consistency among neighbours */ 875 ierr = VecSet(global_vec,0.0);CHKERRQ(ierr); 876 ierr = VecScatterBegin(scatter_ctx,local_vec,global_vec,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 877 ierr = VecScatterEnd(scatter_ctx,local_vec,global_vec,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 878 ierr = VecScatterBegin(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 879 ierr = VecScatterEnd(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 880 ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr); 881 for (i=0;i<graph->nvtxs;i++) { 882 if (PetscRealPart(array[i]) > 0.1) { 883 graph->special_dof[i] = PCBDDCGRAPH_NEUMANN_MARK; 884 } 885 } 886 ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr); 887 } 888 /* Take into account Dirichlet nodes (they overwrite any Neumann boundary node previously set) */ 889 if (dirichlet_is) { 890 ierr = VecSet(local_vec,0.0);CHKERRQ(ierr); 891 ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr); 892 ierr = ISGetLocalSize(dirichlet_is,&is_size);CHKERRQ(ierr); 893 ierr = ISGetIndices(dirichlet_is,(const PetscInt**)&is_indices);CHKERRQ(ierr); 894 for (i=0;i<is_size;i++){ 895 if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */ 896 array[is_indices[i]] = 1.0; 897 } 898 } 899 ierr = ISRestoreIndices(dirichlet_is,(const PetscInt**)&is_indices);CHKERRQ(ierr); 900 ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr); 901 /* Dirichlet nodes: impose consistency among neighbours */ 902 ierr = VecSet(global_vec,0.0);CHKERRQ(ierr); 903 ierr = VecScatterBegin(scatter_ctx,local_vec,global_vec,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 904 ierr = VecScatterEnd(scatter_ctx,local_vec,global_vec,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 905 ierr = VecScatterBegin(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 906 ierr = VecScatterEnd(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 907 ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr); 908 for (i=0;i<graph->nvtxs;i++) { 909 if (PetscRealPart(array[i]) > 0.1) { 910 if (commsize > graph->commsizelimit) { /* dirichlet nodes treated as internal */ 911 ierr = PetscBTSet(graph->touched,i);CHKERRQ(ierr); 912 graph->subset[i] = 0; 913 } 914 graph->special_dof[i] = PCBDDCGRAPH_DIRICHLET_MARK; 915 } 916 } 917 ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr); 918 } 919 /* mark local periodic nodes (if any) and adapt CSR graph (if any) */ 920 if (graph->mirrors) { 921 for (i=0;i<graph->nvtxs;i++) 922 if (graph->mirrors[i]) 923 graph->special_dof[i] = PCBDDCGRAPH_LOCAL_PERIODIC_MARK; 924 925 if (graph->xadj && graph->adjncy) { 926 PetscInt *new_xadj,*new_adjncy; 927 /* sort CSR graph */ 928 for (i=0;i<graph->nvtxs;i++) 929 ierr = PetscSortInt(graph->xadj[i+1]-graph->xadj[i],&graph->adjncy[graph->xadj[i]]);CHKERRQ(ierr); 930 931 /* adapt local CSR graph in case of local periodicity */ 932 k = 0; 933 for (i=0;i<graph->nvtxs;i++) 934 for (j=graph->xadj[i];j<graph->xadj[i+1];j++) 935 k += graph->mirrors[graph->adjncy[j]]; 936 937 ierr = PetscMalloc1(graph->nvtxs+1,&new_xadj);CHKERRQ(ierr); 938 ierr = PetscMalloc1(k+graph->xadj[graph->nvtxs],&new_adjncy);CHKERRQ(ierr); 939 new_xadj[0] = 0; 940 for (i=0;i<graph->nvtxs;i++) { 941 k = graph->xadj[i+1]-graph->xadj[i]; 942 ierr = PetscMemcpy(&new_adjncy[new_xadj[i]],&graph->adjncy[graph->xadj[i]],k*sizeof(PetscInt));CHKERRQ(ierr); 943 new_xadj[i+1] = new_xadj[i]+k; 944 for (j=graph->xadj[i];j<graph->xadj[i+1];j++) { 945 k = graph->mirrors[graph->adjncy[j]]; 946 ierr = PetscMemcpy(&new_adjncy[new_xadj[i+1]],graph->mirrors_set[graph->adjncy[j]],k*sizeof(PetscInt));CHKERRQ(ierr); 947 new_xadj[i+1] += k; 948 } 949 k = new_xadj[i+1]-new_xadj[i]; 950 ierr = PetscSortRemoveDupsInt(&k,&new_adjncy[new_xadj[i]]);CHKERRQ(ierr); 951 new_xadj[i+1] = new_xadj[i]+k; 952 } 953 /* set new CSR into graph */ 954 ierr = PetscFree(graph->xadj);CHKERRQ(ierr); 955 ierr = PetscFree(graph->adjncy);CHKERRQ(ierr); 956 graph->xadj = new_xadj; 957 graph->adjncy = new_adjncy; 958 } 959 } 960 961 /* mark special nodes (if any) -> each will become a single node equivalence class */ 962 if (custom_primal_vertices) { 963 ierr = VecSet(local_vec,0.0);CHKERRQ(ierr); 964 ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr); 965 ierr = ISGetLocalSize(custom_primal_vertices,&is_size);CHKERRQ(ierr); 966 ierr = ISGetIndices(custom_primal_vertices,(const PetscInt**)&is_indices);CHKERRQ(ierr); 967 for (i=0;i<is_size;i++){ 968 if (is_indices[i] > -1 && is_indices[i] < graph->nvtxs) { /* out of bounds indices (if any) are skipped */ 969 array[is_indices[i]] = 1.0; 970 } 971 } 972 ierr = ISRestoreIndices(custom_primal_vertices,(const PetscInt**)&is_indices);CHKERRQ(ierr); 973 ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr); 974 /* special nodes: impose consistency among neighbours */ 975 ierr = VecSet(global_vec,0.0);CHKERRQ(ierr); 976 ierr = VecScatterBegin(scatter_ctx,local_vec,global_vec,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 977 ierr = VecScatterEnd(scatter_ctx,local_vec,global_vec,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 978 ierr = VecScatterBegin(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 979 ierr = VecScatterEnd(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 980 ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr); 981 j = 0; 982 for (i=0;i<graph->nvtxs;i++) { 983 if (PetscRealPart(array[i]) > 0.1 && graph->special_dof[i] != PCBDDCGRAPH_DIRICHLET_MARK) { 984 graph->special_dof[i] = PCBDDCGRAPH_SPECIAL_MARK-j; 985 j++; 986 } 987 } 988 ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr); 989 } 990 991 /* mark interior nodes (if commsize > graph->commsizelimit) as touched and belonging to partition number 0 */ 992 if (commsize > graph->commsizelimit) { 993 for (i=0;i<graph->nvtxs;i++) { 994 if (!graph->count[i]) { 995 ierr = PetscBTSet(graph->touched,i);CHKERRQ(ierr); 996 graph->subset[i] = 0; 997 } 998 } 999 } 1000 1001 /* init graph structure and compute default subsets */ 1002 nodes_touched = 0; 1003 for (i=0;i<graph->nvtxs;i++) { 1004 if (PetscBTLookup(graph->touched,i)) { 1005 nodes_touched++; 1006 } 1007 } 1008 i = 0; 1009 graph->ncc = 0; 1010 total_counts = 0; 1011 1012 /* allocated space for queues */ 1013 if (commsize == graph->commsizelimit) { 1014 ierr = PetscMalloc2(graph->nvtxs+1,&graph->cptr,graph->nvtxs,&graph->queue);CHKERRQ(ierr); 1015 } else { 1016 PetscInt nused = graph->nvtxs - nodes_touched; 1017 ierr = PetscMalloc2(nused+1,&graph->cptr,nused,&graph->queue);CHKERRQ(ierr); 1018 } 1019 1020 while (nodes_touched<graph->nvtxs) { 1021 /* find first untouched node in local ordering */ 1022 while (PetscBTLookup(graph->touched,i)) i++; 1023 ierr = PetscBTSet(graph->touched,i);CHKERRQ(ierr); 1024 graph->subset[i] = graph->ncc+1; 1025 graph->cptr[graph->ncc] = total_counts; 1026 graph->queue[total_counts] = i; 1027 total_counts++; 1028 nodes_touched++; 1029 /* now find all other nodes having the same set of sharing subdomains */ 1030 for (j=i+1;j<graph->nvtxs;j++) { 1031 /* check for same number of sharing subdomains, dof number and same special mark */ 1032 if (!PetscBTLookup(graph->touched,j) && graph->count[i] == graph->count[j] && graph->which_dof[i] == graph->which_dof[j] && graph->special_dof[i] == graph->special_dof[j]) { 1033 /* check for same set of sharing subdomains */ 1034 same_set = PETSC_TRUE; 1035 for (k=0;k<graph->count[j];k++){ 1036 if (graph->neighbours_set[i][k] != graph->neighbours_set[j][k]) { 1037 same_set = PETSC_FALSE; 1038 } 1039 } 1040 /* I found a friend of mine */ 1041 if (same_set) { 1042 ierr = PetscBTSet(graph->touched,j);CHKERRQ(ierr); 1043 graph->subset[j] = graph->ncc+1; 1044 nodes_touched++; 1045 graph->queue[total_counts] = j; 1046 total_counts++; 1047 } 1048 } 1049 } 1050 graph->ncc++; 1051 } 1052 /* set default number of subsets (at this point no info on csr and/or local_subs has been taken into account, so n_subsets = ncc */ 1053 graph->n_subsets = graph->ncc; 1054 ierr = PetscMalloc1(graph->n_subsets,&graph->subset_ncc);CHKERRQ(ierr); 1055 for (i=0;i<graph->n_subsets;i++) { 1056 graph->subset_ncc[i] = 1; 1057 } 1058 /* final pointer */ 1059 graph->cptr[graph->ncc] = total_counts; 1060 1061 /* For consistency reasons (among neighbours), I need to sort (by global ordering) each connected component */ 1062 /* Get a reference node (min index in global ordering) for each subset for tagging messages */ 1063 ierr = PetscMalloc1(graph->ncc,&graph->subset_ref_node);CHKERRQ(ierr); 1064 ierr = PetscMalloc1(graph->cptr[graph->ncc],&queue_global);CHKERRQ(ierr); 1065 ierr = ISLocalToGlobalMappingApply(graph->l2gmap,graph->cptr[graph->ncc],graph->queue,queue_global);CHKERRQ(ierr); 1066 for (j=0;j<graph->ncc;j++) { 1067 ierr = PetscSortIntWithArray(graph->cptr[j+1]-graph->cptr[j],&queue_global[graph->cptr[j]],&graph->queue[graph->cptr[j]]);CHKERRQ(ierr); 1068 graph->subset_ref_node[j] = graph->queue[graph->cptr[j]]; 1069 } 1070 ierr = PetscFree(queue_global);CHKERRQ(ierr); 1071 graph->queue_sorted = PETSC_TRUE; 1072 /* save information on subsets (needed when analyzing the connected components) */ 1073 if (graph->ncc) { 1074 ierr = PetscMalloc2(graph->ncc,&graph->subset_size,graph->ncc,&graph->subset_idxs);CHKERRQ(ierr); 1075 ierr = PetscMalloc1(graph->cptr[graph->ncc],&graph->subset_idxs[0]);CHKERRQ(ierr); 1076 ierr = PetscMemzero(graph->subset_idxs[0],graph->cptr[graph->ncc]*sizeof(PetscInt));CHKERRQ(ierr); 1077 for (j=1;j<graph->ncc;j++) { 1078 graph->subset_size[j-1] = graph->cptr[j] - graph->cptr[j-1]; 1079 graph->subset_idxs[j] = graph->subset_idxs[j-1] + graph->subset_size[j-1]; 1080 } 1081 graph->subset_size[graph->ncc-1] = graph->cptr[graph->ncc] - graph->cptr[graph->ncc-1]; 1082 ierr = PetscMemcpy(graph->subset_idxs[0],graph->queue,graph->cptr[graph->ncc]*sizeof(PetscInt));CHKERRQ(ierr); 1083 } 1084 1085 /* renumber reference nodes */ 1086 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)(graph->l2gmap)),graph->ncc,graph->subset_ref_node,PETSC_COPY_VALUES,&subset_n);CHKERRQ(ierr); 1087 ierr = ISLocalToGlobalMappingApplyIS(graph->l2gmap,subset_n,&subset);CHKERRQ(ierr); 1088 ierr = ISDestroy(&subset_n);CHKERRQ(ierr); 1089 ierr = ISRenumber(subset,NULL,NULL,&subset_n);CHKERRQ(ierr); 1090 ierr = ISDestroy(&subset);CHKERRQ(ierr); 1091 ierr = ISGetLocalSize(subset_n,&k);CHKERRQ(ierr); 1092 if (k != graph->ncc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid size of new subset! %D != %D",k,graph->ncc); 1093 ierr = ISGetIndices(subset_n,&is_indices);CHKERRQ(ierr); 1094 ierr = PetscMemcpy(graph->subset_ref_node,is_indices,graph->ncc*sizeof(PetscInt));CHKERRQ(ierr); 1095 ierr = ISRestoreIndices(subset_n,&is_indices);CHKERRQ(ierr); 1096 ierr = ISDestroy(&subset_n);CHKERRQ(ierr); 1097 1098 /* free workspace */ 1099 ierr = VecDestroy(&local_vec);CHKERRQ(ierr); 1100 ierr = VecDestroy(&local_vec2);CHKERRQ(ierr); 1101 ierr = VecDestroy(&global_vec);CHKERRQ(ierr); 1102 ierr = VecScatterDestroy(&scatter_ctx);CHKERRQ(ierr); 1103 graph->setupcalled = PETSC_TRUE; 1104 PetscFunctionReturn(0); 1105 } 1106 1107 #undef __FUNCT__ 1108 #define __FUNCT__ "PCBDDCGraphResetCSR" 1109 PetscErrorCode PCBDDCGraphResetCSR(PCBDDCGraph graph) 1110 { 1111 PetscErrorCode ierr; 1112 1113 PetscFunctionBegin; 1114 if (graph->freecsr) { 1115 ierr = PetscFree(graph->xadj);CHKERRQ(ierr); 1116 ierr = PetscFree(graph->adjncy);CHKERRQ(ierr); 1117 } else { 1118 graph->xadj = NULL; 1119 graph->adjncy = NULL; 1120 } 1121 graph->freecsr = PETSC_FALSE; 1122 graph->nvtxs_csr = 0; 1123 PetscFunctionReturn(0); 1124 } 1125 1126 #undef __FUNCT__ 1127 #define __FUNCT__ "PCBDDCGraphReset" 1128 PetscErrorCode PCBDDCGraphReset(PCBDDCGraph graph) 1129 { 1130 PetscErrorCode ierr; 1131 1132 PetscFunctionBegin; 1133 graph->setupcalled = PETSC_FALSE; 1134 ierr = ISLocalToGlobalMappingDestroy(&graph->l2gmap);CHKERRQ(ierr); 1135 ierr = PetscFree(graph->subset_ncc);CHKERRQ(ierr); 1136 ierr = PetscFree(graph->subset_ref_node);CHKERRQ(ierr); 1137 if (graph->nvtxs) { 1138 ierr = PetscFree(graph->neighbours_set[0]);CHKERRQ(ierr); 1139 } 1140 ierr = PetscBTDestroy(&graph->touched);CHKERRQ(ierr); 1141 ierr = PetscFree5(graph->count, 1142 graph->neighbours_set, 1143 graph->subset, 1144 graph->which_dof, 1145 graph->special_dof);CHKERRQ(ierr); 1146 ierr = PetscFree2(graph->cptr,graph->queue);CHKERRQ(ierr); 1147 if (graph->mirrors) { 1148 ierr = PetscFree(graph->mirrors_set[0]);CHKERRQ(ierr); 1149 } 1150 ierr = PetscFree2(graph->mirrors,graph->mirrors_set);CHKERRQ(ierr); 1151 if (graph->subset_idxs) { 1152 ierr = PetscFree(graph->subset_idxs[0]);CHKERRQ(ierr); 1153 } 1154 ierr = PetscFree2(graph->subset_size,graph->subset_idxs);CHKERRQ(ierr); 1155 ierr = ISDestroy(&graph->dirdofs);CHKERRQ(ierr); 1156 ierr = ISDestroy(&graph->dirdofsB);CHKERRQ(ierr); 1157 if (graph->n_local_subs) { 1158 ierr = PetscFree(graph->local_subs);CHKERRQ(ierr); 1159 } 1160 graph->has_dirichlet = PETSC_FALSE; 1161 graph->nvtxs = 0; 1162 graph->nvtxs_global = 0; 1163 graph->n_subsets = 0; 1164 graph->custom_minimal_size = 1; 1165 graph->n_local_subs = 0; 1166 PetscFunctionReturn(0); 1167 } 1168 1169 #undef __FUNCT__ 1170 #define __FUNCT__ "PCBDDCGraphInit" 1171 PetscErrorCode PCBDDCGraphInit(PCBDDCGraph graph, ISLocalToGlobalMapping l2gmap, PetscInt N) 1172 { 1173 PetscInt n; 1174 PetscErrorCode ierr; 1175 1176 PetscFunctionBegin; 1177 PetscValidPointer(graph,1); 1178 PetscValidHeaderSpecific(l2gmap,IS_LTOGM_CLASSID,2); 1179 PetscValidLogicalCollectiveInt(l2gmap,N,3); 1180 /* raise an error if already allocated */ 1181 if (graph->nvtxs_global) SETERRQ(PetscObjectComm((PetscObject)l2gmap),PETSC_ERR_PLIB,"BDDCGraph already initialized"); 1182 /* set number of vertices */ 1183 ierr = PetscObjectReference((PetscObject)l2gmap);CHKERRQ(ierr); 1184 graph->l2gmap = l2gmap; 1185 ierr = ISLocalToGlobalMappingGetSize(l2gmap,&n);CHKERRQ(ierr); 1186 graph->nvtxs = n; 1187 graph->nvtxs_global = N; 1188 /* allocate used space */ 1189 ierr = PetscBTCreate(graph->nvtxs,&graph->touched);CHKERRQ(ierr); 1190 ierr = PetscMalloc5(graph->nvtxs,&graph->count, 1191 graph->nvtxs,&graph->neighbours_set, 1192 graph->nvtxs,&graph->subset, 1193 graph->nvtxs,&graph->which_dof, 1194 graph->nvtxs,&graph->special_dof);CHKERRQ(ierr); 1195 /* zeroes memory */ 1196 ierr = PetscMemzero(graph->count,graph->nvtxs*sizeof(PetscInt));CHKERRQ(ierr); 1197 ierr = PetscMemzero(graph->subset,graph->nvtxs*sizeof(PetscInt));CHKERRQ(ierr); 1198 /* use -1 as a default value for which_dof array */ 1199 for (n=0;n<graph->nvtxs;n++) graph->which_dof[n] = -1; 1200 ierr = PetscMemzero(graph->special_dof,graph->nvtxs*sizeof(PetscInt));CHKERRQ(ierr); 1201 /* zeroes first pointer to neighbour set */ 1202 if (graph->nvtxs) { 1203 graph->neighbours_set[0] = 0; 1204 } 1205 /* zeroes workspace for values of ncc */ 1206 graph->subset_ncc = 0; 1207 graph->subset_ref_node = 0; 1208 /* default flag for csr */ 1209 graph->freecsr = PETSC_FALSE; 1210 PetscFunctionReturn(0); 1211 } 1212 1213 #undef __FUNCT__ 1214 #define __FUNCT__ "PCBDDCGraphDestroy" 1215 PetscErrorCode PCBDDCGraphDestroy(PCBDDCGraph* graph) 1216 { 1217 PetscErrorCode ierr; 1218 1219 PetscFunctionBegin; 1220 ierr = PCBDDCGraphReset(*graph);CHKERRQ(ierr); 1221 ierr = PetscFree(*graph);CHKERRQ(ierr); 1222 PetscFunctionReturn(0); 1223 } 1224 1225 #undef __FUNCT__ 1226 #define __FUNCT__ "PCBDDCGraphCreate" 1227 PetscErrorCode PCBDDCGraphCreate(PCBDDCGraph *graph) 1228 { 1229 PCBDDCGraph new_graph; 1230 PetscErrorCode ierr; 1231 1232 PetscFunctionBegin; 1233 ierr = PetscNew(&new_graph);CHKERRQ(ierr); 1234 new_graph->custom_minimal_size = 1; 1235 new_graph->commsizelimit = 1; 1236 *graph = new_graph; 1237 PetscFunctionReturn(0); 1238 } 1239