/[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 2271 - (hide annotations)
Mon Feb 16 05:08:29 2009 UTC (10 years, 1 month ago) by jfenwick
Original Path: trunk/escript/src/DataLazy.cpp
File size: 55741 byte(s)
Merging version 2269 to trunk

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 2271 // #define SIZELIMIT
44 jfenwick 2177 // #define SIZELIMIT if ((m_height>7) || (m_children>127)) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl; resolveToIdentity();}
45 jfenwick 2271 #define SIZELIMIT if ((m_height>7) || (m_children>127)) {resolveToIdentity();}
46 jfenwick 2177
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 2271 // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefts.
352 jfenwick 2157 // 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 2271 p->makeLazyShared();
401 jfenwick 2177 DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
402     makeIdentity(dr);
403 jfenwick 2199 LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
404 jfenwick 1879 }
405 jfenwick 2092 LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
406 jfenwick 1865 }
407    
408 jfenwick 1899
409    
410 jfenwick 1901
411 jfenwick 1886 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
412     : parent(left->getFunctionSpace(),left->getShape()),
413 jfenwick 2066 m_op(op),
414     m_axis_offset(0),
415     m_transpose(0),
416     m_SL(0), m_SM(0), m_SR(0)
417 jfenwick 1886 {
418 jfenwick 2037 if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
419 jfenwick 1886 {
420     throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
421     }
422 jfenwick 2066
423 jfenwick 1886 DataLazy_ptr lleft;
424     if (!left->isLazy())
425     {
426     lleft=DataLazy_ptr(new DataLazy(left));
427     }
428     else
429     {
430     lleft=dynamic_pointer_cast<DataLazy>(left);
431     }
432 jfenwick 1889 m_readytype=lleft->m_readytype;
433 jfenwick 1886 m_left=lleft;
434 jfenwick 2037 m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
435 jfenwick 1886 m_samplesize=getNumDPPSample()*getNoValues();
436 jfenwick 2066 m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
437 jfenwick 2177 m_children=m_left->m_children+1;
438     m_height=m_left->m_height+1;
439     SIZELIMIT
440 jfenwick 1886 }
441    
442    
443 jfenwick 1943 // In this constructor we need to consider interpolation
444 jfenwick 1879 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
445     : parent(resultFS(left,right,op), resultShape(left,right,op)),
446 jfenwick 2066 m_op(op),
447     m_SL(0), m_SM(0), m_SR(0)
448 jfenwick 1879 {
449 jfenwick 2199 LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
450 jfenwick 2037 if ((getOpgroup(op)!=G_BINARY))
451 jfenwick 1886 {
452 jfenwick 1889 throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
453 jfenwick 1886 }
454 jfenwick 1943
455     if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
456     {
457     FunctionSpace fs=getFunctionSpace();
458     Data ltemp(left);
459     Data tmp(ltemp,fs);
460     left=tmp.borrowDataPtr();
461     }
462 jfenwick 2066 if (getFunctionSpace()!=right->getFunctionSpace()) // right needs to be interpolated
463 jfenwick 1943 {
464     Data tmp(Data(right),getFunctionSpace());
465     right=tmp.borrowDataPtr();
466 jfenwick 2199 LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
467 jfenwick 1943 }
468     left->operandCheck(*right);
469    
470 jfenwick 1899 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
471 jfenwick 1879 {
472     m_left=dynamic_pointer_cast<DataLazy>(left);
473 jfenwick 2199 LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
474 jfenwick 1879 }
475     else
476     {
477     m_left=DataLazy_ptr(new DataLazy(left));
478 jfenwick 2199 LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
479 jfenwick 1879 }
480     if (right->isLazy())
481     {
482     m_right=dynamic_pointer_cast<DataLazy>(right);
483 jfenwick 2199 LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
484 jfenwick 1879 }
485     else
486     {
487     m_right=DataLazy_ptr(new DataLazy(right));
488 jfenwick 2199 LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
489 jfenwick 1879 }
490 jfenwick 1889 char lt=m_left->m_readytype;
491     char rt=m_right->m_readytype;
492     if (lt=='E' || rt=='E')
493     {
494     m_readytype='E';
495     }
496     else if (lt=='T' || rt=='T')
497     {
498     m_readytype='T';
499     }
500     else
501     {
502     m_readytype='C';
503     }
504 jfenwick 2066 m_samplesize=getNumDPPSample()*getNoValues();
505     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());
506 jfenwick 1879 m_buffsRequired=calcBuffs(m_left, m_right,m_op);
507 jfenwick 2177 m_children=m_left->m_children+m_right->m_children+2;
508     m_height=max(m_left->m_height,m_right->m_height)+1;
509     SIZELIMIT
510 jfenwick 2092 LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
511 jfenwick 1879 }
512    
513 jfenwick 2066 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
514     : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
515     m_op(op),
516     m_axis_offset(axis_offset),
517     m_transpose(transpose)
518     {
519     if ((getOpgroup(op)!=G_TENSORPROD))
520     {
521     throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
522     }
523     if ((transpose>2) || (transpose<0))
524     {
525     throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
526     }
527     if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
528     {
529     FunctionSpace fs=getFunctionSpace();
530     Data ltemp(left);
531     Data tmp(ltemp,fs);
532     left=tmp.borrowDataPtr();
533     }
534     if (getFunctionSpace()!=right->getFunctionSpace()) // right needs to be interpolated
535     {
536     Data tmp(Data(right),getFunctionSpace());
537     right=tmp.borrowDataPtr();
538     }
539 jfenwick 2195 // left->operandCheck(*right);
540 jfenwick 1879
541 jfenwick 2066 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
542     {
543     m_left=dynamic_pointer_cast<DataLazy>(left);
544     }
545     else
546     {
547     m_left=DataLazy_ptr(new DataLazy(left));
548     }
549     if (right->isLazy())
550     {
551     m_right=dynamic_pointer_cast<DataLazy>(right);
552     }
553     else
554     {
555     m_right=DataLazy_ptr(new DataLazy(right));
556     }
557     char lt=m_left->m_readytype;
558     char rt=m_right->m_readytype;
559     if (lt=='E' || rt=='E')
560     {
561     m_readytype='E';
562     }
563     else if (lt=='T' || rt=='T')
564     {
565     m_readytype='T';
566     }
567     else
568     {
569     m_readytype='C';
570     }
571     m_samplesize=getNumDPPSample()*getNoValues();
572     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());
573     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
574 jfenwick 2177 m_children=m_left->m_children+m_right->m_children+2;
575     m_height=max(m_left->m_height,m_right->m_height)+1;
576     SIZELIMIT
577 jfenwick 2092 LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
578 jfenwick 2066 }
579    
580    
581 jfenwick 2084 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
582 jfenwick 2199 : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
583 jfenwick 2084 m_op(op),
584     m_axis_offset(axis_offset),
585 jfenwick 2147 m_transpose(0),
586     m_tol(0)
587 jfenwick 2084 {
588     if ((getOpgroup(op)!=G_NP1OUT_P))
589     {
590     throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
591     }
592     DataLazy_ptr lleft;
593     if (!left->isLazy())
594     {
595     lleft=DataLazy_ptr(new DataLazy(left));
596     }
597     else
598     {
599     lleft=dynamic_pointer_cast<DataLazy>(left);
600     }
601     m_readytype=lleft->m_readytype;
602     m_left=lleft;
603     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
604     m_samplesize=getNumDPPSample()*getNoValues();
605     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
606 jfenwick 2177 m_children=m_left->m_children+1;
607     m_height=m_left->m_height+1;
608     SIZELIMIT
609 jfenwick 2092 LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
610 jfenwick 2084 }
611    
612 jfenwick 2147 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
613     : parent(left->getFunctionSpace(), left->getShape()),
614     m_op(op),
615     m_axis_offset(0),
616     m_transpose(0),
617     m_tol(tol)
618     {
619     if ((getOpgroup(op)!=G_UNARY_P))
620     {
621     throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
622     }
623     DataLazy_ptr lleft;
624     if (!left->isLazy())
625     {
626     lleft=DataLazy_ptr(new DataLazy(left));
627     }
628     else
629     {
630     lleft=dynamic_pointer_cast<DataLazy>(left);
631     }
632     m_readytype=lleft->m_readytype;
633     m_left=lleft;
634     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
635     m_samplesize=getNumDPPSample()*getNoValues();
636     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
637 jfenwick 2177 m_children=m_left->m_children+1;
638     m_height=m_left->m_height+1;
639     SIZELIMIT
640 jfenwick 2147 LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
641     }
642 jfenwick 2084
643 jfenwick 1865 DataLazy::~DataLazy()
644     {
645     }
646    
647 jfenwick 1879
648     int
649     DataLazy::getBuffsRequired() const
650 jfenwick 1865 {
651 jfenwick 1879 return m_buffsRequired;
652     }
653    
654 jfenwick 1884
655 jfenwick 2066 size_t
656     DataLazy::getMaxSampleSize() const
657     {
658     return m_maxsamplesize;
659     }
660    
661 jfenwick 2271
662    
663     size_t
664     DataLazy::getSampleBufferSize() const
665     {
666     return m_maxsamplesize*(max(1,m_buffsRequired));
667     }
668    
669 jfenwick 1899 /*
670     \brief Evaluates the expression using methods on Data.
671     This does the work for the collapse method.
672     For reasons of efficiency do not call this method on DataExpanded nodes.
673     */
674 jfenwick 1889 DataReady_ptr
675     DataLazy::collapseToReady()
676     {
677     if (m_readytype=='E')
678     { // this is more an efficiency concern than anything else
679     throw DataException("Programmer Error - do not use collapse on Expanded data.");
680     }
681     if (m_op==IDENTITY)
682     {
683     return m_id;
684     }
685     DataReady_ptr pleft=m_left->collapseToReady();
686     Data left(pleft);
687     Data right;
688 jfenwick 2066 if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
689 jfenwick 1889 {
690     right=Data(m_right->collapseToReady());
691     }
692     Data result;
693     switch(m_op)
694     {
695     case ADD:
696     result=left+right;
697     break;
698     case SUB:
699     result=left-right;
700     break;
701     case MUL:
702     result=left*right;
703     break;
704     case DIV:
705     result=left/right;
706     break;
707     case SIN:
708     result=left.sin();
709     break;
710     case COS:
711     result=left.cos();
712     break;
713     case TAN:
714     result=left.tan();
715     break;
716     case ASIN:
717     result=left.asin();
718     break;
719     case ACOS:
720     result=left.acos();
721     break;
722     case ATAN:
723     result=left.atan();
724     break;
725     case SINH:
726     result=left.sinh();
727     break;
728     case COSH:
729     result=left.cosh();
730     break;
731     case TANH:
732     result=left.tanh();
733     break;
734     case ERF:
735     result=left.erf();
736     break;
737     case ASINH:
738     result=left.asinh();
739     break;
740     case ACOSH:
741     result=left.acosh();
742     break;
743     case ATANH:
744     result=left.atanh();
745     break;
746     case LOG10:
747     result=left.log10();
748     break;
749     case LOG:
750     result=left.log();
751     break;
752     case SIGN:
753     result=left.sign();
754     break;
755     case ABS:
756     result=left.abs();
757     break;
758     case NEG:
759     result=left.neg();
760     break;
761     case POS:
762     // it doesn't mean anything for delayed.
763     // it will just trigger a deep copy of the lazy object
764     throw DataException("Programmer error - POS not supported for lazy data.");
765     break;
766     case EXP:
767     result=left.exp();
768     break;
769     case SQRT:
770     result=left.sqrt();
771     break;
772     case RECIP:
773     result=left.oneOver();
774     break;
775     case GZ:
776     result=left.wherePositive();
777     break;
778     case LZ:
779     result=left.whereNegative();
780     break;
781     case GEZ:
782     result=left.whereNonNegative();
783     break;
784     case LEZ:
785     result=left.whereNonPositive();
786     break;
787 jfenwick 2147 case NEZ:
788     result=left.whereNonZero(m_tol);
789     break;
790     case EZ:
791     result=left.whereZero(m_tol);
792     break;
793 jfenwick 2037 case SYM:
794     result=left.symmetric();
795     break;
796     case NSYM:
797     result=left.nonsymmetric();
798     break;
799 jfenwick 2066 case PROD:
800     result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
801     break;
802 jfenwick 2084 case TRANS:
803     result=left.transpose(m_axis_offset);
804     break;
805     case TRACE:
806     result=left.trace(m_axis_offset);
807     break;
808 jfenwick 1889 default:
809 jfenwick 2037 throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
810 jfenwick 1889 }
811     return result.borrowReadyPtr();
812     }
813    
814 jfenwick 1899 /*
815     \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
816     This method uses the original methods on the Data class to evaluate the expressions.
817     For this reason, it should not be used on DataExpanded instances. (To do so would defeat
818     the purpose of using DataLazy in the first place).
819     */
820 jfenwick 1889 void
821     DataLazy::collapse()
822     {
823     if (m_op==IDENTITY)
824     {
825     return;
826     }
827     if (m_readytype=='E')
828     { // this is more an efficiency concern than anything else
829     throw DataException("Programmer Error - do not use collapse on Expanded data.");
830     }
831     m_id=collapseToReady();
832     m_op=IDENTITY;
833     }
834    
835 jfenwick 1899 /*
836 jfenwick 2037 \brief Compute the value of the expression (unary operation) for the given sample.
837 jfenwick 1899 \return Vector which stores the value of the subexpression for the given sample.
838     \param v A vector to store intermediate results.
839     \param offset Index in v to begin storing results.
840     \param sampleNo Sample number to evaluate.
841     \param roffset (output parameter) the offset in the return vector where the result begins.
842    
843     The return value will be an existing vector so do not deallocate it.
844     If the result is stored in v it should be stored at the offset given.
845     Everything from offset to the end of v should be considered available for this method to use.
846     */
847 jfenwick 1898 DataTypes::ValueType*
848     DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
849 jfenwick 1889 {
850     // we assume that any collapsing has been done before we get here
851     // since we only have one argument we don't need to think about only
852     // processing single points.
853     if (m_readytype!='E')
854     {
855     throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
856     }
857 jfenwick 1898 const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
858     const double* left=&((*vleft)[roffset]);
859 jfenwick 1889 double* result=&(v[offset]);
860 jfenwick 1898 roffset=offset;
861 jfenwick 1889 switch (m_op)
862     {
863     case SIN:
864 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
865 jfenwick 1889 break;
866     case COS:
867 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
868 jfenwick 1889 break;
869     case TAN:
870 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
871 jfenwick 1889 break;
872     case ASIN:
873 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
874 jfenwick 1889 break;
875     case ACOS:
876 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
877 jfenwick 1889 break;
878     case ATAN:
879 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
880 jfenwick 1889 break;
881     case SINH:
882 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
883 jfenwick 1889 break;
884     case COSH:
885 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
886 jfenwick 1889 break;
887     case TANH:
888 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
889 jfenwick 1889 break;
890     case ERF:
891 phornby 2050 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
892 jfenwick 1889 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
893     #else
894     tensor_unary_operation(m_samplesize, left, result, ::erf);
895     break;
896     #endif
897     case ASINH:
898 phornby 2050 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
899 jfenwick 1889 tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
900     #else
901     tensor_unary_operation(m_samplesize, left, result, ::asinh);
902     #endif
903     break;
904     case ACOSH:
905 phornby 2050 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
906 jfenwick 1889 tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
907     #else
908     tensor_unary_operation(m_samplesize, left, result, ::acosh);
909     #endif
910     break;
911     case ATANH:
912 phornby 2050 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
913 jfenwick 1889 tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
914     #else
915     tensor_unary_operation(m_samplesize, left, result, ::atanh);
916     #endif
917     break;
918     case LOG10:
919 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
920 jfenwick 1889 break;
921     case LOG:
922 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
923 jfenwick 1889 break;
924     case SIGN:
925     tensor_unary_operation(m_samplesize, left, result, escript::fsign);
926     break;
927     case ABS:
928 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
929 jfenwick 1889 break;
930     case NEG:
931     tensor_unary_operation(m_samplesize, left, result, negate<double>());
932     break;
933     case POS:
934     // it doesn't mean anything for delayed.
935     // it will just trigger a deep copy of the lazy object
936     throw DataException("Programmer error - POS not supported for lazy data.");
937     break;
938     case EXP:
939 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
940 jfenwick 1889 break;
941     case SQRT:
942 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
943 jfenwick 1889 break;
944     case RECIP:
945     tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
946     break;
947     case GZ:
948     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
949     break;
950     case LZ:
951     tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
952     break;
953     case GEZ:
954     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
955     break;
956     case LEZ:
957     tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
958     break;
959 jfenwick 2147 // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
960     case NEZ:
961     tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
962     break;
963     case EZ:
964     tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
965     break;
966 jfenwick 1889
967     default:
968     throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
969     }
970 jfenwick 1898 return &v;
971 jfenwick 1889 }
972    
973    
974 jfenwick 2147
975    
976    
977    
978 jfenwick 2037 /*
979     \brief Compute the value of the expression (unary operation) for the given sample.
980     \return Vector which stores the value of the subexpression for the given sample.
981     \param v A vector to store intermediate results.
982     \param offset Index in v to begin storing results.
983     \param sampleNo Sample number to evaluate.
984     \param roffset (output parameter) the offset in the return vector where the result begins.
985 jfenwick 1898
986 jfenwick 2037 The return value will be an existing vector so do not deallocate it.
987     If the result is stored in v it should be stored at the offset given.
988     Everything from offset to the end of v should be considered available for this method to use.
989     */
990     DataTypes::ValueType*
991     DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
992     {
993     // we assume that any collapsing has been done before we get here
994     // since we only have one argument we don't need to think about only
995     // processing single points.
996     if (m_readytype!='E')
997     {
998     throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
999     }
1000     // since we can't write the result over the input, we need a result offset further along
1001     size_t subroffset=roffset+m_samplesize;
1002 jfenwick 2157 LAZYDEBUG(cerr << "subroffset=" << subroffset << endl;)
1003     const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1004 jfenwick 2037 roffset=offset;
1005 jfenwick 2157 size_t loop=0;
1006     size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1007     size_t step=getNoValues();
1008 jfenwick 2037 switch (m_op)
1009     {
1010     case SYM:
1011 jfenwick 2157 for (loop=0;loop<numsteps;++loop)
1012     {
1013     DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
1014     subroffset+=step;
1015     offset+=step;
1016     }
1017 jfenwick 2037 break;
1018     case NSYM:
1019 jfenwick 2157 for (loop=0;loop<numsteps;++loop)
1020     {
1021     DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
1022     subroffset+=step;
1023     offset+=step;
1024     }
1025 jfenwick 2037 break;
1026     default:
1027     throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1028     }
1029     return &v;
1030     }
1031 jfenwick 1898
1032 jfenwick 2084 /*
1033     \brief Compute the value of the expression (unary operation) for the given sample.
1034     \return Vector which stores the value of the subexpression for the given sample.
1035     \param v A vector to store intermediate results.
1036     \param offset Index in v to begin storing results.
1037     \param sampleNo Sample number to evaluate.
1038     \param roffset (output parameter) the offset in the return vector where the result begins.
1039 jfenwick 1899
1040 jfenwick 2084 The return value will be an existing vector so do not deallocate it.
1041     If the result is stored in v it should be stored at the offset given.
1042     Everything from offset to the end of v should be considered available for this method to use.
1043     */
1044     DataTypes::ValueType*
1045     DataLazy::resolveNP1OUT_P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1046     {
1047     // we assume that any collapsing has been done before we get here
1048     // since we only have one argument we don't need to think about only
1049     // processing single points.
1050     if (m_readytype!='E')
1051     {
1052     throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");
1053     }
1054     // since we can't write the result over the input, we need a result offset further along
1055 jfenwick 2157 size_t subroffset;
1056     const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1057 jfenwick 2166 LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)
1058     LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)
1059 jfenwick 2084 roffset=offset;
1060 jfenwick 2157 size_t loop=0;
1061     size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1062 jfenwick 2166 size_t outstep=getNoValues();
1063     size_t instep=m_left->getNoValues();
1064     LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)
1065 jfenwick 2084 switch (m_op)
1066     {
1067     case TRACE:
1068 jfenwick 2157 for (loop=0;loop<numsteps;++loop)
1069     {
1070 jfenwick 2166 size_t zz=sampleNo*getNumDPPSample()+loop;
1071     if (zz==5800)
1072     {
1073     LAZYDEBUG(cerr << "point=" << zz<< endl;)
1074     LAZYDEBUG(cerr << "Input to trace=" << DataTypes::pointToString(*vleft,m_left->getShape(),subroffset,"") << endl;)
1075     LAZYDEBUG(cerr << "Offset for point=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << " vs ";)
1076     LAZYDEBUG(cerr << subroffset << endl;)
1077     LAZYDEBUG(cerr << "output=" << offset << endl;)
1078     }
1079     DataMaths::trace(*vleft,m_left->getShape(),subroffset, v ,getShape(),offset,m_axis_offset);
1080     if (zz==5800)
1081     {
1082     LAZYDEBUG(cerr << "Result of trace=" << DataTypes::pointToString(v,getShape(),offset,"") << endl;)
1083     }
1084     subroffset+=instep;
1085     offset+=outstep;
1086 jfenwick 2157 }
1087 jfenwick 2084 break;
1088     case TRANS:
1089 jfenwick 2157 for (loop=0;loop<numsteps;++loop)
1090     {
1091     DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
1092 jfenwick 2166 subroffset+=instep;
1093     offset+=outstep;
1094 jfenwick 2157 }
1095 jfenwick 2084 break;
1096     default:
1097     throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1098     }
1099     return &v;
1100     }
1101 jfenwick 2037
1102    
1103 phornby 1987 #define PROC_OP(TYPE,X) \
1104 jfenwick 2152 for (int j=0;j<onumsteps;++j)\
1105     {\
1106     for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1107     { \
1108     LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1109 jfenwick 2157 LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1110 jfenwick 2152 tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1111 jfenwick 2157 LAZYDEBUG(cout << " result= " << resultp[0] << endl;) \
1112 jfenwick 2152 lroffset+=leftstep; \
1113     rroffset+=rightstep; \
1114     }\
1115     lroffset+=oleftstep;\
1116     rroffset+=orightstep;\
1117 jfenwick 1889 }
1118    
1119 jfenwick 1899 /*
1120     \brief Compute the value of the expression (binary operation) for the given sample.
1121     \return Vector which stores the value of the subexpression for the given sample.
1122     \param v A vector to store intermediate results.
1123     \param offset Index in v to begin storing results.
1124     \param sampleNo Sample number to evaluate.
1125     \param roffset (output parameter) the offset in the return vector where the result begins.
1126    
1127     The return value will be an existing vector so do not deallocate it.
1128     If the result is stored in v it should be stored at the offset given.
1129     Everything from offset to the end of v should be considered available for this method to use.
1130     */
1131     // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1132     // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1133     // If both children are expanded, then we can process them in a single operation (we treat
1134     // the whole sample as one big datapoint.
1135     // If one of the children is not expanded, then we need to treat each point in the sample
1136     // individually.
1137 jfenwick 1908 // There is an additional complication when scalar operations are considered.
1138     // For example, 2+Vector.
1139     // In this case each double within the point is treated individually
1140 jfenwick 1898 DataTypes::ValueType*
1141 jfenwick 1899 DataLazy::resolveBinary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1142 jfenwick 1889 {
1143 jfenwick 2092 LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1144 jfenwick 1889
1145 jfenwick 1899 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
1146     // first work out which of the children are expanded
1147 jfenwick 1889 bool leftExp=(m_left->m_readytype=='E');
1148     bool rightExp=(m_right->m_readytype=='E');
1149 jfenwick 2084 if (!leftExp && !rightExp)
1150     {
1151     throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1152     }
1153     bool leftScalar=(m_left->getRank()==0);
1154     bool rightScalar=(m_right->getRank()==0);
1155 jfenwick 2152 if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1156 jfenwick 1908 {
1157 jfenwick 2152 throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1158     }
1159     size_t leftsize=m_left->getNoValues();
1160     size_t rightsize=m_right->getNoValues();
1161     size_t chunksize=1; // how many doubles will be processed in one go
1162     int leftstep=0; // how far should the left offset advance after each step
1163     int rightstep=0;
1164     int numsteps=0; // total number of steps for the inner loop
1165     int oleftstep=0; // the o variables refer to the outer loop
1166     int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1167     int onumsteps=1;
1168    
1169     bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1170     bool RES=(rightExp && rightScalar);
1171     bool LS=(!leftExp && leftScalar); // left is a single scalar
1172     bool RS=(!rightExp && rightScalar);
1173     bool LN=(!leftExp && !leftScalar); // left is a single non-scalar
1174     bool RN=(!rightExp && !rightScalar);
1175     bool LEN=(leftExp && !leftScalar); // left is an expanded non-scalar
1176     bool REN=(rightExp && !rightScalar);
1177    
1178     if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1179     {
1180     chunksize=m_left->getNumDPPSample()*leftsize;
1181     leftstep=0;
1182     rightstep=0;
1183     numsteps=1;
1184     }
1185     else if (LES || RES)
1186     {
1187     chunksize=1;
1188     if (LES) // left is an expanded scalar
1189 jfenwick 2084 {
1190 jfenwick 2152 if (RS)
1191     {
1192     leftstep=1;
1193     rightstep=0;
1194     numsteps=m_left->getNumDPPSample();
1195     }
1196     else // RN or REN
1197     {
1198     leftstep=0;
1199     oleftstep=1;
1200     rightstep=1;
1201 phornby 2171 orightstep=(RN ? -(int)rightsize : 0);
1202 jfenwick 2152 numsteps=rightsize;
1203     onumsteps=m_left->getNumDPPSample();
1204     }
1205 jfenwick 2084 }
1206 jfenwick 2152 else // right is an expanded scalar
1207     {
1208     if (LS)
1209     {
1210     rightstep=1;
1211     leftstep=0;
1212     numsteps=m_right->getNumDPPSample();
1213     }
1214     else
1215     {
1216     rightstep=0;
1217     orightstep=1;
1218     leftstep=1;
1219 jfenwick 2177 oleftstep=(LN ? -(int)leftsize : 0);
1220 jfenwick 2152 numsteps=leftsize;
1221     onumsteps=m_right->getNumDPPSample();
1222     }
1223     }
1224     }
1225     else // this leaves (LEN, RS), (LEN, RN) and their transposes
1226     {
1227     if (LEN) // and Right will be a single value
1228     {
1229     chunksize=rightsize;
1230     leftstep=rightsize;
1231     rightstep=0;
1232     numsteps=m_left->getNumDPPSample();
1233     if (RS)
1234     {
1235     numsteps*=leftsize;
1236     }
1237     }
1238     else // REN
1239     {
1240     chunksize=leftsize;
1241     rightstep=leftsize;
1242     leftstep=0;
1243     numsteps=m_right->getNumDPPSample();
1244     if (LS)
1245     {
1246     numsteps*=rightsize;
1247     }
1248     }
1249     }
1250    
1251     int resultStep=max(leftstep,rightstep); // only one (at most) should be !=0
1252 jfenwick 1899 // Get the values of sub-expressions
1253 jfenwick 2157 const ValueType* left=m_left->resolveSample(v,offset+getMaxSampleSize(),sampleNo,lroffset); // see note on
1254     // calcBufss for why we can't put offset as the 2nd param above
1255     const ValueType* right=m_right->resolveSample(v,offset+2*getMaxSampleSize(),sampleNo,rroffset); // Note
1256 jfenwick 1899 // the right child starts further along.
1257 jfenwick 2152 LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1258     LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1259     LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1260     LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1261     LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1262     LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1263     LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN << endl;)
1264 jfenwick 2195
1265 jfenwick 2199
1266 jfenwick 1899 double* resultp=&(v[offset]); // results are stored at the vector offset we recieved
1267 jfenwick 1889 switch(m_op)
1268     {
1269     case ADD:
1270 phornby 2021 PROC_OP(NO_ARG,plus<double>());
1271 jfenwick 1889 break;
1272 jfenwick 1899 case SUB:
1273 phornby 2021 PROC_OP(NO_ARG,minus<double>());
1274 jfenwick 1899 break;
1275     case MUL:
1276 phornby 2021 PROC_OP(NO_ARG,multiplies<double>());
1277 jfenwick 1899 break;
1278     case DIV:
1279 phornby 2021 PROC_OP(NO_ARG,divides<double>());
1280 jfenwick 1899 break;
1281 jfenwick 1910 case POW:
1282 phornby 1994 PROC_OP(double (double,double),::pow);
1283 jfenwick 1910 break;
1284 jfenwick 1889 default:
1285 jfenwick 1898 throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1286 jfenwick 1889 }
1287 jfenwick 2195 roffset=offset;
1288 jfenwick 1898 return &v;
1289 jfenwick 1889 }
1290    
1291 jfenwick 1898
1292 jfenwick 2152
1293 jfenwick 2066 /*
1294     \brief Compute the value of the expression (tensor product) for the given sample.
1295     \return Vector which stores the value of the subexpression for the given sample.
1296     \param v A vector to store intermediate results.
1297     \param offset Index in v to begin storing results.
1298     \param sampleNo Sample number to evaluate.
1299     \param roffset (output parameter) the offset in the return vector where the result begins.
1300 jfenwick 1898
1301 jfenwick 2066 The return value will be an existing vector so do not deallocate it.
1302     If the result is stored in v it should be stored at the offset given.
1303     Everything from offset to the end of v should be considered available for this method to use.
1304     */
1305     // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1306     // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1307     // unlike the other resolve helpers, we must treat these datapoints separately.
1308     DataTypes::ValueType*
1309     DataLazy::resolveTProd(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1310     {
1311 jfenwick 2195 LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << " to offset " << offset<< endl;)
1312 jfenwick 2066
1313     size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
1314     // first work out which of the children are expanded
1315     bool leftExp=(m_left->m_readytype=='E');
1316     bool rightExp=(m_right->m_readytype=='E');
1317     int steps=getNumDPPSample();
1318 jfenwick 2195 /* int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
1319     int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);*/
1320     int leftStep=(leftExp? m_left->getNoValues() : 0); // do not have scalars as input to this method
1321     int rightStep=(rightExp?m_right->getNoValues() : 0);
1322    
1323     int resultStep=getNoValues();
1324 jfenwick 2066 // Get the values of sub-expressions (leave a gap of one sample for the result).
1325 jfenwick 2199 int gap=offset+m_samplesize;
1326 jfenwick 2195
1327     LAZYDEBUG(cout << "Query left with offset=" << gap << endl;)
1328    
1329 jfenwick 2153 const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);
1330 jfenwick 2199 gap+=m_left->getMaxSampleSize();
1331 jfenwick 2195
1332    
1333     LAZYDEBUG(cout << "Query right with offset=" << gap << endl;)
1334    
1335    
1336     const ValueType* right=m_right->resolveSample(v,gap,sampleNo,rroffset);
1337    
1338     LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) << endl;
1339     cout << getNoValues() << endl;)
1340     LAZYDEBUG(cerr << "Result of left=";)
1341     LAZYDEBUG(cerr << "[" << lroffset << " .. " << lroffset+m_left->getNoValues() << "]" << endl;
1342 jfenwick 2199
1343     for (int i=lroffset, limit=lroffset+(leftExp?m_left->getNoValues()*m_left->getNumDPPSample():m_left->getNoValues());i<limit;++i)
1344 jfenwick 2195 {
1345 jfenwick 2199 cout << "[" << setw(2) << i-lroffset << "] " << setw(10) << (*left)[i] << " ";
1346     if (i%4==0) cout << endl;
1347 jfenwick 2195 })
1348     LAZYDEBUG(cerr << "\nResult of right=" << endl;)
1349 jfenwick 2199 LAZYDEBUG(
1350     for (int i=rroffset, limit=rroffset+(rightExp?m_right->getNoValues()*m_right->getNumDPPSample():m_right->getNoValues());i<limit;++i)
1351 jfenwick 2195 {
1352 jfenwick 2199 cerr << "[" << setw(2)<< i-rroffset << "] " << setw(10) << (*right)[i] << " ";
1353     if (i%4==0) cout << endl;
1354 jfenwick 2195 }
1355     cerr << endl;
1356     )
1357     LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1358 jfenwick 2153 LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1359     LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1360     LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1361     LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1362     LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1363 jfenwick 2195 LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1364    
1365 jfenwick 2066 double* resultp=&(v[offset]); // results are stored at the vector offset we recieved
1366     switch(m_op)
1367     {
1368     case PROD:
1369     for (int i=0;i<steps;++i,resultp+=resultStep)
1370     {
1371 jfenwick 2199
1372 jfenwick 2153 LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)
1373     LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)
1374     LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)
1375 jfenwick 2199
1376 jfenwick 2066 const double *ptr_0 = &((*left)[lroffset]);
1377     const double *ptr_1 = &((*right)[rroffset]);
1378 jfenwick 2199
1379 jfenwick 2153 LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1380     LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1381 jfenwick 2199
1382 jfenwick 2066 matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
1383 jfenwick 2199
1384 jfenwick 2195 LAZYDEBUG(cout << "Results--\n";
1385 jfenwick 2199 {
1386     DataVector dv(getNoValues());
1387 jfenwick 2195 for (int z=0;z<getNoValues();++z)
1388     {
1389 jfenwick 2199 cout << "[" << setw(2) << z<< "] " << setw(10) << resultp[z] << " ";
1390     if (z%4==0) cout << endl;
1391     dv[z]=resultp[z];
1392 jfenwick 2195 }
1393 jfenwick 2199 cout << endl << DataTypes::pointToString(dv,getShape(),0,"RESLT");
1394 jfenwick 2195 cout << "\nWritten to: " << resultp << " resultStep=" << resultStep << endl;
1395 jfenwick 2199 }
1396 jfenwick 2195 )
1397 jfenwick 2066 lroffset+=leftStep;
1398     rroffset+=rightStep;
1399     }
1400     break;
1401     default:
1402     throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1403     }
1404     roffset=offset;
1405     return &v;
1406     }
1407    
1408    
1409    
1410 jfenwick 1899 /*
1411     \brief Compute the value of the expression for the given sample.
1412     \return Vector which stores the value of the subexpression for the given sample.
1413     \param v A vector to store intermediate results.
1414     \param offset Index in v to begin storing results.
1415     \param sampleNo Sample number to evaluate.
1416     \param roffset (output parameter) the offset in the return vector where the result begins.
1417 jfenwick 1898
1418 jfenwick 1899 The return value will be an existing vector so do not deallocate it.
1419     */
1420 jfenwick 1884 // the vector and the offset are a place where the method could write its data if it wishes
1421     // it is not obligated to do so. For example, if it has its own storage already, it can use that.
1422 jfenwick 1885 // Hence the return value to indicate where the data is actually stored.
1423     // Regardless, the storage should be assumed to be used, even if it isn't.
1424 jfenwick 1898
1425     // the roffset is the offset within the returned vector where the data begins
1426     const DataTypes::ValueType*
1427     DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
1428 jfenwick 1879 {
1429 jfenwick 2092 LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1430 jfenwick 1889 // collapse so we have a 'E' node or an IDENTITY for some other type
1431     if (m_readytype!='E' && m_op!=IDENTITY)
1432     {
1433     collapse();
1434     }
1435 jfenwick 1885 if (m_op==IDENTITY)
1436 jfenwick 1865 {
1437 jfenwick 2271 const ValueType& vec=m_id->getVectorRO();
1438 jfenwick 1889 if (m_readytype=='C')
1439     {
1440 jfenwick 1898 roffset=0;
1441 jfenwick 2152 LAZYDEBUG(cout << "Finish sample " << toString() << endl;)
1442 jfenwick 1898 return &(vec);
1443 jfenwick 1889 }
1444 jfenwick 1898 roffset=m_id->getPointOffset(sampleNo, 0);
1445 jfenwick 2152 LAZYDEBUG(cout << "Finish sample " << toString() << endl;)
1446 jfenwick 1898 return &(vec);
1447 jfenwick 1865 }
1448 jfenwick 1889 if (m_readytype!='E')
1449     {
1450     throw DataException("Programmer Error - Collapse did not produce an expanded node.");
1451     }
1452     switch (getOpgroup(m_op))
1453     {
1454 jfenwick 2147 case G_UNARY:
1455     case G_UNARY_P: return resolveUnary(v, offset,sampleNo,roffset);
1456 jfenwick 1898 case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
1457 jfenwick 2066 case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
1458 jfenwick 2084 case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);
1459 jfenwick 2066 case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
1460 jfenwick 1889 default:
1461     throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
1462     }
1463 jfenwick 2152
1464 jfenwick 1889 }
1465    
1466 jfenwick 2271 const DataTypes::ValueType*
1467     DataLazy::resolveSample(BufferGroup& bg, int sampleNo, size_t& roffset)
1468     {
1469     #ifdef _OPENMP
1470     int tid=omp_get_thread_num();
1471     #else
1472     int tid=0;
1473     #endif
1474     return resolveSample(bg.getBuffer(tid),bg.getOffset(tid),sampleNo,roffset);
1475     }
1476    
1477    
1478 jfenwick 2177 // This needs to do the work of the idenity constructor
1479     void
1480     DataLazy::resolveToIdentity()
1481     {
1482     if (m_op==IDENTITY)
1483     return;
1484     DataReady_ptr p=resolve();
1485     makeIdentity(p);
1486     }
1487 jfenwick 1889
1488 jfenwick 2177 void DataLazy::makeIdentity(const DataReady_ptr& p)
1489     {
1490     m_op=IDENTITY;
1491     m_axis_offset=0;
1492     m_transpose=0;
1493     m_SL=m_SM=m_SR=0;
1494     m_children=m_height=0;
1495     m_id=p;
1496     if(p->isConstant()) {m_readytype='C';}
1497     else if(p->isExpanded()) {m_readytype='E';}
1498     else if (p->isTagged()) {m_readytype='T';}
1499     else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1500     m_buffsRequired=1;
1501     m_samplesize=p->getNumDPPSample()*p->getNoValues();
1502     m_maxsamplesize=m_samplesize;
1503     m_left.reset();
1504     m_right.reset();
1505     }
1506    
1507 jfenwick 1899 // To simplify the memory management, all threads operate on one large vector, rather than one each.
1508     // Each sample is evaluated independently and copied into the result DataExpanded.
1509 jfenwick 1879 DataReady_ptr
1510     DataLazy::resolve()
1511     {
1512    
1513 jfenwick 2092 LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)
1514     LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)
1515 jfenwick 1899 if (m_readytype!='E') // if the whole sub-expression is Constant or Tagged, then evaluate it normally
1516 jfenwick 1889 {
1517     collapse();
1518     }
1519 jfenwick 1899 if (m_op==IDENTITY) // So a lazy expression of Constant or Tagged data will be returned here.
1520 jfenwick 1889 {
1521     return m_id;
1522     }
1523     // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
1524 jfenwick 2066 size_t threadbuffersize=m_maxsamplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough
1525 jfenwick 1899 // storage to evaluate its expression
1526 jfenwick 1885 int numthreads=1;
1527     #ifdef _OPENMP
1528 jfenwick 2271 numthreads=omp_get_max_threads();
1529 jfenwick 1885 #endif
1530 jfenwick 1886 ValueType v(numthreads*threadbuffersize);
1531 jfenwick 2092 LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)
1532 jfenwick 1898 DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(), ValueType(getNoValues()));
1533 jfenwick 2271 ValueType& resvec=result->getVectorRW();
1534 jfenwick 1879 DataReady_ptr resptr=DataReady_ptr(result);
1535     int sample;
1536 jfenwick 1898 size_t outoffset; // offset in the output data
1537 jfenwick 1885 int totalsamples=getNumSamples();
1538 jfenwick 1899 const ValueType* res=0; // Vector storing the answer
1539     size_t resoffset=0; // where in the vector to find the answer
1540 jfenwick 2152 LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1541 caltinay 2082 #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
1542 jfenwick 1885 for (sample=0;sample<totalsamples;++sample)
1543 jfenwick 1879 {
1544 jfenwick 2157 if (sample==0) {ENABLEDEBUG}
1545 jfenwick 2166
1546     // if (sample==5800/getNumDPPSample()) {ENABLEDEBUG}
1547 jfenwick 2092 LAZYDEBUG(cout << "################################# " << sample << endl;)
1548 jfenwick 1885 #ifdef _OPENMP
1549 jfenwick 1898 res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
1550 jfenwick 1885 #else
1551 jfenwick 1899 res=resolveSample(v,0,sample,resoffset); // res would normally be v, but not if its a single IDENTITY op.
1552 jfenwick 1885 #endif
1553 jfenwick 2092 LAZYDEBUG(cerr << "-------------------------------- " << endl;)
1554 jfenwick 2157 LAZYDEBUG(cerr<< "Copying sample#" << sample << endl;)
1555 jfenwick 1898 outoffset=result->getPointOffset(sample,0);
1556 jfenwick 2157 LAZYDEBUG(cerr << "offset=" << outoffset << " from offset=" << resoffset << " " << m_samplesize << " doubles" << endl;)
1557 jfenwick 1898 for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset) // copy values into the output vector
1558 jfenwick 1879 {
1559 jfenwick 2166 LAZYDEBUG(cerr << "outoffset=" << outoffset << " resoffset=" << resoffset << " " << (*res)[resoffset]<< endl;)
1560 jfenwick 1898 resvec[outoffset]=(*res)[resoffset];
1561 jfenwick 1879 }
1562 jfenwick 2157 LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)
1563 jfenwick 2092 LAZYDEBUG(cerr << "*********************************" << endl;)
1564 jfenwick 2157 DISABLEDEBUG
1565 jfenwick 1879 }
1566     return resptr;
1567     }
1568    
1569 jfenwick 1865 std::string
1570     DataLazy::toString() const
1571     {
1572 jfenwick 1886 ostringstream oss;
1573     oss << "Lazy Data:";
1574     intoString(oss);
1575     return oss.str();
1576 jfenwick 1865 }
1577    
1578 jfenwick 1899
1579 jfenwick 1886 void
1580     DataLazy::intoString(ostringstream& oss) const
1581     {
1582 jfenwick 2271 // oss << "[" << m_children <<";"<<m_height <<"]";
1583 jfenwick 1886 switch (getOpgroup(m_op))
1584     {
1585 jfenwick 1889 case G_IDENTITY:
1586     if (m_id->isExpanded())
1587     {
1588     oss << "E";
1589     }
1590     else if (m_id->isTagged())
1591     {
1592     oss << "T";
1593     }
1594     else if (m_id->isConstant())
1595     {
1596     oss << "C";
1597     }
1598     else
1599     {
1600     oss << "?";
1601     }
1602 jfenwick 1886 oss << '@' << m_id.get();
1603     break;
1604     case G_BINARY:
1605     oss << '(';
1606     m_left->intoString(oss);
1607     oss << ' ' << opToString(m_op) << ' ';
1608     m_right->intoString(oss);
1609     oss << ')';
1610     break;
1611     case G_UNARY:
1612 jfenwick 2147 case G_UNARY_P:
1613 jfenwick 2037 case G_NP1OUT:
1614 jfenwick 2084 case G_NP1OUT_P:
1615 jfenwick 1886 oss << opToString(m_op) << '(';
1616     m_left->intoString(oss);
1617     oss << ')';
1618     break;
1619 jfenwick 2066 case G_TENSORPROD:
1620     oss << opToString(m_op) << '(';
1621     m_left->intoString(oss);
1622     oss << ", ";
1623     m_right->intoString(oss);
1624     oss << ')';
1625     break;
1626 jfenwick 1886 default:
1627     oss << "UNKNOWN";
1628     }
1629     }
1630    
1631 jfenwick 1865 DataAbstract*
1632     DataLazy::deepCopy()
1633     {
1634 jfenwick 1935 switch (getOpgroup(m_op))
1635 jfenwick 1865 {
1636 jfenwick 1935 case G_IDENTITY: return new DataLazy(m_id->deepCopy()->getPtr());
1637     case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1638     case G_BINARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1639 jfenwick 2066 case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1640     case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1641 jfenwick 1935 default:
1642     throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1643 jfenwick 1865 }
1644     }
1645    
1646    
1647 jfenwick 2066 // There is no single, natural interpretation of getLength on DataLazy.
1648     // Instances of DataReady can look at the size of their vectors.
1649     // For lazy though, it could be the size the data would be if it were resolved;
1650     // or it could be some function of the lengths of the DataReady instances which
1651     // form part of the expression.
1652     // Rather than have people making assumptions, I have disabled the method.
1653 jfenwick 1865 DataTypes::ValueType::size_type
1654     DataLazy::getLength() const
1655     {
1656 jfenwick 2066 throw DataException("getLength() does not make sense for lazy data.");
1657 jfenwick 1865 }
1658    
1659    
1660     DataAbstract*
1661     DataLazy::getSlice(const DataTypes::RegionType& region) const
1662     {
1663 jfenwick 1879 throw DataException("getSlice - not implemented for Lazy objects.");
1664 jfenwick 1865 }
1665    
1666 jfenwick 1935
1667     // To do this we need to rely on our child nodes
1668 jfenwick 1865 DataTypes::ValueType::size_type
1669     DataLazy::getPointOffset(int sampleNo,
1670 jfenwick 1935 int dataPointNo)
1671     {
1672     if (m_op==IDENTITY)
1673     {
1674     return m_id->getPointOffset(sampleNo,dataPointNo);
1675     }
1676     if (m_readytype!='E')
1677     {
1678     collapse();
1679     return m_id->getPointOffset(sampleNo,dataPointNo);
1680     }
1681     // at this point we do not have an identity node and the expression will be Expanded
1682     // so we only need to know which child to ask
1683     if (m_left->m_readytype=='E')
1684     {
1685     return m_left->getPointOffset(sampleNo,dataPointNo);
1686     }
1687     else
1688     {
1689     return m_right->getPointOffset(sampleNo,dataPointNo);
1690     }
1691     }
1692    
1693     // To do this we need to rely on our child nodes
1694     DataTypes::ValueType::size_type
1695     DataLazy::getPointOffset(int sampleNo,
1696 jfenwick 1865 int dataPointNo) const
1697     {
1698 jfenwick 1935 if (m_op==IDENTITY)
1699     {
1700     return m_id->getPointOffset(sampleNo,dataPointNo);
1701     }
1702     if (m_readytype=='E')
1703     {
1704     // at this point we do not have an identity node and the expression will be Expanded
1705     // so we only need to know which child to ask
1706     if (m_left->m_readytype=='E')
1707     {
1708     return m_left->getPointOffset(sampleNo,dataPointNo);
1709     }
1710     else
1711     {
1712     return m_right->getPointOffset(sampleNo,dataPointNo);
1713     }
1714     }
1715     if (m_readytype=='C')
1716     {
1717     return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
1718     }
1719     throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1720 jfenwick 1865 }
1721    
1722 jfenwick 2271
1723     // I have decided to let Data:: handle this issue.
1724 jfenwick 1901 void
1725     DataLazy::setToZero()
1726     {
1727 jfenwick 2271 // DataTypes::ValueType v(getNoValues(),0);
1728     // m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1729     // m_op=IDENTITY;
1730     // m_right.reset();
1731     // m_left.reset();
1732     // m_readytype='C';
1733     // m_buffsRequired=1;
1734    
1735     throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
1736    
1737 jfenwick 1901 }
1738    
1739 jfenwick 2271 bool
1740     DataLazy::actsExpanded() const
1741     {
1742     return (m_readytype=='E');
1743     }
1744    
1745 jfenwick 1865 } // end namespace

  ViewVC Help
Powered by ViewVC 1.1.26