/[escript]/trunk/escript/src/DataLazy.cpp
ViewVC logotype

Annotation of /trunk/escript/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26