/[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 1908 - (hide annotations)
Thu Oct 23 01:35:31 2008 UTC (10 years, 5 months ago) by jfenwick
Original Path: branches/schroedinger/escript/src/DataLazy.cpp
File size: 25778 byte(s)
Branch commit.
Fixed to support scalar op non-scalar operations.

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

  ViewVC Help
Powered by ViewVC 1.1.26