/[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 1886 - (show annotations)
Wed Oct 15 01:34:18 2008 UTC (10 years, 6 months ago) by jfenwick
Original Path: branches/schroedinger/escript/src/DataLazy.cpp
File size: 12266 byte(s)
Branch commit.
Added unary ops up to pos.
toString now prints expression.
Added inlines to UnaryFuncs.h.

Still only supporting DataExpanded.

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

  ViewVC Help
Powered by ViewVC 1.1.26