/[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 1910 - (hide annotations)
Thu Oct 23 03:05:28 2008 UTC (10 years, 6 months ago) by jfenwick
Original Path: branches/schroedinger/escript/src/DataLazy.cpp
File size: 25840 byte(s)
Branch commit.
Support for ** added.

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

  ViewVC Help
Powered by ViewVC 1.1.26