/[escript]/branches/clazy/escriptcore/src/DataLazy.cpp
ViewVC logotype

Contents of /branches/clazy/escriptcore/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6066 - (show annotations)
Tue Mar 15 06:42:55 2016 UTC (3 years ago) by jfenwick
Original Path: trunk/escriptcore/src/DataLazy.cpp
File size: 66516 byte(s)
Code path conversions.

TestDomain can now return a JMPI.
Some tests now use TestDomain (didn't want to undo all that work
from before I found the ASSERT fix).
Moving more code to new path.
Commenting out a bunch of old ops.


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

  ViewVC Help
Powered by ViewVC 1.1.26