/[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 6516 - (show annotations)
Mon Mar 6 03:55:11 2017 UTC (9 months, 1 week ago) by jfenwick
File size: 64262 byte(s)
Properly detect when lazy is complex.

Also, let the code know that anti-hermitian is a valid op

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

  ViewVC Help
Powered by ViewVC 1.1.26