/[escript]/trunk/escriptcore/src/DataLazy.cpp
ViewVC logotype

Contents of /trunk/escriptcore/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6538 - (show annotations)
Tue Mar 14 03:49:20 2017 UTC (2 years, 1 month ago) by jfenwick
File size: 92427 byte(s)
Fix whereZero bugs

Add new opgroup G_UNARY_PR (always real) and put the
whereZero whereNonZero ops in it.

Reorder some declarations to stop the compiler complaining.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2017 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
114 DataLazy_ptr makePromote(DataLazy_ptr p)
115 {
116 if (p->isComplex())
117 {
118 return p;
119 }
120 DataLazy* temp=new DataLazy(p, PROM);
121 return DataLazy_ptr(temp);
122 }
123
124
125 namespace
126 {
127
128
129 // enabling this will print out when ever the maximum stacksize used by resolve increases
130 // it assumes _OPENMP is also in use
131 //#define LAZY_STACK_PROF
132
133
134
135 #ifndef _OPENMP
136 #ifdef LAZY_STACK_PROF
137 #undef LAZY_STACK_PROF
138 #endif
139 #endif
140
141
142 #ifdef LAZY_STACK_PROF
143 std::vector<void*> stackstart(getNumberOfThreads());
144 std::vector<void*> stackend(getNumberOfThreads());
145 size_t maxstackuse=0;
146 #endif
147
148
149 inline int max3(int a, int b, int c)
150 {
151 int t=(a>b?a:b);
152 return (t>c?t:c);
153
154 }
155
156 // return the FunctionSpace of the result of "left op right"
157 FunctionSpace
158 resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
159 {
160 // perhaps this should call interpolate and throw or something?
161 // maybe we need an interpolate node -
162 // that way, if interpolate is required in any other op we can just throw a
163 // programming error exception.
164
165 FunctionSpace l=left->getFunctionSpace();
166 FunctionSpace r=right->getFunctionSpace();
167 if (l!=r)
168 {
169 signed char res=r.getDomain()->preferredInterpolationOnDomain(r.getTypeCode(), l.getTypeCode());
170 if (res==1)
171 {
172 return l;
173 }
174 if (res==-1)
175 {
176 return r;
177 }
178 throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
179 }
180 return l;
181 }
182
183 // return the shape of the result of "left op right"
184 // the shapes resulting from tensor product are more complex to compute so are worked out elsewhere
185 DataTypes::ShapeType
186 resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
187 {
188 if (left->getShape()!=right->getShape())
189 {
190 if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
191 {
192 throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
193 }
194
195 if (left->getRank()==0) // we need to allow scalar * anything
196 {
197 return right->getShape();
198 }
199 if (right->getRank()==0)
200 {
201 return left->getShape();
202 }
203 throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
204 }
205 return left->getShape();
206 }
207
208 // return the shape for "op left"
209
210 DataTypes::ShapeType
211 resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
212 {
213 switch(op)
214 {
215 case TRANS:
216 { // for the scoping of variables
217 const DataTypes::ShapeType& s=left->getShape();
218 DataTypes::ShapeType sh;
219 int rank=left->getRank();
220 if (axis_offset<0 || axis_offset>rank)
221 {
222 stringstream e;
223 e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;
224 throw DataException(e.str());
225 }
226 for (int i=0; i<rank; i++)
227 {
228 int index = (axis_offset+i)%rank;
229 sh.push_back(s[index]); // Append to new shape
230 }
231 return sh;
232 }
233 break;
234 case TRACE:
235 {
236 int rank=left->getRank();
237 if (rank<2)
238 {
239 throw DataException("Trace can only be computed for objects with rank 2 or greater.");
240 }
241 if ((axis_offset>rank-2) || (axis_offset<0))
242 {
243 throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
244 }
245 if (rank==2)
246 {
247 return DataTypes::scalarShape;
248 }
249 else if (rank==3)
250 {
251 DataTypes::ShapeType sh;
252 if (axis_offset==0)
253 {
254 sh.push_back(left->getShape()[2]);
255 }
256 else // offset==1
257 {
258 sh.push_back(left->getShape()[0]);
259 }
260 return sh;
261 }
262 else if (rank==4)
263 {
264 DataTypes::ShapeType sh;
265 const DataTypes::ShapeType& s=left->getShape();
266 if (axis_offset==0)
267 {
268 sh.push_back(s[2]);
269 sh.push_back(s[3]);
270 }
271 else if (axis_offset==1)
272 {
273 sh.push_back(s[0]);
274 sh.push_back(s[3]);
275 }
276 else // offset==2
277 {
278 sh.push_back(s[0]);
279 sh.push_back(s[1]);
280 }
281 return sh;
282 }
283 else // unknown rank
284 {
285 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
286 }
287 }
288 break;
289 default:
290 throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
291 }
292 }
293
294 DataTypes::ShapeType
295 SwapShape(DataAbstract_ptr left, const int axis0, const int axis1)
296 {
297 // This code taken from the Data.cpp swapaxes() method
298 // Some of the checks are probably redundant here
299 int axis0_tmp,axis1_tmp;
300 const DataTypes::ShapeType& s=left->getShape();
301 DataTypes::ShapeType out_shape;
302 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
303 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
304 int rank=left->getRank();
305 if (rank<2) {
306 throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
307 }
308 if (axis0<0 || axis0>rank-1) {
309 stringstream e;
310 e << "Error - Data::swapaxes: axis0 must be between 0 and rank-1=" << (rank-1);
311 throw DataException(e.str());
312 }
313 if (axis1<0 || axis1>rank-1) {
314 stringstream e;
315 e << "Error - Data::swapaxes: axis1 must be between 0 and rank-1=" << (rank-1);
316 throw DataException(e.str());
317 }
318 if (axis0 == axis1) {
319 throw DataException("Error - Data::swapaxes: axis indices must be different.");
320 }
321 if (axis0 > axis1) {
322 axis0_tmp=axis1;
323 axis1_tmp=axis0;
324 } else {
325 axis0_tmp=axis0;
326 axis1_tmp=axis1;
327 }
328 for (int i=0; i<rank; i++) {
329 if (i == axis0_tmp) {
330 out_shape.push_back(s[axis1_tmp]);
331 } else if (i == axis1_tmp) {
332 out_shape.push_back(s[axis0_tmp]);
333 } else {
334 out_shape.push_back(s[i]);
335 }
336 }
337 return out_shape;
338 }
339
340
341 // determine the output shape for the general tensor product operation
342 // the additional parameters return information required later for the product
343 // the majority of this code is copy pasted from C_General_Tensor_Product
344 DataTypes::ShapeType
345 GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
346 {
347
348 // Get rank and shape of inputs
349 int rank0 = left->getRank();
350 int rank1 = right->getRank();
351 const DataTypes::ShapeType& shape0 = left->getShape();
352 const DataTypes::ShapeType& shape1 = right->getShape();
353
354 // Prepare for the loops of the product and verify compatibility of shapes
355 int start0=0, start1=0;
356 if (transpose == 0) {}
357 else if (transpose == 1) { start0 = axis_offset; }
358 else if (transpose == 2) { start1 = rank1-axis_offset; }
359 else { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
360
361 if (rank0<axis_offset)
362 {
363 throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
364 }
365
366 // Adjust the shapes for transpose
367 DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
368 DataTypes::ShapeType tmpShape1(rank1); // than using push_back
369 for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
370 for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
371
372 // Prepare for the loops of the product
373 SL=1, SM=1, SR=1;
374 for (int i=0; i<rank0-axis_offset; i++) {
375 SL *= tmpShape0[i];
376 }
377 for (int i=rank0-axis_offset; i<rank0; i++) {
378 if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
379 throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
380 }
381 SM *= tmpShape0[i];
382 }
383 for (int i=axis_offset; i<rank1; i++) {
384 SR *= tmpShape1[i];
385 }
386
387 // Define the shape of the output (rank of shape is the sum of the loop ranges below)
388 DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
389 { // block to limit the scope of out_index
390 int out_index=0;
391 for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
392 for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
393 }
394
395 if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
396 {
397 ostringstream os;
398 os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
399 throw DataException(os.str());
400 }
401
402 return shape2;
403 }
404
405 } // end anonymous namespace
406
407 void DataLazy::LazyNodeSetup()
408 {
409 #ifdef _OPENMP
410 int numthreads=omp_get_max_threads();
411 if (m_iscompl)
412 {
413 m_samples_c.resize(numthreads*m_samplesize);
414 }
415 else
416 {
417 m_samples_r.resize(numthreads*m_samplesize);
418 }
419 m_sampleids=new int[numthreads];
420 for (int i=0;i<numthreads;++i)
421 {
422 m_sampleids[i]=-1;
423 }
424 #else
425 if (m_iscompl)
426 {
427 m_samples_c.resize(m_samplesize);
428 }
429 else
430 {
431 m_samples_r.resize(m_samplesize);
432 }
433 m_sampleids=new int[1];
434 m_sampleids[0]=-1;
435 #endif // _OPENMP
436 }
437
438
439 // Creates an identity node
440 DataLazy::DataLazy(DataAbstract_ptr p)
441 : parent(p->getFunctionSpace(),p->getShape())
442 ,m_sampleids(0),
443 m_samples_r(1),
444 m_op(IDENTITY),
445 m_opgroup(getOpgroup(m_op))
446 {
447 if (p->isLazy())
448 {
449 // I don't want identity of Lazy.
450 // Question: Why would that be so bad?
451 // Answer: We assume that the child of ID is something we can call getVector on
452 throw DataException("Programmer error - attempt to create identity from a DataLazy.");
453 }
454 else
455 {
456 DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
457 makeIdentity(dr);
458 LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
459 }
460 LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
461 }
462
463 // Wrapping a unary operation
464 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
465 : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),
466 m_op(op),
467 m_opgroup(getOpgroup(m_op)),
468 m_axis_offset(0),
469 m_transpose(0),
470 m_SL(0), m_SM(0), m_SR(0)
471 {
472 ES_opgroup gop=getOpgroup(op);
473 if ((gop!=G_UNARY) && (gop!=G_NP1OUT) && (gop!=G_REDUCTION) && (gop!=G_UNARY_C) && (gop!=G_UNARY_R))
474 {
475 throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
476 }
477 DataLazy_ptr lleft;
478 if (!left->isLazy())
479 {
480 lleft=DataLazy_ptr(new DataLazy(left));
481 }
482 else
483 {
484 lleft=dynamic_pointer_cast<DataLazy>(left);
485 }
486 m_readytype=lleft->m_readytype;
487 m_left=lleft;
488 m_samplesize=getNumDPPSample()*getNoValues();
489 m_children=m_left->m_children+1;
490 m_height=m_left->m_height+1;
491 if (gop==G_UNARY_R)
492 {
493 m_iscompl=false;
494 }
495 else if (gop==G_UNARY_C)
496 {
497 m_iscompl=true;
498 }
499 else
500 {
501 m_iscompl=left->isComplex();
502 }
503 LazyNodeSetup();
504 SIZELIMIT
505 }
506
507
508 // In this constructor we need to consider interpolation and promotion
509 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
510 : parent(resultFS(left,right,op), resultShape(left,right,op)),
511 m_op(op),
512 m_opgroup(getOpgroup(m_op)),
513 m_SL(0), m_SM(0), m_SR(0)
514 {
515 LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
516 if ((getOpgroup(op)!=G_BINARY))
517 {
518 throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
519 }
520
521 if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
522 {
523 FunctionSpace fs=getFunctionSpace();
524 Data ltemp(left);
525 Data tmp(ltemp,fs);
526 left=tmp.borrowDataPtr();
527 }
528 if (getFunctionSpace()!=right->getFunctionSpace()) // right needs to be interpolated
529 {
530 Data tmp(Data(right),getFunctionSpace());
531 right=tmp.borrowDataPtr();
532 LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
533 }
534 left->operandCheck(*right);
535
536 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
537 {
538 m_left=dynamic_pointer_cast<DataLazy>(left);
539 LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
540 }
541 else
542 {
543 m_left=DataLazy_ptr(new DataLazy(left));
544 LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
545 }
546 if (right->isLazy())
547 {
548 m_right=dynamic_pointer_cast<DataLazy>(right);
549 LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
550 }
551 else
552 {
553 m_right=DataLazy_ptr(new DataLazy(right));
554 LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
555 }
556 char lt=m_left->m_readytype;
557 char rt=m_right->m_readytype;
558 if (lt=='E' || rt=='E')
559 {
560 m_readytype='E';
561 }
562 else if (lt=='T' || rt=='T')
563 {
564 m_readytype='T';
565 }
566 else
567 {
568 m_readytype='C';
569 }
570 m_samplesize=getNumDPPSample()*getNoValues();
571 m_children=m_left->m_children+m_right->m_children+2;
572 m_height=max(m_left->m_height,m_right->m_height)+1;
573
574 // now we need to work out if we need to promote anything
575 if (left->isComplex()!=right->isComplex())
576 {
577 if (left->isComplex())
578 {
579 m_right=makePromote(m_right);
580 }
581 else
582 {
583 m_left=makePromote(m_left);
584 }
585 }
586 m_iscompl=m_left->isComplex();
587 LazyNodeSetup();
588 SIZELIMIT
589 LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
590 }
591
592 // need to consider promotion
593 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
594 : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
595 m_op(op),
596 m_opgroup(getOpgroup(m_op)),
597 m_axis_offset(axis_offset),
598 m_transpose(transpose)
599 {
600 if ((getOpgroup(op)!=G_TENSORPROD))
601 {
602 throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
603 }
604 if ((transpose>2) || (transpose<0))
605 {
606 throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
607 }
608 if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
609 {
610 FunctionSpace fs=getFunctionSpace();
611 Data ltemp(left);
612 Data tmp(ltemp,fs);
613 left=tmp.borrowDataPtr();
614 }
615 if (getFunctionSpace()!=right->getFunctionSpace()) // right needs to be interpolated
616 {
617 Data tmp(Data(right),getFunctionSpace());
618 right=tmp.borrowDataPtr();
619 }
620 // left->operandCheck(*right);
621
622 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
623 {
624 m_left=dynamic_pointer_cast<DataLazy>(left);
625 }
626 else
627 {
628 m_left=DataLazy_ptr(new DataLazy(left));
629 }
630 if (right->isLazy())
631 {
632 m_right=dynamic_pointer_cast<DataLazy>(right);
633 }
634 else
635 {
636 m_right=DataLazy_ptr(new DataLazy(right));
637 }
638 char lt=m_left->m_readytype;
639 char rt=m_right->m_readytype;
640 if (lt=='E' || rt=='E')
641 {
642 m_readytype='E';
643 }
644 else if (lt=='T' || rt=='T')
645 {
646 m_readytype='T';
647 }
648 else
649 {
650 m_readytype='C';
651 }
652 m_samplesize=getNumDPPSample()*getNoValues();
653 m_children=m_left->m_children+m_right->m_children+2;
654 m_height=max(m_left->m_height,m_right->m_height)+1;
655
656 // now we need to work out if we need to promote anything
657 if (left->isComplex()!=right->isComplex())
658 {
659 if (left->isComplex())
660 {
661 m_right=makePromote(m_right);
662 }
663 else
664 {
665 m_left=makePromote(m_left);
666 }
667 }
668 m_iscompl=m_left->isComplex();
669 LazyNodeSetup();
670 SIZELIMIT
671 LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
672 }
673
674
675 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
676 : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
677 m_op(op),
678 m_opgroup(getOpgroup(m_op)),
679 m_axis_offset(axis_offset),
680 m_transpose(0),
681 m_tol(0)
682 {
683 if ((getOpgroup(op)!=G_NP1OUT_P))
684 {
685 throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
686 }
687 DataLazy_ptr lleft;
688 if (!left->isLazy())
689 {
690 lleft=DataLazy_ptr(new DataLazy(left));
691 }
692 else
693 {
694 lleft=dynamic_pointer_cast<DataLazy>(left);
695 }
696 m_readytype=lleft->m_readytype;
697 m_left=lleft;
698 m_samplesize=getNumDPPSample()*getNoValues();
699 m_children=m_left->m_children+1;
700 m_height=m_left->m_height+1;
701 m_iscompl=left->isComplex();
702 LazyNodeSetup();
703 SIZELIMIT
704 LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
705 }
706
707 // G_UNARY_P or G_UNARY_PR
708 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
709 : parent(left->getFunctionSpace(), left->getShape()),
710 m_op(op),
711 m_opgroup(getOpgroup(m_op)),
712 m_axis_offset(0),
713 m_transpose(0),
714 m_tol(tol)
715 {
716 if ((m_opgroup!=G_UNARY_P) && (m_opgroup!=G_UNARY_PR))
717 {
718 throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
719 }
720 DataLazy_ptr lleft;
721 if (!left->isLazy())
722 {
723 lleft=DataLazy_ptr(new DataLazy(left));
724 }
725 else
726 {
727 lleft=dynamic_pointer_cast<DataLazy>(left);
728 }
729 m_readytype=lleft->m_readytype;
730 m_left=lleft;
731 m_samplesize=getNumDPPSample()*getNoValues();
732 m_children=m_left->m_children+1;
733 m_height=m_left->m_height+1;
734 if (m_opgroup==G_UNARY_PR)
735 {
736 m_iscompl=false;
737 }
738 else
739 {
740 m_iscompl=left->isComplex();
741 }
742 LazyNodeSetup();
743 SIZELIMIT
744 LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
745 }
746
747
748 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, const int axis0, const int axis1)
749 : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
750 m_op(op),
751 m_opgroup(getOpgroup(m_op)),
752 m_axis_offset(axis0),
753 m_transpose(axis1),
754 m_tol(0)
755 {
756 if ((getOpgroup(op)!=G_NP1OUT_2P))
757 {
758 throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require two integer parameters.");
759 }
760 DataLazy_ptr lleft;
761 if (!left->isLazy())
762 {
763 lleft=DataLazy_ptr(new DataLazy(left));
764 }
765 else
766 {
767 lleft=dynamic_pointer_cast<DataLazy>(left);
768 }
769 m_readytype=lleft->m_readytype;
770 m_left=lleft;
771 m_samplesize=getNumDPPSample()*getNoValues();
772 m_children=m_left->m_children+1;
773 m_height=m_left->m_height+1;
774 m_iscompl=left->isComplex();
775 LazyNodeSetup();
776 SIZELIMIT
777 LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)
778 }
779
780
781
782
783 DataLazy::DataLazy(DataAbstract_ptr mask, DataAbstract_ptr left, DataAbstract_ptr right/*, double tol*/)
784 : parent(left->getFunctionSpace(), left->getShape()),
785 m_op(CONDEVAL),
786 m_opgroup(getOpgroup(m_op)),
787 m_axis_offset(0),
788 m_transpose(0),
789 m_tol(0)
790 {
791 DataLazy_ptr lmask;
792 DataLazy_ptr lleft;
793 DataLazy_ptr lright;
794 if (!mask->isLazy())
795 {
796 lmask=DataLazy_ptr(new DataLazy(mask));
797 }
798 else
799 {
800 lmask=dynamic_pointer_cast<DataLazy>(mask);
801 }
802 if (!left->isLazy())
803 {
804 lleft=DataLazy_ptr(new DataLazy(left));
805 }
806 else
807 {
808 lleft=dynamic_pointer_cast<DataLazy>(left);
809 }
810 if (!right->isLazy())
811 {
812 lright=DataLazy_ptr(new DataLazy(right));
813 }
814 else
815 {
816 lright=dynamic_pointer_cast<DataLazy>(right);
817 }
818 m_readytype=lmask->m_readytype;
819 if ((lleft->m_readytype!=lright->m_readytype) || (lmask->m_readytype!=lleft->m_readytype))
820 {
821 throw DataException("Programmer Error - condEval arguments must have the same readytype");
822 }
823 m_left=lleft;
824 m_right=lright;
825 m_mask=lmask;
826 m_samplesize=getNumDPPSample()*getNoValues();
827 m_children=m_left->m_children+m_right->m_children+m_mask->m_children+1;
828 m_height=max3(m_left->m_height,m_right->m_height,m_mask->m_height)+1;
829
830 // now we need to work out if we need to promote anything
831 if (left->isComplex()!=right->isComplex())
832 {
833 if (left->isComplex())
834 {
835 m_right=makePromote(m_right);
836 }
837 else
838 {
839 m_left=makePromote(m_left);
840 }
841 }
842 m_iscompl=left->isComplex();
843 LazyNodeSetup();
844 SIZELIMIT
845 LAZYDEBUG(cout << "(8)Lazy created with " << m_samplesize << endl;)
846 }
847
848
849
850 DataLazy::~DataLazy()
851 {
852 delete[] m_sampleids;
853 }
854
855
856 /*
857 \brief Evaluates the expression using methods on Data.
858 This does the work for the collapse method.
859 For reasons of efficiency do not call this method on DataExpanded nodes.
860 */
861 DataReady_ptr
862 DataLazy::collapseToReady() const
863 {
864 if (m_readytype=='E')
865 { // this is more an efficiency concern than anything else
866 throw DataException("Programmer Error - do not use collapse on Expanded data.");
867 }
868 if (m_op==IDENTITY)
869 {
870 return m_id;
871 }
872 DataReady_ptr pleft=m_left->collapseToReady();
873 Data left(pleft);
874 Data right;
875 if ((m_opgroup==G_BINARY) || (m_opgroup==G_TENSORPROD))
876 {
877 right=Data(m_right->collapseToReady());
878 }
879 Data result;
880 switch(m_op)
881 {
882 case ADD:
883 result=left+right;
884 break;
885 case SUB:
886 result=left-right;
887 break;
888 case MUL:
889 result=left*right;
890 break;
891 case DIV:
892 result=left/right;
893 break;
894 case SIN:
895 result=left.sin();
896 break;
897 case COS:
898 result=left.cos();
899 break;
900 case TAN:
901 result=left.tan();
902 break;
903 case ASIN:
904 result=left.asin();
905 break;
906 case ACOS:
907 result=left.acos();
908 break;
909 case ATAN:
910 result=left.atan();
911 break;
912 case SINH:
913 result=left.sinh();
914 break;
915 case COSH:
916 result=left.cosh();
917 break;
918 case TANH:
919 result=left.tanh();
920 break;
921 case ERF:
922 result=left.erf();
923 break;
924 case ASINH:
925 result=left.asinh();
926 break;
927 case ACOSH:
928 result=left.acosh();
929 break;
930 case ATANH:
931 result=left.atanh();
932 break;
933 case LOG10:
934 result=left.log10();
935 break;
936 case LOG:
937 result=left.log();
938 break;
939 case SIGN:
940 result=left.sign();
941 break;
942 case ABS:
943 result=left.abs();
944 break;
945 case NEG:
946 result=left.neg();
947 break;
948 case POS:
949 // it doesn't mean anything for delayed.
950 // it will just trigger a deep copy of the lazy object
951 throw DataException("Programmer error - POS not supported for lazy data.");
952 break;
953 case EXP:
954 result=left.exp();
955 break;
956 case SQRT:
957 result=left.sqrt();
958 break;
959 case RECIP:
960 result=left.oneOver();
961 break;
962 case GZ:
963 result=left.wherePositive();
964 break;
965 case LZ:
966 result=left.whereNegative();
967 break;
968 case GEZ:
969 result=left.whereNonNegative();
970 break;
971 case LEZ:
972 result=left.whereNonPositive();
973 break;
974 case NEZ:
975 result=left.whereNonZero(m_tol);
976 break;
977 case EZ:
978 result=left.whereZero(m_tol);
979 break;
980 case SYM:
981 result=left.symmetric();
982 break;
983 case NSYM:
984 result=left.antisymmetric();
985 break;
986 case PROD:
987 result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
988 break;
989 case TRANS:
990 result=left.transpose(m_axis_offset);
991 break;
992 case TRACE:
993 result=left.trace(m_axis_offset);
994 break;
995 case SWAP:
996 result=left.swapaxes(m_axis_offset, m_transpose);
997 break;
998 case MINVAL:
999 result=left.minval();
1000 break;
1001 case MAXVAL:
1002 result=left.minval();
1003 break;
1004 case HER:
1005 result=left.hermitian();
1006 break;
1007 case NHER:
1008 result=left.antihermitian();
1009 break;
1010 case PROM:
1011 result.copy(left);
1012 result.complicate();
1013 break;
1014 break;
1015 default:
1016 throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
1017 }
1018 return result.borrowReadyPtr();
1019 }
1020
1021 /*
1022 \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
1023 This method uses the original methods on the Data class to evaluate the expressions.
1024 For this reason, it should not be used on DataExpanded instances. (To do so would defeat
1025 the purpose of using DataLazy in the first place).
1026 */
1027 void
1028 DataLazy::collapse() const
1029 {
1030 if (m_op==IDENTITY)
1031 {
1032 return;
1033 }
1034 if (m_readytype=='E')
1035 { // this is more an efficiency concern than anything else
1036 throw DataException("Programmer Error - do not use collapse on Expanded data.");
1037 }
1038 m_id=collapseToReady();
1039 m_op=IDENTITY;
1040 }
1041
1042 // The result will be stored in m_samples
1043 // The return value is a pointer to the DataVector, offset is the offset within the return value
1044 const DataTypes::RealVectorType*
1045 DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset) const
1046 {
1047 LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1048 // collapse so we have a 'E' node or an IDENTITY for some other type
1049 if (m_readytype!='E' && m_op!=IDENTITY)
1050 {
1051 collapse();
1052 }
1053 if (m_op==IDENTITY)
1054 {
1055 const RealVectorType& vec=m_id->getVectorRO();
1056 roffset=m_id->getPointOffset(sampleNo, 0);
1057 #ifdef LAZY_STACK_PROF
1058 int x;
1059 if (&x<stackend[omp_get_thread_num()])
1060 {
1061 stackend[omp_get_thread_num()]=&x;
1062 }
1063 #endif
1064 return &(vec);
1065 }
1066 if (m_readytype!='E')
1067 {
1068 throw DataException("Programmer Error - Collapse did not produce an expanded node.");
1069 }
1070 if (m_sampleids[tid]==sampleNo)
1071 {
1072 roffset=tid*m_samplesize;
1073 return &(m_samples_r); // sample is already resolved
1074 }
1075 m_sampleids[tid]=sampleNo;
1076
1077 switch (m_opgroup)
1078 {
1079 // note that G_UNARY_C is not listed here - this is deliberate
1080 // by definition this type can not produce a real_t answer
1081 // so processing it here would be a bug
1082 case G_UNARY:
1083 case G_UNARY_R:
1084 case G_UNARY_PR:
1085 case G_UNARY_P: return resolveNodeUnary(tid, sampleNo, roffset);
1086 case G_BINARY: return resolveNodeBinary(tid, sampleNo, roffset);
1087 case G_NP1OUT: return resolveNodeNP1OUT(tid, sampleNo, roffset);
1088 case G_NP1OUT_P: return resolveNodeNP1OUT_P(tid, sampleNo, roffset);
1089 case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);
1090 case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);
1091 case G_REDUCTION: return resolveNodeReduction(tid, sampleNo, roffset);
1092 case G_CONDEVAL: return resolveNodeCondEval(tid, sampleNo, roffset);
1093 default:
1094 throw DataException("Programmer Error - resolveNodeSample does not know how to process "+opToString(m_op)+" in group "+groupToString(m_opgroup)+".");
1095 }
1096 }
1097
1098 // The result will be stored in m_samples
1099 // The return value is a pointer to the DataVector, offset is the offset within the return value
1100 const DataTypes::CplxVectorType*
1101 DataLazy::resolveNodeSampleCplx(int tid, int sampleNo, size_t& roffset) const
1102 {
1103 LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1104 // collapse so we have a 'E' node or an IDENTITY for some other type
1105 if (m_readytype!='E' && m_op!=IDENTITY)
1106 {
1107 collapse();
1108 }
1109 if (m_op==IDENTITY)
1110 {
1111 const CplxVectorType& vec=m_id->getVectorROC();
1112 roffset=m_id->getPointOffset(sampleNo, 0);
1113 #ifdef LAZY_STACK_PROF
1114 int x;
1115 if (&x<stackend[omp_get_thread_num()])
1116 {
1117 stackend[omp_get_thread_num()]=&x;
1118 }
1119 #endif
1120 return &(vec);
1121 }
1122 if (m_readytype!='E')
1123 {
1124 throw DataException("Programmer Error - Collapse did not produce an expanded node.");
1125 }
1126 if (m_sampleids[tid]==sampleNo)
1127 {
1128 roffset=tid*m_samplesize;
1129 return &(m_samples_c); // sample is already resolved
1130 }
1131 m_sampleids[tid]=sampleNo;
1132
1133 switch (m_opgroup)
1134 {
1135 case G_UNARY:
1136 case G_UNARY_P: return resolveNodeUnaryCplx(tid, sampleNo, roffset);
1137 case G_UNARY_C: return resolveNodeUnary_C(tid, sampleNo, roffset);
1138 case G_BINARY: return resolveNodeBinaryCplx(tid, sampleNo, roffset);
1139 case G_NP1OUT: return resolveNodeNP1OUTCplx(tid, sampleNo, roffset);
1140 case G_NP1OUT_P: return resolveNodeNP1OUT_PCplx(tid, sampleNo, roffset);
1141 case G_TENSORPROD: return resolveNodeTProdCplx(tid, sampleNo, roffset);
1142 case G_NP1OUT_2P: return resolveNodeNP1OUT_2PCplx(tid, sampleNo, roffset);
1143 case G_REDUCTION: return resolveNodeReductionCplx(tid, sampleNo, roffset);
1144 case G_CONDEVAL: return resolveNodeCondEvalCplx(tid, sampleNo, roffset);
1145 default:
1146 throw DataException("Programmer Error - resolveNodeSampleCplx does not know how to process "+opToString(m_op)+".");
1147 }
1148 }
1149
1150 const DataTypes::RealVectorType*
1151 DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset) const
1152 {
1153 // we assume that any collapsing has been done before we get here
1154 // since we only have one argument we don't need to think about only
1155 // processing single points.
1156 // we will also know we won't get identity nodes
1157 if (m_readytype!='E')
1158 {
1159 throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1160 }
1161 if (m_op==IDENTITY)
1162 {
1163 throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1164 }
1165 if (m_op==POS)
1166 {
1167 // this should be prevented earlier
1168 // operation is meaningless for lazy
1169 throw DataException("Programmer error - POS not supported for lazy data.");
1170 }
1171 // The tricky bit here is that we expect a real output, but for G_UNARY_R, the child nodes could be complex
1172 if (((m_opgroup==G_UNARY_R) || (m_opgroup==G_UNARY_PR)) && (m_left->isComplex()))
1173 {
1174 const DataTypes::CplxVectorType* leftres=m_left->resolveNodeSampleCplx(tid, sampleNo, roffset);
1175 const cplx_t* left=&((*leftres)[roffset]);
1176 roffset=m_samplesize*tid;
1177 real_t* result=&(m_samples_r[roffset]);
1178 tensor_unary_array_operation_real(m_samplesize,
1179 left,
1180 result,
1181 m_op,
1182 m_tol);
1183 }
1184 else
1185 {
1186 const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
1187 const real_t* left=&((*leftres)[roffset]);
1188 roffset=m_samplesize*tid;
1189 real_t* result=&(m_samples_r[roffset]);
1190 tensor_unary_array_operation(m_samplesize,
1191 left,
1192 result,
1193 m_op,
1194 m_tol);
1195 }
1196 return &(m_samples_r);
1197 }
1198
1199 const DataTypes::CplxVectorType*
1200 DataLazy::resolveNodeUnaryCplx(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 if (m_op==POS)
1215 {
1216 // this should be prevented earlier
1217 // operation is meaningless for lazy
1218 throw DataException("Programmer error - POS not supported for lazy data.");
1219 }
1220
1221 roffset=m_samplesize*tid;
1222 DataTypes::cplx_t* result=&(m_samples_c[roffset]);
1223 if (m_op==PROM)
1224 {
1225 size_t subroffset=0;
1226 const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1227 const DataTypes::real_t* left=&((*leftres)[subroffset]);
1228
1229 tensor_unary_promote(m_samplesize, left, result);
1230 }
1231 else
1232 {
1233 size_t subroffset=0;
1234 const DataTypes::CplxVectorType* leftres=m_left->resolveNodeSampleCplx(tid, sampleNo, subroffset);
1235 const DataTypes::cplx_t* left=&((*leftres)[subroffset]);
1236 tensor_unary_array_operation(m_samplesize,
1237 left,
1238 result,
1239 m_op,
1240 m_tol);
1241 }
1242 return &(m_samples_c);
1243 }
1244
1245 const DataTypes::RealVectorType*
1246 DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset) const
1247 {
1248 // we assume that any collapsing has been done before we get here
1249 // since we only have one argument we don't need to think about only
1250 // processing single points.
1251 // we will also know we won't get identity nodes
1252 if (m_readytype!='E')
1253 {
1254 throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1255 }
1256 if (m_op==IDENTITY)
1257 {
1258 throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1259 }
1260 size_t loffset=0;
1261 const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);
1262
1263 roffset=m_samplesize*tid;
1264 unsigned int ndpps=getNumDPPSample();
1265 unsigned int psize=DataTypes::noValues(m_left->getShape());
1266 double* result=&(m_samples_r[roffset]);
1267 switch (m_op)
1268 {
1269 case MINVAL:
1270 {
1271 for (unsigned int z=0;z<ndpps;++z)
1272 {
1273 FMin op;
1274 *result=escript::reductionOpVector(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max());
1275 loffset+=psize;
1276 result++;
1277 }
1278 }
1279 break;
1280 case MAXVAL:
1281 {
1282 for (unsigned int z=0;z<ndpps;++z)
1283 {
1284 FMax op;
1285 *result=escript::reductionOpVector(*leftres, m_left->getShape(), loffset, op, numeric_limits<double>::max()*-1);
1286 loffset+=psize;
1287 result++;
1288 }
1289 }
1290 break;
1291 default:
1292 throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1293 }
1294 return &(m_samples_r);
1295 }
1296
1297 const DataTypes::CplxVectorType*
1298 DataLazy::resolveNodeReductionCplx(int tid, int sampleNo, size_t& roffset) const
1299 {
1300 // we assume that any collapsing has been done before we get here
1301 // since we only have one argument we don't need to think about only
1302 // processing single points.
1303 // we will also know we won't get identity nodes
1304 if (m_readytype!='E')
1305 {
1306 throw DataException("Programmer error - resolveReductionCplx should only be called on expanded Data.");
1307 }
1308 if (m_op==IDENTITY)
1309 {
1310 throw DataException("Programmer error - resolveNodeReductionCplx should not be called on identity nodes.");
1311 }
1312 throw DataException("Programmer error - reduction operations MIN and MAX not supported for complex values.");
1313 }
1314
1315
1316 const DataTypes::RealVectorType*
1317 DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset) const
1318 {
1319 // we assume that any collapsing has been done before we get here
1320 // since we only have one argument we don't need to think about only
1321 // processing single points.
1322 if (m_readytype!='E')
1323 {
1324 throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");
1325 }
1326 if (m_op==IDENTITY)
1327 {
1328 throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1329 }
1330 size_t subroffset;
1331 const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1332 roffset=m_samplesize*tid;
1333 size_t loop=0;
1334 size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1335 size_t step=getNoValues();
1336 size_t offset=roffset;
1337 switch (m_op)
1338 {
1339 case SYM:
1340 for (loop=0;loop<numsteps;++loop)
1341 {
1342 escript::symmetric(*leftres,m_left->getShape(),subroffset, m_samples_r, getShape(), offset);
1343 subroffset+=step;
1344 offset+=step;
1345 }
1346 break;
1347 case NSYM:
1348 for (loop=0;loop<numsteps;++loop)
1349 {
1350 escript::antisymmetric(*leftres,m_left->getShape(),subroffset, m_samples_r, getShape(), offset);
1351 subroffset+=step;
1352 offset+=step;
1353 }
1354 break;
1355 default:
1356 throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1357 }
1358 return &m_samples_r;
1359 }
1360
1361
1362 const DataTypes::CplxVectorType*
1363 DataLazy::resolveNodeUnary_C(int tid, int sampleNo, size_t& roffset) const
1364 {
1365 // we assume that any collapsing has been done before we get here
1366 // since we only have one argument we don't need to think about only
1367 // processing single points.
1368 if (m_readytype!='E')
1369 {
1370 throw DataException("Programmer error - resolveNodeUnaryC should only be called on expanded Data.");
1371 }
1372 if (m_op==IDENTITY)
1373 {
1374 throw DataException("Programmer error - resolveNodeUnaryC should not be called on identity nodes.");
1375 }
1376 if (m_op!=PROM)
1377 {
1378 throw DataException("Programmer error - resolveNodeUnaryC can not resolve operator "+opToString(m_op)+".");
1379 }
1380 const DataTypes::RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
1381 const real_t* left=&((*leftres)[roffset]);
1382 roffset=m_samplesize*tid;
1383 cplx_t* result=&(m_samples_c[roffset]);
1384 // rather than calling the tensor_unary_array_operation, I'm just coding this one here
1385 for (int i = 0; i < m_samplesize; ++i) {
1386 result[i] = static_cast<cplx_t>(left[i]);
1387 }
1388 return &(m_samples_c);
1389 }
1390
1391 const DataTypes::CplxVectorType*
1392 DataLazy::resolveNodeNP1OUTCplx(int tid, int sampleNo, size_t& roffset) const
1393 {
1394 // we assume that any collapsing has been done before we get here
1395 // since we only have one argument we don't need to think about only
1396 // processing single points.
1397 if (m_readytype!='E')
1398 {
1399 throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");
1400 }
1401 if (m_op==IDENTITY)
1402 {
1403 throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1404 }
1405 size_t subroffset;
1406 const CplxVectorType* leftres=m_left->resolveNodeSampleCplx(tid, sampleNo, subroffset);
1407 roffset=m_samplesize*tid;
1408 size_t loop=0;
1409 size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1410 size_t step=getNoValues();
1411 size_t offset=roffset;
1412 switch (m_op)
1413 {
1414 case SYM:
1415 for (loop=0;loop<numsteps;++loop)
1416 {
1417 escript::symmetric(*leftres,m_left->getShape(),subroffset, m_samples_c, getShape(), offset);
1418 subroffset+=step;
1419 offset+=step;
1420 }
1421 break;
1422 case NSYM:
1423 for (loop=0;loop<numsteps;++loop)
1424 {
1425 escript::antisymmetric(*leftres,m_left->getShape(),subroffset, m_samples_c, getShape(), offset);
1426 subroffset+=step;
1427 offset+=step;
1428 }
1429 break;
1430 case HER:
1431 for (loop=0;loop<numsteps;++loop)
1432 {
1433 escript::hermitian(*leftres,m_left->getShape(),subroffset, m_samples_c, getShape(), offset);
1434 subroffset+=step;
1435 offset+=step;
1436 }
1437 break;
1438 case NHER:
1439 for (loop=0;loop<numsteps;++loop)
1440 {
1441 escript::antihermitian(*leftres,m_left->getShape(),subroffset, m_samples_c, getShape(), offset);
1442 subroffset+=step;
1443 offset+=step;
1444 }
1445 break;
1446 default:
1447 throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1448 }
1449 return &m_samples_c;
1450 }
1451
1452 const DataTypes::RealVectorType*
1453 DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset) const
1454 {
1455 // we assume that any collapsing has been done before we get here
1456 // since we only have one argument we don't need to think about only
1457 // processing single points.
1458 if (m_readytype!='E')
1459 {
1460 throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");
1461 }
1462 if (m_op==IDENTITY)
1463 {
1464 throw DataException("Programmer error - resolveNodeNP1OUT_P should not be called on identity nodes.");
1465 }
1466 size_t subroffset;
1467 size_t offset;
1468 const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1469 roffset=m_samplesize*tid;
1470 offset=roffset;
1471 size_t loop=0;
1472 size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1473 size_t outstep=getNoValues();
1474 size_t instep=m_left->getNoValues();
1475 switch (m_op)
1476 {
1477 case TRACE:
1478 for (loop=0;loop<numsteps;++loop)
1479 {
1480 escript::trace(*leftres,m_left->getShape(),subroffset, m_samples_r ,getShape(),offset,m_axis_offset);
1481 subroffset+=instep;
1482 offset+=outstep;
1483 }
1484 break;
1485 case TRANS:
1486 for (loop=0;loop<numsteps;++loop)
1487 {
1488 escript::transpose(*leftres,m_left->getShape(),subroffset, m_samples_r, getShape(),offset,m_axis_offset);
1489 subroffset+=instep;
1490 offset+=outstep;
1491 }
1492 break;
1493 default:
1494 throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1495 }
1496 return &m_samples_r;
1497 }
1498
1499 const DataTypes::CplxVectorType*
1500 DataLazy::resolveNodeNP1OUT_PCplx(int tid, int sampleNo, size_t& roffset) const
1501 {
1502 // we assume that any collapsing has been done before we get here
1503 // since we only have one argument we don't need to think about only
1504 // processing single points.
1505 if (m_readytype!='E')
1506 {
1507 throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");
1508 }
1509 if (m_op==IDENTITY)
1510 {
1511 throw DataException("Programmer error - resolveNodeNP1OUT_P should not be called on identity nodes.");
1512 }
1513 size_t subroffset;
1514 size_t offset;
1515 const CplxVectorType* leftres=m_left->resolveNodeSampleCplx(tid, sampleNo, subroffset);
1516 roffset=m_samplesize*tid;
1517 offset=roffset;
1518 size_t loop=0;
1519 size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1520 size_t outstep=getNoValues();
1521 size_t instep=m_left->getNoValues();
1522 switch (m_op)
1523 {
1524 case TRACE:
1525 for (loop=0;loop<numsteps;++loop)
1526 {
1527 escript::trace(*leftres,m_left->getShape(),subroffset, m_samples_c ,getShape(),offset,m_axis_offset);
1528 subroffset+=instep;
1529 offset+=outstep;
1530 }
1531 break;
1532 case TRANS:
1533 for (loop=0;loop<numsteps;++loop)
1534 {
1535 escript::transpose(*leftres,m_left->getShape(),subroffset, m_samples_c, getShape(),offset,m_axis_offset);
1536 subroffset+=instep;
1537 offset+=outstep;
1538 }
1539 break;
1540 default:
1541 throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1542 }
1543 return &m_samples_c;
1544 }
1545
1546 const DataTypes::RealVectorType*
1547 DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset) const
1548 {
1549 if (m_readytype!='E')
1550 {
1551 throw DataException("Programmer error - resolveNodeNP1OUT_2P should only be called on expanded Data.");
1552 }
1553 if (m_op==IDENTITY)
1554 {
1555 throw DataException("Programmer error - resolveNodeNP1OUT_2P should not be called on identity nodes.");
1556 }
1557 size_t subroffset;
1558 size_t offset;
1559 const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1560 roffset=m_samplesize*tid;
1561 offset=roffset;
1562 size_t loop=0;
1563 size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1564 size_t outstep=getNoValues();
1565 size_t instep=m_left->getNoValues();
1566 switch (m_op)
1567 {
1568 case SWAP:
1569 for (loop=0;loop<numsteps;++loop)
1570 {
1571 escript::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples_r, getShape(),offset, m_axis_offset, m_transpose);
1572 subroffset+=instep;
1573 offset+=outstep;
1574 }
1575 break;
1576 default:
1577 throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1578 }
1579 return &m_samples_r;
1580 }
1581
1582
1583 const DataTypes::CplxVectorType*
1584 DataLazy::resolveNodeNP1OUT_2PCplx(int tid, int sampleNo, size_t& roffset) const
1585 {
1586 if (m_readytype!='E')
1587 {
1588 throw DataException("Programmer error - resolveNodeNP1OUT_2P should only be called on expanded Data.");
1589 }
1590 if (m_op==IDENTITY)
1591 {
1592 throw DataException("Programmer error - resolveNodeNP1OUT_2P should not be called on identity nodes.");
1593 }
1594 size_t subroffset;
1595 size_t offset;
1596 const CplxVectorType* leftres=m_left->resolveNodeSampleCplx(tid, sampleNo, subroffset);
1597 roffset=m_samplesize*tid;
1598 offset=roffset;
1599 size_t loop=0;
1600 size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1601 size_t outstep=getNoValues();
1602 size_t instep=m_left->getNoValues();
1603 switch (m_op)
1604 {
1605 case SWAP:
1606 for (loop=0;loop<numsteps;++loop)
1607 {
1608 escript::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples_c, getShape(),offset, m_axis_offset, m_transpose);
1609 subroffset+=instep;
1610 offset+=outstep;
1611 }
1612 break;
1613 default:
1614 throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1615 }
1616 return &m_samples_c;
1617 }
1618
1619 const DataTypes::RealVectorType*
1620 DataLazy::resolveNodeCondEval(int tid, int sampleNo, size_t& roffset) const
1621 {
1622 if (m_readytype!='E')
1623 {
1624 throw DataException("Programmer error - resolveNodeCondEval should only be called on expanded Data.");
1625 }
1626 if (m_op!=CONDEVAL)
1627 {
1628 throw DataException("Programmer error - resolveNodeCondEval should only be called on CONDEVAL nodes.");
1629 }
1630 size_t subroffset;
1631
1632 const RealVectorType* maskres=m_mask->resolveNodeSample(tid, sampleNo, subroffset);
1633 const RealVectorType* srcres=0;
1634 if ((*maskres)[subroffset]>0)
1635 {
1636 srcres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1637 }
1638 else
1639 {
1640 srcres=m_right->resolveNodeSample(tid, sampleNo, subroffset);
1641 }
1642
1643 // Now we need to copy the result
1644
1645 roffset=m_samplesize*tid;
1646 for (int i=0;i<m_samplesize;++i)
1647 {
1648 m_samples_r[roffset+i]=(*srcres)[subroffset+i];
1649 }
1650
1651 return &m_samples_r;
1652 }
1653
1654 const DataTypes::CplxVectorType*
1655 DataLazy::resolveNodeCondEvalCplx(int tid, int sampleNo, size_t& roffset) const
1656 {
1657 if (m_readytype!='E')
1658 {
1659 throw DataException("Programmer error - resolveNodeCondEval should only be called on expanded Data.");
1660 }
1661 if (m_op!=CONDEVAL)
1662 {
1663 throw DataException("Programmer error - resolveNodeCondEval should only be called on CONDEVAL nodes.");
1664 }
1665 size_t subroffset;
1666
1667 const RealVectorType* maskres=m_mask->resolveNodeSample(tid, sampleNo, subroffset);
1668 const CplxVectorType* srcres=0;
1669 if ((*maskres)[subroffset]>0)
1670 {
1671 srcres=m_left->resolveNodeSampleCplx(tid, sampleNo, subroffset);
1672 }
1673 else
1674 {
1675 srcres=m_right->resolveNodeSampleCplx(tid, sampleNo, subroffset);
1676 }
1677
1678 // Now we need to copy the result
1679
1680 roffset=m_samplesize*tid;
1681 for (int i=0;i<m_samplesize;++i)
1682 {
1683 m_samples_c[roffset+i]=(*srcres)[subroffset+i];
1684 }
1685
1686 return &m_samples_c;
1687 }
1688
1689
1690 // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1691 // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1692 // If both children are expanded, then we can process them in a single operation (we treat
1693 // the whole sample as one big datapoint.
1694 // If one of the children is not expanded, then we need to treat each point in the sample
1695 // individually.
1696 // There is an additional complication when scalar operations are considered.
1697 // For example, 2+Vector.
1698 // In this case each double within the point is treated individually
1699 const DataTypes::RealVectorType*
1700 DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset) const
1701 {
1702 LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1703
1704 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
1705 // first work out which of the children are expanded
1706 bool leftExp=(m_left->m_readytype=='E');
1707 bool rightExp=(m_right->m_readytype=='E');
1708 if (!leftExp && !rightExp)
1709 {
1710 throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1711 }
1712 bool leftScalar=(m_left->getRank()==0);
1713 bool rightScalar=(m_right->getRank()==0);
1714 if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1715 {
1716 throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1717 }
1718 size_t leftsize=m_left->getNoValues();
1719 size_t rightsize=m_right->getNoValues();
1720 size_t chunksize=1; // how many doubles will be processed in one go
1721 int leftstep=0; // how far should the left offset advance after each step
1722 int rightstep=0;
1723 int numsteps=0; // total number of steps for the inner loop
1724 int oleftstep=0; // the o variables refer to the outer loop
1725 int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1726 int onumsteps=1;
1727
1728 bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1729 bool RES=(rightExp && rightScalar);
1730 bool LS=(!leftExp && leftScalar); // left is a single scalar
1731 bool RS=(!rightExp && rightScalar);
1732 bool LN=(!leftExp && !leftScalar); // left is a single non-scalar
1733 bool RN=(!rightExp && !rightScalar);
1734 bool LEN=(leftExp && !leftScalar); // left is an expanded non-scalar
1735 bool REN=(rightExp && !rightScalar);
1736
1737 if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1738 {
1739 chunksize=m_left->getNumDPPSample()*leftsize;
1740 leftstep=0;
1741 rightstep=0;
1742 numsteps=1;
1743 }
1744 else if (LES || RES)
1745 {
1746 chunksize=1;
1747 if (LES) // left is an expanded scalar
1748 {
1749 if (RS)
1750 {
1751 leftstep=1;
1752 rightstep=0;
1753 numsteps=m_left->getNumDPPSample();
1754 }
1755 else // RN or REN
1756 {
1757 leftstep=0;
1758 oleftstep=1;
1759 rightstep=1;
1760 orightstep=(RN ? -(int)rightsize : 0);
1761 numsteps=rightsize;
1762 onumsteps=m_left->getNumDPPSample();
1763 }
1764 }
1765 else // right is an expanded scalar
1766 {
1767 if (LS)
1768 {
1769 rightstep=1;
1770 leftstep=0;
1771 numsteps=m_right->getNumDPPSample();
1772 }
1773 else
1774 {
1775 rightstep=0;
1776 orightstep=1;
1777 leftstep=1;
1778 oleftstep=(LN ? -(int)leftsize : 0);
1779 numsteps=leftsize;
1780 onumsteps=m_right->getNumDPPSample();
1781 }
1782 }
1783 }
1784 else // this leaves (LEN, RS), (LEN, RN) and their transposes
1785 {
1786 if (LEN) // and Right will be a single value
1787 {
1788 chunksize=rightsize;
1789 leftstep=rightsize;
1790 rightstep=0;
1791 numsteps=m_left->getNumDPPSample();
1792 if (RS)
1793 {
1794 numsteps*=leftsize;
1795 }
1796 }
1797 else // REN
1798 {
1799 chunksize=leftsize;
1800 rightstep=leftsize;
1801 leftstep=0;
1802 numsteps=m_right->getNumDPPSample();
1803 if (LS)
1804 {
1805 numsteps*=rightsize;
1806 }
1807 }
1808 }
1809
1810 int resultStep=max(leftstep,rightstep); // only one (at most) should be !=0
1811 // Get the values of sub-expressions
1812 const RealVectorType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);
1813 const RealVectorType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
1814 LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1815 LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1816 LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1817 LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1818 LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1819 LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1820 LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN << endl;)
1821
1822 LAZYDEBUG(cout << "Left res["<< lroffset<< "]=" << (*left)[lroffset] << endl;)
1823 LAZYDEBUG(cout << "Right res["<< rroffset<< "]=" << (*right)[rroffset] << endl;)
1824
1825
1826 roffset=m_samplesize*tid;
1827 double* resultp=&(m_samples_r[roffset]); // results are stored at the vector offset we received
1828 switch(m_op)
1829 {
1830 case ADD:
1831 //PROC_OP(NO_ARG,plus<double>());
1832 escript::binaryOpVectorLazyArithmeticHelper<real_t, real_t, real_t>(resultp,
1833 &(*left)[0],
1834 &(*right)[0],
1835 chunksize,
1836 onumsteps,
1837 numsteps,
1838 resultStep,
1839 leftstep,
1840 rightstep,
1841 oleftstep,
1842 orightstep,
1843 lroffset,
1844 rroffset,
1845 escript::ES_optype::ADD);
1846 break;
1847 case SUB:
1848 escript::binaryOpVectorLazyArithmeticHelper<real_t, real_t, real_t>(resultp,
1849 &(*left)[0],
1850 &(*right)[0],
1851 chunksize,
1852 onumsteps,
1853 numsteps,
1854 resultStep,
1855 leftstep,
1856 rightstep,
1857 oleftstep,
1858 orightstep,
1859 lroffset,
1860 rroffset,
1861 escript::ES_optype::SUB);
1862 //PROC_OP(NO_ARG,minus<double>());
1863 break;
1864 case MUL:
1865 //PROC_OP(NO_ARG,multiplies<double>());
1866 escript::binaryOpVectorLazyArithmeticHelper<real_t, real_t, real_t>(resultp,
1867 &(*left)[0],
1868 &(*right)[0],
1869 chunksize,
1870 onumsteps,
1871 numsteps,
1872 resultStep,
1873 leftstep,
1874 rightstep,
1875 oleftstep,
1876 orightstep,
1877 lroffset,
1878 rroffset,
1879 escript::ES_optype::MUL);
1880 break;
1881 case DIV:
1882 //PROC_OP(NO_ARG,divides<double>());
1883 escript::binaryOpVectorLazyArithmeticHelper<real_t, real_t, real_t>(resultp,
1884 &(*left)[0],
1885 &(*right)[0],
1886 chunksize,
1887 onumsteps,
1888 numsteps,
1889 resultStep,
1890 leftstep,
1891 rightstep,
1892 oleftstep,
1893 orightstep,
1894 lroffset,
1895 rroffset,
1896 escript::ES_optype::DIV);
1897 break;
1898 case POW:
1899 //PROC_OP(double (double,double),::pow);
1900 escript::binaryOpVectorLazyArithmeticHelper<real_t, real_t, real_t>(resultp,
1901 &(*left)[0],
1902 &(*right)[0],
1903 chunksize,
1904 onumsteps,
1905 numsteps,
1906 resultStep,
1907 leftstep,
1908 rightstep,
1909 oleftstep,
1910 orightstep,
1911 lroffset,
1912 rroffset,
1913 escript::ES_optype::POW);
1914 break;
1915 default:
1916 throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1917 }
1918 LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples_r[roffset] << endl;)
1919 return &m_samples_r;
1920 }
1921
1922 const DataTypes::CplxVectorType*
1923 DataLazy::resolveNodeBinaryCplx(int tid, int sampleNo, size_t& roffset) const
1924 {
1925 LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1926
1927 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
1928 // first work out which of the children are expanded
1929 bool leftExp=(m_left->m_readytype=='E');
1930 bool rightExp=(m_right->m_readytype=='E');
1931 if (!leftExp && !rightExp)
1932 {
1933 throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1934 }
1935 bool leftScalar=(m_left->getRank()==0);
1936 bool rightScalar=(m_right->getRank()==0);
1937 if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1938 {
1939 throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1940 }
1941 size_t leftsize=m_left->getNoValues();
1942 size_t rightsize=m_right->getNoValues();
1943 size_t chunksize=1; // how many doubles will be processed in one go
1944 int leftstep=0; // how far should the left offset advance after each step
1945 int rightstep=0;
1946 int numsteps=0; // total number of steps for the inner loop
1947 int oleftstep=0; // the o variables refer to the outer loop
1948 int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1949 int onumsteps=1;
1950
1951 bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1952 bool RES=(rightExp && rightScalar);
1953 bool LS=(!leftExp && leftScalar); // left is a single scalar
1954 bool RS=(!rightExp && rightScalar);
1955 bool LN=(!leftExp && !leftScalar); // left is a single non-scalar
1956 bool RN=(!rightExp && !rightScalar);
1957 bool LEN=(leftExp && !leftScalar); // left is an expanded non-scalar
1958 bool REN=(rightExp && !rightScalar);
1959
1960 if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1961 {
1962 chunksize=m_left->getNumDPPSample()*leftsize;
1963 leftstep=0;
1964 rightstep=0;
1965 numsteps=1;
1966 }
1967 else if (LES || RES)
1968 {
1969 chunksize=1;
1970 if (LES) // left is an expanded scalar
1971 {
1972 if (RS)
1973 {
1974 leftstep=1;
1975 rightstep=0;
1976 numsteps=m_left->getNumDPPSample();
1977 }
1978 else // RN or REN
1979 {
1980 leftstep=0;
1981 oleftstep=1;
1982 rightstep=1;
1983 orightstep=(RN ? -(int)rightsize : 0);
1984 numsteps=rightsize;
1985 onumsteps=m_left->getNumDPPSample();
1986 }
1987 }
1988 else // right is an expanded scalar
1989 {
1990 if (LS)
1991 {
1992 rightstep=1;
1993 leftstep=0;
1994 numsteps=m_right->getNumDPPSample();
1995 }
1996 else
1997 {
1998 rightstep=0;
1999 orightstep=1;
2000 leftstep=1;
2001 oleftstep=(LN ? -(int)leftsize : 0);
2002 numsteps=leftsize;
2003 onumsteps=m_right->getNumDPPSample();
2004 }
2005 }
2006 }
2007 else // this leaves (LEN, RS), (LEN, RN) and their transposes
2008 {
2009 if (LEN) // and Right will be a single value
2010 {
2011 chunksize=rightsize;
2012 leftstep=rightsize;
2013 rightstep=0;
2014 numsteps=m_left->getNumDPPSample();
2015 if (RS)
2016 {
2017 numsteps*=leftsize;
2018 }
2019 }
2020 else // REN
2021 {
2022 chunksize=leftsize;
2023 rightstep=leftsize;
2024 leftstep=0;
2025 numsteps=m_right->getNumDPPSample();
2026 if (LS)
2027 {
2028 numsteps*=rightsize;
2029 }
2030 }
2031 }
2032
2033 int resultStep=max(leftstep,rightstep); // only one (at most) should be !=0
2034 // Get the values of sub-expressions
2035 const CplxVectorType* left=m_left->resolveNodeSampleCplx(tid,sampleNo,lroffset);
2036 const CplxVectorType* right=m_right->resolveNodeSampleCplx(tid,sampleNo,rroffset);
2037 LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
2038 LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
2039 LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
2040 LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
2041 LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
2042 LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
2043 LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN << endl;)
2044
2045 LAZYDEBUG(cout << "Left res["<< lroffset<< "]=" << (*left)[lroffset] << endl;)
2046 LAZYDEBUG(cout << "Right res["<< rroffset<< "]=" << (*right)[rroffset] << endl;)
2047
2048
2049 roffset=m_samplesize*tid;
2050 cplx_t* resultp=&(m_samples_c[roffset]); // results are stored at the vector offset we received
2051 switch(m_op)
2052 {
2053 case ADD:
2054 //PROC_OP(NO_ARG,plus<double>());
2055 escript::binaryOpVectorLazyArithmeticHelper<cplx_t, cplx_t, cplx_t>(resultp,
2056 &(*left)[0],
2057 &(*right)[0],
2058 chunksize,
2059 onumsteps,
2060 numsteps,
2061 resultStep,
2062 leftstep,
2063 rightstep,
2064 oleftstep,
2065 orightstep,
2066 lroffset,
2067 rroffset,
2068 escript::ES_optype::ADD);
2069 break;
2070 case SUB:
2071 escript::binaryOpVectorLazyArithmeticHelper<cplx_t, cplx_t, cplx_t>(resultp,
2072 &(*left)[0],
2073 &(*right)[0],
2074 chunksize,
2075 onumsteps,
2076 numsteps,
2077 resultStep,
2078 leftstep,
2079 rightstep,
2080 oleftstep,
2081 orightstep,
2082 lroffset,
2083 rroffset,
2084 escript::ES_optype::SUB);
2085 //PROC_OP(NO_ARG,minus<double>());
2086 break;
2087 case MUL:
2088 //PROC_OP(NO_ARG,multiplies<double>());
2089 escript::binaryOpVectorLazyArithmeticHelper<cplx_t, cplx_t, cplx_t>(resultp,
2090 &(*left)[0],
2091 &(*right)[0],
2092 chunksize,
2093 onumsteps,
2094 numsteps,
2095 resultStep,
2096 leftstep,
2097 rightstep,
2098 oleftstep,
2099 orightstep,
2100 lroffset,
2101 rroffset,
2102 escript::ES_optype::MUL);
2103 break;
2104 case DIV:
2105 //PROC_OP(NO_ARG,divides<double>());
2106 escript::binaryOpVectorLazyArithmeticHelper<cplx_t, cplx_t, cplx_t>(resultp,
2107 &(*left)[0],
2108 &(*right)[0],
2109 chunksize,
2110 onumsteps,
2111 numsteps,
2112 resultStep,
2113 leftstep,
2114 rightstep,
2115 oleftstep,
2116 orightstep,
2117 lroffset,
2118 rroffset,
2119 escript::ES_optype::DIV);
2120 break;
2121 case POW:
2122 //PROC_OP(double (double,double),::pow);
2123 escript::binaryOpVectorLazyArithmeticHelper<cplx_t, cplx_t, cplx_t>(resultp,
2124 &(*left)[0],
2125 &(*right)[0],
2126 chunksize,
2127 onumsteps,
2128 numsteps,
2129 resultStep,
2130 leftstep,
2131 rightstep,
2132 oleftstep,
2133 orightstep,
2134 lroffset,
2135 rroffset,
2136 escript::ES_optype::POW);
2137 break;
2138 default:
2139 throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
2140 }
2141 LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples_c[roffset] << endl;)
2142 return &m_samples_c;
2143 }
2144
2145
2146 // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
2147 // have already been collapsed to IDENTITY. So we must have at least one expanded child.
2148 // unlike the other resolve helpers, we must treat these datapoints separately.
2149 const DataTypes::RealVectorType*
2150 DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset) const
2151 {
2152 LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
2153 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
2154 // first work out which of the children are expanded
2155 bool leftExp=(m_left->m_readytype=='E');
2156 bool rightExp=(m_right->m_readytype=='E');
2157 int steps=getNumDPPSample();
2158 int leftStep=(leftExp? m_left->getNoValues() : 0); // do not have scalars as input to this method
2159 int rightStep=(rightExp?m_right->getNoValues() : 0);
2160
2161 int resultStep=getNoValues();
2162 roffset=m_samplesize*tid;
2163 size_t offset=roffset;
2164 const RealVectorType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
2165 const RealVectorType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
2166 LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) << endl;
2167 cout << getNoValues() << endl;)
2168
2169
2170 LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
2171 LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
2172 LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
2173 LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
2174 LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
2175 LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
2176 LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
2177
2178 double* resultp=&(m_samples_r[offset]); // results are stored at the vector offset we received
2179 switch(m_op)
2180 {
2181 case PROD:
2182 for (int i=0;i<steps;++i,resultp+=resultStep)
2183 {
2184 const double *ptr_0 = &((*left)[lroffset]);
2185 const double *ptr_1 = &((*right)[rroffset]);
2186
2187 LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
2188 LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
2189
2190 matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
2191
2192 lroffset+=leftStep;
2193 rroffset+=rightStep;
2194 }
2195 break;
2196 default:
2197 throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
2198 }
2199 roffset=offset;
2200 return &m_samples_r;
2201 }
2202
2203 const DataTypes::CplxVectorType*
2204 DataLazy::resolveNodeTProdCplx(int tid, int sampleNo, size_t& roffset) const
2205 {
2206 LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
2207
2208 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
2209 // first work out which of the children are expanded
2210 bool leftExp=(m_left->m_readytype=='E');
2211 bool rightExp=(m_right->m_readytype=='E');
2212 int steps=getNumDPPSample();
2213 int leftStep=(leftExp? m_left->getNoValues() : 0); // do not have scalars as input to this method
2214 int rightStep=(rightExp?m_right->getNoValues() : 0);
2215
2216 int resultStep=getNoValues();
2217 roffset=m_samplesize*tid;
2218 size_t offset=roffset;
2219
2220 const CplxVectorType* left=m_left->resolveNodeSampleCplx(tid, sampleNo, lroffset);
2221
2222 const CplxVectorType* right=m_right->resolveNodeSampleCplx(tid, sampleNo, rroffset);
2223
2224 LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) << endl;
2225 cout << getNoValues() << endl;)
2226
2227
2228 LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
2229 LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
2230 LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
2231 LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
2232 LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
2233 LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
2234 LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
2235
2236 cplx_t* resultp=&(m_samples_c[offset]); // results are stored at the vector offset we received
2237 switch(m_op)
2238 {
2239 case PROD:
2240 for (int i=0;i<steps;++i,resultp+=resultStep)
2241 {
2242 const cplx_t *ptr_0 = &((*left)[lroffset]);
2243 const cplx_t *ptr_1 = &((*right)[rroffset]);
2244
2245 LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
2246 LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
2247
2248 matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
2249
2250 lroffset+=leftStep;
2251 rroffset+=rightStep;
2252 }
2253 break;
2254 default:
2255 throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
2256 }
2257 roffset=offset;
2258 return &m_samples_c;
2259 }
2260
2261 const DataTypes::RealVectorType*
2262 DataLazy::resolveSample(int sampleNo, size_t& roffset) const
2263 {
2264 #ifdef _OPENMP
2265 int tid=omp_get_thread_num();
2266 #else
2267 int tid=0;
2268 #endif
2269
2270 #ifdef LAZY_STACK_PROF
2271 stackstart[tid]=&tid;
2272 stackend[tid]=&tid;
2273 const DataTypes::RealVectorType* r=resolveNodeSample(tid, sampleNo, roffset);
2274 size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];
2275 #pragma omp critical
2276 if (d>maxstackuse)
2277 {
2278 cout << "Max resolve Stack use " << d << endl;
2279 maxstackuse=d;
2280 }
2281 return r;
2282 #else
2283 return resolveNodeSample(tid, sampleNo, roffset);
2284 #endif
2285 }
2286
2287
2288 const DataTypes::RealVectorType*
2289 DataLazy::resolveTypedSample(int sampleNo, size_t& roffset, DataTypes::real_t dummy) const
2290 {
2291 #ifdef _OPENMP
2292 int tid=omp_get_thread_num();
2293 #else
2294 int tid=0;
2295 #endif
2296
2297 #ifdef LAZY_STACK_PROF
2298 stackstart[tid]=&tid;
2299 stackend[tid]=&tid;
2300 const DataTypes::RealVectorType* r=resolveNodeSample(tid, sampleNo, roffset);
2301 size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];
2302 #pragma omp critical
2303 if (d>maxstackuse)
2304 {
2305 cout << "Max resolve Stack use " << d << endl;
2306 maxstackuse=d;
2307 }
2308 return r;
2309 #else
2310 return resolveNodeSample(tid, sampleNo, roffset);
2311 #endif
2312 }
2313
2314 const DataTypes::CplxVectorType*
2315 DataLazy::resolveTypedSample(int sampleNo, size_t& roffset, DataTypes::cplx_t dummy) const
2316 {
2317 #ifdef _OPENMP
2318 int tid=omp_get_thread_num();
2319 #else
2320 int tid=0;
2321 #endif
2322
2323 #ifdef LAZY_STACK_PROF
2324 stackstart[tid]=&tid;
2325 stackend[tid]=&tid;
2326 const DataTypes::CplxVectorType* r=resolveNodeSampleCplx(tid, sampleNo, roffset);
2327 size_t d=(size_t)stackstart[tid]-(size_t)stackend[tid];
2328 #pragma omp critical
2329 if (d>maxstackuse)
2330 {
2331 cout << "Max resolve Stack use " << d << endl;
2332 maxstackuse=d;
2333 }
2334 return r;
2335 #else
2336 return resolveNodeSampleCplx(tid, sampleNo, roffset);
2337 #endif
2338 }
2339
2340
2341
2342 // This needs to do the work of the identity constructor
2343 void
2344 DataLazy::resolveToIdentity()
2345 {
2346 if (m_op==IDENTITY)
2347 return;
2348 if (isComplex())
2349 {
2350 DataReady_ptr p=resolveNodeWorkerCplx();
2351 makeIdentity(p);
2352 }
2353 else
2354 {
2355 DataReady_ptr p=resolveNodeWorker();
2356 makeIdentity(p);
2357 }
2358 }
2359
2360 void DataLazy::makeIdentity(const DataReady_ptr& p)
2361 {
2362 m_op=IDENTITY;
2363 m_axis_offset=0;
2364 m_transpose=0;
2365 m_SL=m_SM=m_SR=0;
2366 m_children=m_height=0;
2367 m_id=p;
2368 if(p->isConstant()) {m_readytype='C';}
2369 else if(p->isExpanded()) {m_readytype='E';}
2370 else if (p->isTagged()) {m_readytype='T';}
2371 else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
2372 m_samplesize=p->getNumDPPSample()*p->getNoValues();
2373 m_left.reset();
2374 m_right.reset();
2375 m_iscompl=p->isComplex();
2376 }
2377
2378
2379 DataReady_ptr
2380 DataLazy::resolve()
2381 {
2382 resolveToIdentity();
2383 return m_id;
2384 }
2385
2386
2387 /* This is really a static method but I think that caused problems in windows */
2388 void
2389 DataLazy::resolveGroupWorker(std::vector<DataLazy*>& dats)
2390 {
2391 if (dats.empty())
2392 {
2393 return;
2394 }
2395 vector<DataLazy*> work;
2396 FunctionSpace fs=dats[0]->getFunctionSpace();
2397 bool match=true;
2398 for (int i=dats.size()-1;i>=0;--i)
2399 {
2400 if (dats[i]->m_readytype!='E')
2401 {
2402 dats[i]->collapse();
2403 }
2404 if (dats[i]->m_op!=IDENTITY)
2405 {
2406 work.push_back(dats[i]);
2407 if (fs!=dats[i]->getFunctionSpace())
2408 {
2409 match=false;
2410 }
2411 }
2412 }
2413 if (work.empty())
2414 {
2415 return; // no work to do
2416 }
2417 if (match) // all functionspaces match. Yes I realise this is overly strict
2418 { // it is possible that dats[0] is one of the objects which we discarded and
2419 // all the other functionspaces match.
2420 vector<DataExpanded*> dep;
2421 vector<RealVectorType*> vecs;
2422 for (int i=0;i<work.size();++i)
2423 {
2424 dep.push_back(new DataExpanded(fs,work[i]->getShape(), RealVectorType(work[i]->getNoValues())));
2425 vecs.push_back(&(dep[i]->getVectorRW()));
2426 }
2427 int totalsamples=work[0]->getNumSamples();
2428 const RealVectorType* res=0; // Storage for answer
2429 int sample;
2430 #pragma omp parallel private(sample, res)
2431 {
2432 size_t roffset=0;
2433 #pragma omp for schedule(static)
2434 for (sample=0;sample<totalsamples;++sample)
2435 {
2436 roffset=0;
2437 int j;
2438 for (j=work.size()-1;j>=0;--j)
2439 {
2440 #ifdef _OPENMP
2441 res=work[j]->resolveNodeSample(omp_get_thread_num(),sample,roffset);
2442 #else
2443 res=work[j]->resolveNodeSample(0,sample,roffset);
2444 #endif
2445 RealVectorType::size_type outoffset=dep[j]->getPointOffset(sample,0);
2446 memcpy(&((*vecs[j])[outoffset]),&((*res)[roffset]),work[j]->m_samplesize*sizeof(RealVectorType::ElementType));
2447 }
2448 }
2449 }
2450 // Now we need to load the new results as identity ops into the lazy nodes
2451 for (int i=work.size()-1;i>=0;--i)
2452 {
2453 work[i]->makeIdentity(REFCOUNTNS::dynamic_pointer_cast<DataReady>(dep[i]->getPtr()));
2454 }
2455 }
2456 else // functionspaces do not match
2457 {
2458 for (int i=0;i<work.size();++i)
2459 {
2460 work[i]->resolveToIdentity();
2461 }
2462 }
2463 }
2464
2465
2466
2467 // This version of resolve uses storage in each node to hold results
2468 DataReady_ptr
2469 DataLazy::resolveNodeWorker()
2470 {
2471 if (m_readytype!='E') // if the whole sub-expression is Constant or Tagged, then evaluate it normally
2472 {
2473 collapse();
2474 }
2475 if (m_op==IDENTITY) // So a lazy expression of Constant or Tagged data will be returned here.
2476 {
2477 return m_id;
2478 }
2479 // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
2480 DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(), RealVectorType(getNoValues()));
2481 RealVectorType& resvec=result->getVectorRW();
2482 DataReady_ptr resptr=DataReady_ptr(result);
2483
2484 int sample;
2485 int totalsamples=getNumSamples();
2486 const RealVectorType* res=0; // Storage for answer
2487 LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
2488 #pragma omp parallel private(sample,res)
2489 {
2490 size_t roffset=0;
2491 #ifdef LAZY_STACK_PROF
2492 stackstart[omp_get_thread_num()]=&roffset;
2493 stackend[omp_get_thread_num()]=&roffset;
2494 #endif
2495 #pragma omp for schedule(static)
2496 for (sample=0;sample<totalsamples;++sample)
2497 {
2498 roffset=0;
2499 #ifdef _OPENMP
2500 res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
2501 #else
2502 res=resolveNodeSample(0,sample,roffset);
2503 #endif
2504 LAZYDEBUG(cout << "Sample #" << sample << endl;)
2505 LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
2506 RealVectorType::size_type outoffset=result->getPointOffset(sample,0);
2507 memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(RealVectorType::ElementType));
2508 }
2509 }
2510 #ifdef LAZY_STACK_PROF
2511 for (int i=0;i<getNumberOfThreads();++i)
2512 {
2513 size_t r=((size_t)stackstart[i] - (size_t)stackend[i]);
2514 // cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " << r << endl;
2515 if (r>maxstackuse)
2516 {
2517 maxstackuse=r;
2518 }
2519 }
2520 cout << "Max resolve Stack use=" << maxstackuse << endl;
2521 #endif
2522 return resptr;
2523 }
2524
2525 // This version should only be called on complex lazy nodes
2526 DataReady_ptr
2527 DataLazy::resolveNodeWorkerCplx()
2528 {
2529 if (m_readytype!='E') // if the whole sub-expression is Constant or Tagged, then evaluate it normally
2530 {
2531 collapse();
2532 }
2533 if (m_op==IDENTITY) // So a lazy expression of Constant or Tagged data will be returned here.
2534 {
2535 return m_id;
2536 }
2537 // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
2538 DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(), CplxVectorType(getNoValues()));
2539 CplxVectorType& resvec=result->getVectorRWC();
2540 DataReady_ptr resptr=DataReady_ptr(result);
2541
2542 int sample;
2543 int totalsamples=getNumSamples();
2544 const CplxVectorType* res=0; // Storage for answer
2545 LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
2546 #pragma omp parallel private(sample,res)
2547 {
2548 size_t roffset=0;
2549 #ifdef LAZY_STACK_PROF
2550 stackstart[omp_get_thread_num()]=&roffset;
2551 stackend[omp_get_thread_num()]=&roffset;
2552 #endif
2553 #pragma omp for schedule(static)
2554 for (sample=0;sample<totalsamples;++sample)
2555 {
2556 roffset=0;
2557 #ifdef _OPENMP
2558 res=resolveNodeSampleCplx(omp_get_thread_num(),sample,roffset);
2559 #else
2560 res=resolveNodeSampleCplx(0,sample,roffset);
2561 #endif
2562 LAZYDEBUG(cout << "Sample #" << sample << endl;)
2563 LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
2564 CplxVectorType::size_type outoffset=result->getPointOffset(sample,0);
2565 memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(CplxVectorType::ElementType));
2566 }
2567 }
2568 #ifdef LAZY_STACK_PROF
2569 for (int i=0;i<getNumberOfThreads();++i)
2570 {
2571 size_t r=((size_t)stackstart[i] - (size_t)stackend[i]);
2572 // cout << i << " " << stackstart[i] << " .. " << stackend[i] << " = " << r << endl;
2573 if (r>maxstackuse)
2574 {
2575 maxstackuse=r;
2576 }
2577 }
2578 cout << "Max resolve Stack use=" << maxstackuse << endl;
2579 #endif
2580 return resptr;
2581 }
2582
2583
2584 std::string
2585 DataLazy::toString() const
2586 {
2587 ostringstream oss;
2588 oss << "Lazy Data: [depth=" << m_height<< "] ";
2589 switch (escriptParams.getLazyStrFmt())
2590 {
2591 case 1: // tree format
2592 oss << endl;
2593 intoTreeString(oss,"");
2594 break;
2595 case 2: // just the depth
2596 break;
2597 default:
2598 intoString(oss);
2599 break;
2600 }
2601 return oss.str();
2602 }
2603
2604
2605 void
2606 DataLazy::intoString(ostringstream& oss) const
2607 {
2608 // oss << "[" << m_children <<";"<<m_height <<"]";
2609 switch (m_opgroup)
2610 {
2611 case G_IDENTITY:
2612 if (m_id->isExpanded())
2613 {
2614 oss << "E";
2615 }
2616 else if (m_id->isTagged())
2617 {
2618 oss << "T";
2619 }
2620 else if (m_id->isConstant())
2621 {
2622 oss << "C";
2623 }
2624 else
2625 {
2626 oss << "?";
2627 }
2628 if (m_id->isComplex())
2629 {
2630 oss << "j";
2631 }
2632 oss << '@' << m_id.get();
2633 break;
2634 case G_BINARY:
2635 oss << '(';
2636 m_left->intoString(oss);
2637 oss << ' ' << opToString(m_op) << ' ';
2638 m_right->intoString(oss);
2639 oss << ')';
2640 if (isComplex())
2641 {
2642 oss << "j";
2643 }
2644 break;
2645 case G_UNARY:
2646 case G_UNARY_P:
2647 case G_UNARY_PR:
2648 case G_UNARY_R:
2649 case G_NP1OUT:
2650 case G_NP1OUT_P:
2651 case G_REDUCTION:
2652 oss << opToString(m_op) << '(';
2653 m_left->intoString(oss);
2654 oss << ')';
2655 if (isComplex())
2656 {
2657 oss << "j";
2658 }
2659 break;
2660 case G_TENSORPROD:
2661 oss << opToString(m_op) << '(';
2662 m_left->intoString(oss);
2663 oss << ", ";
2664 m_right->intoString(oss);
2665 oss << ')';
2666 if (isComplex())
2667 {
2668 oss << "j";
2669 }
2670 break;
2671 case G_NP1OUT_2P:
2672 oss << opToString(m_op) << '(';
2673 m_left->intoString(oss);
2674 oss << ", " << m_axis_offset << ", " << m_transpose;
2675 oss << ')';
2676 if (isComplex())
2677 {
2678 oss << "j";
2679 }
2680 break;
2681 case G_CONDEVAL:
2682 oss << opToString(m_op)<< '(' ;
2683 m_mask->intoString(oss);
2684 oss << " ? ";
2685 m_left->intoString(oss);
2686 oss << " : ";
2687 m_right->intoString(oss);
2688 oss << ')';
2689 if (isComplex())
2690 {
2691 oss << "j";
2692 }
2693 break;
2694 case G_UNARY_C:
2695 oss << opToString(m_op) <<'(';
2696 m_left->intoString(oss);
2697 oss << ')';
2698 if (isComplex())
2699 {
2700 oss << "j";
2701 }
2702 break;
2703 default:
2704 oss << "UNKNOWN";
2705 if (isComplex())
2706 {
2707 oss << "j";
2708 }
2709 }
2710 }
2711
2712
2713 void
2714 DataLazy::intoTreeString(ostringstream& oss, string indent) const
2715 {
2716 oss << '[' << m_rank << ':' << setw(3) << m_samplesize << "] " << indent;
2717 ES_opgroup gop=m_opgroup;
2718 switch (gop)
2719 {
2720 case G_IDENTITY:
2721 if (m_id->isExpanded())
2722 {
2723 oss << "E";
2724 }
2725 else if (m_id->isTagged())
2726 {
2727 oss << "T";
2728 }
2729 else if (m_id->isConstant())
2730 {
2731 oss << "C";
2732 }
2733 else
2734 {
2735 oss << "?";
2736 }
2737 if (isComplex())
2738 {
2739 oss << "j";
2740 }
2741 oss << '@' << m_id.get() << endl;
2742 break;
2743 case G_BINARY:
2744 oss << opToString(m_op);
2745 if (isComplex())
2746 {
2747 oss << 'j';
2748 }
2749 oss << endl;
2750 indent+='.';
2751 m_left->intoTreeString(oss, indent);
2752 m_right->intoTreeString(oss, indent);
2753 break;
2754 case G_UNARY:
2755 case G_UNARY_P:
2756 case G_UNARY_PR:
2757 case G_UNARY_R:
2758 case G_NP1OUT:
2759 case G_NP1OUT_P:
2760 case G_REDUCTION:
2761 case G_UNARY_C:
2762 oss << opToString(m_op);
2763 if (isComplex())
2764 {
2765 oss << 'j';
2766 }
2767 oss << endl;
2768 indent+='.';
2769 m_left->intoTreeString(oss, indent);
2770 break;
2771 case G_TENSORPROD:
2772 oss << opToString(m_op);
2773 if (isComplex())
2774 {
2775 oss << 'j';
2776 }
2777 oss << endl;
2778 indent+='.';
2779 m_left->intoTreeString(oss, indent);
2780 m_right->intoTreeString(oss, indent);
2781 break;
2782 case G_NP1OUT_2P:
2783 oss << opToString(m_op);
2784 if (isComplex())
2785 {
2786 oss << 'j';
2787 }
2788 oss << ", " << m_axis_offset << ", " << m_transpose<< endl;
2789 indent+='.';
2790 m_left->intoTreeString(oss, indent);
2791 break;
2792 default:
2793 oss << "UNKNOWN";
2794 if (isComplex())
2795 {
2796 oss << 'j';
2797 }
2798 oss << endl;
2799 }
2800 }
2801
2802
2803 DataAbstract*
2804 DataLazy::deepCopy() const
2805 {
2806 switch (m_opgroup)
2807 {
2808 case G_IDENTITY: return new DataLazy(m_id->deepCopy()->getPtr());
2809 case G_UNARY:
2810 case G_UNARY_R:
2811 case G_UNARY_C:
2812 case G_REDUCTION: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
2813 case G_UNARY_PR:
2814 case G_UNARY_P: return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_tol);
2815 case G_BINARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
2816 case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
2817 case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
2818 case G_NP1OUT_P: return new DataLazy(m_left->deepCopy()->getPtr(),m_op, m_axis_offset);
2819 case G_NP1OUT_2P: return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
2820 case G_CONDEVAL: return new DataLazy(m_mask->deepCopy()->getPtr(), m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr());
2821 default:
2822 throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
2823 }
2824 }
2825
2826 // For this, we don't care what op you were doing because the answer is now zero
2827 DataAbstract*
2828 DataLazy::zeroedCopy() const
2829 {
2830 return new DataLazy(m_id->zeroedCopy()->getPtr());
2831 }
2832
2833 // There is no single, natural interpretation of getLength on DataLazy.
2834 // Instances of DataReady can look at the size of their vectors.
2835 // For lazy though, it could be the size the data would be if it were resolved;
2836 // or it could be some function of the lengths of the DataReady instances which
2837 // form part of the expression.
2838 // Rather than have people making assumptions, I have disabled the method.
2839 DataTypes::RealVectorType::size_type
2840 DataLazy::getLength() const
2841 {
2842 throw DataException("getLength() does not make sense for lazy data.");
2843 }
2844
2845
2846 DataAbstract*
2847 DataLazy::getSlice(const DataTypes::RegionType& region) const
2848 {
2849 throw DataException("getSlice - not implemented for Lazy objects.");
2850 }
2851
2852
2853 // To do this we need to rely on our child nodes
2854 DataTypes::RealVectorType::size_type
2855 DataLazy::getPointOffset(int sampleNo,
2856 int dataPointNo)
2857 {
2858 if (m_op==IDENTITY)
2859 {
2860 return m_id->getPointOffset(sampleNo,dataPointNo);
2861 }
2862 if (m_readytype!='E')
2863 {
2864 collapse();
2865 return m_id->getPointOffset(sampleNo,dataPointNo);
2866 }
2867 // at this point we do not have an identity node and the expression will be Expanded
2868 // so we only need to know which child to ask
2869 if (m_left->m_readytype=='E')
2870 {
2871 return m_left->getPointOffset(sampleNo,dataPointNo);
2872 }
2873 else
2874 {
2875 return m_right->getPointOffset(sampleNo,dataPointNo);
2876 }
2877 }
2878
2879 // To do this we need to rely on our child nodes
2880 DataTypes::RealVectorType::size_type
2881 DataLazy::getPointOffset(int sampleNo,
2882 int dataPointNo) const
2883 {
2884 if (m_op==IDENTITY)
2885 {
2886 return m_id->getPointOffset(sampleNo,dataPointNo);
2887 }
2888 if (m_readytype=='E')
2889 {
2890 // at this point we do not have an identity node and the expression will be Expanded
2891 // so we only need to know which child to ask
2892 if (m_left->m_readytype=='E')
2893 {
2894 return m_left->getPointOffset(sampleNo,dataPointNo);
2895 }
2896 else
2897 {
2898 return m_right->getPointOffset(sampleNo,dataPointNo);
2899 }
2900 }
2901 if (m_readytype=='C')
2902 {
2903 return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
2904 }
2905 throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
2906 }
2907
2908
2909 // I have decided to let Data:: handle this issue.
2910 void
2911 DataLazy::setToZero()
2912 {
2913 // DataTypes::RealVectorType v(getNoValues(),0);
2914 // m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
2915 // m_op=IDENTITY;
2916 // m_right.reset();
2917 // m_left.reset();
2918 // m_readytype='C';
2919 // m_buffsRequired=1;
2920
2921 (void)privdebug; // to stop the compiler complaining about unused privdebug
2922 throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
2923 }
2924
2925 bool
2926 DataLazy::actsExpanded() const
2927 {
2928 return (m_readytype=='E');
2929 }
2930
2931 } // end namespace
2932

  ViewVC Help
Powered by ViewVC 1.1.26