/[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 2770 - (hide annotations)
Wed Nov 25 01:24:51 2009 UTC (9 years, 5 months ago) by jfenwick
Original Path: trunk/escript/src/DataLazy.cpp
File size: 54086 byte(s)
Removed buffer implementation of Lazy.
Removed supporting Alloc/Free Sample buffer calls.
1 jfenwick 1865
2     /*******************************************************
3     *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 jfenwick 1865 * 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 2472 #include "EscriptParams.h"
32    
33 jfenwick 2199 #include <iomanip> // for some fancy formatting in debug
34    
35 jfenwick 2157 // #define LAZYDEBUG(X) if (privdebug){X;}
36 jfenwick 2092 #define LAZYDEBUG(X)
37 jfenwick 2157 namespace
38     {
39     bool privdebug=false;
40 jfenwick 2092
41 jfenwick 2157 #define ENABLEDEBUG privdebug=true;
42     #define DISABLEDEBUG privdebug=false;
43     }
44    
45 jfenwick 2501 // #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl;resolveToIdentity();}
46 jfenwick 2177
47 jfenwick 2472 #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {resolveToIdentity();}
48    
49 jfenwick 2501
50 jfenwick 1899 /*
51     How does DataLazy work?
52     ~~~~~~~~~~~~~~~~~~~~~~~
53    
54     Each instance represents a single operation on one or two other DataLazy instances. These arguments are normally
55     denoted left and right.
56    
57     A special operation, IDENTITY, stores an instance of DataReady in the m_id member.
58     This means that all "internal" nodes in the structure are instances of DataLazy.
59    
60     Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
61     Note that IDENITY is not considered a unary operation.
62    
63     I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
64     It must however form a DAG (directed acyclic graph).
65     I will refer to individual DataLazy objects with the structure as nodes.
66    
67     Each node also stores:
68     - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
69     evaluated.
70     - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to
71     evaluate the expression.
72     - m_samplesize ~ the number of doubles stored in a sample.
73    
74     When a new node is created, the above values are computed based on the values in the child nodes.
75     Eg: if left requires 4 samples and right requires 6 then left+right requires 7 samples.
76    
77     The resolve method, which produces a DataReady from a DataLazy, does the following:
78     1) Create a DataReady to hold the new result.
79     2) Allocate a vector (v) big enough to hold m_buffsrequired samples.
80     3) For each sample, call resolveSample with v, to get its values and copy them into the result object.
81    
82     (In the case of OMP, multiple samples are resolved in parallel so the vector needs to be larger.)
83    
84     resolveSample returns a Vector* and an offset within that vector where the result is stored.
85     Normally, this would be v, but for identity nodes their internal vector is returned instead.
86    
87     The convention that I use, is that the resolve methods should store their results starting at the offset they are passed.
88    
89     For expressions which evaluate to Constant or Tagged, there is a different evaluation method.
90     The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.
91 jfenwick 2037
92 jfenwick 2147 To add a new operator you need to do the following (plus anything I might have forgotten - adding a new group for example):
93 jfenwick 2037 1) Add to the ES_optype.
94     2) determine what opgroup your operation belongs to (X)
95     3) add a string for the op to the end of ES_opstrings
96     4) increase ES_opcount
97     5) add an entry (X) to opgroups
98     6) add an entry to the switch in collapseToReady
99     7) add an entry to resolveX
100 jfenwick 1899 */
101    
102    
103 jfenwick 1865 using namespace std;
104 jfenwick 1868 using namespace boost;
105 jfenwick 1865
106     namespace escript
107     {
108    
109     namespace
110     {
111    
112 jfenwick 1886 enum ES_opgroup
113     {
114     G_UNKNOWN,
115     G_IDENTITY,
116 jfenwick 1899 G_BINARY, // pointwise operations with two arguments
117 jfenwick 2037 G_UNARY, // pointwise operations with one argument
118 jfenwick 2147 G_UNARY_P, // pointwise operations with one argument, requiring a parameter
119 jfenwick 2066 G_NP1OUT, // non-pointwise op with one output
120 jfenwick 2084 G_NP1OUT_P, // non-pointwise op with one output requiring a parameter
121 jfenwick 2496 G_TENSORPROD, // general tensor product
122 jfenwick 2721 G_NP1OUT_2P, // non-pointwise op with one output requiring two params
123     G_REDUCTION // non-pointwise unary op with a scalar output
124 jfenwick 1886 };
125    
126    
127    
128    
129 jfenwick 1910 string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
130     "sin","cos","tan",
131 jfenwick 1886 "asin","acos","atan","sinh","cosh","tanh","erf",
132     "asinh","acosh","atanh",
133 jfenwick 1888 "log10","log","sign","abs","neg","pos","exp","sqrt",
134 jfenwick 2147 "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
135 jfenwick 2066 "symmetric","nonsymmetric",
136 jfenwick 2084 "prod",
137 jfenwick 2496 "transpose", "trace",
138 jfenwick 2721 "swapaxes",
139     "minval", "maxval"};
140     int ES_opcount=43;
141 jfenwick 1910 ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
142     G_UNARY,G_UNARY,G_UNARY, //10
143     G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 17
144     G_UNARY,G_UNARY,G_UNARY, // 20
145 jfenwick 2147 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 28
146     G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P, // 35
147 jfenwick 2066 G_NP1OUT,G_NP1OUT,
148 jfenwick 2084 G_TENSORPROD,
149 jfenwick 2496 G_NP1OUT_P, G_NP1OUT_P,
150 jfenwick 2721 G_NP1OUT_2P,
151     G_REDUCTION, G_REDUCTION};
152 jfenwick 1886 inline
153     ES_opgroup
154     getOpgroup(ES_optype op)
155     {
156     return opgroups[op];
157     }
158    
159 jfenwick 1865 // return the FunctionSpace of the result of "left op right"
160     FunctionSpace
161     resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
162     {
163     // perhaps this should call interpolate and throw or something?
164     // maybe we need an interpolate node -
165     // that way, if interpolate is required in any other op we can just throw a
166     // programming error exception.
167 jfenwick 1879
168 jfenwick 1943 FunctionSpace l=left->getFunctionSpace();
169     FunctionSpace r=right->getFunctionSpace();
170     if (l!=r)
171     {
172     if (r.probeInterpolation(l))
173     {
174     return l;
175     }
176     if (l.probeInterpolation(r))
177     {
178     return r;
179     }
180     throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
181     }
182     return l;
183 jfenwick 1865 }
184    
185     // return the shape of the result of "left op right"
186 jfenwick 2066 // the shapes resulting from tensor product are more complex to compute so are worked out elsewhere
187 jfenwick 1865 DataTypes::ShapeType
188     resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
189     {
190 jfenwick 1879 if (left->getShape()!=right->getShape())
191     {
192 jfenwick 2037 if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
193 jfenwick 1908 {
194     throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
195     }
196 jfenwick 2721
197 jfenwick 1908 if (left->getRank()==0) // we need to allow scalar * anything
198     {
199     return right->getShape();
200     }
201     if (right->getRank()==0)
202     {
203     return left->getShape();
204     }
205     throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
206 jfenwick 1879 }
207     return left->getShape();
208 jfenwick 1865 }
209    
210 jfenwick 2084 // return the shape for "op left"
211    
212     DataTypes::ShapeType
213 jfenwick 2199 resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
214 jfenwick 2084 {
215     switch(op)
216     {
217     case TRANS:
218 jfenwick 2199 { // for the scoping of variables
219     const DataTypes::ShapeType& s=left->getShape();
220     DataTypes::ShapeType sh;
221     int rank=left->getRank();
222     if (axis_offset<0 || axis_offset>rank)
223     {
224     throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
225     }
226     for (int i=0; i<rank; i++)
227     {
228     int index = (axis_offset+i)%rank;
229     sh.push_back(s[index]); // Append to new shape
230     }
231     return sh;
232     }
233 jfenwick 2084 break;
234     case TRACE:
235 jfenwick 2199 {
236     int rank=left->getRank();
237     if (rank<2)
238     {
239     throw DataException("Trace can only be computed for objects with rank 2 or greater.");
240     }
241     if ((axis_offset>rank-2) || (axis_offset<0))
242     {
243     throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
244     }
245     if (rank==2)
246     {
247     return DataTypes::scalarShape;
248     }
249     else if (rank==3)
250     {
251     DataTypes::ShapeType sh;
252     if (axis_offset==0)
253     {
254     sh.push_back(left->getShape()[2]);
255     }
256     else // offset==1
257     {
258     sh.push_back(left->getShape()[0]);
259     }
260     return sh;
261     }
262     else if (rank==4)
263     {
264     DataTypes::ShapeType sh;
265     const DataTypes::ShapeType& s=left->getShape();
266     if (axis_offset==0)
267     {
268     sh.push_back(s[2]);
269     sh.push_back(s[3]);
270     }
271     else if (axis_offset==1)
272     {
273     sh.push_back(s[0]);
274     sh.push_back(s[3]);
275     }
276     else // offset==2
277     {
278     sh.push_back(s[0]);
279     sh.push_back(s[1]);
280     }
281     return sh;
282     }
283     else // unknown rank
284     {
285     throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
286     }
287     }
288 jfenwick 2084 break;
289     default:
290     throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
291     }
292     }
293    
294 jfenwick 2496 DataTypes::ShapeType
295     SwapShape(DataAbstract_ptr left, const int axis0, const int axis1)
296     {
297     // This code taken from the Data.cpp swapaxes() method
298     // Some of the checks are probably redundant here
299     int axis0_tmp,axis1_tmp;
300     const DataTypes::ShapeType& s=left->getShape();
301     DataTypes::ShapeType out_shape;
302     // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
303     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
304     int rank=left->getRank();
305     if (rank<2) {
306     throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
307     }
308     if (axis0<0 || axis0>rank-1) {
309     throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
310     }
311     if (axis1<0 || axis1>rank-1) {
312     throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
313     }
314     if (axis0 == axis1) {
315     throw DataException("Error - Data::swapaxes: axis indices must be different.");
316     }
317     if (axis0 > axis1) {
318     axis0_tmp=axis1;
319     axis1_tmp=axis0;
320     } else {
321     axis0_tmp=axis0;
322     axis1_tmp=axis1;
323     }
324     for (int i=0; i<rank; i++) {
325     if (i == axis0_tmp) {
326     out_shape.push_back(s[axis1_tmp]);
327     } else if (i == axis1_tmp) {
328     out_shape.push_back(s[axis0_tmp]);
329     } else {
330     out_shape.push_back(s[i]);
331     }
332     }
333     return out_shape;
334     }
335    
336    
337 jfenwick 2066 // determine the output shape for the general tensor product operation
338     // the additional parameters return information required later for the product
339     // the majority of this code is copy pasted from C_General_Tensor_Product
340     DataTypes::ShapeType
341     GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
342 jfenwick 1865 {
343 jfenwick 2066
344     // Get rank and shape of inputs
345     int rank0 = left->getRank();
346     int rank1 = right->getRank();
347     const DataTypes::ShapeType& shape0 = left->getShape();
348     const DataTypes::ShapeType& shape1 = right->getShape();
349    
350     // Prepare for the loops of the product and verify compatibility of shapes
351     int start0=0, start1=0;
352     if (transpose == 0) {}
353     else if (transpose == 1) { start0 = axis_offset; }
354     else if (transpose == 2) { start1 = rank1-axis_offset; }
355     else { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
356    
357 jfenwick 2085 if (rank0<axis_offset)
358     {
359     throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
360     }
361 jfenwick 2066
362     // Adjust the shapes for transpose
363     DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
364     DataTypes::ShapeType tmpShape1(rank1); // than using push_back
365     for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
366     for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
367    
368     // Prepare for the loops of the product
369     SL=1, SM=1, SR=1;
370     for (int i=0; i<rank0-axis_offset; i++) {
371     SL *= tmpShape0[i];
372     }
373     for (int i=rank0-axis_offset; i<rank0; i++) {
374     if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
375     throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
376     }
377     SM *= tmpShape0[i];
378     }
379     for (int i=axis_offset; i<rank1; i++) {
380     SR *= tmpShape1[i];
381     }
382    
383     // Define the shape of the output (rank of shape is the sum of the loop ranges below)
384     DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
385     { // block to limit the scope of out_index
386     int out_index=0;
387     for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
388     for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
389     }
390 jfenwick 2086
391     if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
392     {
393     ostringstream os;
394     os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
395     throw DataException(os.str());
396     }
397    
398 jfenwick 2066 return shape2;
399 jfenwick 1865 }
400    
401     } // end anonymous namespace
402    
403    
404 jfenwick 1899
405     // Return a string representing the operation
406 jfenwick 1865 const std::string&
407     opToString(ES_optype op)
408     {
409     if (op<0 || op>=ES_opcount)
410     {
411     op=UNKNOWNOP;
412     }
413     return ES_opstrings[op];
414     }
415    
416 jfenwick 2500 void DataLazy::LazyNodeSetup()
417     {
418     #ifdef _OPENMP
419     int numthreads=omp_get_max_threads();
420     m_samples.resize(numthreads*m_samplesize);
421     m_sampleids=new int[numthreads];
422     for (int i=0;i<numthreads;++i)
423     {
424     m_sampleids[i]=-1;
425     }
426     #else
427     m_samples.resize(m_samplesize);
428     m_sampleids=new int[1];
429     m_sampleids[0]=-1;
430     #endif // _OPENMP
431     }
432 jfenwick 1865
433 jfenwick 2500
434     // Creates an identity node
435 jfenwick 1865 DataLazy::DataLazy(DataAbstract_ptr p)
436 jfenwick 2177 : parent(p->getFunctionSpace(),p->getShape())
437 jfenwick 2500 ,m_sampleids(0),
438     m_samples(1)
439 jfenwick 1865 {
440 jfenwick 1879 if (p->isLazy())
441     {
442     // I don't want identity of Lazy.
443     // Question: Why would that be so bad?
444     // Answer: We assume that the child of ID is something we can call getVector on
445     throw DataException("Programmer error - attempt to create identity from a DataLazy.");
446     }
447     else
448     {
449 jfenwick 2271 p->makeLazyShared();
450 jfenwick 2177 DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
451     makeIdentity(dr);
452 jfenwick 2199 LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
453 jfenwick 1879 }
454 jfenwick 2092 LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
455 jfenwick 1865 }
456    
457 jfenwick 1886 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
458 jfenwick 2721 : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),
459 jfenwick 2066 m_op(op),
460     m_axis_offset(0),
461     m_transpose(0),
462     m_SL(0), m_SM(0), m_SR(0)
463 jfenwick 1886 {
464 jfenwick 2721 if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT) && (getOpgroup(op)!=G_REDUCTION))
465 jfenwick 1886 {
466     throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
467     }
468 jfenwick 2066
469 jfenwick 1886 DataLazy_ptr lleft;
470     if (!left->isLazy())
471     {
472     lleft=DataLazy_ptr(new DataLazy(left));
473     }
474     else
475     {
476     lleft=dynamic_pointer_cast<DataLazy>(left);
477     }
478 jfenwick 1889 m_readytype=lleft->m_readytype;
479 jfenwick 1886 m_left=lleft;
480     m_samplesize=getNumDPPSample()*getNoValues();
481 jfenwick 2177 m_children=m_left->m_children+1;
482     m_height=m_left->m_height+1;
483 jfenwick 2500 LazyNodeSetup();
484 jfenwick 2177 SIZELIMIT
485 jfenwick 1886 }
486    
487    
488 jfenwick 1943 // In this constructor we need to consider interpolation
489 jfenwick 1879 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
490     : parent(resultFS(left,right,op), resultShape(left,right,op)),
491 jfenwick 2066 m_op(op),
492     m_SL(0), m_SM(0), m_SR(0)
493 jfenwick 1879 {
494 jfenwick 2199 LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
495 jfenwick 2037 if ((getOpgroup(op)!=G_BINARY))
496 jfenwick 1886 {
497 jfenwick 1889 throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
498 jfenwick 1886 }
499 jfenwick 1943
500     if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
501     {
502     FunctionSpace fs=getFunctionSpace();
503     Data ltemp(left);
504     Data tmp(ltemp,fs);
505     left=tmp.borrowDataPtr();
506     }
507 jfenwick 2066 if (getFunctionSpace()!=right->getFunctionSpace()) // right needs to be interpolated
508 jfenwick 1943 {
509     Data tmp(Data(right),getFunctionSpace());
510     right=tmp.borrowDataPtr();
511 jfenwick 2199 LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
512 jfenwick 1943 }
513     left->operandCheck(*right);
514    
515 jfenwick 1899 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
516 jfenwick 1879 {
517     m_left=dynamic_pointer_cast<DataLazy>(left);
518 jfenwick 2199 LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
519 jfenwick 1879 }
520     else
521     {
522     m_left=DataLazy_ptr(new DataLazy(left));
523 jfenwick 2199 LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
524 jfenwick 1879 }
525     if (right->isLazy())
526     {
527     m_right=dynamic_pointer_cast<DataLazy>(right);
528 jfenwick 2199 LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
529 jfenwick 1879 }
530     else
531     {
532     m_right=DataLazy_ptr(new DataLazy(right));
533 jfenwick 2199 LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
534 jfenwick 1879 }
535 jfenwick 1889 char lt=m_left->m_readytype;
536     char rt=m_right->m_readytype;
537     if (lt=='E' || rt=='E')
538     {
539     m_readytype='E';
540     }
541     else if (lt=='T' || rt=='T')
542     {
543     m_readytype='T';
544     }
545     else
546     {
547     m_readytype='C';
548     }
549 jfenwick 2066 m_samplesize=getNumDPPSample()*getNoValues();
550 jfenwick 2177 m_children=m_left->m_children+m_right->m_children+2;
551     m_height=max(m_left->m_height,m_right->m_height)+1;
552 jfenwick 2500 LazyNodeSetup();
553 jfenwick 2177 SIZELIMIT
554 jfenwick 2092 LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
555 jfenwick 1879 }
556    
557 jfenwick 2066 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
558     : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
559     m_op(op),
560     m_axis_offset(axis_offset),
561     m_transpose(transpose)
562     {
563     if ((getOpgroup(op)!=G_TENSORPROD))
564     {
565     throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
566     }
567     if ((transpose>2) || (transpose<0))
568     {
569     throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
570     }
571     if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
572     {
573     FunctionSpace fs=getFunctionSpace();
574     Data ltemp(left);
575     Data tmp(ltemp,fs);
576     left=tmp.borrowDataPtr();
577     }
578     if (getFunctionSpace()!=right->getFunctionSpace()) // right needs to be interpolated
579     {
580     Data tmp(Data(right),getFunctionSpace());
581     right=tmp.borrowDataPtr();
582     }
583 jfenwick 2195 // left->operandCheck(*right);
584 jfenwick 1879
585 jfenwick 2066 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
586     {
587     m_left=dynamic_pointer_cast<DataLazy>(left);
588     }
589     else
590     {
591     m_left=DataLazy_ptr(new DataLazy(left));
592     }
593     if (right->isLazy())
594     {
595     m_right=dynamic_pointer_cast<DataLazy>(right);
596     }
597     else
598     {
599     m_right=DataLazy_ptr(new DataLazy(right));
600     }
601     char lt=m_left->m_readytype;
602     char rt=m_right->m_readytype;
603     if (lt=='E' || rt=='E')
604     {
605     m_readytype='E';
606     }
607     else if (lt=='T' || rt=='T')
608     {
609     m_readytype='T';
610     }
611     else
612     {
613     m_readytype='C';
614     }
615     m_samplesize=getNumDPPSample()*getNoValues();
616 jfenwick 2177 m_children=m_left->m_children+m_right->m_children+2;
617     m_height=max(m_left->m_height,m_right->m_height)+1;
618 jfenwick 2500 LazyNodeSetup();
619 jfenwick 2177 SIZELIMIT
620 jfenwick 2092 LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
621 jfenwick 2066 }
622    
623    
624 jfenwick 2084 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
625 jfenwick 2199 : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
626 jfenwick 2084 m_op(op),
627     m_axis_offset(axis_offset),
628 jfenwick 2147 m_transpose(0),
629     m_tol(0)
630 jfenwick 2084 {
631     if ((getOpgroup(op)!=G_NP1OUT_P))
632     {
633     throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
634     }
635     DataLazy_ptr lleft;
636     if (!left->isLazy())
637     {
638     lleft=DataLazy_ptr(new DataLazy(left));
639     }
640     else
641     {
642     lleft=dynamic_pointer_cast<DataLazy>(left);
643     }
644     m_readytype=lleft->m_readytype;
645     m_left=lleft;
646     m_samplesize=getNumDPPSample()*getNoValues();
647 jfenwick 2177 m_children=m_left->m_children+1;
648     m_height=m_left->m_height+1;
649 jfenwick 2500 LazyNodeSetup();
650 jfenwick 2177 SIZELIMIT
651 jfenwick 2092 LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
652 jfenwick 2084 }
653    
654 jfenwick 2147 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
655     : parent(left->getFunctionSpace(), left->getShape()),
656     m_op(op),
657     m_axis_offset(0),
658     m_transpose(0),
659     m_tol(tol)
660     {
661     if ((getOpgroup(op)!=G_UNARY_P))
662     {
663     throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
664     }
665     DataLazy_ptr lleft;
666     if (!left->isLazy())
667     {
668     lleft=DataLazy_ptr(new DataLazy(left));
669     }
670     else
671     {
672     lleft=dynamic_pointer_cast<DataLazy>(left);
673     }
674     m_readytype=lleft->m_readytype;
675     m_left=lleft;
676     m_samplesize=getNumDPPSample()*getNoValues();
677 jfenwick 2177 m_children=m_left->m_children+1;
678     m_height=m_left->m_height+1;
679 jfenwick 2500 LazyNodeSetup();
680 jfenwick 2177 SIZELIMIT
681 jfenwick 2147 LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
682     }
683 jfenwick 2084
684 jfenwick 2496
685     DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, const int axis0, const int axis1)
686     : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
687     m_op(op),
688     m_axis_offset(axis0),
689     m_transpose(axis1),
690     m_tol(0)
691     {
692     if ((getOpgroup(op)!=G_NP1OUT_2P))
693     {
694     throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require two integer parameters.");
695     }
696     DataLazy_ptr lleft;
697     if (!left->isLazy())
698     {
699     lleft=DataLazy_ptr(new DataLazy(left));
700     }
701     else
702     {
703     lleft=dynamic_pointer_cast<DataLazy>(left);
704     }
705     m_readytype=lleft->m_readytype;
706     m_left=lleft;
707     m_samplesize=getNumDPPSample()*getNoValues();
708     m_children=m_left->m_children+1;
709     m_height=m_left->m_height+1;
710 jfenwick 2500 LazyNodeSetup();
711 jfenwick 2496 SIZELIMIT
712     LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)
713     }
714    
715 jfenwick 1865 DataLazy::~DataLazy()
716     {
717 jfenwick 2500 delete[] m_sampleids;
718 jfenwick 1865 }
719    
720 jfenwick 1879
721 jfenwick 1899 /*
722     \brief Evaluates the expression using methods on Data.
723     This does the work for the collapse method.
724     For reasons of efficiency do not call this method on DataExpanded nodes.
725     */
726 jfenwick 1889 DataReady_ptr
727     DataLazy::collapseToReady()
728     {
729     if (m_readytype=='E')
730     { // this is more an efficiency concern than anything else
731     throw DataException("Programmer Error - do not use collapse on Expanded data.");
732     }
733     if (m_op==IDENTITY)
734     {
735     return m_id;
736     }
737     DataReady_ptr pleft=m_left->collapseToReady();
738     Data left(pleft);
739     Data right;
740 jfenwick 2066 if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
741 jfenwick 1889 {
742     right=Data(m_right->collapseToReady());
743     }
744     Data result;
745     switch(m_op)
746     {
747     case ADD:
748     result=left+right;
749     break;
750     case SUB:
751     result=left-right;
752     break;
753     case MUL:
754     result=left*right;
755     break;
756     case DIV:
757     result=left/right;
758     break;
759     case SIN:
760     result=left.sin();
761     break;
762     case COS:
763     result=left.cos();
764     break;
765     case TAN:
766     result=left.tan();
767     break;
768     case ASIN:
769     result=left.asin();
770     break;
771     case ACOS:
772     result=left.acos();
773     break;
774     case ATAN:
775     result=left.atan();
776     break;
777     case SINH:
778     result=left.sinh();
779     break;
780     case COSH:
781     result=left.cosh();
782     break;
783     case TANH:
784     result=left.tanh();
785     break;
786     case ERF:
787     result=left.erf();
788     break;
789     case ASINH:
790     result=left.asinh();
791     break;
792     case ACOSH:
793     result=left.acosh();
794     break;
795     case ATANH:
796     result=left.atanh();
797     break;
798     case LOG10:
799     result=left.log10();
800     break;
801     case LOG:
802     result=left.log();
803     break;
804     case SIGN:
805     result=left.sign();
806     break;
807     case ABS:
808     result=left.abs();
809     break;
810     case NEG:
811     result=left.neg();
812     break;
813     case POS:
814     // it doesn't mean anything for delayed.
815     // it will just trigger a deep copy of the lazy object
816     throw DataException("Programmer error - POS not supported for lazy data.");
817     break;
818     case EXP:
819     result=left.exp();
820     break;
821     case SQRT:
822     result=left.sqrt();
823     break;
824     case RECIP:
825     result=left.oneOver();
826     break;
827     case GZ:
828     result=left.wherePositive();
829     break;
830     case LZ:
831     result=left.whereNegative();
832     break;
833     case GEZ:
834     result=left.whereNonNegative();
835     break;
836     case LEZ:
837     result=left.whereNonPositive();
838     break;
839 jfenwick 2147 case NEZ:
840     result=left.whereNonZero(m_tol);
841     break;
842     case EZ:
843     result=left.whereZero(m_tol);
844     break;
845 jfenwick 2037 case SYM:
846     result=left.symmetric();
847     break;
848     case NSYM:
849     result=left.nonsymmetric();
850     break;
851 jfenwick 2066 case PROD:
852     result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
853     break;
854 jfenwick 2084 case TRANS:
855     result=left.transpose(m_axis_offset);
856     break;
857     case TRACE:
858     result=left.trace(m_axis_offset);
859     break;
860 jfenwick 2496 case SWAP:
861     result=left.swapaxes(m_axis_offset, m_transpose);
862     break;
863 jfenwick 2721 case MINVAL:
864     result=left.minval();
865     break;
866     case MAXVAL:
867     result=left.minval();
868     break;
869 jfenwick 1889 default:
870 jfenwick 2037 throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
871 jfenwick 1889 }
872     return result.borrowReadyPtr();
873     }
874    
875 jfenwick 1899 /*
876     \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
877     This method uses the original methods on the Data class to evaluate the expressions.
878     For this reason, it should not be used on DataExpanded instances. (To do so would defeat
879     the purpose of using DataLazy in the first place).
880     */
881 jfenwick 1889 void
882     DataLazy::collapse()
883     {
884     if (m_op==IDENTITY)
885     {
886     return;
887     }
888     if (m_readytype=='E')
889     { // this is more an efficiency concern than anything else
890     throw DataException("Programmer Error - do not use collapse on Expanded data.");
891     }
892     m_id=collapseToReady();
893     m_op=IDENTITY;
894     }
895    
896 jfenwick 1899
897 jfenwick 1889
898    
899    
900 jfenwick 2147
901 phornby 1987 #define PROC_OP(TYPE,X) \
902 jfenwick 2152 for (int j=0;j<onumsteps;++j)\
903     {\
904     for (int i=0;i<numsteps;++i,resultp+=resultStep) \
905     { \
906     LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
907 jfenwick 2157 LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
908 jfenwick 2152 tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
909 jfenwick 2157 LAZYDEBUG(cout << " result= " << resultp[0] << endl;) \
910 jfenwick 2152 lroffset+=leftstep; \
911     rroffset+=rightstep; \
912     }\
913     lroffset+=oleftstep;\
914     rroffset+=orightstep;\
915 jfenwick 1889 }
916    
917 jfenwick 1899
918 jfenwick 2500 // The result will be stored in m_samples
919     // The return value is a pointer to the DataVector, offset is the offset within the return value
920     const DataTypes::ValueType*
921     DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset)
922     {
923     LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
924     // collapse so we have a 'E' node or an IDENTITY for some other type
925     if (m_readytype!='E' && m_op!=IDENTITY)
926     {
927     collapse();
928     }
929     if (m_op==IDENTITY)
930     {
931     const ValueType& vec=m_id->getVectorRO();
932     roffset=m_id->getPointOffset(sampleNo, 0);
933     return &(vec);
934     }
935     if (m_readytype!='E')
936     {
937     throw DataException("Programmer Error - Collapse did not produce an expanded node.");
938     }
939 jfenwick 2501 if (m_sampleids[tid]==sampleNo)
940     {
941     roffset=tid*m_samplesize;
942     return &(m_samples); // sample is already resolved
943     }
944     m_sampleids[tid]=sampleNo;
945 jfenwick 2500 switch (getOpgroup(m_op))
946     {
947     case G_UNARY:
948     case G_UNARY_P: return resolveNodeUnary(tid, sampleNo, roffset);
949     case G_BINARY: return resolveNodeBinary(tid, sampleNo, roffset);
950     case G_NP1OUT: return resolveNodeNP1OUT(tid, sampleNo, roffset);
951     case G_NP1OUT_P: return resolveNodeNP1OUT_P(tid, sampleNo, roffset);
952     case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);
953     case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);
954 jfenwick 2721 case G_REDUCTION: return resolveNodeReduction(tid, sampleNo, roffset);
955 jfenwick 2500 default:
956     throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
957     }
958     }
959    
960     const DataTypes::ValueType*
961     DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset)
962     {
963     // we assume that any collapsing has been done before we get here
964     // since we only have one argument we don't need to think about only
965     // processing single points.
966     // we will also know we won't get identity nodes
967     if (m_readytype!='E')
968     {
969     throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
970     }
971     if (m_op==IDENTITY)
972     {
973     throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
974     }
975     const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
976     const double* left=&((*leftres)[roffset]);
977     roffset=m_samplesize*tid;
978     double* result=&(m_samples[roffset]);
979     switch (m_op)
980     {
981     case SIN:
982     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
983     break;
984     case COS:
985     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
986     break;
987     case TAN:
988     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
989     break;
990     case ASIN:
991     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
992     break;
993     case ACOS:
994     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
995     break;
996     case ATAN:
997     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
998     break;
999     case SINH:
1000     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
1001     break;
1002     case COSH:
1003     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
1004     break;
1005     case TANH:
1006     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
1007     break;
1008     case ERF:
1009     #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1010     throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1011     #else
1012     tensor_unary_operation(m_samplesize, left, result, ::erf);
1013     break;
1014     #endif
1015     case ASINH:
1016     #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1017     tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
1018     #else
1019     tensor_unary_operation(m_samplesize, left, result, ::asinh);
1020     #endif
1021     break;
1022     case ACOSH:
1023     #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1024     tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
1025     #else
1026     tensor_unary_operation(m_samplesize, left, result, ::acosh);
1027     #endif
1028     break;
1029     case ATANH:
1030     #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1031     tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
1032     #else
1033     tensor_unary_operation(m_samplesize, left, result, ::atanh);
1034     #endif
1035     break;
1036     case LOG10:
1037     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
1038     break;
1039     case LOG:
1040     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
1041     break;
1042     case SIGN:
1043     tensor_unary_operation(m_samplesize, left, result, escript::fsign);
1044     break;
1045     case ABS:
1046     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
1047     break;
1048     case NEG:
1049     tensor_unary_operation(m_samplesize, left, result, negate<double>());
1050     break;
1051     case POS:
1052     // it doesn't mean anything for delayed.
1053     // it will just trigger a deep copy of the lazy object
1054     throw DataException("Programmer error - POS not supported for lazy data.");
1055     break;
1056     case EXP:
1057     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
1058     break;
1059     case SQRT:
1060     tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
1061     break;
1062     case RECIP:
1063     tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
1064     break;
1065     case GZ:
1066     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
1067     break;
1068     case LZ:
1069     tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
1070     break;
1071     case GEZ:
1072     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
1073     break;
1074     case LEZ:
1075     tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
1076     break;
1077     // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
1078     case NEZ:
1079     tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
1080     break;
1081     case EZ:
1082     tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
1083     break;
1084    
1085     default:
1086     throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1087     }
1088     return &(m_samples);
1089     }
1090    
1091    
1092     const DataTypes::ValueType*
1093 jfenwick 2721 DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset)
1094     {
1095     // we assume that any collapsing has been done before we get here
1096     // since we only have one argument we don't need to think about only
1097     // processing single points.
1098     // we will also know we won't get identity nodes
1099     if (m_readytype!='E')
1100     {
1101     throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1102     }
1103     if (m_op==IDENTITY)
1104     {
1105     throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1106     }
1107     size_t loffset=0;
1108     const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);
1109    
1110     roffset=m_samplesize*tid;
1111 jfenwick 2734 unsigned int ndpps=getNumDPPSample();
1112     unsigned int psize=DataTypes::noValues(getShape());
1113 jfenwick 2721 double* result=&(m_samples[roffset]);
1114     switch (m_op)
1115     {
1116     case MINVAL:
1117     {
1118 jfenwick 2734 for (unsigned int z=0;z<ndpps;++z)
1119     {
1120     FMin op;
1121     *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());
1122     loffset+=psize;
1123     result++;
1124     }
1125 jfenwick 2721 }
1126     break;
1127     case MAXVAL:
1128     {
1129 jfenwick 2734 for (unsigned int z=0;z<ndpps;++z)
1130     {
1131 jfenwick 2721 FMax op;
1132     *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);
1133 jfenwick 2734 loffset+=psize;
1134     result++;
1135     }
1136 jfenwick 2721 }
1137     break;
1138     default:
1139     throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1140     }
1141     return &(m_samples);
1142     }
1143    
1144     const DataTypes::ValueType*
1145 jfenwick 2500 DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset)
1146     {
1147     // we assume that any collapsing has been done before we get here
1148     // since we only have one argument we don't need to think about only
1149     // processing single points.
1150     if (m_readytype!='E')
1151     {
1152     throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");
1153     }
1154     if (m_op==IDENTITY)
1155     {
1156     throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1157     }
1158     size_t subroffset;
1159     const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1160     roffset=m_samplesize*tid;
1161     size_t loop=0;
1162     size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1163     size_t step=getNoValues();
1164     size_t offset=roffset;
1165     switch (m_op)
1166     {
1167     case SYM:
1168     for (loop=0;loop<numsteps;++loop)
1169     {
1170     DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1171     subroffset+=step;
1172     offset+=step;
1173     }
1174     break;
1175     case NSYM:
1176     for (loop=0;loop<numsteps;++loop)
1177     {
1178     DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1179     subroffset+=step;
1180     offset+=step;
1181     }
1182     break;
1183     default:
1184     throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1185     }
1186     return &m_samples;
1187     }
1188    
1189     const DataTypes::ValueType*
1190     DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset)
1191     {
1192     // we assume that any collapsing has been done before we get here
1193     // since we only have one argument we don't need to think about only
1194     // processing single points.
1195     if (m_readytype!='E')
1196     {
1197     throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");
1198     }
1199     if (m_op==IDENTITY)
1200     {
1201     throw DataException("Programmer error - resolveNodeNP1OUT_P should not be called on identity nodes.");
1202     }
1203     size_t subroffset;
1204     size_t offset;
1205     const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1206     roffset=m_samplesize*tid;
1207     offset=roffset;
1208     size_t loop=0;
1209     size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1210     size_t outstep=getNoValues();
1211     size_t instep=m_left->getNoValues();
1212     switch (m_op)
1213     {
1214     case TRACE:
1215     for (loop=0;loop<numsteps;++loop)
1216     {
1217     DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);
1218     subroffset+=instep;
1219     offset+=outstep;
1220     }
1221     break;
1222     case TRANS:
1223     for (loop=0;loop<numsteps;++loop)
1224     {
1225     DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);
1226     subroffset+=instep;
1227     offset+=outstep;
1228     }
1229     break;
1230     default:
1231     throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1232     }
1233     return &m_samples;
1234     }
1235    
1236    
1237     const DataTypes::ValueType*
1238     DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset)
1239     {
1240     if (m_readytype!='E')
1241     {
1242     throw DataException("Programmer error - resolveNodeNP1OUT_2P should only be called on expanded Data.");
1243     }
1244     if (m_op==IDENTITY)
1245     {
1246     throw DataException("Programmer error - resolveNodeNP1OUT_2P should not be called on identity nodes.");
1247     }
1248     size_t subroffset;
1249     size_t offset;
1250     const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1251     roffset=m_samplesize*tid;
1252     offset=roffset;
1253     size_t loop=0;
1254     size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1255     size_t outstep=getNoValues();
1256     size_t instep=m_left->getNoValues();
1257     switch (m_op)
1258     {
1259     case SWAP:
1260     for (loop=0;loop<numsteps;++loop)
1261     {
1262     DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);
1263     subroffset+=instep;
1264     offset+=outstep;
1265     }
1266     break;
1267     default:
1268     throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1269     }
1270     return &m_samples;
1271     }
1272    
1273    
1274    
1275     // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1276     // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1277     // If both children are expanded, then we can process them in a single operation (we treat
1278     // the whole sample as one big datapoint.
1279     // If one of the children is not expanded, then we need to treat each point in the sample
1280     // individually.
1281     // There is an additional complication when scalar operations are considered.
1282     // For example, 2+Vector.
1283     // In this case each double within the point is treated individually
1284     const DataTypes::ValueType*
1285     DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset)
1286     {
1287     LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1288    
1289     size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
1290     // first work out which of the children are expanded
1291     bool leftExp=(m_left->m_readytype=='E');
1292     bool rightExp=(m_right->m_readytype=='E');
1293     if (!leftExp && !rightExp)
1294     {
1295     throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1296     }
1297     bool leftScalar=(m_left->getRank()==0);
1298     bool rightScalar=(m_right->getRank()==0);
1299     if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1300     {
1301     throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1302     }
1303     size_t leftsize=m_left->getNoValues();
1304     size_t rightsize=m_right->getNoValues();
1305     size_t chunksize=1; // how many doubles will be processed in one go
1306     int leftstep=0; // how far should the left offset advance after each step
1307     int rightstep=0;
1308     int numsteps=0; // total number of steps for the inner loop
1309     int oleftstep=0; // the o variables refer to the outer loop
1310     int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1311     int onumsteps=1;
1312    
1313     bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1314     bool RES=(rightExp && rightScalar);
1315     bool LS=(!leftExp && leftScalar); // left is a single scalar
1316     bool RS=(!rightExp && rightScalar);
1317     bool LN=(!leftExp && !leftScalar); // left is a single non-scalar
1318     bool RN=(!rightExp && !rightScalar);
1319     bool LEN=(leftExp && !leftScalar); // left is an expanded non-scalar
1320     bool REN=(rightExp && !rightScalar);
1321    
1322     if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1323     {
1324     chunksize=m_left->getNumDPPSample()*leftsize;
1325     leftstep=0;
1326     rightstep=0;
1327     numsteps=1;
1328     }
1329     else if (LES || RES)
1330     {
1331     chunksize=1;
1332     if (LES) // left is an expanded scalar
1333     {
1334     if (RS)
1335     {
1336     leftstep=1;
1337     rightstep=0;
1338     numsteps=m_left->getNumDPPSample();
1339     }
1340     else // RN or REN
1341     {
1342     leftstep=0;
1343     oleftstep=1;
1344     rightstep=1;
1345     orightstep=(RN ? -(int)rightsize : 0);
1346     numsteps=rightsize;
1347     onumsteps=m_left->getNumDPPSample();
1348     }
1349     }
1350     else // right is an expanded scalar
1351     {
1352     if (LS)
1353     {
1354     rightstep=1;
1355     leftstep=0;
1356     numsteps=m_right->getNumDPPSample();
1357     }
1358     else
1359     {
1360     rightstep=0;
1361     orightstep=1;
1362     leftstep=1;
1363     oleftstep=(LN ? -(int)leftsize : 0);
1364     numsteps=leftsize;
1365     onumsteps=m_right->getNumDPPSample();
1366     }
1367     }
1368     }
1369     else // this leaves (LEN, RS), (LEN, RN) and their transposes
1370     {
1371     if (LEN) // and Right will be a single value
1372     {
1373     chunksize=rightsize;
1374     leftstep=rightsize;
1375     rightstep=0;
1376     numsteps=m_left->getNumDPPSample();
1377     if (RS)
1378     {
1379     numsteps*=leftsize;
1380     }
1381     }
1382     else // REN
1383     {
1384     chunksize=leftsize;
1385     rightstep=leftsize;
1386     leftstep=0;
1387     numsteps=m_right->getNumDPPSample();
1388     if (LS)
1389     {
1390     numsteps*=rightsize;
1391     }
1392     }
1393     }
1394    
1395     int resultStep=max(leftstep,rightstep); // only one (at most) should be !=0
1396     // Get the values of sub-expressions
1397     const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);
1398     const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
1399     LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1400     LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1401     LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1402     LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1403     LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1404     LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1405     LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN << endl;)
1406    
1407     LAZYDEBUG(cout << "Left res["<< lroffset<< "]=" << (*left)[lroffset] << endl;)
1408     LAZYDEBUG(cout << "Right res["<< rroffset<< "]=" << (*right)[rroffset] << endl;)
1409    
1410    
1411     roffset=m_samplesize*tid;
1412     double* resultp=&(m_samples[roffset]); // results are stored at the vector offset we recieved
1413     switch(m_op)
1414     {
1415     case ADD:
1416     PROC_OP(NO_ARG,plus<double>());
1417     break;
1418     case SUB:
1419     PROC_OP(NO_ARG,minus<double>());
1420     break;
1421     case MUL:
1422     PROC_OP(NO_ARG,multiplies<double>());
1423     break;
1424     case DIV:
1425     PROC_OP(NO_ARG,divides<double>());
1426     break;
1427     case POW:
1428     PROC_OP(double (double,double),::pow);
1429     break;
1430     default:
1431     throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1432     }
1433     LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
1434     return &m_samples;
1435     }
1436    
1437    
1438     // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1439     // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1440     // unlike the other resolve helpers, we must treat these datapoints separately.
1441     const DataTypes::ValueType*
1442     DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset)
1443     {
1444     LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1445    
1446     size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
1447     // first work out which of the children are expanded
1448     bool leftExp=(m_left->m_readytype=='E');
1449     bool rightExp=(m_right->m_readytype=='E');
1450     int steps=getNumDPPSample();
1451     int leftStep=(leftExp? m_left->getNoValues() : 0); // do not have scalars as input to this method
1452     int rightStep=(rightExp?m_right->getNoValues() : 0);
1453    
1454     int resultStep=getNoValues();
1455     roffset=m_samplesize*tid;
1456     size_t offset=roffset;
1457    
1458     const ValueType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
1459    
1460     const ValueType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
1461    
1462     LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) << endl;
1463     cout << getNoValues() << endl;)
1464    
1465    
1466     LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1467     LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1468     LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1469     LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1470     LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1471     LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1472     LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1473    
1474     double* resultp=&(m_samples[offset]); // results are stored at the vector offset we recieved
1475     switch(m_op)
1476     {
1477     case PROD:
1478     for (int i=0;i<steps;++i,resultp+=resultStep)
1479     {
1480     const double *ptr_0 = &((*left)[lroffset]);
1481     const double *ptr_1 = &((*right)[rroffset]);
1482    
1483     LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1484     LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1485    
1486     matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
1487    
1488     lroffset+=leftStep;
1489     rroffset+=rightStep;
1490     }
1491     break;
1492     default:
1493     throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1494     }
1495     roffset=offset;
1496     return &m_samples;
1497     }
1498    
1499 jfenwick 1898
1500     const DataTypes::ValueType*
1501 jfenwick 2770 DataLazy::resolveSample(int sampleNo, size_t& roffset)
1502 jfenwick 1879 {
1503 jfenwick 2271 #ifdef _OPENMP
1504     int tid=omp_get_thread_num();
1505     #else
1506     int tid=0;
1507     #endif
1508 jfenwick 2521 return resolveNodeSample(tid, sampleNo, roffset);
1509 jfenwick 2271 }
1510    
1511    
1512 jfenwick 2497 // This needs to do the work of the identity constructor
1513 jfenwick 2177 void
1514     DataLazy::resolveToIdentity()
1515     {
1516     if (m_op==IDENTITY)
1517     return;
1518 jfenwick 2500 DataReady_ptr p=resolveNodeWorker();
1519 jfenwick 2177 makeIdentity(p);
1520     }
1521 jfenwick 1889
1522 jfenwick 2177 void DataLazy::makeIdentity(const DataReady_ptr& p)
1523     {
1524     m_op=IDENTITY;
1525     m_axis_offset=0;
1526     m_transpose=0;
1527     m_SL=m_SM=m_SR=0;
1528     m_children=m_height=0;
1529     m_id=p;
1530     if(p->isConstant()) {m_readytype='C';}
1531     else if(p->isExpanded()) {m_readytype='E';}
1532     else if (p->isTagged()) {m_readytype='T';}
1533     else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1534     m_samplesize=p->getNumDPPSample()*p->getNoValues();
1535     m_left.reset();
1536     m_right.reset();
1537     }
1538    
1539 jfenwick 2497
1540     DataReady_ptr
1541     DataLazy::resolve()
1542     {
1543     resolveToIdentity();
1544     return m_id;
1545     }
1546    
1547 jfenwick 2500 // This version of resolve uses storage in each node to hold results
1548     DataReady_ptr
1549     DataLazy::resolveNodeWorker()
1550     {
1551     if (m_readytype!='E') // if the whole sub-expression is Constant or Tagged, then evaluate it normally
1552     {
1553     collapse();
1554     }
1555     if (m_op==IDENTITY) // So a lazy expression of Constant or Tagged data will be returned here.
1556     {
1557     return m_id;
1558     }
1559     // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
1560     DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(), ValueType(getNoValues()));
1561     ValueType& resvec=result->getVectorRW();
1562     DataReady_ptr resptr=DataReady_ptr(result);
1563    
1564     int sample;
1565     int totalsamples=getNumSamples();
1566     const ValueType* res=0; // Storage for answer
1567     LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1568     #pragma omp parallel for private(sample,res) schedule(static)
1569     for (sample=0;sample<totalsamples;++sample)
1570     {
1571     size_t roffset=0;
1572     #ifdef _OPENMP
1573     res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
1574     #else
1575     res=resolveNodeSample(0,sample,roffset);
1576     #endif
1577     LAZYDEBUG(cout << "Sample #" << sample << endl;)
1578     LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
1579     DataVector::size_type outoffset=result->getPointOffset(sample,0);
1580     memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(DataVector::ElementType));
1581     }
1582     return resptr;
1583     }
1584    
1585 jfenwick 1865 std::string
1586     DataLazy::toString() const
1587     {
1588 jfenwick 1886 ostringstream oss;
1589     oss << "Lazy Data:";
1590 jfenwick 2737 if (escriptParams.getPRINT_LAZY_TREE()==0)
1591     {
1592     intoString(oss);
1593     }
1594     else
1595     {
1596     oss << endl;
1597     intoTreeString(oss,"");
1598     }
1599 jfenwick 1886 return oss.str();
1600 jfenwick 1865 }
1601    
1602 jfenwick 1899
1603 jfenwick 1886 void
1604     DataLazy::intoString(ostringstream& oss) const
1605     {
1606 jfenwick 2271 // oss << "[" << m_children <<";"<<m_height <<"]";
1607 jfenwick 1886 switch (getOpgroup(m_op))
1608     {
1609 jfenwick 1889 case G_IDENTITY:
1610     if (m_id->isExpanded())
1611     {
1612     oss << "E";
1613     }
1614     else if (m_id->isTagged())
1615     {
1616     oss << "T";
1617     }
1618     else if (m_id->isConstant())
1619     {
1620     oss << "C";
1621     }
1622     else
1623     {
1624     oss << "?";
1625     }
1626 jfenwick 1886 oss << '@' << m_id.get();
1627     break;
1628     case G_BINARY:
1629     oss << '(';
1630     m_left->intoString(oss);
1631     oss << ' ' << opToString(m_op) << ' ';
1632     m_right->intoString(oss);
1633     oss << ')';
1634     break;
1635     case G_UNARY:
1636 jfenwick 2147 case G_UNARY_P:
1637 jfenwick 2037 case G_NP1OUT:
1638 jfenwick 2084 case G_NP1OUT_P:
1639 jfenwick 2721 case G_REDUCTION:
1640 jfenwick 1886 oss << opToString(m_op) << '(';
1641     m_left->intoString(oss);
1642     oss << ')';
1643     break;
1644 jfenwick 2066 case G_TENSORPROD:
1645     oss << opToString(m_op) << '(';
1646     m_left->intoString(oss);
1647     oss << ", ";
1648     m_right->intoString(oss);
1649     oss << ')';
1650     break;
1651 jfenwick 2496 case G_NP1OUT_2P:
1652     oss << opToString(m_op) << '(';
1653     m_left->intoString(oss);
1654     oss << ", " << m_axis_offset << ", " << m_transpose;
1655     oss << ')';
1656     break;
1657 jfenwick 1886 default:
1658     oss << "UNKNOWN";
1659     }
1660     }
1661    
1662 jfenwick 2737
1663     void
1664     DataLazy::intoTreeString(ostringstream& oss, string indent) const
1665     {
1666     oss << '[' << m_rank << ':' << setw(3) << m_samplesize << "] " << indent;
1667     switch (getOpgroup(m_op))
1668     {
1669     case G_IDENTITY:
1670     if (m_id->isExpanded())
1671     {
1672     oss << "E";
1673     }
1674     else if (m_id->isTagged())
1675     {
1676     oss << "T";
1677     }
1678     else if (m_id->isConstant())
1679     {
1680     oss << "C";
1681     }
1682     else
1683     {
1684     oss << "?";
1685     }
1686     oss << '@' << m_id.get() << endl;
1687     break;
1688     case G_BINARY:
1689     oss << opToString(m_op) << endl;
1690     indent+='.';
1691     m_left->intoTreeString(oss, indent);
1692     m_right->intoTreeString(oss, indent);
1693     break;
1694     case G_UNARY:
1695     case G_UNARY_P:
1696     case G_NP1OUT:
1697     case G_NP1OUT_P:
1698     case G_REDUCTION:
1699     oss << opToString(m_op) << endl;
1700     indent+='.';
1701     m_left->intoTreeString(oss, indent);
1702     break;
1703     case G_TENSORPROD:
1704     oss << opToString(m_op) << endl;
1705     indent+='.';
1706     m_left->intoTreeString(oss, indent);
1707     m_right->intoTreeString(oss, indent);
1708     break;
1709     case G_NP1OUT_2P:
1710     oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;
1711     indent+='.';
1712     m_left->intoTreeString(oss, indent);
1713     break;
1714     default:
1715     oss << "UNKNOWN";
1716     }
1717     }
1718    
1719    
1720 jfenwick 1865 DataAbstract*
1721     DataLazy::deepCopy()
1722     {
1723 jfenwick 1935 switch (getOpgroup(m_op))
1724 jfenwick 1865 {
1725 jfenwick 1935 case G_IDENTITY: return new DataLazy(m_id->deepCopy()->getPtr());
1726 jfenwick 2721 case G_UNARY:
1727     case G_REDUCTION: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1728     case G_UNARY_P: return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_tol);
1729 jfenwick 1935 case G_BINARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1730 jfenwick 2066 case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1731     case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1732 jfenwick 2721 case G_NP1OUT_P: return new DataLazy(m_left->deepCopy()->getPtr(),m_op, m_axis_offset);
1733     case G_NP1OUT_2P: return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1734 jfenwick 1935 default:
1735     throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1736 jfenwick 1865 }
1737     }
1738    
1739    
1740 jfenwick 2721
1741 jfenwick 2066 // There is no single, natural interpretation of getLength on DataLazy.
1742     // Instances of DataReady can look at the size of their vectors.
1743     // For lazy though, it could be the size the data would be if it were resolved;
1744     // or it could be some function of the lengths of the DataReady instances which
1745     // form part of the expression.
1746     // Rather than have people making assumptions, I have disabled the method.
1747 jfenwick 1865 DataTypes::ValueType::size_type
1748     DataLazy::getLength() const
1749     {
1750 jfenwick 2066 throw DataException("getLength() does not make sense for lazy data.");
1751 jfenwick 1865 }
1752    
1753    
1754     DataAbstract*
1755     DataLazy::getSlice(const DataTypes::RegionType& region) const
1756     {
1757 jfenwick 1879 throw DataException("getSlice - not implemented for Lazy objects.");
1758 jfenwick 1865 }
1759    
1760 jfenwick 1935
1761     // To do this we need to rely on our child nodes
1762 jfenwick 1865 DataTypes::ValueType::size_type
1763     DataLazy::getPointOffset(int sampleNo,
1764 jfenwick 1935 int dataPointNo)
1765     {
1766     if (m_op==IDENTITY)
1767     {
1768     return m_id->getPointOffset(sampleNo,dataPointNo);
1769     }
1770     if (m_readytype!='E')
1771     {
1772     collapse();
1773     return m_id->getPointOffset(sampleNo,dataPointNo);
1774     }
1775     // at this point we do not have an identity node and the expression will be Expanded
1776     // so we only need to know which child to ask
1777     if (m_left->m_readytype=='E')
1778     {
1779     return m_left->getPointOffset(sampleNo,dataPointNo);
1780     }
1781     else
1782     {
1783     return m_right->getPointOffset(sampleNo,dataPointNo);
1784     }
1785     }
1786    
1787     // To do this we need to rely on our child nodes
1788     DataTypes::ValueType::size_type
1789     DataLazy::getPointOffset(int sampleNo,
1790 jfenwick 1865 int dataPointNo) const
1791     {
1792 jfenwick 1935 if (m_op==IDENTITY)
1793     {
1794     return m_id->getPointOffset(sampleNo,dataPointNo);
1795     }
1796     if (m_readytype=='E')
1797     {
1798     // at this point we do not have an identity node and the expression will be Expanded
1799     // so we only need to know which child to ask
1800     if (m_left->m_readytype=='E')
1801     {
1802     return m_left->getPointOffset(sampleNo,dataPointNo);
1803     }
1804     else
1805     {
1806     return m_right->getPointOffset(sampleNo,dataPointNo);
1807     }
1808     }
1809     if (m_readytype=='C')
1810     {
1811     return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
1812     }
1813     throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1814 jfenwick 1865 }
1815    
1816 jfenwick 2271
1817     // I have decided to let Data:: handle this issue.
1818 jfenwick 1901 void
1819     DataLazy::setToZero()
1820     {
1821 jfenwick 2271 // DataTypes::ValueType v(getNoValues(),0);
1822     // m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1823     // m_op=IDENTITY;
1824     // m_right.reset();
1825     // m_left.reset();
1826     // m_readytype='C';
1827     // m_buffsRequired=1;
1828    
1829 jfenwick 2514 privdebug=privdebug; // to stop the compiler complaining about unused privdebug
1830 jfenwick 2271 throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
1831 jfenwick 1901 }
1832    
1833 jfenwick 2271 bool
1834     DataLazy::actsExpanded() const
1835     {
1836     return (m_readytype=='E');
1837     }
1838    
1839 jfenwick 1865 } // end namespace

  ViewVC Help
Powered by ViewVC 1.1.26