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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4696 - (show annotations)
Wed Feb 19 07:29:50 2014 UTC (5 years, 2 months ago) by jfenwick
File size: 11960 byte(s)
Correcting python libname for sage.
Fixing bug where the number of values in the shape was not considered for buffer and message size.

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

  ViewVC Help
Powered by ViewVC 1.1.26