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

  ViewVC Help
Powered by ViewVC 1.1.26