/[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 2086 - (hide annotations)
Mon Nov 24 02:38:50 2008 UTC (10 years, 4 months ago) by jfenwick
Original Path: trunk/escript/src/DataLazy.cpp
File size: 43128 byte(s)
Added checks in C_GeneralTensorProduct (Data:: and Delayed forms) as 
well as the DataAbstract Constructor to prevent Objects with Rank>4 
being created.

Moved the relevant #define into systemdep.

Removed some comments.

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

  ViewVC Help
Powered by ViewVC 1.1.26