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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2003 - (hide annotations)
Sun Nov 9 23:57:05 2008 UTC (10 years, 10 months ago) by jfenwick
Original Path: branches/schroedinger_upto2002/escript/src/DataLazy.cpp
File size: 28137 byte(s)
Merging from trunk again

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 1950 #include "Utils.h"
30 jfenwick 1865
31 jfenwick 1899 /*
32     How does DataLazy work?
33     ~~~~~~~~~~~~~~~~~~~~~~~
34    
35     Each instance represents a single operation on one or two other DataLazy instances. These arguments are normally
36     denoted left and right.
37    
38     A special operation, IDENTITY, stores an instance of DataReady in the m_id member.
39     This means that all "internal" nodes in the structure are instances of DataLazy.
40    
41     Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
42     Note that IDENITY is not considered a unary operation.
43    
44     I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
45     It must however form a DAG (directed acyclic graph).
46     I will refer to individual DataLazy objects with the structure as nodes.
47    
48     Each node also stores:
49     - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
50     evaluated.
51     - m_length ~ how many values would be stored in the answer if the expression was evaluated.
52     - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to
53     evaluate the expression.
54     - m_samplesize ~ the number of doubles stored in a sample.
55    
56     When a new node is created, the above values are computed based on the values in the child nodes.
57     Eg: if left requires 4 samples and right requires 6 then left+right requires 7 samples.
58    
59     The resolve method, which produces a DataReady from a DataLazy, does the following:
60     1) Create a DataReady to hold the new result.
61     2) Allocate a vector (v) big enough to hold m_buffsrequired samples.
62     3) For each sample, call resolveSample with v, to get its values and copy them into the result object.
63    
64     (In the case of OMP, multiple samples are resolved in parallel so the vector needs to be larger.)
65    
66     resolveSample returns a Vector* and an offset within that vector where the result is stored.
67     Normally, this would be v, but for identity nodes their internal vector is returned instead.
68    
69     The convention that I use, is that the resolve methods should store their results starting at the offset they are passed.
70    
71     For expressions which evaluate to Constant or Tagged, there is a different evaluation method.
72     The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.
73     */
74    
75    
76 jfenwick 1865 using namespace std;
77 jfenwick 1868 using namespace boost;
78 jfenwick 1865
79     namespace escript
80     {
81    
82     namespace
83     {
84    
85 jfenwick 1886 enum ES_opgroup
86     {
87     G_UNKNOWN,
88     G_IDENTITY,
89 jfenwick 1899 G_BINARY, // pointwise operations with two arguments
90     G_UNARY // pointwise operations with one argument
91 jfenwick 1886 };
92    
93    
94    
95    
96 jfenwick 1910 string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
97     "sin","cos","tan",
98 jfenwick 1886 "asin","acos","atan","sinh","cosh","tanh","erf",
99     "asinh","acosh","atanh",
100 jfenwick 1888 "log10","log","sign","abs","neg","pos","exp","sqrt",
101     "1/","where>0","where<0","where>=0","where<=0"};
102 jfenwick 1926 int ES_opcount=33;
103 jfenwick 1910 ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
104     G_UNARY,G_UNARY,G_UNARY, //10
105     G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 17
106     G_UNARY,G_UNARY,G_UNARY, // 20
107     G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 28
108 jfenwick 1888 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};
109 jfenwick 1886 inline
110     ES_opgroup
111     getOpgroup(ES_optype op)
112     {
113     return opgroups[op];
114     }
115    
116 jfenwick 1865 // return the FunctionSpace of the result of "left op right"
117     FunctionSpace
118     resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
119     {
120     // perhaps this should call interpolate and throw or something?
121     // maybe we need an interpolate node -
122     // that way, if interpolate is required in any other op we can just throw a
123     // programming error exception.
124 jfenwick 1879
125 jfenwick 1943 FunctionSpace l=left->getFunctionSpace();
126     FunctionSpace r=right->getFunctionSpace();
127     if (l!=r)
128     {
129     if (r.probeInterpolation(l))
130     {
131     return l;
132     }
133     if (l.probeInterpolation(r))
134     {
135     return r;
136     }
137     throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
138     }
139     return l;
140 jfenwick 1865 }
141    
142     // return the shape of the result of "left op right"
143     DataTypes::ShapeType
144     resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
145     {
146 jfenwick 1879 if (left->getShape()!=right->getShape())
147     {
148 jfenwick 1908 if (getOpgroup(op)!=G_BINARY)
149     {
150     throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
151     }
152     if (left->getRank()==0) // we need to allow scalar * anything
153     {
154     return right->getShape();
155     }
156     if (right->getRank()==0)
157     {
158     return left->getShape();
159     }
160     throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
161 jfenwick 1879 }
162     return left->getShape();
163 jfenwick 1865 }
164    
165 jfenwick 1899 // determine the number of points in the result of "left op right"
166 jfenwick 1865 size_t
167     resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
168     {
169 jfenwick 1886 switch (getOpgroup(op))
170 jfenwick 1865 {
171 jfenwick 1886 case G_BINARY: return left->getLength();
172     case G_UNARY: return left->getLength();
173 jfenwick 1865 default:
174     throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
175     }
176     }
177    
178 jfenwick 1899 // determine the number of samples requires to evaluate an expression combining left and right
179 jfenwick 1879 int
180     calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
181     {
182 jfenwick 1886 switch(getOpgroup(op))
183 jfenwick 1879 {
184 jfenwick 1886 case G_IDENTITY: return 1;
185     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
186     case G_UNARY: return max(left->getBuffsRequired(),1);
187 jfenwick 1879 default:
188     throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
189     }
190     }
191    
192 jfenwick 1899
193 jfenwick 1865 } // end anonymous namespace
194    
195    
196 jfenwick 1899
197     // Return a string representing the operation
198 jfenwick 1865 const std::string&
199     opToString(ES_optype op)
200     {
201     if (op<0 || op>=ES_opcount)
202     {
203     op=UNKNOWNOP;
204     }
205     return ES_opstrings[op];
206     }
207    
208    
209     DataLazy::DataLazy(DataAbstract_ptr p)
210     : parent(p->getFunctionSpace(),p->getShape()),
211     m_op(IDENTITY)
212     {
213 jfenwick 1879 if (p->isLazy())
214     {
215     // I don't want identity of Lazy.
216     // Question: Why would that be so bad?
217     // Answer: We assume that the child of ID is something we can call getVector on
218     throw DataException("Programmer error - attempt to create identity from a DataLazy.");
219     }
220     else
221     {
222     m_id=dynamic_pointer_cast<DataReady>(p);
223 jfenwick 1889 if(p->isConstant()) {m_readytype='C';}
224     else if(p->isExpanded()) {m_readytype='E';}
225     else if (p->isTagged()) {m_readytype='T';}
226     else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}
227 jfenwick 1879 }
228     m_length=p->getLength();
229 jfenwick 1886 m_buffsRequired=1;
230 jfenwick 1879 m_samplesize=getNumDPPSample()*getNoValues();
231     cout << "(1)Lazy created with " << m_samplesize << endl;
232 jfenwick 1865 }
233    
234 jfenwick 1899
235    
236 jfenwick 1901
237 jfenwick 1886 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
238     : parent(left->getFunctionSpace(),left->getShape()),
239     m_op(op)
240     {
241     if (getOpgroup(op)!=G_UNARY)
242     {
243     throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
244     }
245     DataLazy_ptr lleft;
246     if (!left->isLazy())
247     {
248     lleft=DataLazy_ptr(new DataLazy(left));
249     }
250     else
251     {
252     lleft=dynamic_pointer_cast<DataLazy>(left);
253     }
254 jfenwick 1889 m_readytype=lleft->m_readytype;
255 jfenwick 1886 m_length=left->getLength();
256     m_left=lleft;
257     m_buffsRequired=1;
258     m_samplesize=getNumDPPSample()*getNoValues();
259     }
260    
261    
262 jfenwick 1943 // In this constructor we need to consider interpolation
263 jfenwick 1879 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
264     : parent(resultFS(left,right,op), resultShape(left,right,op)),
265     m_op(op)
266     {
267 jfenwick 1886 if (getOpgroup(op)!=G_BINARY)
268     {
269 jfenwick 1889 throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
270 jfenwick 1886 }
271 jfenwick 1943
272     if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
273     {
274     FunctionSpace fs=getFunctionSpace();
275     Data ltemp(left);
276     Data tmp(ltemp,fs);
277     left=tmp.borrowDataPtr();
278     }
279     if (getFunctionSpace()!=right->getFunctionSpace()) // left needs to be interpolated
280     {
281     Data tmp(Data(right),getFunctionSpace());
282     right=tmp.borrowDataPtr();
283     }
284     left->operandCheck(*right);
285    
286 jfenwick 1899 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
287 jfenwick 1879 {
288     m_left=dynamic_pointer_cast<DataLazy>(left);
289     }
290     else
291     {
292     m_left=DataLazy_ptr(new DataLazy(left));
293     }
294     if (right->isLazy())
295     {
296     m_right=dynamic_pointer_cast<DataLazy>(right);
297     }
298     else
299     {
300     m_right=DataLazy_ptr(new DataLazy(right));
301     }
302 jfenwick 1889 char lt=m_left->m_readytype;
303     char rt=m_right->m_readytype;
304     if (lt=='E' || rt=='E')
305     {
306     m_readytype='E';
307     }
308     else if (lt=='T' || rt=='T')
309     {
310     m_readytype='T';
311     }
312     else
313     {
314     m_readytype='C';
315     }
316 jfenwick 1879 m_length=resultLength(m_left,m_right,m_op);
317 jfenwick 1899 m_samplesize=getNumDPPSample()*getNoValues();
318 jfenwick 1879 m_buffsRequired=calcBuffs(m_left, m_right,m_op);
319     cout << "(3)Lazy created with " << m_samplesize << endl;
320     }
321    
322    
323 jfenwick 1865 DataLazy::~DataLazy()
324     {
325     }
326    
327 jfenwick 1879
328     int
329     DataLazy::getBuffsRequired() const
330 jfenwick 1865 {
331 jfenwick 1879 return m_buffsRequired;
332     }
333    
334 jfenwick 1884
335 jfenwick 1899 /*
336     \brief Evaluates the expression using methods on Data.
337     This does the work for the collapse method.
338     For reasons of efficiency do not call this method on DataExpanded nodes.
339     */
340 jfenwick 1889 DataReady_ptr
341     DataLazy::collapseToReady()
342     {
343     if (m_readytype=='E')
344     { // this is more an efficiency concern than anything else
345     throw DataException("Programmer Error - do not use collapse on Expanded data.");
346     }
347     if (m_op==IDENTITY)
348     {
349     return m_id;
350     }
351     DataReady_ptr pleft=m_left->collapseToReady();
352     Data left(pleft);
353     Data right;
354     if (getOpgroup(m_op)==G_BINARY)
355     {
356     right=Data(m_right->collapseToReady());
357     }
358     Data result;
359     switch(m_op)
360     {
361     case ADD:
362     result=left+right;
363     break;
364     case SUB:
365     result=left-right;
366     break;
367     case MUL:
368     result=left*right;
369     break;
370     case DIV:
371     result=left/right;
372     break;
373     case SIN:
374     result=left.sin();
375     break;
376     case COS:
377     result=left.cos();
378     break;
379     case TAN:
380     result=left.tan();
381     break;
382     case ASIN:
383     result=left.asin();
384     break;
385     case ACOS:
386     result=left.acos();
387     break;
388     case ATAN:
389     result=left.atan();
390     break;
391     case SINH:
392     result=left.sinh();
393     break;
394     case COSH:
395     result=left.cosh();
396     break;
397     case TANH:
398     result=left.tanh();
399     break;
400     case ERF:
401     result=left.erf();
402     break;
403     case ASINH:
404     result=left.asinh();
405     break;
406     case ACOSH:
407     result=left.acosh();
408     break;
409     case ATANH:
410     result=left.atanh();
411     break;
412     case LOG10:
413     result=left.log10();
414     break;
415     case LOG:
416     result=left.log();
417     break;
418     case SIGN:
419     result=left.sign();
420     break;
421     case ABS:
422     result=left.abs();
423     break;
424     case NEG:
425     result=left.neg();
426     break;
427     case POS:
428     // it doesn't mean anything for delayed.
429     // it will just trigger a deep copy of the lazy object
430     throw DataException("Programmer error - POS not supported for lazy data.");
431     break;
432     case EXP:
433     result=left.exp();
434     break;
435     case SQRT:
436     result=left.sqrt();
437     break;
438     case RECIP:
439     result=left.oneOver();
440     break;
441     case GZ:
442     result=left.wherePositive();
443     break;
444     case LZ:
445     result=left.whereNegative();
446     break;
447     case GEZ:
448     result=left.whereNonNegative();
449     break;
450     case LEZ:
451     result=left.whereNonPositive();
452     break;
453     default:
454     throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
455     }
456     return result.borrowReadyPtr();
457     }
458    
459 jfenwick 1899 /*
460     \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
461     This method uses the original methods on the Data class to evaluate the expressions.
462     For this reason, it should not be used on DataExpanded instances. (To do so would defeat
463     the purpose of using DataLazy in the first place).
464     */
465 jfenwick 1889 void
466     DataLazy::collapse()
467     {
468     if (m_op==IDENTITY)
469     {
470     return;
471     }
472     if (m_readytype=='E')
473     { // this is more an efficiency concern than anything else
474     throw DataException("Programmer Error - do not use collapse on Expanded data.");
475     }
476     m_id=collapseToReady();
477     m_op=IDENTITY;
478     }
479    
480 jfenwick 1899 /*
481     \brief Compute the value of the expression (binary operation) for the given sample.
482     \return Vector which stores the value of the subexpression for the given sample.
483     \param v A vector to store intermediate results.
484     \param offset Index in v to begin storing results.
485     \param sampleNo Sample number to evaluate.
486     \param roffset (output parameter) the offset in the return vector where the result begins.
487    
488     The return value will be an existing vector so do not deallocate it.
489     If the result is stored in v it should be stored at the offset given.
490     Everything from offset to the end of v should be considered available for this method to use.
491     */
492 jfenwick 1898 DataTypes::ValueType*
493     DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
494 jfenwick 1889 {
495     // we assume that any collapsing has been done before we get here
496     // since we only have one argument we don't need to think about only
497     // processing single points.
498     if (m_readytype!='E')
499     {
500     throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
501     }
502 jfenwick 1898 const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
503     const double* left=&((*vleft)[roffset]);
504 jfenwick 1889 double* result=&(v[offset]);
505 jfenwick 1898 roffset=offset;
506 jfenwick 1889 switch (m_op)
507     {
508     case SIN:
509 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
510 jfenwick 1889 break;
511     case COS:
512 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
513 jfenwick 1889 break;
514     case TAN:
515 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
516 jfenwick 1889 break;
517     case ASIN:
518 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
519 jfenwick 1889 break;
520     case ACOS:
521 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
522 jfenwick 1889 break;
523     case ATAN:
524 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
525 jfenwick 1889 break;
526     case SINH:
527 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
528 jfenwick 1889 break;
529     case COSH:
530 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
531 jfenwick 1889 break;
532     case TANH:
533 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
534 jfenwick 1889 break;
535     case ERF:
536     #ifdef _WIN32
537     throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
538     #else
539     tensor_unary_operation(m_samplesize, left, result, ::erf);
540     break;
541     #endif
542     case ASINH:
543     #ifdef _WIN32
544     tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
545     #else
546     tensor_unary_operation(m_samplesize, left, result, ::asinh);
547     #endif
548     break;
549     case ACOSH:
550     #ifdef _WIN32
551     tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
552     #else
553     tensor_unary_operation(m_samplesize, left, result, ::acosh);
554     #endif
555     break;
556     case ATANH:
557     #ifdef _WIN32
558     tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
559     #else
560     tensor_unary_operation(m_samplesize, left, result, ::atanh);
561     #endif
562     break;
563     case LOG10:
564 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
565 jfenwick 1889 break;
566     case LOG:
567 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
568 jfenwick 1889 break;
569     case SIGN:
570     tensor_unary_operation(m_samplesize, left, result, escript::fsign);
571     break;
572     case ABS:
573 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
574 jfenwick 1889 break;
575     case NEG:
576     tensor_unary_operation(m_samplesize, left, result, negate<double>());
577     break;
578     case POS:
579     // it doesn't mean anything for delayed.
580     // it will just trigger a deep copy of the lazy object
581     throw DataException("Programmer error - POS not supported for lazy data.");
582     break;
583     case EXP:
584 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
585 jfenwick 1889 break;
586     case SQRT:
587 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
588 jfenwick 1889 break;
589     case RECIP:
590     tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
591     break;
592     case GZ:
593     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
594     break;
595     case LZ:
596     tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
597     break;
598     case GEZ:
599     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
600     break;
601     case LEZ:
602     tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
603     break;
604    
605     default:
606     throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
607     }
608 jfenwick 1898 return &v;
609 jfenwick 1889 }
610    
611    
612 jfenwick 1898
613    
614 jfenwick 1899
615 phornby 1987 #define PROC_OP(TYPE,X) \
616 jfenwick 1908 for (int i=0;i<steps;++i,resultp+=resultStep) \
617 jfenwick 1889 { \
618 phornby 1994 tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
619 jfenwick 1899 lroffset+=leftStep; \
620     rroffset+=rightStep; \
621 jfenwick 1889 }
622    
623 jfenwick 1899 /*
624     \brief Compute the value of the expression (binary operation) for the given sample.
625     \return Vector which stores the value of the subexpression for the given sample.
626     \param v A vector to store intermediate results.
627     \param offset Index in v to begin storing results.
628     \param sampleNo Sample number to evaluate.
629     \param roffset (output parameter) the offset in the return vector where the result begins.
630    
631     The return value will be an existing vector so do not deallocate it.
632     If the result is stored in v it should be stored at the offset given.
633     Everything from offset to the end of v should be considered available for this method to use.
634     */
635     // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
636     // have already been collapsed to IDENTITY. So we must have at least one expanded child.
637     // If both children are expanded, then we can process them in a single operation (we treat
638     // the whole sample as one big datapoint.
639     // If one of the children is not expanded, then we need to treat each point in the sample
640     // individually.
641 jfenwick 1908 // There is an additional complication when scalar operations are considered.
642     // For example, 2+Vector.
643     // In this case each double within the point is treated individually
644 jfenwick 1898 DataTypes::ValueType*
645 jfenwick 1899 DataLazy::resolveBinary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
646 jfenwick 1889 {
647     cout << "Resolve binary: " << toString() << endl;
648    
649 jfenwick 1899 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
650     // first work out which of the children are expanded
651 jfenwick 1889 bool leftExp=(m_left->m_readytype=='E');
652     bool rightExp=(m_right->m_readytype=='E');
653 jfenwick 1899 bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
654 jfenwick 1908 int steps=(bigloops?1:getNumDPPSample());
655 jfenwick 1899 size_t chunksize=(bigloops? m_samplesize : getNoValues()); // if bigloops, pretend the whole sample is a datapoint
656 jfenwick 1908 if (m_left->getRank()!=m_right->getRank()) // need to deal with scalar * ? ops
657     {
658     EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");
659     steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());
660     chunksize=1; // for scalar
661     }
662     int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
663     int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
664     int resultStep=max(leftStep,rightStep); // only one (at most) should be !=0
665 jfenwick 1899 // Get the values of sub-expressions
666 jfenwick 1898 const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
667 jfenwick 1899 const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
668     // the right child starts further along.
669     double* resultp=&(v[offset]); // results are stored at the vector offset we recieved
670 jfenwick 1889 switch(m_op)
671     {
672     case ADD:
673 phornby 1987 PROC_OP(/**/,plus<double>());
674 jfenwick 1889 break;
675 jfenwick 1899 case SUB:
676 phornby 1987 PROC_OP(/**/,minus<double>());
677 jfenwick 1899 break;
678     case MUL:
679 phornby 1987 PROC_OP(/**/,multiplies<double>());
680 jfenwick 1899 break;
681     case DIV:
682 phornby 1987 PROC_OP(/**/,divides<double>());
683 jfenwick 1899 break;
684 jfenwick 1910 case POW:
685 phornby 1994 PROC_OP(double (double,double),::pow);
686 jfenwick 1910 break;
687 jfenwick 1889 default:
688 jfenwick 1898 throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
689 jfenwick 1889 }
690 jfenwick 1899 roffset=offset;
691 jfenwick 1898 return &v;
692 jfenwick 1889 }
693    
694 jfenwick 1898
695    
696 jfenwick 1899 /*
697     \brief Compute the value of the expression for the given sample.
698     \return Vector which stores the value of the subexpression for the given sample.
699     \param v A vector to store intermediate results.
700     \param offset Index in v to begin storing results.
701     \param sampleNo Sample number to evaluate.
702     \param roffset (output parameter) the offset in the return vector where the result begins.
703 jfenwick 1898
704 jfenwick 1899 The return value will be an existing vector so do not deallocate it.
705     */
706 jfenwick 1884 // the vector and the offset are a place where the method could write its data if it wishes
707     // it is not obligated to do so. For example, if it has its own storage already, it can use that.
708 jfenwick 1885 // Hence the return value to indicate where the data is actually stored.
709     // Regardless, the storage should be assumed to be used, even if it isn't.
710 jfenwick 1898
711     // the roffset is the offset within the returned vector where the data begins
712     const DataTypes::ValueType*
713     DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
714 jfenwick 1879 {
715 jfenwick 1889 cout << "Resolve sample " << toString() << endl;
716     // collapse so we have a 'E' node or an IDENTITY for some other type
717     if (m_readytype!='E' && m_op!=IDENTITY)
718     {
719     collapse();
720     }
721 jfenwick 1885 if (m_op==IDENTITY)
722 jfenwick 1865 {
723 jfenwick 1879 const ValueType& vec=m_id->getVector();
724 jfenwick 1889 if (m_readytype=='C')
725     {
726 jfenwick 1898 roffset=0;
727     return &(vec);
728 jfenwick 1889 }
729 jfenwick 1898 roffset=m_id->getPointOffset(sampleNo, 0);
730     return &(vec);
731 jfenwick 1865 }
732 jfenwick 1889 if (m_readytype!='E')
733     {
734     throw DataException("Programmer Error - Collapse did not produce an expanded node.");
735     }
736     switch (getOpgroup(m_op))
737     {
738 jfenwick 1898 case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);
739     case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
740 jfenwick 1889 default:
741     throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
742     }
743     }
744    
745    
746 jfenwick 1899 // To simplify the memory management, all threads operate on one large vector, rather than one each.
747     // Each sample is evaluated independently and copied into the result DataExpanded.
748 jfenwick 1879 DataReady_ptr
749     DataLazy::resolve()
750     {
751    
752     cout << "Sample size=" << m_samplesize << endl;
753     cout << "Buffers=" << m_buffsRequired << endl;
754    
755 jfenwick 1899 if (m_readytype!='E') // if the whole sub-expression is Constant or Tagged, then evaluate it normally
756 jfenwick 1889 {
757     collapse();
758     }
759 jfenwick 1899 if (m_op==IDENTITY) // So a lazy expression of Constant or Tagged data will be returned here.
760 jfenwick 1889 {
761     return m_id;
762     }
763     // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
764 jfenwick 1899 size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough
765     // storage to evaluate its expression
766 jfenwick 1885 int numthreads=1;
767     #ifdef _OPENMP
768 jfenwick 1889 numthreads=getNumberOfThreads();
769 jfenwick 1885 int threadnum=0;
770     #endif
771 jfenwick 1886 ValueType v(numthreads*threadbuffersize);
772 jfenwick 1879 cout << "Buffer created with size=" << v.size() << endl;
773 jfenwick 1898 DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(), ValueType(getNoValues()));
774 jfenwick 1879 ValueType& resvec=result->getVector();
775     DataReady_ptr resptr=DataReady_ptr(result);
776     int sample;
777 jfenwick 1898 size_t outoffset; // offset in the output data
778 jfenwick 1885 int totalsamples=getNumSamples();
779 jfenwick 1899 const ValueType* res=0; // Vector storing the answer
780     size_t resoffset=0; // where in the vector to find the answer
781 jfenwick 1898 #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)
782 jfenwick 1885 for (sample=0;sample<totalsamples;++sample)
783 jfenwick 1879 {
784 jfenwick 1898 cout << "################################# " << sample << endl;
785 jfenwick 1885 #ifdef _OPENMP
786 jfenwick 1898 res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
787 jfenwick 1885 #else
788 jfenwick 1899 res=resolveSample(v,0,sample,resoffset); // res would normally be v, but not if its a single IDENTITY op.
789 jfenwick 1885 #endif
790 jfenwick 1898 cerr << "-------------------------------- " << endl;
791     outoffset=result->getPointOffset(sample,0);
792     cerr << "offset=" << outoffset << endl;
793     for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset) // copy values into the output vector
794 jfenwick 1879 {
795 jfenwick 1898 resvec[outoffset]=(*res)[resoffset];
796 jfenwick 1879 }
797 jfenwick 1898 cerr << "*********************************" << endl;
798 jfenwick 1879 }
799     return resptr;
800     }
801    
802 jfenwick 1865 std::string
803     DataLazy::toString() const
804     {
805 jfenwick 1886 ostringstream oss;
806     oss << "Lazy Data:";
807     intoString(oss);
808     return oss.str();
809 jfenwick 1865 }
810    
811 jfenwick 1899
812 jfenwick 1886 void
813     DataLazy::intoString(ostringstream& oss) const
814     {
815     switch (getOpgroup(m_op))
816     {
817 jfenwick 1889 case G_IDENTITY:
818     if (m_id->isExpanded())
819     {
820     oss << "E";
821     }
822     else if (m_id->isTagged())
823     {
824     oss << "T";
825     }
826     else if (m_id->isConstant())
827     {
828     oss << "C";
829     }
830     else
831     {
832     oss << "?";
833     }
834 jfenwick 1886 oss << '@' << m_id.get();
835     break;
836     case G_BINARY:
837     oss << '(';
838     m_left->intoString(oss);
839     oss << ' ' << opToString(m_op) << ' ';
840     m_right->intoString(oss);
841     oss << ')';
842     break;
843     case G_UNARY:
844     oss << opToString(m_op) << '(';
845     m_left->intoString(oss);
846     oss << ')';
847     break;
848     default:
849     oss << "UNKNOWN";
850     }
851     }
852    
853 jfenwick 1865 DataAbstract*
854     DataLazy::deepCopy()
855     {
856 jfenwick 1935 switch (getOpgroup(m_op))
857 jfenwick 1865 {
858 jfenwick 1935 case G_IDENTITY: return new DataLazy(m_id->deepCopy()->getPtr());
859     case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
860     case G_BINARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
861     default:
862     throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
863 jfenwick 1865 }
864     }
865    
866    
867     DataTypes::ValueType::size_type
868     DataLazy::getLength() const
869     {
870 jfenwick 1879 return m_length;
871 jfenwick 1865 }
872    
873    
874     DataAbstract*
875     DataLazy::getSlice(const DataTypes::RegionType& region) const
876     {
877 jfenwick 1879 throw DataException("getSlice - not implemented for Lazy objects.");
878 jfenwick 1865 }
879    
880 jfenwick 1935
881     // To do this we need to rely on our child nodes
882 jfenwick 1865 DataTypes::ValueType::size_type
883     DataLazy::getPointOffset(int sampleNo,
884 jfenwick 1935 int dataPointNo)
885     {
886     if (m_op==IDENTITY)
887     {
888     return m_id->getPointOffset(sampleNo,dataPointNo);
889     }
890     if (m_readytype!='E')
891     {
892     collapse();
893     return m_id->getPointOffset(sampleNo,dataPointNo);
894     }
895     // at this point we do not have an identity node and the expression will be Expanded
896     // so we only need to know which child to ask
897     if (m_left->m_readytype=='E')
898     {
899     return m_left->getPointOffset(sampleNo,dataPointNo);
900     }
901     else
902     {
903     return m_right->getPointOffset(sampleNo,dataPointNo);
904     }
905     }
906    
907     // To do this we need to rely on our child nodes
908     DataTypes::ValueType::size_type
909     DataLazy::getPointOffset(int sampleNo,
910 jfenwick 1865 int dataPointNo) const
911     {
912 jfenwick 1935 if (m_op==IDENTITY)
913     {
914     return m_id->getPointOffset(sampleNo,dataPointNo);
915     }
916     if (m_readytype=='E')
917     {
918     // at this point we do not have an identity node and the expression will be Expanded
919     // so we only need to know which child to ask
920     if (m_left->m_readytype=='E')
921     {
922     return m_left->getPointOffset(sampleNo,dataPointNo);
923     }
924     else
925     {
926     return m_right->getPointOffset(sampleNo,dataPointNo);
927     }
928     }
929     if (m_readytype=='C')
930     {
931     return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
932     }
933     throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
934 jfenwick 1865 }
935    
936 jfenwick 1901 // It would seem that DataTagged will need to be treated differently since even after setting all tags
937     // to zero, all the tags from all the DataTags would be in the result.
938     // However since they all have the same value (0) whether they are there or not should not matter.
939     // So I have decided that for all types this method will create a constant 0.
940     // It can be promoted up as required.
941     // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
942     // but we can deal with that if it arrises.
943     void
944     DataLazy::setToZero()
945     {
946     DataTypes::ValueType v(getNoValues(),0);
947     m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
948     m_op=IDENTITY;
949     m_right.reset();
950     m_left.reset();
951     m_readytype='C';
952     m_buffsRequired=1;
953     }
954    
955 jfenwick 1865 } // end namespace

  ViewVC Help
Powered by ViewVC 1.1.26