/[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 5972 - (hide annotations)
Wed Feb 24 04:05:30 2016 UTC (3 years, 1 month ago) by caltinay
Original Path: trunk/escriptcore/src/DataLazy.cpp
File size: 67893 byte(s)
Major rework of our exceptions. We now have specific
AssertException
NotImplementedError
ValueError
which translate to the corresponding python exception type.
I have gone through a few places and replaced things but not everywhere.


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

  ViewVC Help
Powered by ViewVC 1.1.26