/[escript]/branches/arrexp_2137_win_merge/escript/src/DataLazy.cpp
ViewVC logotype

Contents of /branches/arrexp_2137_win_merge/escript/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26