/[escript]/branches/schroedinger/escript/src/DataLazy.cpp
ViewVC logotype

Annotation of /branches/schroedinger/escript/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26