/[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 2066 - (hide annotations)
Thu Nov 20 05:31:33 2008 UTC (10 years, 4 months ago) by jfenwick
Original Path: trunk/escript/src/DataLazy.cpp
File size: 39119 byte(s)
Fixed Data::toString to look at the amount of data actually stored rather than the number of points in the domain.

Added support for GeneralTensorProduct to LazyData
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 1950 #include "Utils.h"
30 jfenwick 1865
31 jfenwick 1899 /*
32     How does DataLazy work?
33     ~~~~~~~~~~~~~~~~~~~~~~~
34    
35     Each instance represents a single operation on one or two other DataLazy instances. These arguments are normally
36     denoted left and right.
37    
38     A special operation, IDENTITY, stores an instance of DataReady in the m_id member.
39     This means that all "internal" nodes in the structure are instances of DataLazy.
40    
41     Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
42     Note that IDENITY is not considered a unary operation.
43    
44     I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
45     It must however form a DAG (directed acyclic graph).
46     I will refer to individual DataLazy objects with the structure as nodes.
47    
48     Each node also stores:
49     - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
50     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 jfenwick 2037
73     To add a new operator you need to do the following (plus anything I might have forgotten):
74     1) Add to the ES_optype.
75     2) determine what opgroup your operation belongs to (X)
76     3) add a string for the op to the end of ES_opstrings
77     4) increase ES_opcount
78     5) add an entry (X) to opgroups
79     6) add an entry to the switch in collapseToReady
80     7) add an entry to resolveX
81 jfenwick 1899 */
82    
83    
84 jfenwick 1865 using namespace std;
85 jfenwick 1868 using namespace boost;
86 jfenwick 1865
87     namespace escript
88     {
89    
90     namespace
91     {
92    
93 jfenwick 1886 enum ES_opgroup
94     {
95     G_UNKNOWN,
96     G_IDENTITY,
97 jfenwick 1899 G_BINARY, // pointwise operations with two arguments
98 jfenwick 2037 G_UNARY, // pointwise operations with one argument
99 jfenwick 2066 G_NP1OUT, // non-pointwise op with one output
100     G_TENSORPROD // general tensor product
101 jfenwick 1886 };
102    
103    
104    
105    
106 jfenwick 1910 string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
107     "sin","cos","tan",
108 jfenwick 1886 "asin","acos","atan","sinh","cosh","tanh","erf",
109     "asinh","acosh","atanh",
110 jfenwick 1888 "log10","log","sign","abs","neg","pos","exp","sqrt",
111 jfenwick 2037 "1/","where>0","where<0","where>=0","where<=0",
112 jfenwick 2066 "symmetric","nonsymmetric",
113     "prod"};
114     int ES_opcount=36;
115 jfenwick 1910 ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
116     G_UNARY,G_UNARY,G_UNARY, //10
117     G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 17
118     G_UNARY,G_UNARY,G_UNARY, // 20
119     G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 28
120 jfenwick 2037 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 33
121 jfenwick 2066 G_NP1OUT,G_NP1OUT,
122     G_TENSORPROD};
123 jfenwick 1886 inline
124     ES_opgroup
125     getOpgroup(ES_optype op)
126     {
127     return opgroups[op];
128     }
129    
130 jfenwick 1865 // return the FunctionSpace of the result of "left op right"
131     FunctionSpace
132     resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
133     {
134     // perhaps this should call interpolate and throw or something?
135     // maybe we need an interpolate node -
136     // that way, if interpolate is required in any other op we can just throw a
137     // programming error exception.
138 jfenwick 1879
139 jfenwick 1943 FunctionSpace l=left->getFunctionSpace();
140     FunctionSpace r=right->getFunctionSpace();
141     if (l!=r)
142     {
143     if (r.probeInterpolation(l))
144     {
145     return l;
146     }
147     if (l.probeInterpolation(r))
148     {
149     return r;
150     }
151     throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
152     }
153     return l;
154 jfenwick 1865 }
155    
156     // return the shape of the result of "left op right"
157 jfenwick 2066 // the shapes resulting from tensor product are more complex to compute so are worked out elsewhere
158 jfenwick 1865 DataTypes::ShapeType
159     resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
160     {
161 jfenwick 1879 if (left->getShape()!=right->getShape())
162     {
163 jfenwick 2037 if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
164 jfenwick 1908 {
165     throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
166     }
167     if (left->getRank()==0) // we need to allow scalar * anything
168     {
169     return right->getShape();
170     }
171     if (right->getRank()==0)
172     {
173     return left->getShape();
174     }
175     throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
176 jfenwick 1879 }
177     return left->getShape();
178 jfenwick 1865 }
179    
180 jfenwick 2066 // determine the output shape for the general tensor product operation
181     // the additional parameters return information required later for the product
182     // the majority of this code is copy pasted from C_General_Tensor_Product
183     DataTypes::ShapeType
184     GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
185 jfenwick 1865 {
186 jfenwick 2066
187     // Get rank and shape of inputs
188     int rank0 = left->getRank();
189     int rank1 = right->getRank();
190     const DataTypes::ShapeType& shape0 = left->getShape();
191     const DataTypes::ShapeType& shape1 = right->getShape();
192    
193     // Prepare for the loops of the product and verify compatibility of shapes
194     int start0=0, start1=0;
195     if (transpose == 0) {}
196     else if (transpose == 1) { start0 = axis_offset; }
197     else if (transpose == 2) { start1 = rank1-axis_offset; }
198     else { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
199    
200    
201     // Adjust the shapes for transpose
202     DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
203     DataTypes::ShapeType tmpShape1(rank1); // than using push_back
204     for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
205     for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
206    
207     // Prepare for the loops of the product
208     SL=1, SM=1, SR=1;
209     for (int i=0; i<rank0-axis_offset; i++) {
210     SL *= tmpShape0[i];
211     }
212     for (int i=rank0-axis_offset; i<rank0; i++) {
213     if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
214     throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
215     }
216     SM *= tmpShape0[i];
217     }
218     for (int i=axis_offset; i<rank1; i++) {
219     SR *= tmpShape1[i];
220     }
221    
222     // Define the shape of the output (rank of shape is the sum of the loop ranges below)
223     DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
224     { // block to limit the scope of out_index
225     int out_index=0;
226     for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
227     for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
228     }
229     return shape2;
230 jfenwick 1865 }
231    
232 jfenwick 2066
233     // determine the number of points in the result of "left op right"
234     // note that determining the resultLength for G_TENSORPROD is more complex and will not be processed here
235     // size_t
236     // resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
237     // {
238     // switch (getOpgroup(op))
239     // {
240     // case G_BINARY: return left->getLength();
241     // case G_UNARY: return left->getLength();
242     // case G_NP1OUT: return left->getLength();
243     // default:
244     // throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
245     // }
246     // }
247    
248 jfenwick 1899 // determine the number of samples requires to evaluate an expression combining left and right
249 jfenwick 2037 // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
250 jfenwick 2066 // The same goes for G_TENSORPROD
251 jfenwick 1879 int
252     calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
253     {
254 jfenwick 1886 switch(getOpgroup(op))
255 jfenwick 1879 {
256 jfenwick 1886 case G_IDENTITY: return 1;
257     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
258     case G_UNARY: return max(left->getBuffsRequired(),1);
259 jfenwick 2037 case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
260 jfenwick 2066 case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
261 jfenwick 1879 default:
262     throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
263     }
264     }
265    
266 jfenwick 1899
267 jfenwick 1865 } // end anonymous namespace
268    
269    
270 jfenwick 1899
271     // Return a string representing the operation
272 jfenwick 1865 const std::string&
273     opToString(ES_optype op)
274     {
275     if (op<0 || op>=ES_opcount)
276     {
277     op=UNKNOWNOP;
278     }
279     return ES_opstrings[op];
280     }
281    
282    
283     DataLazy::DataLazy(DataAbstract_ptr p)
284     : parent(p->getFunctionSpace(),p->getShape()),
285 jfenwick 2066 m_op(IDENTITY),
286     m_axis_offset(0),
287     m_transpose(0),
288     m_SL(0), m_SM(0), m_SR(0)
289 jfenwick 1865 {
290 jfenwick 1879 if (p->isLazy())
291     {
292     // I don't want identity of Lazy.
293     // Question: Why would that be so bad?
294     // Answer: We assume that the child of ID is something we can call getVector on
295     throw DataException("Programmer error - attempt to create identity from a DataLazy.");
296     }
297     else
298     {
299     m_id=dynamic_pointer_cast<DataReady>(p);
300 jfenwick 1889 if(p->isConstant()) {m_readytype='C';}
301     else if(p->isExpanded()) {m_readytype='E';}
302     else if (p->isTagged()) {m_readytype='T';}
303     else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}
304 jfenwick 1879 }
305 jfenwick 1886 m_buffsRequired=1;
306 jfenwick 1879 m_samplesize=getNumDPPSample()*getNoValues();
307 jfenwick 2066 m_maxsamplesize=m_samplesize;
308 jfenwick 1879 cout << "(1)Lazy created with " << m_samplesize << endl;
309 jfenwick 1865 }
310    
311 jfenwick 1899
312    
313 jfenwick 1901
314 jfenwick 1886 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
315     : parent(left->getFunctionSpace(),left->getShape()),
316 jfenwick 2066 m_op(op),
317     m_axis_offset(0),
318     m_transpose(0),
319     m_SL(0), m_SM(0), m_SR(0)
320 jfenwick 1886 {
321 jfenwick 2037 if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
322 jfenwick 1886 {
323     throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
324     }
325 jfenwick 2066
326 jfenwick 1886 DataLazy_ptr lleft;
327     if (!left->isLazy())
328     {
329     lleft=DataLazy_ptr(new DataLazy(left));
330     }
331     else
332     {
333     lleft=dynamic_pointer_cast<DataLazy>(left);
334     }
335 jfenwick 1889 m_readytype=lleft->m_readytype;
336 jfenwick 1886 m_left=lleft;
337 jfenwick 2037 m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
338 jfenwick 1886 m_samplesize=getNumDPPSample()*getNoValues();
339 jfenwick 2066 m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
340 jfenwick 1886 }
341    
342    
343 jfenwick 1943 // In this constructor we need to consider interpolation
344 jfenwick 1879 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
345     : parent(resultFS(left,right,op), resultShape(left,right,op)),
346 jfenwick 2066 m_op(op),
347     m_SL(0), m_SM(0), m_SR(0)
348 jfenwick 1879 {
349 jfenwick 2037 if ((getOpgroup(op)!=G_BINARY))
350 jfenwick 1886 {
351 jfenwick 1889 throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
352 jfenwick 1886 }
353 jfenwick 1943
354     if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
355     {
356     FunctionSpace fs=getFunctionSpace();
357     Data ltemp(left);
358     Data tmp(ltemp,fs);
359     left=tmp.borrowDataPtr();
360     }
361 jfenwick 2066 if (getFunctionSpace()!=right->getFunctionSpace()) // right needs to be interpolated
362 jfenwick 1943 {
363     Data tmp(Data(right),getFunctionSpace());
364     right=tmp.borrowDataPtr();
365     }
366     left->operandCheck(*right);
367    
368 jfenwick 1899 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
369 jfenwick 1879 {
370     m_left=dynamic_pointer_cast<DataLazy>(left);
371     }
372     else
373     {
374     m_left=DataLazy_ptr(new DataLazy(left));
375     }
376     if (right->isLazy())
377     {
378     m_right=dynamic_pointer_cast<DataLazy>(right);
379     }
380     else
381     {
382     m_right=DataLazy_ptr(new DataLazy(right));
383     }
384 jfenwick 1889 char lt=m_left->m_readytype;
385     char rt=m_right->m_readytype;
386     if (lt=='E' || rt=='E')
387     {
388     m_readytype='E';
389     }
390     else if (lt=='T' || rt=='T')
391     {
392     m_readytype='T';
393     }
394     else
395     {
396     m_readytype='C';
397     }
398 jfenwick 2066 m_samplesize=getNumDPPSample()*getNoValues();
399     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());
400 jfenwick 1879 m_buffsRequired=calcBuffs(m_left, m_right,m_op);
401     cout << "(3)Lazy created with " << m_samplesize << endl;
402     }
403    
404 jfenwick 2066 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
405     : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
406     m_op(op),
407     m_axis_offset(axis_offset),
408     m_transpose(transpose)
409     {
410     if ((getOpgroup(op)!=G_TENSORPROD))
411     {
412     throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
413     }
414     if ((transpose>2) || (transpose<0))
415     {
416     throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
417     }
418     if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
419     {
420     FunctionSpace fs=getFunctionSpace();
421     Data ltemp(left);
422     Data tmp(ltemp,fs);
423     left=tmp.borrowDataPtr();
424     }
425     if (getFunctionSpace()!=right->getFunctionSpace()) // right needs to be interpolated
426     {
427     Data tmp(Data(right),getFunctionSpace());
428     right=tmp.borrowDataPtr();
429     }
430     left->operandCheck(*right);
431 jfenwick 1879
432 jfenwick 2066 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
433     {
434     m_left=dynamic_pointer_cast<DataLazy>(left);
435     }
436     else
437     {
438     m_left=DataLazy_ptr(new DataLazy(left));
439     }
440     if (right->isLazy())
441     {
442     m_right=dynamic_pointer_cast<DataLazy>(right);
443     }
444     else
445     {
446     m_right=DataLazy_ptr(new DataLazy(right));
447     }
448     char lt=m_left->m_readytype;
449     char rt=m_right->m_readytype;
450     if (lt=='E' || rt=='E')
451     {
452     m_readytype='E';
453     }
454     else if (lt=='T' || rt=='T')
455     {
456     m_readytype='T';
457     }
458     else
459     {
460     m_readytype='C';
461     }
462     m_samplesize=getNumDPPSample()*getNoValues();
463     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());
464     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
465     cout << "(4)Lazy created with " << m_samplesize << endl;
466     }
467    
468    
469 jfenwick 1865 DataLazy::~DataLazy()
470     {
471     }
472    
473 jfenwick 1879
474     int
475     DataLazy::getBuffsRequired() const
476 jfenwick 1865 {
477 jfenwick 1879 return m_buffsRequired;
478     }
479    
480 jfenwick 1884
481 jfenwick 2066 size_t
482     DataLazy::getMaxSampleSize() const
483     {
484     return m_maxsamplesize;
485     }
486    
487 jfenwick 1899 /*
488     \brief Evaluates the expression using methods on Data.
489     This does the work for the collapse method.
490     For reasons of efficiency do not call this method on DataExpanded nodes.
491     */
492 jfenwick 1889 DataReady_ptr
493     DataLazy::collapseToReady()
494     {
495     if (m_readytype=='E')
496     { // this is more an efficiency concern than anything else
497     throw DataException("Programmer Error - do not use collapse on Expanded data.");
498     }
499     if (m_op==IDENTITY)
500     {
501     return m_id;
502     }
503     DataReady_ptr pleft=m_left->collapseToReady();
504     Data left(pleft);
505     Data right;
506 jfenwick 2066 if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
507 jfenwick 1889 {
508     right=Data(m_right->collapseToReady());
509     }
510     Data result;
511     switch(m_op)
512     {
513     case ADD:
514     result=left+right;
515     break;
516     case SUB:
517     result=left-right;
518     break;
519     case MUL:
520     result=left*right;
521     break;
522     case DIV:
523     result=left/right;
524     break;
525     case SIN:
526     result=left.sin();
527     break;
528     case COS:
529     result=left.cos();
530     break;
531     case TAN:
532     result=left.tan();
533     break;
534     case ASIN:
535     result=left.asin();
536     break;
537     case ACOS:
538     result=left.acos();
539     break;
540     case ATAN:
541     result=left.atan();
542     break;
543     case SINH:
544     result=left.sinh();
545     break;
546     case COSH:
547     result=left.cosh();
548     break;
549     case TANH:
550     result=left.tanh();
551     break;
552     case ERF:
553     result=left.erf();
554     break;
555     case ASINH:
556     result=left.asinh();
557     break;
558     case ACOSH:
559     result=left.acosh();
560     break;
561     case ATANH:
562     result=left.atanh();
563     break;
564     case LOG10:
565     result=left.log10();
566     break;
567     case LOG:
568     result=left.log();
569     break;
570     case SIGN:
571     result=left.sign();
572     break;
573     case ABS:
574     result=left.abs();
575     break;
576     case NEG:
577     result=left.neg();
578     break;
579     case POS:
580     // it doesn't mean anything for delayed.
581     // it will just trigger a deep copy of the lazy object
582     throw DataException("Programmer error - POS not supported for lazy data.");
583     break;
584     case EXP:
585     result=left.exp();
586     break;
587     case SQRT:
588     result=left.sqrt();
589     break;
590     case RECIP:
591     result=left.oneOver();
592     break;
593     case GZ:
594     result=left.wherePositive();
595     break;
596     case LZ:
597     result=left.whereNegative();
598     break;
599     case GEZ:
600     result=left.whereNonNegative();
601     break;
602     case LEZ:
603     result=left.whereNonPositive();
604     break;
605 jfenwick 2037 case SYM:
606     result=left.symmetric();
607     break;
608     case NSYM:
609     result=left.nonsymmetric();
610     break;
611 jfenwick 2066 case PROD:
612     result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
613     break;
614 jfenwick 1889 default:
615 jfenwick 2037 throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
616 jfenwick 1889 }
617     return result.borrowReadyPtr();
618     }
619    
620 jfenwick 1899 /*
621     \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
622     This method uses the original methods on the Data class to evaluate the expressions.
623     For this reason, it should not be used on DataExpanded instances. (To do so would defeat
624     the purpose of using DataLazy in the first place).
625     */
626 jfenwick 1889 void
627     DataLazy::collapse()
628     {
629     if (m_op==IDENTITY)
630     {
631     return;
632     }
633     if (m_readytype=='E')
634     { // this is more an efficiency concern than anything else
635     throw DataException("Programmer Error - do not use collapse on Expanded data.");
636     }
637     m_id=collapseToReady();
638     m_op=IDENTITY;
639     }
640    
641 jfenwick 1899 /*
642 jfenwick 2037 \brief Compute the value of the expression (unary operation) for the given sample.
643 jfenwick 1899 \return Vector which stores the value of the subexpression for the given sample.
644     \param v A vector to store intermediate results.
645     \param offset Index in v to begin storing results.
646     \param sampleNo Sample number to evaluate.
647     \param roffset (output parameter) the offset in the return vector where the result begins.
648    
649     The return value will be an existing vector so do not deallocate it.
650     If the result is stored in v it should be stored at the offset given.
651     Everything from offset to the end of v should be considered available for this method to use.
652     */
653 jfenwick 1898 DataTypes::ValueType*
654     DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
655 jfenwick 1889 {
656     // we assume that any collapsing has been done before we get here
657     // since we only have one argument we don't need to think about only
658     // processing single points.
659     if (m_readytype!='E')
660     {
661     throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
662     }
663 jfenwick 1898 const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
664     const double* left=&((*vleft)[roffset]);
665 jfenwick 1889 double* result=&(v[offset]);
666 jfenwick 1898 roffset=offset;
667 jfenwick 1889 switch (m_op)
668     {
669     case SIN:
670 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
671 jfenwick 1889 break;
672     case COS:
673 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
674 jfenwick 1889 break;
675     case TAN:
676 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
677 jfenwick 1889 break;
678     case ASIN:
679 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
680 jfenwick 1889 break;
681     case ACOS:
682 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
683 jfenwick 1889 break;
684     case ATAN:
685 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
686 jfenwick 1889 break;
687     case SINH:
688 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
689 jfenwick 1889 break;
690     case COSH:
691 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
692 jfenwick 1889 break;
693     case TANH:
694 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
695 jfenwick 1889 break;
696     case ERF:
697 phornby 2050 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
698 jfenwick 1889 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
699     #else
700     tensor_unary_operation(m_samplesize, left, result, ::erf);
701     break;
702     #endif
703     case ASINH:
704 phornby 2050 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
705 jfenwick 1889 tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
706     #else
707     tensor_unary_operation(m_samplesize, left, result, ::asinh);
708     #endif
709     break;
710     case ACOSH:
711 phornby 2050 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
712 jfenwick 1889 tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
713     #else
714     tensor_unary_operation(m_samplesize, left, result, ::acosh);
715     #endif
716     break;
717     case ATANH:
718 phornby 2050 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
719 jfenwick 1889 tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
720     #else
721     tensor_unary_operation(m_samplesize, left, result, ::atanh);
722     #endif
723     break;
724     case LOG10:
725 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
726 jfenwick 1889 break;
727     case LOG:
728 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
729 jfenwick 1889 break;
730     case SIGN:
731     tensor_unary_operation(m_samplesize, left, result, escript::fsign);
732     break;
733     case ABS:
734 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
735 jfenwick 1889 break;
736     case NEG:
737     tensor_unary_operation(m_samplesize, left, result, negate<double>());
738     break;
739     case POS:
740     // it doesn't mean anything for delayed.
741     // it will just trigger a deep copy of the lazy object
742     throw DataException("Programmer error - POS not supported for lazy data.");
743     break;
744     case EXP:
745 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
746 jfenwick 1889 break;
747     case SQRT:
748 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
749 jfenwick 1889 break;
750     case RECIP:
751     tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
752     break;
753     case GZ:
754     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
755     break;
756     case LZ:
757     tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
758     break;
759     case GEZ:
760     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
761     break;
762     case LEZ:
763     tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
764     break;
765    
766     default:
767     throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
768     }
769 jfenwick 1898 return &v;
770 jfenwick 1889 }
771    
772    
773 jfenwick 2037 /*
774     \brief Compute the value of the expression (unary operation) for the given sample.
775     \return Vector which stores the value of the subexpression for the given sample.
776     \param v A vector to store intermediate results.
777     \param offset Index in v to begin storing results.
778     \param sampleNo Sample number to evaluate.
779     \param roffset (output parameter) the offset in the return vector where the result begins.
780 jfenwick 1898
781 jfenwick 2037 The return value will be an existing vector so do not deallocate it.
782     If the result is stored in v it should be stored at the offset given.
783     Everything from offset to the end of v should be considered available for this method to use.
784     */
785     DataTypes::ValueType*
786     DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
787     {
788     // we assume that any collapsing has been done before we get here
789     // since we only have one argument we don't need to think about only
790     // processing single points.
791     if (m_readytype!='E')
792     {
793     throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
794     }
795     // since we can't write the result over the input, we need a result offset further along
796     size_t subroffset=roffset+m_samplesize;
797     const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);
798     roffset=offset;
799     switch (m_op)
800     {
801     case SYM:
802     DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
803     break;
804     case NSYM:
805     DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
806     break;
807     default:
808     throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
809     }
810     return &v;
811     }
812 jfenwick 1898
813 jfenwick 1899
814 jfenwick 2037
815    
816 phornby 1987 #define PROC_OP(TYPE,X) \
817 jfenwick 1908 for (int i=0;i<steps;++i,resultp+=resultStep) \
818 jfenwick 1889 { \
819 phornby 1994 tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
820 jfenwick 1899 lroffset+=leftStep; \
821     rroffset+=rightStep; \
822 jfenwick 1889 }
823    
824 jfenwick 1899 /*
825     \brief Compute the value of the expression (binary operation) for the given sample.
826     \return Vector which stores the value of the subexpression for the given sample.
827     \param v A vector to store intermediate results.
828     \param offset Index in v to begin storing results.
829     \param sampleNo Sample number to evaluate.
830     \param roffset (output parameter) the offset in the return vector where the result begins.
831    
832     The return value will be an existing vector so do not deallocate it.
833     If the result is stored in v it should be stored at the offset given.
834     Everything from offset to the end of v should be considered available for this method to use.
835     */
836     // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
837     // have already been collapsed to IDENTITY. So we must have at least one expanded child.
838     // If both children are expanded, then we can process them in a single operation (we treat
839     // the whole sample as one big datapoint.
840     // If one of the children is not expanded, then we need to treat each point in the sample
841     // individually.
842 jfenwick 1908 // There is an additional complication when scalar operations are considered.
843     // For example, 2+Vector.
844     // In this case each double within the point is treated individually
845 jfenwick 1898 DataTypes::ValueType*
846 jfenwick 1899 DataLazy::resolveBinary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
847 jfenwick 1889 {
848     cout << "Resolve binary: " << toString() << endl;
849    
850 jfenwick 1899 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
851     // first work out which of the children are expanded
852 jfenwick 1889 bool leftExp=(m_left->m_readytype=='E');
853     bool rightExp=(m_right->m_readytype=='E');
854 jfenwick 1899 bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
855 jfenwick 1908 int steps=(bigloops?1:getNumDPPSample());
856 jfenwick 1899 size_t chunksize=(bigloops? m_samplesize : getNoValues()); // if bigloops, pretend the whole sample is a datapoint
857 jfenwick 1908 if (m_left->getRank()!=m_right->getRank()) // need to deal with scalar * ? ops
858     {
859     EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");
860     steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());
861     chunksize=1; // for scalar
862     }
863     int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
864     int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
865     int resultStep=max(leftStep,rightStep); // only one (at most) should be !=0
866 jfenwick 1899 // Get the values of sub-expressions
867 jfenwick 1898 const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
868 jfenwick 1899 const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
869     // the right child starts further along.
870     double* resultp=&(v[offset]); // results are stored at the vector offset we recieved
871 jfenwick 1889 switch(m_op)
872     {
873     case ADD:
874 phornby 2021 PROC_OP(NO_ARG,plus<double>());
875 jfenwick 1889 break;
876 jfenwick 1899 case SUB:
877 phornby 2021 PROC_OP(NO_ARG,minus<double>());
878 jfenwick 1899 break;
879     case MUL:
880 phornby 2021 PROC_OP(NO_ARG,multiplies<double>());
881 jfenwick 1899 break;
882     case DIV:
883 phornby 2021 PROC_OP(NO_ARG,divides<double>());
884 jfenwick 1899 break;
885 jfenwick 1910 case POW:
886 phornby 1994 PROC_OP(double (double,double),::pow);
887 jfenwick 1910 break;
888 jfenwick 1889 default:
889 jfenwick 1898 throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
890 jfenwick 1889 }
891 jfenwick 1899 roffset=offset;
892 jfenwick 1898 return &v;
893 jfenwick 1889 }
894    
895 jfenwick 1898
896 jfenwick 2066 /*
897     \brief Compute the value of the expression (tensor product) for the given sample.
898     \return Vector which stores the value of the subexpression for the given sample.
899     \param v A vector to store intermediate results.
900     \param offset Index in v to begin storing results.
901     \param sampleNo Sample number to evaluate.
902     \param roffset (output parameter) the offset in the return vector where the result begins.
903 jfenwick 1898
904 jfenwick 2066 The return value will be an existing vector so do not deallocate it.
905     If the result is stored in v it should be stored at the offset given.
906     Everything from offset to the end of v should be considered available for this method to use.
907     */
908     // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
909     // have already been collapsed to IDENTITY. So we must have at least one expanded child.
910     // unlike the other resolve helpers, we must treat these datapoints separately.
911     DataTypes::ValueType*
912     DataLazy::resolveTProd(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
913     {
914     cout << "Resolve TensorProduct: " << toString() << endl;
915    
916     size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
917     // first work out which of the children are expanded
918     bool leftExp=(m_left->m_readytype=='E');
919     bool rightExp=(m_right->m_readytype=='E');
920     int steps=getNumDPPSample();
921     int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
922     int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
923     int resultStep=max(leftStep,rightStep); // only one (at most) should be !=0
924     // Get the values of sub-expressions (leave a gap of one sample for the result).
925     const ValueType* left=m_left->resolveSample(v,offset+m_samplesize,sampleNo,lroffset);
926     const ValueType* right=m_right->resolveSample(v,offset+2*m_samplesize,sampleNo,rroffset);
927     double* resultp=&(v[offset]); // results are stored at the vector offset we recieved
928     switch(m_op)
929     {
930     case PROD:
931     for (int i=0;i<steps;++i,resultp+=resultStep)
932     {
933     const double *ptr_0 = &((*left)[lroffset]);
934     const double *ptr_1 = &((*right)[rroffset]);
935     matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
936     lroffset+=leftStep;
937     rroffset+=rightStep;
938     }
939     break;
940     default:
941     throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
942     }
943     roffset=offset;
944     return &v;
945     }
946    
947    
948    
949 jfenwick 1899 /*
950     \brief Compute the value of the expression for the given sample.
951     \return Vector which stores the value of the subexpression for the given sample.
952     \param v A vector to store intermediate results.
953     \param offset Index in v to begin storing results.
954     \param sampleNo Sample number to evaluate.
955     \param roffset (output parameter) the offset in the return vector where the result begins.
956 jfenwick 1898
957 jfenwick 1899 The return value will be an existing vector so do not deallocate it.
958     */
959 jfenwick 1884 // the vector and the offset are a place where the method could write its data if it wishes
960     // it is not obligated to do so. For example, if it has its own storage already, it can use that.
961 jfenwick 1885 // Hence the return value to indicate where the data is actually stored.
962     // Regardless, the storage should be assumed to be used, even if it isn't.
963 jfenwick 1898
964     // the roffset is the offset within the returned vector where the data begins
965     const DataTypes::ValueType*
966     DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
967 jfenwick 1879 {
968 jfenwick 1889 cout << "Resolve sample " << toString() << endl;
969     // collapse so we have a 'E' node or an IDENTITY for some other type
970     if (m_readytype!='E' && m_op!=IDENTITY)
971     {
972     collapse();
973     }
974 jfenwick 1885 if (m_op==IDENTITY)
975 jfenwick 1865 {
976 jfenwick 1879 const ValueType& vec=m_id->getVector();
977 jfenwick 1889 if (m_readytype=='C')
978     {
979 jfenwick 1898 roffset=0;
980     return &(vec);
981 jfenwick 1889 }
982 jfenwick 1898 roffset=m_id->getPointOffset(sampleNo, 0);
983     return &(vec);
984 jfenwick 1865 }
985 jfenwick 1889 if (m_readytype!='E')
986     {
987     throw DataException("Programmer Error - Collapse did not produce an expanded node.");
988     }
989     switch (getOpgroup(m_op))
990     {
991 jfenwick 1898 case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);
992     case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
993 jfenwick 2066 case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
994     case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
995 jfenwick 1889 default:
996     throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
997     }
998     }
999    
1000    
1001 jfenwick 1899 // To simplify the memory management, all threads operate on one large vector, rather than one each.
1002     // Each sample is evaluated independently and copied into the result DataExpanded.
1003 jfenwick 1879 DataReady_ptr
1004     DataLazy::resolve()
1005     {
1006    
1007     cout << "Sample size=" << m_samplesize << endl;
1008     cout << "Buffers=" << m_buffsRequired << endl;
1009    
1010 jfenwick 1899 if (m_readytype!='E') // if the whole sub-expression is Constant or Tagged, then evaluate it normally
1011 jfenwick 1889 {
1012     collapse();
1013     }
1014 jfenwick 1899 if (m_op==IDENTITY) // So a lazy expression of Constant or Tagged data will be returned here.
1015 jfenwick 1889 {
1016     return m_id;
1017     }
1018     // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
1019 jfenwick 2066 size_t threadbuffersize=m_maxsamplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough
1020 jfenwick 1899 // storage to evaluate its expression
1021 jfenwick 1885 int numthreads=1;
1022     #ifdef _OPENMP
1023 jfenwick 1889 numthreads=getNumberOfThreads();
1024 jfenwick 1885 int threadnum=0;
1025     #endif
1026 jfenwick 1886 ValueType v(numthreads*threadbuffersize);
1027 jfenwick 1879 cout << "Buffer created with size=" << v.size() << endl;
1028 jfenwick 1898 DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(), ValueType(getNoValues()));
1029 jfenwick 1879 ValueType& resvec=result->getVector();
1030     DataReady_ptr resptr=DataReady_ptr(result);
1031     int sample;
1032 jfenwick 1898 size_t outoffset; // offset in the output data
1033 jfenwick 1885 int totalsamples=getNumSamples();
1034 jfenwick 1899 const ValueType* res=0; // Vector storing the answer
1035     size_t resoffset=0; // where in the vector to find the answer
1036 jfenwick 1898 #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)
1037 jfenwick 1885 for (sample=0;sample<totalsamples;++sample)
1038 jfenwick 1879 {
1039 jfenwick 1898 cout << "################################# " << sample << endl;
1040 jfenwick 1885 #ifdef _OPENMP
1041 jfenwick 1898 res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
1042 jfenwick 1885 #else
1043 jfenwick 1899 res=resolveSample(v,0,sample,resoffset); // res would normally be v, but not if its a single IDENTITY op.
1044 jfenwick 1885 #endif
1045 jfenwick 1898 cerr << "-------------------------------- " << endl;
1046     outoffset=result->getPointOffset(sample,0);
1047     cerr << "offset=" << outoffset << endl;
1048     for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset) // copy values into the output vector
1049 jfenwick 1879 {
1050 jfenwick 1898 resvec[outoffset]=(*res)[resoffset];
1051 jfenwick 1879 }
1052 jfenwick 1898 cerr << "*********************************" << endl;
1053 jfenwick 1879 }
1054     return resptr;
1055     }
1056    
1057 jfenwick 1865 std::string
1058     DataLazy::toString() const
1059     {
1060 jfenwick 1886 ostringstream oss;
1061     oss << "Lazy Data:";
1062     intoString(oss);
1063     return oss.str();
1064 jfenwick 1865 }
1065    
1066 jfenwick 1899
1067 jfenwick 1886 void
1068     DataLazy::intoString(ostringstream& oss) const
1069     {
1070     switch (getOpgroup(m_op))
1071     {
1072 jfenwick 1889 case G_IDENTITY:
1073     if (m_id->isExpanded())
1074     {
1075     oss << "E";
1076     }
1077     else if (m_id->isTagged())
1078     {
1079     oss << "T";
1080     }
1081     else if (m_id->isConstant())
1082     {
1083     oss << "C";
1084     }
1085     else
1086     {
1087     oss << "?";
1088     }
1089 jfenwick 1886 oss << '@' << m_id.get();
1090     break;
1091     case G_BINARY:
1092     oss << '(';
1093     m_left->intoString(oss);
1094     oss << ' ' << opToString(m_op) << ' ';
1095     m_right->intoString(oss);
1096     oss << ')';
1097     break;
1098     case G_UNARY:
1099 jfenwick 2037 case G_NP1OUT:
1100 jfenwick 1886 oss << opToString(m_op) << '(';
1101     m_left->intoString(oss);
1102     oss << ')';
1103     break;
1104 jfenwick 2066 case G_TENSORPROD:
1105     oss << opToString(m_op) << '(';
1106     m_left->intoString(oss);
1107     oss << ", ";
1108     m_right->intoString(oss);
1109     oss << ')';
1110     break;
1111 jfenwick 1886 default:
1112     oss << "UNKNOWN";
1113     }
1114     }
1115    
1116 jfenwick 1865 DataAbstract*
1117     DataLazy::deepCopy()
1118     {
1119 jfenwick 1935 switch (getOpgroup(m_op))
1120 jfenwick 1865 {
1121 jfenwick 1935 case G_IDENTITY: return new DataLazy(m_id->deepCopy()->getPtr());
1122     case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1123     case G_BINARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1124 jfenwick 2066 case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1125     case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1126 jfenwick 1935 default:
1127     throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1128 jfenwick 1865 }
1129     }
1130    
1131    
1132 jfenwick 2066 // There is no single, natural interpretation of getLength on DataLazy.
1133     // Instances of DataReady can look at the size of their vectors.
1134     // For lazy though, it could be the size the data would be if it were resolved;
1135     // or it could be some function of the lengths of the DataReady instances which
1136     // form part of the expression.
1137     // Rather than have people making assumptions, I have disabled the method.
1138 jfenwick 1865 DataTypes::ValueType::size_type
1139     DataLazy::getLength() const
1140     {
1141 jfenwick 2066 throw DataException("getLength() does not make sense for lazy data.");
1142 jfenwick 1865 }
1143    
1144    
1145     DataAbstract*
1146     DataLazy::getSlice(const DataTypes::RegionType& region) const
1147     {
1148 jfenwick 1879 throw DataException("getSlice - not implemented for Lazy objects.");
1149 jfenwick 1865 }
1150    
1151 jfenwick 1935
1152     // To do this we need to rely on our child nodes
1153 jfenwick 1865 DataTypes::ValueType::size_type
1154     DataLazy::getPointOffset(int sampleNo,
1155 jfenwick 1935 int dataPointNo)
1156     {
1157     if (m_op==IDENTITY)
1158     {
1159     return m_id->getPointOffset(sampleNo,dataPointNo);
1160     }
1161     if (m_readytype!='E')
1162     {
1163     collapse();
1164     return m_id->getPointOffset(sampleNo,dataPointNo);
1165     }
1166     // at this point we do not have an identity node and the expression will be Expanded
1167     // so we only need to know which child to ask
1168     if (m_left->m_readytype=='E')
1169     {
1170     return m_left->getPointOffset(sampleNo,dataPointNo);
1171     }
1172     else
1173     {
1174     return m_right->getPointOffset(sampleNo,dataPointNo);
1175     }
1176     }
1177    
1178     // To do this we need to rely on our child nodes
1179     DataTypes::ValueType::size_type
1180     DataLazy::getPointOffset(int sampleNo,
1181 jfenwick 1865 int dataPointNo) const
1182     {
1183 jfenwick 1935 if (m_op==IDENTITY)
1184     {
1185     return m_id->getPointOffset(sampleNo,dataPointNo);
1186     }
1187     if (m_readytype=='E')
1188     {
1189     // at this point we do not have an identity node and the expression will be Expanded
1190     // so we only need to know which child to ask
1191     if (m_left->m_readytype=='E')
1192     {
1193     return m_left->getPointOffset(sampleNo,dataPointNo);
1194     }
1195     else
1196     {
1197     return m_right->getPointOffset(sampleNo,dataPointNo);
1198     }
1199     }
1200     if (m_readytype=='C')
1201     {
1202     return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
1203     }
1204     throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1205 jfenwick 1865 }
1206    
1207 jfenwick 1901 // It would seem that DataTagged will need to be treated differently since even after setting all tags
1208     // to zero, all the tags from all the DataTags would be in the result.
1209     // However since they all have the same value (0) whether they are there or not should not matter.
1210     // So I have decided that for all types this method will create a constant 0.
1211     // It can be promoted up as required.
1212     // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
1213     // but we can deal with that if it arrises.
1214     void
1215     DataLazy::setToZero()
1216     {
1217     DataTypes::ValueType v(getNoValues(),0);
1218     m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1219     m_op=IDENTITY;
1220     m_right.reset();
1221     m_left.reset();
1222     m_readytype='C';
1223     m_buffsRequired=1;
1224     }
1225    
1226 jfenwick 1865 } // end namespace

  ViewVC Help
Powered by ViewVC 1.1.26