/[escript]/branches/doubleplusgood/finley/src/ElementFile_distributeByRankOfDOF.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/finley/src/ElementFile_distributeByRankOfDOF.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4332 - (show annotations)
Thu Mar 21 04:21:14 2013 UTC (5 years, 11 months ago) by jfenwick
File size: 10192 byte(s)
like sand though the hourglass
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 /************************************************************************************/
18
19 /* Finley: ElementFile: this will redistribute the Elements including overlap by rank */
20
21 /************************************************************************************/
22
23 #include "ElementFile.h"
24 #ifdef _OPENMP
25 #include <omp.h>
26 #endif
27
28 /************************************************************************************/
29
30 void Finley_ElementFile_distributeByRankOfDOF(Finley_ElementFile* self, Esys_MPI_rank* mpiRankOfDOF, index_t* Id) {
31 size_t size_size;
32 Esys_MPI_rank myRank, p, *Owner_buffer=NULL, loc_proc_mask_max;
33 dim_t e, j, i, size, *send_count=NULL, *recv_count=NULL, *newOwner=NULL, *loc_proc_mask=NULL, *loc_send_count=NULL,
34 newNumElements, numElementsInBuffer, numNodes, NN;
35 index_t *send_offset=NULL, *recv_offset=NULL, *Id_buffer=NULL, *Tag_buffer=NULL, *Nodes_buffer=NULL, k;
36 bool_t *proc_mask=NULL;
37 #ifdef ESYS_MPI
38 dim_t numRequests=0;
39 MPI_Request* mpi_requests=NULL;
40 MPI_Status* mpi_stati=NULL;
41 #endif
42 if (self==NULL) return;
43 myRank=self->MPIInfo->rank;
44 size=self->MPIInfo->size;
45 size_size=size*sizeof(dim_t);
46 numNodes=self->numNodes;
47 NN=self->numNodes;
48 if (size>1) {
49 #ifdef ESYS_MPI
50 mpi_requests=new MPI_Request[8*size];
51 mpi_stati=new MPI_Status[8*size];
52 Finley_checkPtr(mpi_requests);
53 Finley_checkPtr(mpi_stati);
54 #endif
55
56 /* count the number of elements that have to be sent to each processor (send_count)
57 and define a new element owner as the processor with the largest number of DOFs and the smallest id */
58 send_count=new dim_t[size];
59 recv_count=new dim_t[size];
60 newOwner=new Esys_MPI_rank[self->numElements];
61 if ( !( Finley_checkPtr(send_count) || Finley_checkPtr(recv_count) || Finley_checkPtr(newOwner) ) ) {
62 memset(send_count, 0, size_size);
63 #pragma omp parallel private(p,loc_proc_mask,loc_send_count)
64 {
65 loc_proc_mask=new dim_t[size];
66 loc_send_count=new dim_t[size];
67 memset(loc_send_count, 0, size_size);
68 #pragma omp for private(e,j,loc_proc_mask_max) schedule(static)
69 for (e=0;e<self->numElements;e++) {
70 if (self->Owner[e] == myRank) {
71 newOwner[e]=myRank;
72 memset(loc_proc_mask, 0, size_size);
73 for(j=0;j<numNodes;j++) {
74 p=mpiRankOfDOF[self->Nodes[INDEX2(j,e,NN)]];
75 loc_proc_mask[p]++;
76 }
77 loc_proc_mask_max=0;
78 for (p=0;p<size;++p) {
79 if (loc_proc_mask[p]>0) loc_send_count[p]++;
80 if (loc_proc_mask[p]>loc_proc_mask_max) {
81 newOwner[e]=p;
82 loc_proc_mask_max=loc_proc_mask[p];
83 }
84 }
85 } else {
86 newOwner[e]=-1;
87 }
88 }
89 #pragma omp critical
90 {
91 for (p=0;p<size;++p) send_count[p]+=loc_send_count[p];
92 }
93 delete[] loc_proc_mask;
94 delete[] loc_send_count;
95 }
96 #ifdef ESYS_MPI
97 MPI_Alltoall(send_count,1,MPI_INT,recv_count,1,MPI_INT,self->MPIInfo->comm);
98 #else
99 for (p=0;p<size;++p) recv_count[p]=send_count[p];
100 #endif
101 /* get the new number of elements for this processor */
102 newNumElements=0;
103 for (p=0;p<size;++p) newNumElements+=recv_count[p];
104
105 /* get the new number of elements for this processor */
106 numElementsInBuffer=0;
107 for (p=0;p<size;++p) numElementsInBuffer+=send_count[p];
108 /* allocate buffers */
109 Id_buffer=new index_t[numElementsInBuffer];
110 Tag_buffer=new index_t[numElementsInBuffer];
111 Owner_buffer=new Esys_MPI_rank[numElementsInBuffer];
112 Nodes_buffer=new index_t[numElementsInBuffer*NN];
113 send_offset=new index_t[size];
114 recv_offset=new index_t[size];
115 proc_mask=new bool_t[size];
116 if ( !( Finley_checkPtr(Id_buffer) || Finley_checkPtr(Tag_buffer) || Finley_checkPtr(Owner_buffer) ||
117 Finley_checkPtr(Nodes_buffer) || Finley_checkPtr(send_offset) || Finley_checkPtr(recv_offset) ||
118 Finley_checkPtr(proc_mask) )) {
119
120 /* calculate the offsets for the processor buffers */
121 recv_offset[0]=0;
122 for (p=0;p<size-1;++p) recv_offset[p+1]=recv_offset[p]+recv_count[p];
123 send_offset[0]=0;
124 for (p=0;p<size-1;++p) send_offset[p+1]=send_offset[p]+send_count[p];
125
126 memset(send_count, 0, size_size);
127 /* copy element into buffers. proc_mask makes sure that an element is copied once only for each processor */
128 for (e=0;e<self->numElements;e++) {
129 if (self->Owner[e] == myRank) {
130 memset(proc_mask, TRUE, size_size);
131 for(j=0;j<numNodes;j++) {
132 p=mpiRankOfDOF[self->Nodes[INDEX2(j,e,NN)]];
133 if (proc_mask[p]) {
134 k=send_offset[p]+send_count[p];
135 Id_buffer[k]=self->Id[e];
136 Tag_buffer[k]=self->Tag[e];
137 Owner_buffer[k]=newOwner[e];
138 for (i=0;i<numNodes;i++) Nodes_buffer[INDEX2(i,k,NN)]=Id[self->Nodes[INDEX2(i,e,NN)]];
139 send_count[p]++;
140 proc_mask[p]=FALSE;
141 }
142 }
143 }
144 }
145 /* allocate new tables */
146 Finley_ElementFile_allocTable(self,newNumElements);
147
148 /* start to receive new elements */
149 for (p=0;p<size;++p) {
150 if (recv_count[p]>0) {
151 #ifdef ESYS_MPI
152 MPI_Irecv(&(self->Id[recv_offset[p]]), recv_count[p],
153 MPI_INT, p, self->MPIInfo->msg_tag_counter+myRank,
154 self->MPIInfo->comm, &mpi_requests[numRequests]);
155 numRequests++;
156 MPI_Irecv(&(self->Tag[recv_offset[p]]), recv_count[p],
157 MPI_INT, p, self->MPIInfo->msg_tag_counter+size+myRank,
158 self->MPIInfo->comm, &mpi_requests[numRequests]);
159 numRequests++;
160 MPI_Irecv(&(self->Owner[recv_offset[p]]), recv_count[p],
161 MPI_INT, p, self->MPIInfo->msg_tag_counter+2*size+myRank,
162 self->MPIInfo->comm, &mpi_requests[numRequests]);
163 numRequests++;
164 MPI_Irecv(&(self->Nodes[recv_offset[p]*NN]), recv_count[p]*NN,
165 MPI_INT, p, self->MPIInfo->msg_tag_counter+3*size+myRank,
166 self->MPIInfo->comm, &mpi_requests[numRequests]);
167 numRequests++;
168 #endif
169 }
170 }
171 /* now the buffers can be sent away */
172 for (p=0;p<size;++p) {
173 if (send_count[p]>0) {
174 #ifdef ESYS_MPI
175 MPI_Issend(&(Id_buffer[send_offset[p]]), send_count[p],
176 MPI_INT, p, self->MPIInfo->msg_tag_counter+p,
177 self->MPIInfo->comm, &mpi_requests[numRequests]);
178 numRequests++;
179 MPI_Issend(&(Tag_buffer[send_offset[p]]), send_count[p],
180 MPI_INT, p, self->MPIInfo->msg_tag_counter+size+p,
181 self->MPIInfo->comm, &mpi_requests[numRequests]);
182 numRequests++;
183 MPI_Issend(&(Owner_buffer[send_offset[p]]), send_count[p],
184 MPI_INT, p, self->MPIInfo->msg_tag_counter+2*size+p,
185 self->MPIInfo->comm, &mpi_requests[numRequests]);
186 numRequests++;
187 MPI_Issend(&(Nodes_buffer[send_offset[p]*NN]), send_count[p]*NN,
188 MPI_INT, p, self->MPIInfo->msg_tag_counter+3*size+p,
189 self->MPIInfo->comm, &mpi_requests[numRequests]);
190 numRequests++;
191 #endif
192
193 }
194 }
195 self->MPIInfo->msg_tag_counter+=4*size;
196 /* wait for the requests to be finalized */
197 #ifdef ESYS_MPI
198 MPI_Waitall(numRequests,mpi_requests,mpi_stati);
199 #endif
200 }
201 /* clear buffer */
202 delete[] Id_buffer;
203 delete[] Tag_buffer;
204 delete[] Owner_buffer;
205 delete[] Nodes_buffer;
206 delete[] send_offset;
207 delete[] recv_offset;
208 delete[] proc_mask;
209 }
210 #ifdef ESYS_MPI
211 delete[] mpi_requests;
212 delete[] mpi_stati;
213 #endif
214 delete[] send_count;
215 delete[] recv_count;
216 delete[] newOwner;
217 } else {
218 #pragma omp for private(e,i) schedule(static)
219 for (e=0;e<self->numElements;e++) {
220 self->Owner[e]=myRank;
221 for (i=0;i<numNodes;i++) self->Nodes[INDEX2(i,e,NN)]=Id[self->Nodes[INDEX2(i,e,NN)]];
222 }
223 }
224 return;
225 }
226

  ViewVC Help
Powered by ViewVC 1.1.26