/[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 1943 - (hide annotations)
Wed Oct 29 04:05:14 2008 UTC (10 years, 5 months ago) by jfenwick
Original Path: branches/schroedinger/escript/src/DataLazy.cpp
File size: 27785 byte(s)
Branch commit
Added interpolation.

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

  ViewVC Help
Powered by ViewVC 1.1.26