/[escript]/trunk/escript/src/DataLazy.cpp
ViewVC logotype

Annotation of /trunk/escript/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

1 jfenwick 1865
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 jfenwick 1885 #ifdef _OPENMP
23     #include <omp.h>
24     #endif
25 jfenwick 1865 #include "FunctionSpace.h"
26     #include "DataTypes.h"
27 jfenwick 1868 #include "Data.h"
28 jfenwick 1886 #include "UnaryFuncs.h" // for escript::fsign
29 jfenwick 1865
30     using namespace std;
31 jfenwick 1868 using namespace boost;
32 jfenwick 1865
33     namespace escript
34     {
35    
36     const std::string&
37     opToString(ES_optype op);
38    
39     namespace
40     {
41    
42 jfenwick 1886
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 jfenwick 1888 "log10","log","sign","abs","neg","pos","exp","sqrt",
59     "1/","where>0","where<0","where>=0","where<=0"};
60     int ES_opcount=32;
61 jfenwick 1886 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 jfenwick 1888 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 jfenwick 1886 inline
67     ES_opgroup
68     getOpgroup(ES_optype op)
69     {
70     return opgroups[op];
71     }
72    
73 jfenwick 1865 // 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 jfenwick 1879
82    
83     if (left->getFunctionSpace()!=right->getFunctionSpace())
84     {
85     throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");
86     }
87 jfenwick 1865 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 jfenwick 1879 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 jfenwick 1865 }
100    
101     size_t
102     resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
103     {
104 jfenwick 1886 switch (getOpgroup(op))
105 jfenwick 1865 {
106 jfenwick 1886 case G_BINARY: return left->getLength();
107     case G_UNARY: return left->getLength();
108 jfenwick 1865 default:
109     throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
110     }
111     }
112    
113 jfenwick 1879 int
114     calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
115     {
116 jfenwick 1886 switch(getOpgroup(op))
117 jfenwick 1879 {
118 jfenwick 1886 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 jfenwick 1879 default:
122     throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
123     }
124     }
125    
126 jfenwick 1865
127 jfenwick 1879
128 jfenwick 1865 } // 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 jfenwick 1879 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 jfenwick 1886 m_buffsRequired=1;
160 jfenwick 1879 m_samplesize=getNumDPPSample()*getNoValues();
161     cout << "(1)Lazy created with " << m_samplesize << endl;
162 jfenwick 1865 }
163    
164 jfenwick 1886 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 jfenwick 1879 DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)
189 jfenwick 1865 : parent(resultFS(left,right,op), resultShape(left,right,op)),
190     m_left(left),
191     m_right(right),
192     m_op(op)
193     {
194 jfenwick 1886 if (getOpgroup(op)!=G_BINARY)
195     {
196     throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
197     }
198 jfenwick 1879 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 jfenwick 1865 }
203    
204 jfenwick 1879 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 jfenwick 1886 if (getOpgroup(op)!=G_BINARY)
209     {
210     throw DataException("Programmer error - constructor DataLazy(left, op) will only process BINARY operations.");
211     }
212 jfenwick 1879 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 jfenwick 1865 DataLazy::~DataLazy()
237     {
238     }
239    
240 jfenwick 1879
241     int
242     DataLazy::getBuffsRequired() const
243 jfenwick 1865 {
244 jfenwick 1879 return m_buffsRequired;
245     }
246    
247 jfenwick 1884
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 jfenwick 1885 // 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 jfenwick 1884 const double*
253 jfenwick 1879 DataLazy::resolveSample(ValueType& v,int sampleNo, size_t offset ) const
254     {
255 jfenwick 1885 if (m_op==IDENTITY)
256 jfenwick 1865 {
257 jfenwick 1879 const ValueType& vec=m_id->getVector();
258 jfenwick 1884 return &(vec[m_id->getPointOffset(sampleNo, 0)]);
259 jfenwick 1865 }
260 jfenwick 1879 size_t rightoffset=offset+m_samplesize;
261 jfenwick 1884 const double* left=m_left->resolveSample(v,sampleNo,offset);
262 jfenwick 1886 const double* right=0;
263     if (getOpgroup(m_op)==G_BINARY)
264     {
265     right=m_right->resolveSample(v,sampleNo,rightoffset);
266     }
267 jfenwick 1884 double* result=&(v[offset]);
268 jfenwick 1865 {
269 jfenwick 1879 switch(m_op)
270     {
271     case ADD: // since these are pointwise ops, pretend each sample is one point
272 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, plus<double>());
273 jfenwick 1879 break;
274     case SUB:
275 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, minus<double>());
276 jfenwick 1879 break;
277     case MUL:
278 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());
279 jfenwick 1879 break;
280     case DIV:
281 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, divides<double>());
282 jfenwick 1879 break;
283 jfenwick 1886 // 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 jfenwick 1888 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 jfenwick 1868 default:
382     throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
383 jfenwick 1879 }
384 jfenwick 1868 }
385 jfenwick 1884 return result;
386 jfenwick 1865 }
387    
388 jfenwick 1879 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 jfenwick 1885 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 jfenwick 1886 ValueType v(numthreads*threadbuffersize);
404 jfenwick 1879 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 jfenwick 1885 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 jfenwick 1879 {
415 jfenwick 1885 #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 jfenwick 1879 {
423 jfenwick 1885 resvec[resoffset]=res[i];
424 jfenwick 1879 }
425     }
426     return resptr;
427     }
428    
429 jfenwick 1865 std::string
430     DataLazy::toString() const
431     {
432 jfenwick 1886 ostringstream oss;
433     oss << "Lazy Data:";
434     intoString(oss);
435     return oss.str();
436 jfenwick 1865 }
437    
438 jfenwick 1886 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 jfenwick 1865 // 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 jfenwick 1879 return m_length;
480 jfenwick 1865 }
481    
482    
483     DataAbstract*
484     DataLazy::getSlice(const DataTypes::RegionType& region) const
485     {
486 jfenwick 1879 throw DataException("getSlice - not implemented for Lazy objects.");
487 jfenwick 1865 }
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