/[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 2770 - (show annotations)
Wed Nov 25 01:24:51 2009 UTC (9 years, 4 months ago) by jfenwick
Original Path: trunk/escript/src/DataLazy.cpp
File size: 54086 byte(s)
Removed buffer implementation of Lazy.
Removed supporting Alloc/Free Sample buffer calls.
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 #include "DataLazy.h"
16 #ifdef USE_NETCDF
17 #include <netcdfcpp.h>
18 #endif
19 #ifdef PASO_MPI
20 #include <mpi.h>
21 #endif
22 #ifdef _OPENMP
23 #include <omp.h>
24 #endif
25 #include "FunctionSpace.h"
26 #include "DataTypes.h"
27 #include "Data.h"
28 #include "UnaryFuncs.h" // for escript::fsign
29 #include "Utils.h"
30
31 #include "EscriptParams.h"
32
33 #include <iomanip> // for some fancy formatting in debug
34
35 // #define LAZYDEBUG(X) if (privdebug){X;}
36 #define LAZYDEBUG(X)
37 namespace
38 {
39 bool privdebug=false;
40
41 #define ENABLEDEBUG privdebug=true;
42 #define DISABLEDEBUG privdebug=false;
43 }
44
45 // #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl;resolveToIdentity();}
46
47 #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {resolveToIdentity();}
48
49
50 /*
51 How does DataLazy work?
52 ~~~~~~~~~~~~~~~~~~~~~~~
53
54 Each instance represents a single operation on one or two other DataLazy instances. These arguments are normally
55 denoted left and right.
56
57 A special operation, IDENTITY, stores an instance of DataReady in the m_id member.
58 This means that all "internal" nodes in the structure are instances of DataLazy.
59
60 Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
61 Note that IDENITY is not considered a unary operation.
62
63 I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
64 It must however form a DAG (directed acyclic graph).
65 I will refer to individual DataLazy objects with the structure as nodes.
66
67 Each node also stores:
68 - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
69 evaluated.
70 - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to
71 evaluate the expression.
72 - m_samplesize ~ the number of doubles stored in a sample.
73
74 When a new node is created, the above values are computed based on the values in the child nodes.
75 Eg: if left requires 4 samples and right requires 6 then left+right requires 7 samples.
76
77 The resolve method, which produces a DataReady from a DataLazy, does the following:
78 1) Create a DataReady to hold the new result.
79 2) Allocate a vector (v) big enough to hold m_buffsrequired samples.
80 3) For each sample, call resolveSample with v, to get its values and copy them into the result object.
81
82 (In the case of OMP, multiple samples are resolved in parallel so the vector needs to be larger.)
83
84 resolveSample returns a Vector* and an offset within that vector where the result is stored.
85 Normally, this would be v, but for identity nodes their internal vector is returned instead.
86
87 The convention that I use, is that the resolve methods should store their results starting at the offset they are passed.
88
89 For expressions which evaluate to Constant or Tagged, there is a different evaluation method.
90 The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.
91
92 To add a new operator you need to do the following (plus anything I might have forgotten - adding a new group for example):
93 1) Add to the ES_optype.
94 2) determine what opgroup your operation belongs to (X)
95 3) add a string for the op to the end of ES_opstrings
96 4) increase ES_opcount
97 5) add an entry (X) to opgroups
98 6) add an entry to the switch in collapseToReady
99 7) add an entry to resolveX
100 */
101
102
103 using namespace std;
104 using namespace boost;
105
106 namespace escript
107 {
108
109 namespace
110 {
111
112 enum ES_opgroup
113 {
114 G_UNKNOWN,
115 G_IDENTITY,
116 G_BINARY, // pointwise operations with two arguments
117 G_UNARY, // pointwise operations with one argument
118 G_UNARY_P, // pointwise operations with one argument, requiring a parameter
119 G_NP1OUT, // non-pointwise op with one output
120 G_NP1OUT_P, // non-pointwise op with one output requiring a parameter
121 G_TENSORPROD, // general tensor product
122 G_NP1OUT_2P, // non-pointwise op with one output requiring two params
123 G_REDUCTION // non-pointwise unary op with a scalar output
124 };
125
126
127
128
129 string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
130 "sin","cos","tan",
131 "asin","acos","atan","sinh","cosh","tanh","erf",
132 "asinh","acosh","atanh",
133 "log10","log","sign","abs","neg","pos","exp","sqrt",
134 "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
135 "symmetric","nonsymmetric",
136 "prod",
137 "transpose", "trace",
138 "swapaxes",
139 "minval", "maxval"};
140 int ES_opcount=43;
141 ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
142 G_UNARY,G_UNARY,G_UNARY, //10
143 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 17
144 G_UNARY,G_UNARY,G_UNARY, // 20
145 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 28
146 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P, // 35
147 G_NP1OUT,G_NP1OUT,
148 G_TENSORPROD,
149 G_NP1OUT_P, G_NP1OUT_P,
150 G_NP1OUT_2P,
151 G_REDUCTION, G_REDUCTION};
152 inline
153 ES_opgroup
154 getOpgroup(ES_optype op)
155 {
156 return opgroups[op];
157 }
158
159 // return the FunctionSpace of the result of "left op right"
160 FunctionSpace
161 resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
162 {
163 // perhaps this should call interpolate and throw or something?
164 // maybe we need an interpolate node -
165 // that way, if interpolate is required in any other op we can just throw a
166 // programming error exception.
167
168 FunctionSpace l=left->getFunctionSpace();
169 FunctionSpace r=right->getFunctionSpace();
170 if (l!=r)
171 {
172 if (r.probeInterpolation(l))
173 {
174 return l;
175 }
176 if (l.probeInterpolation(r))
177 {
178 return r;
179 }
180 throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
181 }
182 return l;
183 }
184
185 // return the shape of the result of "left op right"
186 // the shapes resulting from tensor product are more complex to compute so are worked out elsewhere
187 DataTypes::ShapeType
188 resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
189 {
190 if (left->getShape()!=right->getShape())
191 {
192 if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
193 {
194 throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
195 }
196
197 if (left->getRank()==0) // we need to allow scalar * anything
198 {
199 return right->getShape();
200 }
201 if (right->getRank()==0)
202 {
203 return left->getShape();
204 }
205 throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
206 }
207 return left->getShape();
208 }
209
210 // return the shape for "op left"
211
212 DataTypes::ShapeType
213 resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
214 {
215 switch(op)
216 {
217 case TRANS:
218 { // for the scoping of variables
219 const DataTypes::ShapeType& s=left->getShape();
220 DataTypes::ShapeType sh;
221 int rank=left->getRank();
222 if (axis_offset<0 || axis_offset>rank)
223 {
224 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
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 throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
310 }
311 if (axis1<0 || axis1>rank-1) {
312 throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
313 }
314 if (axis0 == axis1) {
315 throw DataException("Error - Data::swapaxes: axis indices must be different.");
316 }
317 if (axis0 > axis1) {
318 axis0_tmp=axis1;
319 axis1_tmp=axis0;
320 } else {
321 axis0_tmp=axis0;
322 axis1_tmp=axis1;
323 }
324 for (int i=0; i<rank; i++) {
325 if (i == axis0_tmp) {
326 out_shape.push_back(s[axis1_tmp]);
327 } else if (i == axis1_tmp) {
328 out_shape.push_back(s[axis0_tmp]);
329 } else {
330 out_shape.push_back(s[i]);
331 }
332 }
333 return out_shape;
334 }
335
336
337 // determine the output shape for the general tensor product operation
338 // the additional parameters return information required later for the product
339 // the majority of this code is copy pasted from C_General_Tensor_Product
340 DataTypes::ShapeType
341 GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
342 {
343
344 // Get rank and shape of inputs
345 int rank0 = left->getRank();
346 int rank1 = right->getRank();
347 const DataTypes::ShapeType& shape0 = left->getShape();
348 const DataTypes::ShapeType& shape1 = right->getShape();
349
350 // Prepare for the loops of the product and verify compatibility of shapes
351 int start0=0, start1=0;
352 if (transpose == 0) {}
353 else if (transpose == 1) { start0 = axis_offset; }
354 else if (transpose == 2) { start1 = rank1-axis_offset; }
355 else { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
356
357 if (rank0<axis_offset)
358 {
359 throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
360 }
361
362 // Adjust the shapes for transpose
363 DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
364 DataTypes::ShapeType tmpShape1(rank1); // than using push_back
365 for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
366 for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
367
368 // Prepare for the loops of the product
369 SL=1, SM=1, SR=1;
370 for (int i=0; i<rank0-axis_offset; i++) {
371 SL *= tmpShape0[i];
372 }
373 for (int i=rank0-axis_offset; i<rank0; i++) {
374 if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
375 throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
376 }
377 SM *= tmpShape0[i];
378 }
379 for (int i=axis_offset; i<rank1; i++) {
380 SR *= tmpShape1[i];
381 }
382
383 // Define the shape of the output (rank of shape is the sum of the loop ranges below)
384 DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
385 { // block to limit the scope of out_index
386 int out_index=0;
387 for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
388 for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
389 }
390
391 if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
392 {
393 ostringstream os;
394 os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
395 throw DataException(os.str());
396 }
397
398 return shape2;
399 }
400
401 } // end anonymous namespace
402
403
404
405 // Return a string representing the operation
406 const std::string&
407 opToString(ES_optype op)
408 {
409 if (op<0 || op>=ES_opcount)
410 {
411 op=UNKNOWNOP;
412 }
413 return ES_opstrings[op];
414 }
415
416 void DataLazy::LazyNodeSetup()
417 {
418 #ifdef _OPENMP
419 int numthreads=omp_get_max_threads();
420 m_samples.resize(numthreads*m_samplesize);
421 m_sampleids=new int[numthreads];
422 for (int i=0;i<numthreads;++i)
423 {
424 m_sampleids[i]=-1;
425 }
426 #else
427 m_samples.resize(m_samplesize);
428 m_sampleids=new int[1];
429 m_sampleids[0]=-1;
430 #endif // _OPENMP
431 }
432
433
434 // Creates an identity node
435 DataLazy::DataLazy(DataAbstract_ptr p)
436 : parent(p->getFunctionSpace(),p->getShape())
437 ,m_sampleids(0),
438 m_samples(1)
439 {
440 if (p->isLazy())
441 {
442 // I don't want identity of Lazy.
443 // Question: Why would that be so bad?
444 // Answer: We assume that the child of ID is something we can call getVector on
445 throw DataException("Programmer error - attempt to create identity from a DataLazy.");
446 }
447 else
448 {
449 p->makeLazyShared();
450 DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
451 makeIdentity(dr);
452 LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
453 }
454 LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
455 }
456
457 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
458 : parent(left->getFunctionSpace(),(getOpgroup(op)!=G_REDUCTION)?left->getShape():DataTypes::scalarShape),
459 m_op(op),
460 m_axis_offset(0),
461 m_transpose(0),
462 m_SL(0), m_SM(0), m_SR(0)
463 {
464 if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT) && (getOpgroup(op)!=G_REDUCTION))
465 {
466 throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
467 }
468
469 DataLazy_ptr lleft;
470 if (!left->isLazy())
471 {
472 lleft=DataLazy_ptr(new DataLazy(left));
473 }
474 else
475 {
476 lleft=dynamic_pointer_cast<DataLazy>(left);
477 }
478 m_readytype=lleft->m_readytype;
479 m_left=lleft;
480 m_samplesize=getNumDPPSample()*getNoValues();
481 m_children=m_left->m_children+1;
482 m_height=m_left->m_height+1;
483 LazyNodeSetup();
484 SIZELIMIT
485 }
486
487
488 // In this constructor we need to consider interpolation
489 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
490 : parent(resultFS(left,right,op), resultShape(left,right,op)),
491 m_op(op),
492 m_SL(0), m_SM(0), m_SR(0)
493 {
494 LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
495 if ((getOpgroup(op)!=G_BINARY))
496 {
497 throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
498 }
499
500 if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
501 {
502 FunctionSpace fs=getFunctionSpace();
503 Data ltemp(left);
504 Data tmp(ltemp,fs);
505 left=tmp.borrowDataPtr();
506 }
507 if (getFunctionSpace()!=right->getFunctionSpace()) // right needs to be interpolated
508 {
509 Data tmp(Data(right),getFunctionSpace());
510 right=tmp.borrowDataPtr();
511 LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
512 }
513 left->operandCheck(*right);
514
515 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
516 {
517 m_left=dynamic_pointer_cast<DataLazy>(left);
518 LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
519 }
520 else
521 {
522 m_left=DataLazy_ptr(new DataLazy(left));
523 LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
524 }
525 if (right->isLazy())
526 {
527 m_right=dynamic_pointer_cast<DataLazy>(right);
528 LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
529 }
530 else
531 {
532 m_right=DataLazy_ptr(new DataLazy(right));
533 LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
534 }
535 char lt=m_left->m_readytype;
536 char rt=m_right->m_readytype;
537 if (lt=='E' || rt=='E')
538 {
539 m_readytype='E';
540 }
541 else if (lt=='T' || rt=='T')
542 {
543 m_readytype='T';
544 }
545 else
546 {
547 m_readytype='C';
548 }
549 m_samplesize=getNumDPPSample()*getNoValues();
550 m_children=m_left->m_children+m_right->m_children+2;
551 m_height=max(m_left->m_height,m_right->m_height)+1;
552 LazyNodeSetup();
553 SIZELIMIT
554 LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
555 }
556
557 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
558 : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
559 m_op(op),
560 m_axis_offset(axis_offset),
561 m_transpose(transpose)
562 {
563 if ((getOpgroup(op)!=G_TENSORPROD))
564 {
565 throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
566 }
567 if ((transpose>2) || (transpose<0))
568 {
569 throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
570 }
571 if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
572 {
573 FunctionSpace fs=getFunctionSpace();
574 Data ltemp(left);
575 Data tmp(ltemp,fs);
576 left=tmp.borrowDataPtr();
577 }
578 if (getFunctionSpace()!=right->getFunctionSpace()) // right needs to be interpolated
579 {
580 Data tmp(Data(right),getFunctionSpace());
581 right=tmp.borrowDataPtr();
582 }
583 // left->operandCheck(*right);
584
585 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
586 {
587 m_left=dynamic_pointer_cast<DataLazy>(left);
588 }
589 else
590 {
591 m_left=DataLazy_ptr(new DataLazy(left));
592 }
593 if (right->isLazy())
594 {
595 m_right=dynamic_pointer_cast<DataLazy>(right);
596 }
597 else
598 {
599 m_right=DataLazy_ptr(new DataLazy(right));
600 }
601 char lt=m_left->m_readytype;
602 char rt=m_right->m_readytype;
603 if (lt=='E' || rt=='E')
604 {
605 m_readytype='E';
606 }
607 else if (lt=='T' || rt=='T')
608 {
609 m_readytype='T';
610 }
611 else
612 {
613 m_readytype='C';
614 }
615 m_samplesize=getNumDPPSample()*getNoValues();
616 m_children=m_left->m_children+m_right->m_children+2;
617 m_height=max(m_left->m_height,m_right->m_height)+1;
618 LazyNodeSetup();
619 SIZELIMIT
620 LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
621 }
622
623
624 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
625 : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
626 m_op(op),
627 m_axis_offset(axis_offset),
628 m_transpose(0),
629 m_tol(0)
630 {
631 if ((getOpgroup(op)!=G_NP1OUT_P))
632 {
633 throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
634 }
635 DataLazy_ptr lleft;
636 if (!left->isLazy())
637 {
638 lleft=DataLazy_ptr(new DataLazy(left));
639 }
640 else
641 {
642 lleft=dynamic_pointer_cast<DataLazy>(left);
643 }
644 m_readytype=lleft->m_readytype;
645 m_left=lleft;
646 m_samplesize=getNumDPPSample()*getNoValues();
647 m_children=m_left->m_children+1;
648 m_height=m_left->m_height+1;
649 LazyNodeSetup();
650 SIZELIMIT
651 LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
652 }
653
654 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
655 : parent(left->getFunctionSpace(), left->getShape()),
656 m_op(op),
657 m_axis_offset(0),
658 m_transpose(0),
659 m_tol(tol)
660 {
661 if ((getOpgroup(op)!=G_UNARY_P))
662 {
663 throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
664 }
665 DataLazy_ptr lleft;
666 if (!left->isLazy())
667 {
668 lleft=DataLazy_ptr(new DataLazy(left));
669 }
670 else
671 {
672 lleft=dynamic_pointer_cast<DataLazy>(left);
673 }
674 m_readytype=lleft->m_readytype;
675 m_left=lleft;
676 m_samplesize=getNumDPPSample()*getNoValues();
677 m_children=m_left->m_children+1;
678 m_height=m_left->m_height+1;
679 LazyNodeSetup();
680 SIZELIMIT
681 LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
682 }
683
684
685 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, const int axis0, const int axis1)
686 : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
687 m_op(op),
688 m_axis_offset(axis0),
689 m_transpose(axis1),
690 m_tol(0)
691 {
692 if ((getOpgroup(op)!=G_NP1OUT_2P))
693 {
694 throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require two integer parameters.");
695 }
696 DataLazy_ptr lleft;
697 if (!left->isLazy())
698 {
699 lleft=DataLazy_ptr(new DataLazy(left));
700 }
701 else
702 {
703 lleft=dynamic_pointer_cast<DataLazy>(left);
704 }
705 m_readytype=lleft->m_readytype;
706 m_left=lleft;
707 m_samplesize=getNumDPPSample()*getNoValues();
708 m_children=m_left->m_children+1;
709 m_height=m_left->m_height+1;
710 LazyNodeSetup();
711 SIZELIMIT
712 LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)
713 }
714
715 DataLazy::~DataLazy()
716 {
717 delete[] m_sampleids;
718 }
719
720
721 /*
722 \brief Evaluates the expression using methods on Data.
723 This does the work for the collapse method.
724 For reasons of efficiency do not call this method on DataExpanded nodes.
725 */
726 DataReady_ptr
727 DataLazy::collapseToReady()
728 {
729 if (m_readytype=='E')
730 { // this is more an efficiency concern than anything else
731 throw DataException("Programmer Error - do not use collapse on Expanded data.");
732 }
733 if (m_op==IDENTITY)
734 {
735 return m_id;
736 }
737 DataReady_ptr pleft=m_left->collapseToReady();
738 Data left(pleft);
739 Data right;
740 if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
741 {
742 right=Data(m_right->collapseToReady());
743 }
744 Data result;
745 switch(m_op)
746 {
747 case ADD:
748 result=left+right;
749 break;
750 case SUB:
751 result=left-right;
752 break;
753 case MUL:
754 result=left*right;
755 break;
756 case DIV:
757 result=left/right;
758 break;
759 case SIN:
760 result=left.sin();
761 break;
762 case COS:
763 result=left.cos();
764 break;
765 case TAN:
766 result=left.tan();
767 break;
768 case ASIN:
769 result=left.asin();
770 break;
771 case ACOS:
772 result=left.acos();
773 break;
774 case ATAN:
775 result=left.atan();
776 break;
777 case SINH:
778 result=left.sinh();
779 break;
780 case COSH:
781 result=left.cosh();
782 break;
783 case TANH:
784 result=left.tanh();
785 break;
786 case ERF:
787 result=left.erf();
788 break;
789 case ASINH:
790 result=left.asinh();
791 break;
792 case ACOSH:
793 result=left.acosh();
794 break;
795 case ATANH:
796 result=left.atanh();
797 break;
798 case LOG10:
799 result=left.log10();
800 break;
801 case LOG:
802 result=left.log();
803 break;
804 case SIGN:
805 result=left.sign();
806 break;
807 case ABS:
808 result=left.abs();
809 break;
810 case NEG:
811 result=left.neg();
812 break;
813 case POS:
814 // it doesn't mean anything for delayed.
815 // it will just trigger a deep copy of the lazy object
816 throw DataException("Programmer error - POS not supported for lazy data.");
817 break;
818 case EXP:
819 result=left.exp();
820 break;
821 case SQRT:
822 result=left.sqrt();
823 break;
824 case RECIP:
825 result=left.oneOver();
826 break;
827 case GZ:
828 result=left.wherePositive();
829 break;
830 case LZ:
831 result=left.whereNegative();
832 break;
833 case GEZ:
834 result=left.whereNonNegative();
835 break;
836 case LEZ:
837 result=left.whereNonPositive();
838 break;
839 case NEZ:
840 result=left.whereNonZero(m_tol);
841 break;
842 case EZ:
843 result=left.whereZero(m_tol);
844 break;
845 case SYM:
846 result=left.symmetric();
847 break;
848 case NSYM:
849 result=left.nonsymmetric();
850 break;
851 case PROD:
852 result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
853 break;
854 case TRANS:
855 result=left.transpose(m_axis_offset);
856 break;
857 case TRACE:
858 result=left.trace(m_axis_offset);
859 break;
860 case SWAP:
861 result=left.swapaxes(m_axis_offset, m_transpose);
862 break;
863 case MINVAL:
864 result=left.minval();
865 break;
866 case MAXVAL:
867 result=left.minval();
868 break;
869 default:
870 throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
871 }
872 return result.borrowReadyPtr();
873 }
874
875 /*
876 \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
877 This method uses the original methods on the Data class to evaluate the expressions.
878 For this reason, it should not be used on DataExpanded instances. (To do so would defeat
879 the purpose of using DataLazy in the first place).
880 */
881 void
882 DataLazy::collapse()
883 {
884 if (m_op==IDENTITY)
885 {
886 return;
887 }
888 if (m_readytype=='E')
889 { // this is more an efficiency concern than anything else
890 throw DataException("Programmer Error - do not use collapse on Expanded data.");
891 }
892 m_id=collapseToReady();
893 m_op=IDENTITY;
894 }
895
896
897
898
899
900
901 #define PROC_OP(TYPE,X) \
902 for (int j=0;j<onumsteps;++j)\
903 {\
904 for (int i=0;i<numsteps;++i,resultp+=resultStep) \
905 { \
906 LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
907 LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
908 tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
909 LAZYDEBUG(cout << " result= " << resultp[0] << endl;) \
910 lroffset+=leftstep; \
911 rroffset+=rightstep; \
912 }\
913 lroffset+=oleftstep;\
914 rroffset+=orightstep;\
915 }
916
917
918 // The result will be stored in m_samples
919 // The return value is a pointer to the DataVector, offset is the offset within the return value
920 const DataTypes::ValueType*
921 DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset)
922 {
923 LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
924 // collapse so we have a 'E' node or an IDENTITY for some other type
925 if (m_readytype!='E' && m_op!=IDENTITY)
926 {
927 collapse();
928 }
929 if (m_op==IDENTITY)
930 {
931 const ValueType& vec=m_id->getVectorRO();
932 roffset=m_id->getPointOffset(sampleNo, 0);
933 return &(vec);
934 }
935 if (m_readytype!='E')
936 {
937 throw DataException("Programmer Error - Collapse did not produce an expanded node.");
938 }
939 if (m_sampleids[tid]==sampleNo)
940 {
941 roffset=tid*m_samplesize;
942 return &(m_samples); // sample is already resolved
943 }
944 m_sampleids[tid]=sampleNo;
945 switch (getOpgroup(m_op))
946 {
947 case G_UNARY:
948 case G_UNARY_P: return resolveNodeUnary(tid, sampleNo, roffset);
949 case G_BINARY: return resolveNodeBinary(tid, sampleNo, roffset);
950 case G_NP1OUT: return resolveNodeNP1OUT(tid, sampleNo, roffset);
951 case G_NP1OUT_P: return resolveNodeNP1OUT_P(tid, sampleNo, roffset);
952 case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);
953 case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);
954 case G_REDUCTION: return resolveNodeReduction(tid, sampleNo, roffset);
955 default:
956 throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
957 }
958 }
959
960 const DataTypes::ValueType*
961 DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset)
962 {
963 // we assume that any collapsing has been done before we get here
964 // since we only have one argument we don't need to think about only
965 // processing single points.
966 // we will also know we won't get identity nodes
967 if (m_readytype!='E')
968 {
969 throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
970 }
971 if (m_op==IDENTITY)
972 {
973 throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
974 }
975 const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
976 const double* left=&((*leftres)[roffset]);
977 roffset=m_samplesize*tid;
978 double* result=&(m_samples[roffset]);
979 switch (m_op)
980 {
981 case SIN:
982 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
983 break;
984 case COS:
985 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
986 break;
987 case TAN:
988 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
989 break;
990 case ASIN:
991 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
992 break;
993 case ACOS:
994 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
995 break;
996 case ATAN:
997 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
998 break;
999 case SINH:
1000 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
1001 break;
1002 case COSH:
1003 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
1004 break;
1005 case TANH:
1006 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
1007 break;
1008 case ERF:
1009 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1010 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1011 #else
1012 tensor_unary_operation(m_samplesize, left, result, ::erf);
1013 break;
1014 #endif
1015 case ASINH:
1016 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1017 tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
1018 #else
1019 tensor_unary_operation(m_samplesize, left, result, ::asinh);
1020 #endif
1021 break;
1022 case ACOSH:
1023 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1024 tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
1025 #else
1026 tensor_unary_operation(m_samplesize, left, result, ::acosh);
1027 #endif
1028 break;
1029 case ATANH:
1030 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1031 tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
1032 #else
1033 tensor_unary_operation(m_samplesize, left, result, ::atanh);
1034 #endif
1035 break;
1036 case LOG10:
1037 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
1038 break;
1039 case LOG:
1040 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
1041 break;
1042 case SIGN:
1043 tensor_unary_operation(m_samplesize, left, result, escript::fsign);
1044 break;
1045 case ABS:
1046 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
1047 break;
1048 case NEG:
1049 tensor_unary_operation(m_samplesize, left, result, negate<double>());
1050 break;
1051 case POS:
1052 // it doesn't mean anything for delayed.
1053 // it will just trigger a deep copy of the lazy object
1054 throw DataException("Programmer error - POS not supported for lazy data.");
1055 break;
1056 case EXP:
1057 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
1058 break;
1059 case SQRT:
1060 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
1061 break;
1062 case RECIP:
1063 tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
1064 break;
1065 case GZ:
1066 tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
1067 break;
1068 case LZ:
1069 tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
1070 break;
1071 case GEZ:
1072 tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
1073 break;
1074 case LEZ:
1075 tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
1076 break;
1077 // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
1078 case NEZ:
1079 tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
1080 break;
1081 case EZ:
1082 tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
1083 break;
1084
1085 default:
1086 throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1087 }
1088 return &(m_samples);
1089 }
1090
1091
1092 const DataTypes::ValueType*
1093 DataLazy::resolveNodeReduction(int tid, int sampleNo, size_t& roffset)
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::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, loffset);
1109
1110 roffset=m_samplesize*tid;
1111 unsigned int ndpps=getNumDPPSample();
1112 unsigned int psize=DataTypes::noValues(getShape());
1113 double* result=&(m_samples[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=DataMaths::reductionOp(*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=DataMaths::reductionOp(*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);
1142 }
1143
1144 const DataTypes::ValueType*
1145 DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset)
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 ValueType* 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 DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1171 subroffset+=step;
1172 offset+=step;
1173 }
1174 break;
1175 case NSYM:
1176 for (loop=0;loop<numsteps;++loop)
1177 {
1178 DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, 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;
1187 }
1188
1189 const DataTypes::ValueType*
1190 DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset)
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 ValueType* 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 DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,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 DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, 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;
1234 }
1235
1236
1237 const DataTypes::ValueType*
1238 DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset)
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 ValueType* 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 DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, 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;
1271 }
1272
1273
1274
1275 // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1276 // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1277 // If both children are expanded, then we can process them in a single operation (we treat
1278 // the whole sample as one big datapoint.
1279 // If one of the children is not expanded, then we need to treat each point in the sample
1280 // individually.
1281 // There is an additional complication when scalar operations are considered.
1282 // For example, 2+Vector.
1283 // In this case each double within the point is treated individually
1284 const DataTypes::ValueType*
1285 DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset)
1286 {
1287 LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1288
1289 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
1290 // first work out which of the children are expanded
1291 bool leftExp=(m_left->m_readytype=='E');
1292 bool rightExp=(m_right->m_readytype=='E');
1293 if (!leftExp && !rightExp)
1294 {
1295 throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1296 }
1297 bool leftScalar=(m_left->getRank()==0);
1298 bool rightScalar=(m_right->getRank()==0);
1299 if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1300 {
1301 throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1302 }
1303 size_t leftsize=m_left->getNoValues();
1304 size_t rightsize=m_right->getNoValues();
1305 size_t chunksize=1; // how many doubles will be processed in one go
1306 int leftstep=0; // how far should the left offset advance after each step
1307 int rightstep=0;
1308 int numsteps=0; // total number of steps for the inner loop
1309 int oleftstep=0; // the o variables refer to the outer loop
1310 int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1311 int onumsteps=1;
1312
1313 bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1314 bool RES=(rightExp && rightScalar);
1315 bool LS=(!leftExp && leftScalar); // left is a single scalar
1316 bool RS=(!rightExp && rightScalar);
1317 bool LN=(!leftExp && !leftScalar); // left is a single non-scalar
1318 bool RN=(!rightExp && !rightScalar);
1319 bool LEN=(leftExp && !leftScalar); // left is an expanded non-scalar
1320 bool REN=(rightExp && !rightScalar);
1321
1322 if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1323 {
1324 chunksize=m_left->getNumDPPSample()*leftsize;
1325 leftstep=0;
1326 rightstep=0;
1327 numsteps=1;
1328 }
1329 else if (LES || RES)
1330 {
1331 chunksize=1;
1332 if (LES) // left is an expanded scalar
1333 {
1334 if (RS)
1335 {
1336 leftstep=1;
1337 rightstep=0;
1338 numsteps=m_left->getNumDPPSample();
1339 }
1340 else // RN or REN
1341 {
1342 leftstep=0;
1343 oleftstep=1;
1344 rightstep=1;
1345 orightstep=(RN ? -(int)rightsize : 0);
1346 numsteps=rightsize;
1347 onumsteps=m_left->getNumDPPSample();
1348 }
1349 }
1350 else // right is an expanded scalar
1351 {
1352 if (LS)
1353 {
1354 rightstep=1;
1355 leftstep=0;
1356 numsteps=m_right->getNumDPPSample();
1357 }
1358 else
1359 {
1360 rightstep=0;
1361 orightstep=1;
1362 leftstep=1;
1363 oleftstep=(LN ? -(int)leftsize : 0);
1364 numsteps=leftsize;
1365 onumsteps=m_right->getNumDPPSample();
1366 }
1367 }
1368 }
1369 else // this leaves (LEN, RS), (LEN, RN) and their transposes
1370 {
1371 if (LEN) // and Right will be a single value
1372 {
1373 chunksize=rightsize;
1374 leftstep=rightsize;
1375 rightstep=0;
1376 numsteps=m_left->getNumDPPSample();
1377 if (RS)
1378 {
1379 numsteps*=leftsize;
1380 }
1381 }
1382 else // REN
1383 {
1384 chunksize=leftsize;
1385 rightstep=leftsize;
1386 leftstep=0;
1387 numsteps=m_right->getNumDPPSample();
1388 if (LS)
1389 {
1390 numsteps*=rightsize;
1391 }
1392 }
1393 }
1394
1395 int resultStep=max(leftstep,rightstep); // only one (at most) should be !=0
1396 // Get the values of sub-expressions
1397 const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);
1398 const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
1399 LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1400 LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1401 LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1402 LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1403 LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1404 LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1405 LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN << endl;)
1406
1407 LAZYDEBUG(cout << "Left res["<< lroffset<< "]=" << (*left)[lroffset] << endl;)
1408 LAZYDEBUG(cout << "Right res["<< rroffset<< "]=" << (*right)[rroffset] << endl;)
1409
1410
1411 roffset=m_samplesize*tid;
1412 double* resultp=&(m_samples[roffset]); // results are stored at the vector offset we recieved
1413 switch(m_op)
1414 {
1415 case ADD:
1416 PROC_OP(NO_ARG,plus<double>());
1417 break;
1418 case SUB:
1419 PROC_OP(NO_ARG,minus<double>());
1420 break;
1421 case MUL:
1422 PROC_OP(NO_ARG,multiplies<double>());
1423 break;
1424 case DIV:
1425 PROC_OP(NO_ARG,divides<double>());
1426 break;
1427 case POW:
1428 PROC_OP(double (double,double),::pow);
1429 break;
1430 default:
1431 throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1432 }
1433 LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
1434 return &m_samples;
1435 }
1436
1437
1438 // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1439 // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1440 // unlike the other resolve helpers, we must treat these datapoints separately.
1441 const DataTypes::ValueType*
1442 DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset)
1443 {
1444 LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1445
1446 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
1447 // first work out which of the children are expanded
1448 bool leftExp=(m_left->m_readytype=='E');
1449 bool rightExp=(m_right->m_readytype=='E');
1450 int steps=getNumDPPSample();
1451 int leftStep=(leftExp? m_left->getNoValues() : 0); // do not have scalars as input to this method
1452 int rightStep=(rightExp?m_right->getNoValues() : 0);
1453
1454 int resultStep=getNoValues();
1455 roffset=m_samplesize*tid;
1456 size_t offset=roffset;
1457
1458 const ValueType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
1459
1460 const ValueType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
1461
1462 LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) << endl;
1463 cout << getNoValues() << endl;)
1464
1465
1466 LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1467 LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1468 LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1469 LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1470 LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1471 LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1472 LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1473
1474 double* resultp=&(m_samples[offset]); // results are stored at the vector offset we recieved
1475 switch(m_op)
1476 {
1477 case PROD:
1478 for (int i=0;i<steps;++i,resultp+=resultStep)
1479 {
1480 const double *ptr_0 = &((*left)[lroffset]);
1481 const double *ptr_1 = &((*right)[rroffset]);
1482
1483 LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1484 LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1485
1486 matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
1487
1488 lroffset+=leftStep;
1489 rroffset+=rightStep;
1490 }
1491 break;
1492 default:
1493 throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1494 }
1495 roffset=offset;
1496 return &m_samples;
1497 }
1498
1499
1500 const DataTypes::ValueType*
1501 DataLazy::resolveSample(int sampleNo, size_t& roffset)
1502 {
1503 #ifdef _OPENMP
1504 int tid=omp_get_thread_num();
1505 #else
1506 int tid=0;
1507 #endif
1508 return resolveNodeSample(tid, sampleNo, roffset);
1509 }
1510
1511
1512 // This needs to do the work of the identity constructor
1513 void
1514 DataLazy::resolveToIdentity()
1515 {
1516 if (m_op==IDENTITY)
1517 return;
1518 DataReady_ptr p=resolveNodeWorker();
1519 makeIdentity(p);
1520 }
1521
1522 void DataLazy::makeIdentity(const DataReady_ptr& p)
1523 {
1524 m_op=IDENTITY;
1525 m_axis_offset=0;
1526 m_transpose=0;
1527 m_SL=m_SM=m_SR=0;
1528 m_children=m_height=0;
1529 m_id=p;
1530 if(p->isConstant()) {m_readytype='C';}
1531 else if(p->isExpanded()) {m_readytype='E';}
1532 else if (p->isTagged()) {m_readytype='T';}
1533 else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1534 m_samplesize=p->getNumDPPSample()*p->getNoValues();
1535 m_left.reset();
1536 m_right.reset();
1537 }
1538
1539
1540 DataReady_ptr
1541 DataLazy::resolve()
1542 {
1543 resolveToIdentity();
1544 return m_id;
1545 }
1546
1547 // This version of resolve uses storage in each node to hold results
1548 DataReady_ptr
1549 DataLazy::resolveNodeWorker()
1550 {
1551 if (m_readytype!='E') // if the whole sub-expression is Constant or Tagged, then evaluate it normally
1552 {
1553 collapse();
1554 }
1555 if (m_op==IDENTITY) // So a lazy expression of Constant or Tagged data will be returned here.
1556 {
1557 return m_id;
1558 }
1559 // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
1560 DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(), ValueType(getNoValues()));
1561 ValueType& resvec=result->getVectorRW();
1562 DataReady_ptr resptr=DataReady_ptr(result);
1563
1564 int sample;
1565 int totalsamples=getNumSamples();
1566 const ValueType* res=0; // Storage for answer
1567 LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1568 #pragma omp parallel for private(sample,res) schedule(static)
1569 for (sample=0;sample<totalsamples;++sample)
1570 {
1571 size_t roffset=0;
1572 #ifdef _OPENMP
1573 res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
1574 #else
1575 res=resolveNodeSample(0,sample,roffset);
1576 #endif
1577 LAZYDEBUG(cout << "Sample #" << sample << endl;)
1578 LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
1579 DataVector::size_type outoffset=result->getPointOffset(sample,0);
1580 memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(DataVector::ElementType));
1581 }
1582 return resptr;
1583 }
1584
1585 std::string
1586 DataLazy::toString() const
1587 {
1588 ostringstream oss;
1589 oss << "Lazy Data:";
1590 if (escriptParams.getPRINT_LAZY_TREE()==0)
1591 {
1592 intoString(oss);
1593 }
1594 else
1595 {
1596 oss << endl;
1597 intoTreeString(oss,"");
1598 }
1599 return oss.str();
1600 }
1601
1602
1603 void
1604 DataLazy::intoString(ostringstream& oss) const
1605 {
1606 // oss << "[" << m_children <<";"<<m_height <<"]";
1607 switch (getOpgroup(m_op))
1608 {
1609 case G_IDENTITY:
1610 if (m_id->isExpanded())
1611 {
1612 oss << "E";
1613 }
1614 else if (m_id->isTagged())
1615 {
1616 oss << "T";
1617 }
1618 else if (m_id->isConstant())
1619 {
1620 oss << "C";
1621 }
1622 else
1623 {
1624 oss << "?";
1625 }
1626 oss << '@' << m_id.get();
1627 break;
1628 case G_BINARY:
1629 oss << '(';
1630 m_left->intoString(oss);
1631 oss << ' ' << opToString(m_op) << ' ';
1632 m_right->intoString(oss);
1633 oss << ')';
1634 break;
1635 case G_UNARY:
1636 case G_UNARY_P:
1637 case G_NP1OUT:
1638 case G_NP1OUT_P:
1639 case G_REDUCTION:
1640 oss << opToString(m_op) << '(';
1641 m_left->intoString(oss);
1642 oss << ')';
1643 break;
1644 case G_TENSORPROD:
1645 oss << opToString(m_op) << '(';
1646 m_left->intoString(oss);
1647 oss << ", ";
1648 m_right->intoString(oss);
1649 oss << ')';
1650 break;
1651 case G_NP1OUT_2P:
1652 oss << opToString(m_op) << '(';
1653 m_left->intoString(oss);
1654 oss << ", " << m_axis_offset << ", " << m_transpose;
1655 oss << ')';
1656 break;
1657 default:
1658 oss << "UNKNOWN";
1659 }
1660 }
1661
1662
1663 void
1664 DataLazy::intoTreeString(ostringstream& oss, string indent) const
1665 {
1666 oss << '[' << m_rank << ':' << setw(3) << m_samplesize << "] " << indent;
1667 switch (getOpgroup(m_op))
1668 {
1669 case G_IDENTITY:
1670 if (m_id->isExpanded())
1671 {
1672 oss << "E";
1673 }
1674 else if (m_id->isTagged())
1675 {
1676 oss << "T";
1677 }
1678 else if (m_id->isConstant())
1679 {
1680 oss << "C";
1681 }
1682 else
1683 {
1684 oss << "?";
1685 }
1686 oss << '@' << m_id.get() << endl;
1687 break;
1688 case G_BINARY:
1689 oss << opToString(m_op) << endl;
1690 indent+='.';
1691 m_left->intoTreeString(oss, indent);
1692 m_right->intoTreeString(oss, indent);
1693 break;
1694 case G_UNARY:
1695 case G_UNARY_P:
1696 case G_NP1OUT:
1697 case G_NP1OUT_P:
1698 case G_REDUCTION:
1699 oss << opToString(m_op) << endl;
1700 indent+='.';
1701 m_left->intoTreeString(oss, indent);
1702 break;
1703 case G_TENSORPROD:
1704 oss << opToString(m_op) << endl;
1705 indent+='.';
1706 m_left->intoTreeString(oss, indent);
1707 m_right->intoTreeString(oss, indent);
1708 break;
1709 case G_NP1OUT_2P:
1710 oss << opToString(m_op) << ", " << m_axis_offset << ", " << m_transpose<< endl;
1711 indent+='.';
1712 m_left->intoTreeString(oss, indent);
1713 break;
1714 default:
1715 oss << "UNKNOWN";
1716 }
1717 }
1718
1719
1720 DataAbstract*
1721 DataLazy::deepCopy()
1722 {
1723 switch (getOpgroup(m_op))
1724 {
1725 case G_IDENTITY: return new DataLazy(m_id->deepCopy()->getPtr());
1726 case G_UNARY:
1727 case G_REDUCTION: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1728 case G_UNARY_P: return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_tol);
1729 case G_BINARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1730 case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1731 case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1732 case G_NP1OUT_P: return new DataLazy(m_left->deepCopy()->getPtr(),m_op, m_axis_offset);
1733 case G_NP1OUT_2P: return new DataLazy(m_left->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1734 default:
1735 throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1736 }
1737 }
1738
1739
1740
1741 // There is no single, natural interpretation of getLength on DataLazy.
1742 // Instances of DataReady can look at the size of their vectors.
1743 // For lazy though, it could be the size the data would be if it were resolved;
1744 // or it could be some function of the lengths of the DataReady instances which
1745 // form part of the expression.
1746 // Rather than have people making assumptions, I have disabled the method.
1747 DataTypes::ValueType::size_type
1748 DataLazy::getLength() const
1749 {
1750 throw DataException("getLength() does not make sense for lazy data.");
1751 }
1752
1753
1754 DataAbstract*
1755 DataLazy::getSlice(const DataTypes::RegionType& region) const
1756 {
1757 throw DataException("getSlice - not implemented for Lazy objects.");
1758 }
1759
1760
1761 // To do this we need to rely on our child nodes
1762 DataTypes::ValueType::size_type
1763 DataLazy::getPointOffset(int sampleNo,
1764 int dataPointNo)
1765 {
1766 if (m_op==IDENTITY)
1767 {
1768 return m_id->getPointOffset(sampleNo,dataPointNo);
1769 }
1770 if (m_readytype!='E')
1771 {
1772 collapse();
1773 return m_id->getPointOffset(sampleNo,dataPointNo);
1774 }
1775 // at this point we do not have an identity node and the expression will be Expanded
1776 // so we only need to know which child to ask
1777 if (m_left->m_readytype=='E')
1778 {
1779 return m_left->getPointOffset(sampleNo,dataPointNo);
1780 }
1781 else
1782 {
1783 return m_right->getPointOffset(sampleNo,dataPointNo);
1784 }
1785 }
1786
1787 // To do this we need to rely on our child nodes
1788 DataTypes::ValueType::size_type
1789 DataLazy::getPointOffset(int sampleNo,
1790 int dataPointNo) const
1791 {
1792 if (m_op==IDENTITY)
1793 {
1794 return m_id->getPointOffset(sampleNo,dataPointNo);
1795 }
1796 if (m_readytype=='E')
1797 {
1798 // at this point we do not have an identity node and the expression will be Expanded
1799 // so we only need to know which child to ask
1800 if (m_left->m_readytype=='E')
1801 {
1802 return m_left->getPointOffset(sampleNo,dataPointNo);
1803 }
1804 else
1805 {
1806 return m_right->getPointOffset(sampleNo,dataPointNo);
1807 }
1808 }
1809 if (m_readytype=='C')
1810 {
1811 return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
1812 }
1813 throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1814 }
1815
1816
1817 // I have decided to let Data:: handle this issue.
1818 void
1819 DataLazy::setToZero()
1820 {
1821 // DataTypes::ValueType v(getNoValues(),0);
1822 // m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1823 // m_op=IDENTITY;
1824 // m_right.reset();
1825 // m_left.reset();
1826 // m_readytype='C';
1827 // m_buffsRequired=1;
1828
1829 privdebug=privdebug; // to stop the compiler complaining about unused privdebug
1830 throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
1831 }
1832
1833 bool
1834 DataLazy::actsExpanded() const
1835 {
1836 return (m_readytype=='E');
1837 }
1838
1839 } // end namespace

  ViewVC Help
Powered by ViewVC 1.1.26