/[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 6042 - (show annotations)
Wed Mar 9 04:30:36 2016 UTC (3 years ago) by jfenwick
Original Path: trunk/escriptcore/src/DataLazy.cpp
File size: 68126 byte(s)
Preliminary to lazy cleanup


Lazy is not templated for complex as yet but hopefully
we can relieve the dependency on the briarpatch of
the current methods.


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
995
996
997
998
999 #define PROC_OP(TYPE,X) \
1000 for (int j=0;j<onumsteps;++j)\
1001 {\
1002 for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1003 { \
1004 LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1005 LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1006 tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1007 LAZYDEBUG(cout << " result= " << resultp[0] << endl;) \
1008 lroffset+=leftstep; \
1009 rroffset+=rightstep; \
1010 }\
1011 lroffset+=oleftstep;\
1012 rroffset+=orightstep;\
1013 }
1014
1015
1016 // The result will be stored in m_samples
1017 // The return value is a pointer to the DataVector, offset is the offset within the return value
1018 const DataTypes::RealVectorType*
1019 DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset) const
1020 {
1021 LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1022 // collapse so we have a 'E' node or an IDENTITY for some other type
1023 if (m_readytype!='E' && m_op!=IDENTITY)
1024 {
1025 collapse();
1026 }
1027 if (m_op==IDENTITY)
1028 {
1029 const RealVectorType& vec=m_id->getVectorRO();
1030 roffset=m_id->getPointOffset(sampleNo, 0);
1031 #ifdef LAZY_STACK_PROF
1032 int x;
1033 if (&x<stackend[omp_get_thread_num()])
1034 {
1035 stackend[omp_get_thread_num()]=&x;
1036 }
1037 #endif
1038 return &(vec);
1039 }
1040 if (m_readytype!='E')
1041 {
1042 throw DataException("Programmer Error - Collapse did not produce an expanded node.");
1043 }
1044 if (m_sampleids[tid]==sampleNo)
1045 {
1046 roffset=tid*m_samplesize;
1047 return &(m_samples); // sample is already resolved
1048 }
1049 m_sampleids[tid]=sampleNo;
1050
1051 switch (getOpgroup(m_op))
1052 {
1053 case G_UNARY:
1054 case G_UNARY_P: return resolveNodeUnary(tid, sampleNo, roffset);
1055 case G_BINARY: return resolveNodeBinary(tid, sampleNo, roffset);
1056 case G_NP1OUT: return resolveNodeNP1OUT(tid, sampleNo, roffset);
1057 case G_NP1OUT_P: return resolveNodeNP1OUT_P(tid, sampleNo, roffset);
1058 case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);
1059 case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);
1060 case G_REDUCTION: return resolveNodeReduction(tid, sampleNo, roffset);
1061 case G_CONDEVAL: return resolveNodeCondEval(tid, sampleNo, roffset);
1062 default:
1063 throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
1064 }
1065 }
1066
1067 const DataTypes::RealVectorType*
1068 DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset) const
1069 {
1070 // we assume that any collapsing has been done before we get here
1071 // since we only have one argument we don't need to think about only
1072 // processing single points.
1073 // we will also know we won't get identity nodes
1074 if (m_readytype!='E')
1075 {
1076 throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1077 }
1078 if (m_op==IDENTITY)
1079 {
1080 throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1081 }
1082 const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
1083 const double* left=&((*leftres)[roffset]);
1084 roffset=m_samplesize*tid;
1085 double* result=&(m_samples[roffset]);
1086 switch (m_op)
1087 {
1088 case SIN:
1089 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
1090 break;
1091 case COS:
1092 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
1093 break;
1094 case TAN:
1095 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
1096 break;
1097 case ASIN:
1098 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
1099 break;
1100 case ACOS:
1101 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
1102 break;
1103 case ATAN:
1104 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
1105 break;
1106 case SINH:
1107 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
1108 break;
1109 case COSH:
1110 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
1111 break;
1112 case TANH:
1113 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
1114 break;
1115 case ERF:
1116 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1117 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1118 #else
1119 tensor_unary_operation(m_samplesize, left, result, ::erf);
1120 break;
1121 #endif
1122 case ASINH:
1123 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1124 tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
1125 #else
1126 tensor_unary_operation(m_samplesize, left, result, ::asinh);
1127 #endif
1128 break;
1129 case ACOSH:
1130 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1131 tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
1132 #else
1133 tensor_unary_operation(m_samplesize, left, result, ::acosh);
1134 #endif
1135 break;
1136 case ATANH:
1137 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1138 tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
1139 #else
1140 tensor_unary_operation(m_samplesize, left, result, ::atanh);
1141 #endif
1142 break;
1143 case LOG10:
1144 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
1145 break;
1146 case LOG:
1147 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
1148 break;
1149 case SIGN:
1150 tensor_unary_operation(m_samplesize, left, result, escript::fsign);
1151 break;
1152 case ABS:
1153 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
1154 break;
1155 case NEG:
1156 tensor_unary_operation(m_samplesize, left, result, negate<double>());
1157 break;
1158 case POS:
1159 // it doesn't mean anything for delayed.
1160 // it will just trigger a deep copy of the lazy object
1161 throw DataException("Programmer error - POS not supported for lazy data.");
1162 break;
1163 case EXP:
1164 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
1165 break;
1166 case SQRT:
1167 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
1168 break;
1169 case RECIP:
1170 tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
1171 break;
1172 case GZ:
1173 tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
1174 break;
1175 case LZ:
1176 tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
1177 break;
1178 case GEZ:
1179 tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
1180 break;
1181 case LEZ:
1182 tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
1183 break;
1184 // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
1185 case NEZ:
1186 tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
1187 break;
1188 case EZ:
1189 tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
1190 break;
1191
1192 default:
1193 throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1194 }
1195 return &(m_samples);
1196 }
1197
1198
1199 const DataTypes::RealVectorType*
1200 DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset) const
1201 {
1202 // we assume that any collapsing has been done before we get here
1203 // since we only have one argument we don't need to think about only
1204 // processing single points.
1205 // we will also know we won't get identity nodes
1206 if (m_readytype!='E')
1207 {
1208 throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1209 }
1210 if (m_op==IDENTITY)
1211 {
1212 throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1213 }
1214 size_t loffset=0;
1215 const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);
1216
1217 roffset=m_samplesize*tid;
1218 unsigned int ndpps=getNumDPPSample();
1219 unsigned int psize=DataTypes::noValues(m_left->getShape());
1220 double* result=&(m_samples[roffset]);
1221 switch (m_op)
1222 {
1223 case MINVAL:
1224 {
1225 for (unsigned int z=0;z<ndpps;++z)
1226 {
1227 FMin op;
1228 *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());
1229 loffset+=psize;
1230 result++;
1231 }
1232 }
1233 break;
1234 case MAXVAL:
1235 {
1236 for (unsigned int z=0;z<ndpps;++z)
1237 {
1238 FMax op;
1239 *result=DataMaths::reductionOp(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);
1240 loffset+=psize;
1241 result++;
1242 }
1243 }
1244 break;
1245 default:
1246 throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1247 }
1248 return &(m_samples);
1249 }
1250
1251 const DataTypes::RealVectorType*
1252 DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset) const
1253 {
1254 // we assume that any collapsing has been done before we get here
1255 // since we only have one argument we don't need to think about only
1256 // processing single points.
1257 if (m_readytype!='E')
1258 {
1259 throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");
1260 }
1261 if (m_op==IDENTITY)
1262 {
1263 throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1264 }
1265 size_t subroffset;
1266 const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1267 roffset=m_samplesize*tid;
1268 size_t loop=0;
1269 size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1270 size_t step=getNoValues();
1271 size_t offset=roffset;
1272 switch (m_op)
1273 {
1274 case SYM:
1275 for (loop=0;loop<numsteps;++loop)
1276 {
1277 DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1278 subroffset+=step;
1279 offset+=step;
1280 }
1281 break;
1282 case NSYM:
1283 for (loop=0;loop<numsteps;++loop)
1284 {
1285 DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1286 subroffset+=step;
1287 offset+=step;
1288 }
1289 break;
1290 default:
1291 throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1292 }
1293 return &m_samples;
1294 }
1295
1296 const DataTypes::RealVectorType*
1297 DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset) const
1298 {
1299 // we assume that any collapsing has been done before we get here
1300 // since we only have one argument we don't need to think about only
1301 // processing single points.
1302 if (m_readytype!='E')
1303 {
1304 throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");
1305 }
1306 if (m_op==IDENTITY)
1307 {
1308 throw DataException("Programmer error - resolveNodeNP1OUT_P should not be called on identity nodes.");
1309 }
1310 size_t subroffset;
1311 size_t offset;
1312 const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1313 roffset=m_samplesize*tid;
1314 offset=roffset;
1315 size_t loop=0;
1316 size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1317 size_t outstep=getNoValues();
1318 size_t instep=m_left->getNoValues();
1319 switch (m_op)
1320 {
1321 case TRACE:
1322 for (loop=0;loop<numsteps;++loop)
1323 {
1324 DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);
1325 subroffset+=instep;
1326 offset+=outstep;
1327 }
1328 break;
1329 case TRANS:
1330 for (loop=0;loop<numsteps;++loop)
1331 {
1332 DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);
1333 subroffset+=instep;
1334 offset+=outstep;
1335 }
1336 break;
1337 default:
1338 throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1339 }
1340 return &m_samples;
1341 }
1342
1343
1344 const DataTypes::RealVectorType*
1345 DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset) const
1346 {
1347 if (m_readytype!='E')
1348 {
1349 throw DataException("Programmer error - resolveNodeNP1OUT_2P should only be called on expanded Data.");
1350 }
1351 if (m_op==IDENTITY)
1352 {
1353 throw DataException("Programmer error - resolveNodeNP1OUT_2P should not be called on identity nodes.");
1354 }
1355 size_t subroffset;
1356 size_t offset;
1357 const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1358 roffset=m_samplesize*tid;
1359 offset=roffset;
1360 size_t loop=0;
1361 size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1362 size_t outstep=getNoValues();
1363 size_t instep=m_left->getNoValues();
1364 switch (m_op)
1365 {
1366 case SWAP:
1367 for (loop=0;loop<numsteps;++loop)
1368 {
1369 DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);
1370 subroffset+=instep;
1371 offset+=outstep;
1372 }
1373 break;
1374 default:
1375 throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1376 }
1377 return &m_samples;
1378 }
1379
1380 const DataTypes::RealVectorType*
1381 DataLazy::resolveNodeCondEval(int tid, int sampleNo, size_t& roffset) const
1382 {
1383 if (m_readytype!='E')
1384 {
1385 throw DataException("Programmer error - resolveNodeCondEval should only be called on expanded Data.");
1386 }
1387 if (m_op!=CONDEVAL)
1388 {
1389 throw DataException("Programmer error - resolveNodeCondEval should only be called on CONDEVAL nodes.");
1390 }
1391 size_t subroffset;
1392
1393 const RealVectorType* maskres=m_mask->resolveNodeSample(tid, sampleNo, subroffset);
1394 const RealVectorType* srcres=0;
1395 if ((*maskres)[subroffset]>0)
1396 {
1397 srcres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1398 }
1399 else
1400 {
1401 srcres=m_right->resolveNodeSample(tid, sampleNo, subroffset);
1402 }
1403
1404 // Now we need to copy the result
1405
1406 roffset=m_samplesize*tid;
1407 for (int i=0;i<m_samplesize;++i)
1408 {
1409 m_samples[roffset+i]=(*srcres)[subroffset+i];
1410 }
1411
1412 return &m_samples;
1413 }
1414
1415 // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1416 // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1417 // If both children are expanded, then we can process them in a single operation (we treat
1418 // the whole sample as one big datapoint.
1419 // If one of the children is not expanded, then we need to treat each point in the sample
1420 // individually.
1421 // There is an additional complication when scalar operations are considered.
1422 // For example, 2+Vector.
1423 // In this case each double within the point is treated individually
1424 const DataTypes::RealVectorType*
1425 DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset) const
1426 {
1427 LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1428
1429 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
1430 // first work out which of the children are expanded
1431 bool leftExp=(m_left->m_readytype=='E');
1432 bool rightExp=(m_right->m_readytype=='E');
1433 if (!leftExp && !rightExp)
1434 {
1435 throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1436 }
1437 bool leftScalar=(m_left->getRank()==0);
1438 bool rightScalar=(m_right->getRank()==0);
1439 if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1440 {
1441 throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1442 }
1443 size_t leftsize=m_left->getNoValues();
1444 size_t rightsize=m_right->getNoValues();
1445 size_t chunksize=1; // how many doubles will be processed in one go
1446 int leftstep=0; // how far should the left offset advance after each step
1447 int rightstep=0;
1448 int numsteps=0; // total number of steps for the inner loop
1449 int oleftstep=0; // the o variables refer to the outer loop
1450 int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1451 int onumsteps=1;
1452
1453 bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1454 bool RES=(rightExp && rightScalar);
1455 bool LS=(!leftExp && leftScalar); // left is a single scalar
1456 bool RS=(!rightExp && rightScalar);
1457 bool LN=(!leftExp && !leftScalar); // left is a single non-scalar
1458 bool RN=(!rightExp && !rightScalar);
1459 bool LEN=(leftExp && !leftScalar); // left is an expanded non-scalar
1460 bool REN=(rightExp && !rightScalar);
1461
1462 if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1463 {
1464 chunksize=m_left->getNumDPPSample()*leftsize;
1465 leftstep=0;
1466 rightstep=0;
1467 numsteps=1;
1468 }
1469 else if (LES || RES)
1470 {
1471 chunksize=1;
1472 if (LES) // left is an expanded scalar
1473 {
1474 if (RS)
1475 {
1476 leftstep=1;
1477 rightstep=0;
1478 numsteps=m_left->getNumDPPSample();
1479 }
1480 else // RN or REN
1481 {
1482 leftstep=0;
1483 oleftstep=1;
1484 rightstep=1;
1485 orightstep=(RN ? -(int)rightsize : 0);
1486 numsteps=rightsize;
1487 onumsteps=m_left->getNumDPPSample();
1488 }
1489 }
1490 else // right is an expanded scalar
1491 {
1492 if (LS)
1493 {
1494 rightstep=1;
1495 leftstep=0;
1496 numsteps=m_right->getNumDPPSample();
1497 }
1498 else
1499 {
1500 rightstep=0;
1501 orightstep=1;
1502 leftstep=1;
1503 oleftstep=(LN ? -(int)leftsize : 0);
1504 numsteps=leftsize;
1505 onumsteps=m_right->getNumDPPSample();
1506 }
1507 }
1508 }
1509 else // this leaves (LEN, RS), (LEN, RN) and their transposes
1510 {
1511 if (LEN) // and Right will be a single value
1512 {
1513 chunksize=rightsize;
1514 leftstep=rightsize;
1515 rightstep=0;
1516 numsteps=m_left->getNumDPPSample();
1517 if (RS)
1518 {
1519 numsteps*=leftsize;
1520 }
1521 }
1522 else // REN
1523 {
1524 chunksize=leftsize;
1525 rightstep=leftsize;
1526 leftstep=0;
1527 numsteps=m_right->getNumDPPSample();
1528 if (LS)
1529 {
1530 numsteps*=rightsize;
1531 }
1532 }
1533 }
1534
1535 int resultStep=max(leftstep,rightstep); // only one (at most) should be !=0
1536 // Get the values of sub-expressions
1537 const RealVectorType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);
1538 const RealVectorType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
1539 LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1540 LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1541 LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1542 LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1543 LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1544 LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1545 LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN << endl;)
1546
1547 LAZYDEBUG(cout << "Left res["<< lroffset<< "]=" << (*left)[lroffset] << endl;)
1548 LAZYDEBUG(cout << "Right res["<< rroffset<< "]=" << (*right)[rroffset] << endl;)
1549
1550
1551 roffset=m_samplesize*tid;
1552 double* resultp=&(m_samples[roffset]); // results are stored at the vector offset we received
1553 switch(m_op)
1554 {
1555 case ADD:
1556 //PROC_OP(NO_ARG,plus<double>());
1557 DataMaths::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1558 &(*left)[0],
1559 &(*right)[0],
1560 chunksize,
1561 onumsteps,
1562 numsteps,
1563 resultStep,
1564 leftstep,
1565 rightstep,
1566 oleftstep,
1567 orightstep,
1568 escript::ESFunction::PLUSF);
1569 break;
1570 case SUB:
1571 PROC_OP(NO_ARG,minus<double>());
1572 break;
1573 case MUL:
1574 PROC_OP(NO_ARG,multiplies<double>());
1575 break;
1576 case DIV:
1577 PROC_OP(NO_ARG,divides<double>());
1578 break;
1579 case POW:
1580 PROC_OP(double (double,double),::pow);
1581 break;
1582 default:
1583 throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1584 }
1585 LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
1586 return &m_samples;
1587 }
1588
1589
1590 // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1591 // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1592 // unlike the other resolve helpers, we must treat these datapoints separately.
1593 const DataTypes::RealVectorType*
1594 DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset) const
1595 {
1596 LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1597
1598 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
1599 // first work out which of the children are expanded
1600 bool leftExp=(m_left->m_readytype=='E');
1601 bool rightExp=(m_right->m_readytype=='E');
1602 int steps=getNumDPPSample();
1603 int leftStep=(leftExp? m_left->getNoValues() : 0); // do not have scalars as input to this method
1604 int rightStep=(rightExp?m_right->getNoValues() : 0);
1605
1606 int resultStep=getNoValues();
1607 roffset=m_samplesize*tid;
1608 size_t offset=roffset;
1609
1610 const RealVectorType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
1611
1612 const RealVectorType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
1613
1614 LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) << endl;
1615 cout << getNoValues() << endl;)
1616
1617
1618 LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1619 LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1620 LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1621 LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1622 LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1623 LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1624 LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1625
1626 double* resultp=&(m_samples[offset]); // results are stored at the vector offset we received
1627 switch(m_op)
1628 {
1629 case PROD:
1630 for (int i=0;i<steps;++i,resultp+=resultStep)
1631 {
1632 const double *ptr_0 = &((*left)[lroffset]);
1633 const double *ptr_1 = &((*right)[rroffset]);
1634
1635 LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1636 LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1637
1638 matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
1639
1640 lroffset+=leftStep;
1641 rroffset+=rightStep;
1642 }
1643 break;
1644 default:
1645 throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1646 }
1647 roffset=offset;
1648 return &m_samples;
1649 }
1650
1651
1652 const DataTypes::RealVectorType*
1653 DataLazy::resolveSample(int sampleNo, size_t& roffset) const
1654 {
1655 #ifdef _OPENMP
1656 int tid=omp_get_thread_num();
1657 #else
1658 int tid=0;
1659 #endif
1660
1661 #ifdef LAZY_STACK_PROF
1662 stackstart[tid]=&tid;
1663 stackend[tid]=&tid;
1664 const DataTypes::RealVectorType* r=resolveNodeSample(tid, sampleNo, roffset);
1665 size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];
1666 #pragma omp critical
1667 if (d>maxstackuse)
1668 {
1669 cout << "Max resolve Stack use " << d << endl;
1670 maxstackuse=d;
1671 }
1672 return r;
1673 #else
1674 return resolveNodeSample(tid, sampleNo, roffset);
1675 #endif
1676 }
1677
1678
1679 // This needs to do the work of the identity constructor
1680 void
1681 DataLazy::resolveToIdentity()
1682 {
1683 if (m_op==IDENTITY)
1684 return;
1685 DataReady_ptr p=resolveNodeWorker();
1686 makeIdentity(p);
1687 }
1688
1689 void DataLazy::makeIdentity(const DataReady_ptr& p)
1690 {
1691 m_op=IDENTITY;
1692 m_axis_offset=0;
1693 m_transpose=0;
1694 m_SL=m_SM=m_SR=0;
1695 m_children=m_height=0;
1696 m_id=p;
1697 if(p->isConstant()) {m_readytype='C';}
1698 else if(p->isExpanded()) {m_readytype='E';}
1699 else if (p->isTagged()) {m_readytype='T';}
1700 else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1701 m_samplesize=p->getNumDPPSample()*p->getNoValues();
1702 m_left.reset();
1703 m_right.reset();
1704 }
1705
1706
1707 DataReady_ptr
1708 DataLazy::resolve()
1709 {
1710 resolveToIdentity();
1711 return m_id;
1712 }
1713
1714
1715 /* This is really a static method but I think that caused problems in windows */
1716 void
1717 DataLazy::resolveGroupWorker(std::vector<DataLazy*>& dats)
1718 {
1719 if (dats.empty())
1720 {
1721 return;
1722 }
1723 vector<DataLazy*> work;
1724 FunctionSpace fs=dats[0]->getFunctionSpace();
1725 bool match=true;
1726 for (int i=dats.size()-1;i>=0;--i)
1727 {
1728 if (dats[i]->m_readytype!='E')
1729 {
1730 dats[i]->collapse();
1731 }
1732 if (dats[i]->m_op!=IDENTITY)
1733 {
1734 work.push_back(dats[i]);
1735 if (fs!=dats[i]->getFunctionSpace())
1736 {
1737 match=false;
1738 }
1739 }
1740 }
1741 if (work.empty())
1742 {
1743 return; // no work to do
1744 }
1745 if (match) // all functionspaces match. Yes I realise this is overly strict
1746 { // it is possible that dats[0] is one of the objects which we discarded and
1747 // all the other functionspaces match.
1748 vector<DataExpanded*> dep;
1749 vector<RealVectorType*> vecs;
1750 for (int i=0;i<work.size();++i)
1751 {
1752 dep.push_back(new DataExpanded(fs,work[i]->getShape(), RealVectorType(work[i]->getNoValues())));
1753 vecs.push_back(&(dep[i]->getVectorRW()));
1754 }
1755 int totalsamples=work[0]->getNumSamples();
1756 const RealVectorType* res=0; // Storage for answer
1757 int sample;
1758 #pragma omp parallel private(sample, res)
1759 {
1760 size_t roffset=0;
1761 #pragma omp for schedule(static)
1762 for (sample=0;sample<totalsamples;++sample)
1763 {
1764 roffset=0;
1765 int j;
1766 for (j=work.size()-1;j>=0;--j)
1767 {
1768 #ifdef _OPENMP
1769 res=work[j]->resolveNodeSample(omp_get_thread_num(),sample,roffset);
1770 #else
1771 res=work[j]->resolveNodeSample(0,sample,roffset);
1772 #endif
1773 RealVectorType::size_type outoffset=dep[j]->getPointOffset(sample,0);
1774 memcpy(&((*vecs[j])[outoffset]),&((*res)[roffset]),work[j]->m_samplesize*sizeof(RealVectorType::ElementType));
1775 }
1776 }
1777 }
1778 // Now we need to load the new results as identity ops into the lazy nodes
1779 for (int i=work.size()-1;i>=0;--i)
1780 {
1781 work[i]->makeIdentity(REFCOUNTNS::dynamic_pointer_cast<DataReady>(dep[i]->getPtr()));
1782 }
1783 }
1784 else // functionspaces do not match
1785 {
1786 for (int i=0;i<work.size();++i)
1787 {
1788 work[i]->resolveToIdentity();
1789 }
1790 }
1791 }
1792
1793
1794
1795 // This version of resolve uses storage in each node to hold results
1796 DataReady_ptr
1797 DataLazy::resolveNodeWorker()
1798 {
1799 if (m_readytype!='E') // if the whole sub-expression is Constant or Tagged, then evaluate it normally
1800 {
1801 collapse();
1802 }
1803 if (m_op==IDENTITY) // So a lazy expression of Constant or Tagged data will be returned here.
1804 {
1805 return m_id;
1806 }
1807 // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
1808 DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(), RealVectorType(getNoValues()));
1809 RealVectorType& resvec=result->getVectorRW();
1810 DataReady_ptr resptr=DataReady_ptr(result);
1811
1812 int sample;
1813 int totalsamples=getNumSamples();
1814 const RealVectorType* res=0; // Storage for answer
1815 LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1816 #pragma omp parallel private(sample,res)
1817 {
1818 size_t roffset=0;
1819 #ifdef LAZY_STACK_PROF
1820 stackstart[omp_get_thread_num()]=&roffset;
1821 stackend[omp_get_thread_num()]=&roffset;
1822 #endif
1823 #pragma omp for schedule(static)
1824 for (sample=0;sample<totalsamples;++sample)
1825 {
1826 roffset=0;
1827 #ifdef _OPENMP
1828 res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
1829 #else
1830 res=resolveNodeSample(0,sample,roffset);
1831 #endif
1832 LAZYDEBUG(cout << "Sample #" << sample << endl;)
1833 LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
1834 RealVectorType::size_type outoffset=result->getPointOffset(sample,0);
1835 memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(RealVectorType::ElementType));
1836 }
1837 }
1838 #ifdef LAZY_STACK_PROF
1839 for (int i=0;i<getNumberOfThreads();++i)
1840 {
1841 size_t r=((size_t)stackstart[i] - (size_t)stackend[i]);
1842 // cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " << r << endl;
1843 if (r>maxstackuse)
1844 {
1845 maxstackuse=r;
1846 }
1847 }
1848 cout << "Max resolve Stack use=" << maxstackuse << endl;
1849 #endif
1850 return resptr;
1851 }
1852
1853 std::string
1854 DataLazy::toString() const
1855 {
1856 ostringstream oss;
1857 oss << "Lazy Data: [depth=" << m_height<< "] ";
1858 switch (escriptParams.getLAZY_STR_FMT())
1859 {
1860 case 1: // tree format
1861 oss << endl;
1862 intoTreeString(oss,"");
1863 break;
1864 case 2: // just the depth
1865 break;
1866 default:
1867 intoString(oss);
1868 break;
1869 }
1870 return oss.str();
1871 }
1872
1873
1874 void
1875 DataLazy::intoString(ostringstream& oss) const
1876 {
1877 // oss << "[" << m_children <<";"<<m_height <<"]";
1878 switch (getOpgroup(m_op))
1879 {
1880 case G_IDENTITY:
1881 if (m_id->isExpanded())
1882 {
1883 oss << "E";
1884 }
1885 else if (m_id->isTagged())
1886 {
1887 oss << "T";
1888 }
1889 else if (m_id->isConstant())
1890 {
1891 oss << "C";
1892 }
1893 else
1894 {
1895 oss << "?";
1896 }
1897 oss << '@' << m_id.get();
1898 break;
1899 case G_BINARY:
1900 oss << '(';
1901 m_left->intoString(oss);
1902 oss << ' ' << opToString(m_op) << ' ';
1903 m_right->intoString(oss);
1904 oss << ')';
1905 break;
1906 case G_UNARY:
1907 case G_UNARY_P:
1908 case G_NP1OUT:
1909 case G_NP1OUT_P:
1910 case G_REDUCTION:
1911 oss << opToString(m_op) << '(';
1912 m_left->intoString(oss);
1913 oss << ')';
1914 break;
1915 case G_TENSORPROD:
1916 oss << opToString(m_op) << '(';
1917 m_left->intoString(oss);
1918 oss << ", ";
1919 m_right->intoString(oss);
1920 oss << ')';
1921 break;
1922 case G_NP1OUT_2P:
1923 oss << opToString(m_op) << '(';
1924 m_left->intoString(oss);
1925 oss << ", " << m_axis_offset << ", " << m_transpose;
1926 oss << ')';
1927 break;
1928 case G_CONDEVAL:
1929 oss << opToString(m_op)<< '(' ;
1930 m_mask->intoString(oss);
1931 oss << " ? ";
1932 m_left->intoString(oss);
1933 oss << " : ";
1934 m_right->intoString(oss);
1935 oss << ')';
1936 break;
1937 default:
1938 oss << "UNKNOWN";
1939 }
1940 }
1941
1942
1943 void
1944 DataLazy::intoTreeString(ostringstream& oss, string indent) const
1945 {
1946 oss << '[' << m_rank << ':' << setw(3) << m_samplesize << "] " << indent;
1947 switch (getOpgroup(m_op))
1948 {
1949 case G_IDENTITY:
1950 if (m_id->isExpanded())
1951 {
1952 oss << "E";
1953 }
1954 else if (m_id->isTagged())
1955 {
1956 oss << "T";
1957 }
1958 else if (m_id->isConstant())
1959 {
1960 oss << "C";
1961 }
1962 else
1963 {
1964 oss << "?";
1965 }
1966 oss << '@' << m_id.get() << endl;
1967 break;
1968 case G_BINARY:
1969 oss << opToString(m_op) << endl;
1970 indent+='.';
1971 m_left->intoTreeString(oss, indent);
1972 m_right->intoTreeString(oss, indent);
1973 break;
1974 case G_UNARY:
1975 case G_UNARY_P:
1976 case G_NP1OUT:
1977 case G_NP1OUT_P:
1978 case G_REDUCTION:
1979 oss << opToString(m_op) << endl;
1980 indent+='.';
1981 m_left->intoTreeString(oss, indent);
1982 break;
1983 case G_TENSORPROD:
1984 oss << opToString(m_op) << endl;
1985 indent+='.';
1986 m_left->intoTreeString(oss, indent);
1987 m_right->intoTreeString(oss, indent);
1988 break;
1989 case G_NP1OUT_2P:
1990 oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;
1991 indent+='.';
1992 m_left->intoTreeString(oss, indent);
1993 break;
1994 default:
1995 oss << "UNKNOWN";
1996 }
1997 }
1998
1999
2000 DataAbstract*
2001 DataLazy::deepCopy() const
2002 {
2003 switch (getOpgroup(m_op))
2004 {
2005 case G_IDENTITY: return new DataLazy(m_id->deepCopy()->getPtr());
2006 case G_UNARY:
2007 case G_REDUCTION: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
2008 case G_UNARY_P: return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_tol);
2009 case G_BINARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
2010 case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
2011 case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
2012 case G_NP1OUT_P: return new DataLazy(m_left->deepCopy()->getPtr(),m_op, m_axis_offset);
2013 case G_NP1OUT_2P: return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
2014 default:
2015 throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
2016 }
2017 }
2018
2019
2020
2021 // There is no single, natural interpretation of getLength on DataLazy.
2022 // Instances of DataReady can look at the size of their vectors.
2023 // For lazy though, it could be the size the data would be if it were resolved;
2024 // or it could be some function of the lengths of the DataReady instances which
2025 // form part of the expression.
2026 // Rather than have people making assumptions, I have disabled the method.
2027 DataTypes::RealVectorType::size_type
2028 DataLazy::getLength() const
2029 {
2030 throw DataException("getLength() does not make sense for lazy data.");
2031 }
2032
2033
2034 DataAbstract*
2035 DataLazy::getSlice(const DataTypes::RegionType& region) const
2036 {
2037 throw DataException("getSlice - not implemented for Lazy objects.");
2038 }
2039
2040
2041 // To do this we need to rely on our child nodes
2042 DataTypes::RealVectorType::size_type
2043 DataLazy::getPointOffset(int sampleNo,
2044 int dataPointNo)
2045 {
2046 if (m_op==IDENTITY)
2047 {
2048 return m_id->getPointOffset(sampleNo,dataPointNo);
2049 }
2050 if (m_readytype!='E')
2051 {
2052 collapse();
2053 return m_id->getPointOffset(sampleNo,dataPointNo);
2054 }
2055 // at this point we do not have an identity node and the expression will be Expanded
2056 // so we only need to know which child to ask
2057 if (m_left->m_readytype=='E')
2058 {
2059 return m_left->getPointOffset(sampleNo,dataPointNo);
2060 }
2061 else
2062 {
2063 return m_right->getPointOffset(sampleNo,dataPointNo);
2064 }
2065 }
2066
2067 // To do this we need to rely on our child nodes
2068 DataTypes::RealVectorType::size_type
2069 DataLazy::getPointOffset(int sampleNo,
2070 int dataPointNo) const
2071 {
2072 if (m_op==IDENTITY)
2073 {
2074 return m_id->getPointOffset(sampleNo,dataPointNo);
2075 }
2076 if (m_readytype=='E')
2077 {
2078 // at this point we do not have an identity node and the expression will be Expanded
2079 // so we only need to know which child to ask
2080 if (m_left->m_readytype=='E')
2081 {
2082 return m_left->getPointOffset(sampleNo,dataPointNo);
2083 }
2084 else
2085 {
2086 return m_right->getPointOffset(sampleNo,dataPointNo);
2087 }
2088 }
2089 if (m_readytype=='C')
2090 {
2091 return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
2092 }
2093 throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
2094 }
2095
2096
2097 // I have decided to let Data:: handle this issue.
2098 void
2099 DataLazy::setToZero()
2100 {
2101 // DataTypes::RealVectorType v(getNoValues(),0);
2102 // m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
2103 // m_op=IDENTITY;
2104 // m_right.reset();
2105 // m_left.reset();
2106 // m_readytype='C';
2107 // m_buffsRequired=1;
2108
2109 (void)privdebug; // to stop the compiler complaining about unused privdebug
2110 throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
2111 }
2112
2113 bool
2114 DataLazy::actsExpanded() const
2115 {
2116 return (m_readytype=='E');
2117 }
2118
2119 } // end namespace
2120

  ViewVC Help
Powered by ViewVC 1.1.26