/[escript]/trunk/ripley/src/blocktools.cpp
ViewVC logotype

Contents of /trunk/ripley/src/blocktools.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6000 - (show annotations)
Tue Mar 1 00:24:43 2016 UTC (3 years, 1 month ago) by caltinay
File size: 11969 byte(s)
a few more include rearrangements.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2014-2016 by The 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 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include "blocktools.h"
18
19 #include <cstring> // for memset
20 #include <iostream> // for the debug method
21
22 using namespace std;
23
24 BlockGrid::BlockGrid(coord_t x, coord_t y, coord_t z)
25 : xmax(x), ymax(y), zmax(z)
26 {}
27
28 neighbourID_t BlockGrid::getNID(coord_t x, coord_t y, coord_t z) const
29 {
30 return x+y*(xmax+1)+z*(xmax+1)*(ymax+1);
31 }
32
33 // generate all incoming com messages for this block.
34 // for each subblock (27 of them), there may be an x, y, z direction to search in
35 void generateInNeighbours(coord_t blockx, coord_t blocky, coord_t blockz, messvec& v);
36
37
38 // generate all outgoing com messages for this block
39 void generateOutNeighbours(coord_t blockx, coord_t blocky, coord_t blockz, messvec& v);
40
41 // generate all incoming com messages for this block.
42 // for each subblock (27 of them), there may be an x, y, z direction to search in
43 void BlockGrid::generateInNeighbours(coord_t blockx, coord_t blocky, coord_t blockz, messvec& v)
44 {
45 neighbourID_t myid=getNID(blockx, blocky, blockz);
46 unsigned char deltax=0;
47 unsigned char deltay=0;
48 unsigned char deltaz=blockz?1:0; // do not look to a lower layer if we are on the bottom layer of _blocks_.
49 for (unsigned char z=0;z<3;++z)
50 {
51 deltay=blocky?1:0; // do not look at a lower row of blocks if there isn't one
52 for (unsigned char y=0;y<3;++y)
53 {
54 deltax=blockx?1:0;
55 for (unsigned char x=0;x<3;++x)
56 {
57 // now we will have a set of possible directions to look in (delta?==0 means don't attempt to
58 // change this component
59 if (deltax+deltay+deltaz) // if we have an import to do
60 {
61 coord_t srcx=blockx-deltax;
62 coord_t srcy=blocky-deltay;
63 coord_t srcz=blockz-deltaz;
64 message m;
65 m.sourceID=getNID(srcx, srcy, srcz);
66 m.destID=myid;
67 m.tag=getTag(x,y,z,deltax==1, deltay==1, deltaz==1);
68 m.srcbuffid=getSrcBuffID(x,y,z,deltax==1, deltay==1, deltaz==1);
69 m.destbuffid=x+y*3+z*9;
70 v.push_back(m);
71 }
72 deltax=0; // we are no longer on the left hand edge
73 }
74 deltay=0; // only y=0 imports
75 }
76 deltaz=0; // since only the bottom sublayer looks to the row below
77 }
78 }
79
80
81 // generate all outgoing com messages for this block
82 void BlockGrid::generateOutNeighbours(coord_t blockx, coord_t blocky, coord_t blockz, messvec& v)
83 {
84 messvec vv;
85 neighbourID_t myid=getNID(blockx, blocky, blockz);
86 for (unsigned char z=0;z<2;++z)
87 {
88 if (z && (blockz==zmax)) // don't look up if there is no up
89 {
90 break;
91 }
92 for (unsigned char y=0;y<2;++y)
93 {
94 if (y && (blocky==ymax))
95 {
96 break;
97 }
98 for (unsigned char x=0;x<2;++x)
99 {
100 if (x && (blockx==xmax))
101 {
102 break;
103 }
104 if (x+y+z>0) // don't look at your own neighbours
105 {
106 generateInNeighbours(blockx+x, blocky+y, blockz+z, vv);
107 }
108 }
109 }
110 }
111 // Now we reverse the direction of the coms (since we want Out not In
112 for (int i=0;i<vv.size();++i)
113 {
114 if (vv[i].sourceID==myid)
115 {
116 v.push_back(vv[i]);
117 }
118 }
119 }
120
121
122 double* Block::getOutBuffer(unsigned char subx, unsigned char suby, unsigned char subz)
123 {
124 unsigned char bid=subx+suby*3+subz*3*3; // (bid is "blockid")
125 if (bid==13) // there is no buffer for block 13
126 {
127 return 0; // don't ask for this buffer because refusal may offend
128 }
129 return outbuffptr[bid];
130 }
131
132 double* Block::getInBuffer(unsigned char subx, unsigned char suby, unsigned char subz)
133 {
134 unsigned char bid=subx+suby*3+subz*3*3; // (bid is "blockid")
135 if (bid==13) // there is no buffer for block 13
136 {
137 return 0; // don't ask for this buffer because refusal may offend
138 }
139 return inbuffptr[bid];
140 }
141
142 size_t Block::getBuffSize(unsigned char subx, unsigned char suby, unsigned char subz)
143 {
144 unsigned char bid=subx+suby*3+subz*3*3; // (bid is "blockid")
145 if (bid==13) // there is no buffer for block 13
146 {
147 return 0;
148 }
149 return dims[bid][0]*dims[bid][1]*dims[bid][2]*dpsize;
150 }
151
152 double* Block::getOutBuffer(unsigned char bid)
153 {
154 if (bid==13) // there is no buffer for block 13
155 {
156 return 0; // don't ask for this buffer because refusal may offend
157 }
158 return outbuffptr[bid];
159 }
160
161 double* Block::getInBuffer(unsigned char bid)
162 {
163 if (bid==13) // there is no buffer for block 13
164 {
165 return 0; // don't ask for this buffer because refusal may offend
166 }
167 return inbuffptr[bid];
168 }
169
170 size_t Block::getBuffSize(unsigned char bid)
171 {
172 if (bid==13) // there is no buffer for block 13
173 {
174 return 0;
175 }
176 return dims[bid][0]*dims[bid][1]*dims[bid][2]*dpsize;
177 }
178
179
180
181 Block::~Block()
182 {
183 delete[] inbuff;
184 delete[] outbuff;
185 }
186
187 void Block::populateDimsTable()
188 {
189 for (int i=0;i<27;++i)
190 {
191 for (int j=0;j<3;++j)
192 {
193 dims[i][j]=inset;
194 }
195 }
196 for (int i=1;i<27;i+=3)
197 {
198 dims[i][0]=xmidlen;
199 }
200 for (int l=0;l<3;++l)
201 {
202 dims[3+l*9][1]=ymidlen;
203 dims[4+l*9][1]=ymidlen;
204 dims[5+l*9][1]=ymidlen;
205 }
206 for (int i=9;i<18;++i)
207 {
208 dims[i][2]=zmidlen;
209 }
210 }
211
212
213 // gives the relative start location within a flat
214 void Block::populateOffsetTable(size_t inset, size_t xmidlen, size_t ymidlen, size_t zmidlen)
215 {
216
217 size_t cur=0;
218 for (int i=0;i<27;++i)
219 {
220 flatoffsets[i]=cur;
221 cur+=dims[i][0]*dims[i][1]*dims[i][2]*dpsize;
222 }
223 for (int i=0;i<13;++i)
224 {
225 buffoffsets[i]=flatoffsets[i];
226 }
227 buffoffsets[13]=0;
228 for (int i=14;i<27;++i)
229 {
230 buffoffsets[i]=flatoffsets[i]-flatoffsets[14]+flatoffsets[13];
231 }
232
233 }
234
235
236 void Block::createBuffArrays(double* startaddress, double* buffptr[27], size_t inset, size_t xmidlen, size_t ymidlen, size_t zmidlen)
237 {
238 buffptr[0]=startaddress;
239 for (int i=0;i<27;++i)
240 {
241 buffptr[i]=startaddress+buffoffsets[i];
242 }
243 buffptr[13]=0; // since the middle should never be sent anywhere
244 }
245
246 void Block::setUsed(unsigned char buffid)
247 {
248 used[buffid]=true;
249 }
250
251 void Block::copyAllToBuffer(double* src)
252 {
253 for (unsigned char i=0;i<13;++i)
254 {
255 copyToBuffer(i, src);
256
257 }
258 for (unsigned char i=14;i<27;++i)
259 {
260 copyToBuffer(i, src);
261 }
262 }
263
264 void Block::copyUsedFromBuffer(double* dest)
265 {
266 for (unsigned char i=0;i<27;++i)
267 {
268 if (used[i])
269 {
270 copyFromBuffer(i, dest);
271 }
272 }
273 }
274
275 // s? specifiy the size (in points) of each dimension
276 // maxb? gives the largest block number in each dimension in the overall grid (number from zero)
277 Block::Block(size_t sx, size_t sy, size_t sz, size_t inset, size_t xmidlen,
278 size_t ymidlen, size_t zmidlen, unsigned int dpp) : dpsize(dpp)
279 {
280 this->sx=sx;
281 this->sy=sy;
282 this->sz=sz;
283 this->inset=inset;
284 this->xmidlen=xmidlen;
285 this->ymidlen=ymidlen;
286 this->zmidlen=zmidlen;
287 populateDimsTable();
288
289 size_t totalbuff=0;
290 for (int i=0;i<27;++i)
291 {
292 used[i]=false;
293 totalbuff+=dims[i][0]*dims[i][1]*dims[i][2];
294 }
295 totalbuff-=dims[13][0]*dims[13][1]*dims[13][2];
296
297 totalbuff*=dpsize; // to account for non-scalars
298 inbuff=new double[totalbuff];
299 outbuff=new double[totalbuff];
300 memset(inbuff,0,totalbuff*sizeof(double));
301 memset(outbuff,0,totalbuff*sizeof(double));
302 populateOffsetTable(inset, xmidlen, ymidlen, zmidlen);
303 createBuffArrays(inbuff, inbuffptr, inset, xmidlen, ymidlen, zmidlen);
304 createBuffArrays(outbuff, outbuffptr, inset, xmidlen, ymidlen, zmidlen);
305
306 }
307
308 // where does the subblock specified start in a source array
309 size_t Block::startOffset(unsigned char subx, unsigned char suby, unsigned char subz)
310 {
311 size_t off=0;
312 off+=((subx==0) ? 0 :((subx==1)?inset:(inset+xmidlen) ));
313 // now the y component
314 size_t ystep=((suby==0) ? 0 :((suby==1)?inset:(inset+ymidlen) ));
315 size_t zstep=((subz==0) ? 0 :((subz==1)?inset:(inset+zmidlen) ));
316 off+=ystep*(2*inset+xmidlen);
317 off+=zstep*(2*inset+xmidlen)*(2*inset+ymidlen);
318 return off*dpsize;
319 }
320
321
322 // This is only intended for debugging, so I have not made it fast
323 void Block::displayBlock(unsigned char subx, unsigned char suby, unsigned char subz, bool out)
324 {
325
326 unsigned char bid=subx+suby*3+subz*3*3; // (bid is "blockid")
327 double* b=out?outbuffptr[bid]:inbuffptr[bid];
328 for (int z=0;z<dims[bid][2];++z)
329 {
330 std::cout << std::endl << "Row " << z << std::endl;
331
332 for (int y=0;y<dims[bid][1];++y)
333 {
334 for (int x=0;x<dims[bid][0];++x)
335 {
336 if (dpsize==1)
337 {
338 std::cout << b[x+y*dims[bid][0]+z*dims[bid][1]*dims[bid][0]] << ' ';
339 }
340 else
341 {
342 std::cout << '(';
343 for (int i=0;i<dpsize;++i)
344 {
345 std::cout << b[(x+y*dims[bid][0]+z*dims[bid][1]*dims[bid][0])*dpsize+i] << ", ";
346 }
347 std::cout << ") ";
348 }
349 }
350 std::cout << std::endl;
351 }
352 }
353 }
354
355
356 // Copy a 3d region from a flat array into a buffer
357 void Block::copyToBuffer(unsigned char bid, double* src)
358 {
359 if (bid==13) // there is no buffer for block 13
360 {
361 return;
362 }
363 // where does the strided content start?
364 double* start=src+startOffset(bid%3, bid%9/3, bid/9);
365 double* dest=outbuffptr[bid];
366
367 size_t zlim=dims[bid][2]; // how big is the block
368 size_t ylim=dims[bid][1];
369 size_t xlim=dims[bid][0];
370 size_t totaly=(2*inset+ymidlen);
371 for (size_t z=0;z<zlim;++z)
372 {
373 for (size_t y=0;y<ylim;++y)
374 {
375 memcpy(dest, start, xlim*sizeof(double)*dpsize);
376 dest+=xlim*dpsize;
377 start+=(2*inset+xmidlen)*dpsize;
378 }
379 // we are at the end of the slab so we need to jump to the next level up
380 start+=((totaly-ylim)*(2*inset+xmidlen))*dpsize;
381 }
382 }
383
384
385 // Copy a 3d region from a buffer into a flat array
386 void Block::copyFromBuffer(unsigned char bid, double* dest)
387 {
388 if (bid==13) // there is no buffer for block 13
389 {
390 return;
391 }
392 double* start=dest+startOffset(bid%3, bid%9/3, bid/9);
393 double* src=inbuffptr[bid];
394 size_t zlim=dims[bid][2]; // how big is the block
395 size_t ylim=dims[bid][1];
396 size_t xlim=dims[bid][0];
397 size_t totaly=(2*inset+ymidlen);
398 for (size_t z=0;z<zlim;++z)
399 {
400 for (size_t y=0;y<ylim;++y)
401 {
402 memcpy(start, src, xlim*sizeof(double)*dpsize);
403 src+=xlim*dpsize;
404 start+=(2*inset+xmidlen)*dpsize;
405 }
406 // we are at the end of the slab so we need to jump to the next level up
407 start+=(totaly-ylim)*(2*inset+xmidlen)*dpsize;
408 }
409 }
410
411 // Returns the MPI message tag to use for a transfer between the two subblocks
412 int getTag(unsigned char sourcex, unsigned char sourcey, unsigned char sourcez, unsigned char targetx, unsigned char targety, unsigned char targetz)
413 {
414 return sourcex*100000+sourcey*10000+sourcez*1000+targetx*100+targety*10+targetz;
415 }
416
417 // computes the tag based on the destination and the direction it comes from
418 // the booleans indicate whether a negative shift in that direction is required
419 int getTag(unsigned char destx, unsigned char desty, unsigned char destz, bool deltax, bool deltay, bool deltaz)
420 {
421 unsigned char sourcex=deltax?2:destx;
422 unsigned char sourcey=deltay?2:desty;
423 unsigned char sourcez=deltaz?2:destz;
424
425 return sourcex*100000+sourcey*10000+sourcez*1000+destx*100+desty*10+destz;
426 }
427
428
429 // the booleans indicate whether a negative shift in that direction is required
430 unsigned char getSrcBuffID(unsigned char destx, unsigned char desty, unsigned char destz, bool deltax, bool deltay, bool deltaz)
431 {
432 unsigned char sourcex=deltax?2:destx;
433 unsigned char sourcey=deltay?2:desty;
434 unsigned char sourcez=deltaz?2:destz;
435
436 return sourcex+sourcey*3+sourcez*9;
437 }
438

  ViewVC Help
Powered by ViewVC 1.1.26