/[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 2199 - (hide annotations)
Thu Jan 8 06:10:52 2009 UTC (10 years, 3 months ago) by jfenwick
Original Path: trunk/escript/src/DataLazy.cpp
File size: 55552 byte(s)
Misc fixes:
Added some svn:ignore properties for output files that were cluttering things up.

Lazy fixes:
Fixed shape calculations for TRACE and TRANSPOSE for rank>2.
Adjusted unit test accordingly.

As a Temporary change to DataC.cpp to test for lazy data in DataC's expanded check.
This is wrong but would only affect people using lazy data.
The proper fix will come when the numarray removal code moves over from the branch.

Made tensor product AUTOLAZY capable.
Fixed some bugs resolving tensor products (incorrect offsets in buffers).
Macro'd some stray couts.

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

  ViewVC Help
Powered by ViewVC 1.1.26