/[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 2824 - (hide annotations)
Wed Dec 16 01:22:56 2009 UTC (9 years, 3 months ago) by jfenwick
Original Path: trunk/escript/src/DataLazy.cpp
File size: 57953 byte(s)
ResolveGroup fixed so it actually resolves groups.
Children, when you iterate backwards remember its i>=0 not i>0

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

  ViewVC Help
Powered by ViewVC 1.1.26