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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1899 - (hide annotations)
Mon Oct 20 05:13:24 2008 UTC (10 years, 7 months ago) by jfenwick
Original Path: branches/schroedinger/escript/src/DataLazy.cpp
File size: 24053 byte(s)
Branch commit.
Added some doco to DataLazy.
Made Data::integrate aware of DataLazy.


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

  ViewVC Help
Powered by ViewVC 1.1.26