/[escript]/branches/schroedinger/escript/src/DataLazy.cpp
ViewVC logotype

Annotation of /branches/schroedinger/escript/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1886 - (hide annotations)
Wed Oct 15 01:34:18 2008 UTC (11 years, 1 month ago) by jfenwick
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 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     "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 jfenwick 1865 // 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 jfenwick 1879
81    
82     if (left->getFunctionSpace()!=right->getFunctionSpace())
83     {
84     throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");
85     }
86 jfenwick 1865 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 jfenwick 1879 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 jfenwick 1865 }
99    
100     size_t
101     resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
102     {
103 jfenwick 1886 switch (getOpgroup(op))
104 jfenwick 1865 {
105 jfenwick 1886 case G_BINARY: return left->getLength();
106     case G_UNARY: return left->getLength();
107 jfenwick 1865 default:
108     throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
109     }
110     }
111    
112 jfenwick 1879 int
113     calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
114     {
115 jfenwick 1886 switch(getOpgroup(op))
116 jfenwick 1879 {
117 jfenwick 1886 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 jfenwick 1879 default:
121     throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
122     }
123     }
124    
125 jfenwick 1865
126 jfenwick 1879
127 jfenwick 1865 } // 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 jfenwick 1879 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 jfenwick 1886 m_buffsRequired=1;
159 jfenwick 1879 m_samplesize=getNumDPPSample()*getNoValues();
160     cout << "(1)Lazy created with " << m_samplesize << endl;
161 jfenwick 1865 }
162    
163 jfenwick 1886 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 jfenwick 1879 DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)
188 jfenwick 1865 : parent(resultFS(left,right,op), resultShape(left,right,op)),
189     m_left(left),
190     m_right(right),
191     m_op(op)
192     {
193 jfenwick 1886 if (getOpgroup(op)!=G_BINARY)
194     {
195     throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
196     }
197 jfenwick 1879 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 jfenwick 1865 }
202    
203 jfenwick 1879 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 jfenwick 1886 if (getOpgroup(op)!=G_BINARY)
208     {
209     throw DataException("Programmer error - constructor DataLazy(left, op) will only process BINARY operations.");
210     }
211 jfenwick 1879 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 jfenwick 1865 DataLazy::~DataLazy()
236     {
237     }
238    
239 jfenwick 1879
240     int
241     DataLazy::getBuffsRequired() const
242 jfenwick 1865 {
243 jfenwick 1879 return m_buffsRequired;
244     }
245    
246 jfenwick 1884
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 jfenwick 1885 // 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 jfenwick 1884 const double*
252 jfenwick 1879 DataLazy::resolveSample(ValueType& v,int sampleNo, size_t offset ) const
253     {
254 jfenwick 1885 if (m_op==IDENTITY)
255 jfenwick 1865 {
256 jfenwick 1879 const ValueType& vec=m_id->getVector();
257 jfenwick 1884 return &(vec[m_id->getPointOffset(sampleNo, 0)]);
258 jfenwick 1865 }
259 jfenwick 1879 size_t rightoffset=offset+m_samplesize;
260 jfenwick 1884 const double* left=m_left->resolveSample(v,sampleNo,offset);
261 jfenwick 1886 const double* right=0;
262     if (getOpgroup(m_op)==G_BINARY)
263     {
264     right=m_right->resolveSample(v,sampleNo,rightoffset);
265     }
266 jfenwick 1884 double* result=&(v[offset]);
267 jfenwick 1865 {
268 jfenwick 1879 switch(m_op)
269     {
270     case ADD: // since these are pointwise ops, pretend each sample is one point
271 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, plus<double>());
272 jfenwick 1879 break;
273     case SUB:
274 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, minus<double>());
275 jfenwick 1879 break;
276     case MUL:
277 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());
278 jfenwick 1879 break;
279     case DIV:
280 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, divides<double>());
281 jfenwick 1879 break;
282 jfenwick 1886 // 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 jfenwick 1868 default:
359     throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
360 jfenwick 1879 }
361 jfenwick 1868 }
362 jfenwick 1884 return result;
363 jfenwick 1865 }
364    
365 jfenwick 1879 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 jfenwick 1885 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 jfenwick 1886 ValueType v(numthreads*threadbuffersize);
381 jfenwick 1879 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 jfenwick 1885 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 jfenwick 1879 {
392 jfenwick 1885 #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 jfenwick 1879 {
400 jfenwick 1885 resvec[resoffset]=res[i];
401 jfenwick 1879 }
402     }
403     return resptr;
404     }
405    
406 jfenwick 1865 std::string
407     DataLazy::toString() const
408     {
409 jfenwick 1886 ostringstream oss;
410     oss << "Lazy Data:";
411     intoString(oss);
412     return oss.str();
413 jfenwick 1865 }
414    
415 jfenwick 1886 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 jfenwick 1865 // 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 jfenwick 1879 return m_length;
457 jfenwick 1865 }
458    
459    
460     DataAbstract*
461     DataLazy::getSlice(const DataTypes::RegionType& region) const
462     {
463 jfenwick 1879 throw DataException("getSlice - not implemented for Lazy objects.");
464 jfenwick 1865 }
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