/[escript]/branches/clazy/escriptcore/src/DataLazy.cpp
ViewVC logotype

Contents of /branches/clazy/escriptcore/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1888 - (show annotations)
Wed Oct 15 04:00:21 2008 UTC (10 years, 6 months ago) by jfenwick
Original Path: branches/schroedinger/escript/src/DataLazy.cpp
File size: 13096 byte(s)
Branch commit.
Added more binary ops.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 #include "DataLazy.h"
16 #ifdef USE_NETCDF
17 #include <netcdfcpp.h>
18 #endif
19 #ifdef PASO_MPI
20 #include <mpi.h>
21 #endif
22 #ifdef _OPENMP
23 #include <omp.h>
24 #endif
25 #include "FunctionSpace.h"
26 #include "DataTypes.h"
27 #include "Data.h"
28 #include "UnaryFuncs.h" // for escript::fsign
29
30 using namespace std;
31 using namespace boost;
32
33 namespace escript
34 {
35
36 const std::string&
37 opToString(ES_optype op);
38
39 namespace
40 {
41
42
43
44 enum ES_opgroup
45 {
46 G_UNKNOWN,
47 G_IDENTITY,
48 G_BINARY,
49 G_UNARY
50 };
51
52
53
54
55 string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","sin","cos","tan",
56 "asin","acos","atan","sinh","cosh","tanh","erf",
57 "asinh","acosh","atanh",
58 "log10","log","sign","abs","neg","pos","exp","sqrt",
59 "1/","where>0","where<0","where>=0","where<=0"};
60 int ES_opcount=32;
61 ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY,G_UNARY,G_UNARY,G_UNARY, //9
62 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 16
63 G_UNARY,G_UNARY,G_UNARY, // 19
64 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 27
65 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};
66 inline
67 ES_opgroup
68 getOpgroup(ES_optype op)
69 {
70 return opgroups[op];
71 }
72
73 // return the FunctionSpace of the result of "left op right"
74 FunctionSpace
75 resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
76 {
77 // perhaps this should call interpolate and throw or something?
78 // maybe we need an interpolate node -
79 // that way, if interpolate is required in any other op we can just throw a
80 // programming error exception.
81
82
83 if (left->getFunctionSpace()!=right->getFunctionSpace())
84 {
85 throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");
86 }
87 return left->getFunctionSpace();
88 }
89
90 // return the shape of the result of "left op right"
91 DataTypes::ShapeType
92 resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
93 {
94 if (left->getShape()!=right->getShape())
95 {
96 throw DataException("Shapes not the same - shapes must match for lazy data.");
97 }
98 return left->getShape();
99 }
100
101 size_t
102 resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
103 {
104 switch (getOpgroup(op))
105 {
106 case G_BINARY: return left->getLength();
107 case G_UNARY: return left->getLength();
108 default:
109 throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
110 }
111 }
112
113 int
114 calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
115 {
116 switch(getOpgroup(op))
117 {
118 case G_IDENTITY: return 1;
119 case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
120 case G_UNARY: return max(left->getBuffsRequired(),1);
121 default:
122 throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
123 }
124 }
125
126
127
128 } // end anonymous namespace
129
130
131 const std::string&
132 opToString(ES_optype op)
133 {
134 if (op<0 || op>=ES_opcount)
135 {
136 op=UNKNOWNOP;
137 }
138 return ES_opstrings[op];
139 }
140
141
142 DataLazy::DataLazy(DataAbstract_ptr p)
143 : parent(p->getFunctionSpace(),p->getShape()),
144 m_op(IDENTITY)
145 {
146 if (p->isLazy())
147 {
148 // TODO: fix this. We could make the new node a copy of p?
149 // I don't want identity of Lazy.
150 // Question: Why would that be so bad?
151 // Answer: We assume that the child of ID is something we can call getVector on
152 throw DataException("Programmer error - attempt to create identity from a DataLazy.");
153 }
154 else
155 {
156 m_id=dynamic_pointer_cast<DataReady>(p);
157 }
158 m_length=p->getLength();
159 m_buffsRequired=1;
160 m_samplesize=getNumDPPSample()*getNoValues();
161 cout << "(1)Lazy created with " << m_samplesize << endl;
162 }
163
164 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
165 : parent(left->getFunctionSpace(),left->getShape()),
166 m_op(op)
167 {
168 if (getOpgroup(op)!=G_UNARY)
169 {
170 throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
171 }
172 DataLazy_ptr lleft;
173 if (!left->isLazy())
174 {
175 lleft=DataLazy_ptr(new DataLazy(left));
176 }
177 else
178 {
179 lleft=dynamic_pointer_cast<DataLazy>(left);
180 }
181 m_length=left->getLength();
182 m_left=lleft;
183 m_buffsRequired=1;
184 m_samplesize=getNumDPPSample()*getNoValues();
185 }
186
187
188 DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)
189 : parent(resultFS(left,right,op), resultShape(left,right,op)),
190 m_left(left),
191 m_right(right),
192 m_op(op)
193 {
194 if (getOpgroup(op)!=G_BINARY)
195 {
196 throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
197 }
198 m_length=resultLength(m_left,m_right,m_op);
199 m_samplesize=getNumDPPSample()*getNoValues();
200 m_buffsRequired=calcBuffs(m_left, m_right, m_op);
201 cout << "(2)Lazy created with " << m_samplesize << endl;
202 }
203
204 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
205 : parent(resultFS(left,right,op), resultShape(left,right,op)),
206 m_op(op)
207 {
208 if (getOpgroup(op)!=G_BINARY)
209 {
210 throw DataException("Programmer error - constructor DataLazy(left, op) will only process BINARY operations.");
211 }
212 if (left->isLazy())
213 {
214 m_left=dynamic_pointer_cast<DataLazy>(left);
215 }
216 else
217 {
218 m_left=DataLazy_ptr(new DataLazy(left));
219 }
220 if (right->isLazy())
221 {
222 m_right=dynamic_pointer_cast<DataLazy>(right);
223 }
224 else
225 {
226 m_right=DataLazy_ptr(new DataLazy(right));
227 }
228
229 m_length=resultLength(m_left,m_right,m_op);
230 m_samplesize=getNumDPPSample()*getNoValues();
231 m_buffsRequired=calcBuffs(m_left, m_right,m_op);
232 cout << "(3)Lazy created with " << m_samplesize << endl;
233 }
234
235
236 DataLazy::~DataLazy()
237 {
238 }
239
240
241 int
242 DataLazy::getBuffsRequired() const
243 {
244 return m_buffsRequired;
245 }
246
247
248 // the vector and the offset are a place where the method could write its data if it wishes
249 // it is not obligated to do so. For example, if it has its own storage already, it can use that.
250 // Hence the return value to indicate where the data is actually stored.
251 // Regardless, the storage should be assumed to be used, even if it isn't.
252 const double*
253 DataLazy::resolveSample(ValueType& v,int sampleNo, size_t offset ) const
254 {
255 if (m_op==IDENTITY)
256 {
257 const ValueType& vec=m_id->getVector();
258 return &(vec[m_id->getPointOffset(sampleNo, 0)]);
259 }
260 size_t rightoffset=offset+m_samplesize;
261 const double* left=m_left->resolveSample(v,sampleNo,offset);
262 const double* right=0;
263 if (getOpgroup(m_op)==G_BINARY)
264 {
265 right=m_right->resolveSample(v,sampleNo,rightoffset);
266 }
267 double* result=&(v[offset]);
268 {
269 switch(m_op)
270 {
271 case ADD: // since these are pointwise ops, pretend each sample is one point
272 tensor_binary_operation(m_samplesize, left, right, result, plus<double>());
273 break;
274 case SUB:
275 tensor_binary_operation(m_samplesize, left, right, result, minus<double>());
276 break;
277 case MUL:
278 tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());
279 break;
280 case DIV:
281 tensor_binary_operation(m_samplesize, left, right, result, divides<double>());
282 break;
283 // unary ops
284 case SIN:
285 tensor_unary_operation(m_samplesize, left, result, ::sin);
286 break;
287 case COS:
288 tensor_unary_operation(m_samplesize, left, result, ::cos);
289 break;
290 case TAN:
291 tensor_unary_operation(m_samplesize, left, result, ::tan);
292 break;
293 case ASIN:
294 tensor_unary_operation(m_samplesize, left, result, ::asin);
295 break;
296 case ACOS:
297 tensor_unary_operation(m_samplesize, left, result, ::acos);
298 break;
299 case ATAN:
300 tensor_unary_operation(m_samplesize, left, result, ::atan);
301 break;
302 case SINH:
303 tensor_unary_operation(m_samplesize, left, result, ::sinh);
304 break;
305 case COSH:
306 tensor_unary_operation(m_samplesize, left, result, ::cosh);
307 break;
308 case TANH:
309 tensor_unary_operation(m_samplesize, left, result, ::tanh);
310 break;
311 case ERF:
312 #ifdef _WIN32
313 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
314 #else
315 tensor_unary_operation(m_samplesize, left, result, ::erf);
316 break;
317 #endif
318 case ASINH:
319 #ifdef _WIN32
320 tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
321 #else
322 tensor_unary_operation(m_samplesize, left, result, ::asinh);
323 #endif
324 break;
325 case ACOSH:
326 #ifdef _WIN32
327 tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
328 #else
329 tensor_unary_operation(m_samplesize, left, result, ::acosh);
330 #endif
331 break;
332 case ATANH:
333 #ifdef _WIN32
334 tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
335 #else
336 tensor_unary_operation(m_samplesize, left, result, ::atanh);
337 #endif
338 break;
339 case LOG10:
340 tensor_unary_operation(m_samplesize, left, result, ::log10);
341 break;
342 case LOG:
343 tensor_unary_operation(m_samplesize, left, result, ::log);
344 break;
345 case SIGN:
346 tensor_unary_operation(m_samplesize, left, result, escript::fsign);
347 break;
348 case ABS:
349 tensor_unary_operation(m_samplesize, left, result, ::fabs);
350 break;
351 case NEG:
352 tensor_unary_operation(m_samplesize, left, result, negate<double>());
353 break;
354 case POS:
355 // it doesn't mean anything for delayed.
356 // it will just trigger a deep copy of the lazy object
357 throw DataException("Programmer error - POS not supported for lazy data.");
358 break;
359 case EXP:
360 tensor_unary_operation(m_samplesize, left, result, ::exp);
361 break;
362 case SQRT:
363 tensor_unary_operation(m_samplesize, left, result, ::sqrt);
364 break;
365 case RECIP:
366 tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
367 break;
368 case GZ:
369 tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
370 break;
371 case LZ:
372 tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
373 break;
374 case GEZ:
375 tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
376 break;
377 case LEZ:
378 tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
379 break;
380
381 default:
382 throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
383 }
384 }
385 return result;
386 }
387
388 DataReady_ptr
389 DataLazy::resolve()
390 {
391 // This is broken! We need to have a buffer per thread!
392 // so the allocation of v needs to move inside the loop somehow
393
394 cout << "Sample size=" << m_samplesize << endl;
395 cout << "Buffers=" << m_buffsRequired << endl;
396
397 size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);
398 int numthreads=1;
399 #ifdef _OPENMP
400 numthreads=omp_get_max_threads();
401 int threadnum=0;
402 #endif
403 ValueType v(numthreads*threadbuffersize);
404 cout << "Buffer created with size=" << v.size() << endl;
405 ValueType dummy(getNoValues());
406 DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),dummy);
407 ValueType& resvec=result->getVector();
408 DataReady_ptr resptr=DataReady_ptr(result);
409 int sample;
410 int resoffset;
411 int totalsamples=getNumSamples();
412 #pragma omp parallel for private(sample,resoffset,threadnum) schedule(static)
413 for (sample=0;sample<totalsamples;++sample)
414 {
415 #ifdef _OPENMP
416 const double* res=resolveSample(v,sample,threadbuffersize*omp_get_thread_num());
417 #else
418 const double* res=resolveSample(v,sample,0); // this would normally be v, but not if its a single IDENTITY op.
419 #endif
420 resoffset=result->getPointOffset(sample,0);
421 for (int i=0;i<m_samplesize;++i,++resoffset) // copy values into the output vector
422 {
423 resvec[resoffset]=res[i];
424 }
425 }
426 return resptr;
427 }
428
429 std::string
430 DataLazy::toString() const
431 {
432 ostringstream oss;
433 oss << "Lazy Data:";
434 intoString(oss);
435 return oss.str();
436 }
437
438 void
439 DataLazy::intoString(ostringstream& oss) const
440 {
441 switch (getOpgroup(m_op))
442 {
443 case G_IDENTITY:
444 oss << '@' << m_id.get();
445 break;
446 case G_BINARY:
447 oss << '(';
448 m_left->intoString(oss);
449 oss << ' ' << opToString(m_op) << ' ';
450 m_right->intoString(oss);
451 oss << ')';
452 break;
453 case G_UNARY:
454 oss << opToString(m_op) << '(';
455 m_left->intoString(oss);
456 oss << ')';
457 break;
458 default:
459 oss << "UNKNOWN";
460 }
461 }
462
463 // Note that in this case, deepCopy does not make copies of the leaves.
464 // Hopefully copy on write (or whatever we end up using) will take care of this.
465 DataAbstract*
466 DataLazy::deepCopy()
467 {
468 if (m_op==IDENTITY)
469 {
470 return new DataLazy(m_left); // we don't need to copy the child here
471 }
472 return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
473 }
474
475
476 DataTypes::ValueType::size_type
477 DataLazy::getLength() const
478 {
479 return m_length;
480 }
481
482
483 DataAbstract*
484 DataLazy::getSlice(const DataTypes::RegionType& region) const
485 {
486 throw DataException("getSlice - not implemented for Lazy objects.");
487 }
488
489 DataTypes::ValueType::size_type
490 DataLazy::getPointOffset(int sampleNo,
491 int dataPointNo) const
492 {
493 throw DataException("getPointOffset - not implemented for Lazy objects - yet.");
494 }
495
496 } // end namespace

  ViewVC Help
Powered by ViewVC 1.1.26