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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1950 - (hide annotations)
Thu Oct 30 00:59:34 2008 UTC (10 years, 5 months ago) by jfenwick
Original Path: branches/schroedinger_upto1946/escript/src/DataLazy.cpp
File size: 27804 byte(s)
Branch commit
Brought schroedinger merge upto trunk@1946.

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

  ViewVC Help
Powered by ViewVC 1.1.26