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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2037 - (show annotations)
Thu Nov 13 06:17:12 2008 UTC (10 years, 10 months ago) by jfenwick
File size: 30790 byte(s)
Fixed some warnings in the unit tests.
Added support for symmetric and nonsymmetric operations on LazyData.
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 /*
32 How does DataLazy work?
33 ~~~~~~~~~~~~~~~~~~~~~~~
34
35 Each instance represents a single operation on one or two other DataLazy instances. These arguments are normally
36 denoted left and right.
37
38 A special operation, IDENTITY, stores an instance of DataReady in the m_id member.
39 This means that all "internal" nodes in the structure are instances of DataLazy.
40
41 Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
42 Note that IDENITY is not considered a unary operation.
43
44 I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
45 It must however form a DAG (directed acyclic graph).
46 I will refer to individual DataLazy objects with the structure as nodes.
47
48 Each node also stores:
49 - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
50 evaluated.
51 - m_length ~ how many values would be stored in the answer if the expression was evaluated.
52 - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to
53 evaluate the expression.
54 - m_samplesize ~ the number of doubles stored in a sample.
55
56 When a new node is created, the above values are computed based on the values in the child nodes.
57 Eg: if left requires 4 samples and right requires 6 then left+right requires 7 samples.
58
59 The resolve method, which produces a DataReady from a DataLazy, does the following:
60 1) Create a DataReady to hold the new result.
61 2) Allocate a vector (v) big enough to hold m_buffsrequired samples.
62 3) For each sample, call resolveSample with v, to get its values and copy them into the result object.
63
64 (In the case of OMP, multiple samples are resolved in parallel so the vector needs to be larger.)
65
66 resolveSample returns a Vector* and an offset within that vector where the result is stored.
67 Normally, this would be v, but for identity nodes their internal vector is returned instead.
68
69 The convention that I use, is that the resolve methods should store their results starting at the offset they are passed.
70
71 For expressions which evaluate to Constant or Tagged, there is a different evaluation method.
72 The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.
73
74 To add a new operator you need to do the following (plus anything I might have forgotten):
75 1) Add to the ES_optype.
76 2) determine what opgroup your operation belongs to (X)
77 3) add a string for the op to the end of ES_opstrings
78 4) increase ES_opcount
79 5) add an entry (X) to opgroups
80 6) add an entry to the switch in collapseToReady
81 7) add an entry to resolveX
82 */
83
84
85 using namespace std;
86 using namespace boost;
87
88 namespace escript
89 {
90
91 namespace
92 {
93
94 enum ES_opgroup
95 {
96 G_UNKNOWN,
97 G_IDENTITY,
98 G_BINARY, // pointwise operations with two arguments
99 G_UNARY, // pointwise operations with one argument
100 G_NP1OUT // non-pointwise op with one output
101 };
102
103
104
105
106 string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
107 "sin","cos","tan",
108 "asin","acos","atan","sinh","cosh","tanh","erf",
109 "asinh","acosh","atanh",
110 "log10","log","sign","abs","neg","pos","exp","sqrt",
111 "1/","where>0","where<0","where>=0","where<=0",
112 "symmetric","nonsymmetric"};
113 int ES_opcount=35;
114 ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
115 G_UNARY,G_UNARY,G_UNARY, //10
116 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 17
117 G_UNARY,G_UNARY,G_UNARY, // 20
118 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 28
119 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 33
120 G_NP1OUT,G_NP1OUT};
121 inline
122 ES_opgroup
123 getOpgroup(ES_optype op)
124 {
125 return opgroups[op];
126 }
127
128 // return the FunctionSpace of the result of "left op right"
129 FunctionSpace
130 resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
131 {
132 // perhaps this should call interpolate and throw or something?
133 // maybe we need an interpolate node -
134 // that way, if interpolate is required in any other op we can just throw a
135 // programming error exception.
136
137 FunctionSpace l=left->getFunctionSpace();
138 FunctionSpace r=right->getFunctionSpace();
139 if (l!=r)
140 {
141 if (r.probeInterpolation(l))
142 {
143 return l;
144 }
145 if (l.probeInterpolation(r))
146 {
147 return r;
148 }
149 throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
150 }
151 return l;
152 }
153
154 // return the shape of the result of "left op right"
155 DataTypes::ShapeType
156 resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
157 {
158 if (left->getShape()!=right->getShape())
159 {
160 if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
161 {
162 throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
163 }
164 if (left->getRank()==0) // we need to allow scalar * anything
165 {
166 return right->getShape();
167 }
168 if (right->getRank()==0)
169 {
170 return left->getShape();
171 }
172 throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
173 }
174 return left->getShape();
175 }
176
177 // determine the number of points in the result of "left op right"
178 size_t
179 resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
180 {
181 switch (getOpgroup(op))
182 {
183 case G_BINARY: return left->getLength();
184 case G_UNARY: return left->getLength();
185 case G_NP1OUT: return left->getLength();
186 default:
187 throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
188 }
189 }
190
191 // determine the number of samples requires to evaluate an expression combining left and right
192 // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
193 int
194 calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
195 {
196 switch(getOpgroup(op))
197 {
198 case G_IDENTITY: return 1;
199 case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
200 case G_UNARY: return max(left->getBuffsRequired(),1);
201 case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
202 default:
203 throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
204 }
205 }
206
207
208 } // end anonymous namespace
209
210
211
212 // Return a string representing the operation
213 const std::string&
214 opToString(ES_optype op)
215 {
216 if (op<0 || op>=ES_opcount)
217 {
218 op=UNKNOWNOP;
219 }
220 return ES_opstrings[op];
221 }
222
223
224 DataLazy::DataLazy(DataAbstract_ptr p)
225 : parent(p->getFunctionSpace(),p->getShape()),
226 m_op(IDENTITY)
227 {
228 if (p->isLazy())
229 {
230 // I don't want identity of Lazy.
231 // Question: Why would that be so bad?
232 // Answer: We assume that the child of ID is something we can call getVector on
233 throw DataException("Programmer error - attempt to create identity from a DataLazy.");
234 }
235 else
236 {
237 m_id=dynamic_pointer_cast<DataReady>(p);
238 if(p->isConstant()) {m_readytype='C';}
239 else if(p->isExpanded()) {m_readytype='E';}
240 else if (p->isTagged()) {m_readytype='T';}
241 else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}
242 }
243 m_length=p->getLength();
244 m_buffsRequired=1;
245 m_samplesize=getNumDPPSample()*getNoValues();
246 cout << "(1)Lazy created with " << m_samplesize << endl;
247 }
248
249
250
251
252 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
253 : parent(left->getFunctionSpace(),left->getShape()),
254 m_op(op)
255 {
256 if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
257 {
258 throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
259 }
260 DataLazy_ptr lleft;
261 if (!left->isLazy())
262 {
263 lleft=DataLazy_ptr(new DataLazy(left));
264 }
265 else
266 {
267 lleft=dynamic_pointer_cast<DataLazy>(left);
268 }
269 m_readytype=lleft->m_readytype;
270 m_length=left->getLength();
271 m_left=lleft;
272 m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
273 m_samplesize=getNumDPPSample()*getNoValues();
274 }
275
276
277 // In this constructor we need to consider interpolation
278 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
279 : parent(resultFS(left,right,op), resultShape(left,right,op)),
280 m_op(op)
281 {
282 if ((getOpgroup(op)!=G_BINARY))
283 {
284 cout << opToString(op) << endl;
285 throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
286 }
287
288 if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
289 {
290 FunctionSpace fs=getFunctionSpace();
291 Data ltemp(left);
292 Data tmp(ltemp,fs);
293 left=tmp.borrowDataPtr();
294 }
295 if (getFunctionSpace()!=right->getFunctionSpace()) // left needs to be interpolated
296 {
297 Data tmp(Data(right),getFunctionSpace());
298 right=tmp.borrowDataPtr();
299 }
300 left->operandCheck(*right);
301
302 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
303 {
304 m_left=dynamic_pointer_cast<DataLazy>(left);
305 }
306 else
307 {
308 m_left=DataLazy_ptr(new DataLazy(left));
309 }
310 if (right->isLazy())
311 {
312 m_right=dynamic_pointer_cast<DataLazy>(right);
313 }
314 else
315 {
316 m_right=DataLazy_ptr(new DataLazy(right));
317 }
318 char lt=m_left->m_readytype;
319 char rt=m_right->m_readytype;
320 if (lt=='E' || rt=='E')
321 {
322 m_readytype='E';
323 }
324 else if (lt=='T' || rt=='T')
325 {
326 m_readytype='T';
327 }
328 else
329 {
330 m_readytype='C';
331 }
332 m_length=resultLength(m_left,m_right,m_op);
333 m_samplesize=getNumDPPSample()*getNoValues();
334 m_buffsRequired=calcBuffs(m_left, m_right,m_op);
335 cout << "(3)Lazy created with " << m_samplesize << endl;
336 }
337
338
339 DataLazy::~DataLazy()
340 {
341 }
342
343
344 int
345 DataLazy::getBuffsRequired() const
346 {
347 return m_buffsRequired;
348 }
349
350
351 /*
352 \brief Evaluates the expression using methods on Data.
353 This does the work for the collapse method.
354 For reasons of efficiency do not call this method on DataExpanded nodes.
355 */
356 DataReady_ptr
357 DataLazy::collapseToReady()
358 {
359 if (m_readytype=='E')
360 { // this is more an efficiency concern than anything else
361 throw DataException("Programmer Error - do not use collapse on Expanded data.");
362 }
363 if (m_op==IDENTITY)
364 {
365 return m_id;
366 }
367 DataReady_ptr pleft=m_left->collapseToReady();
368 Data left(pleft);
369 Data right;
370 if (getOpgroup(m_op)==G_BINARY)
371 {
372 right=Data(m_right->collapseToReady());
373 }
374 Data result;
375 switch(m_op)
376 {
377 case ADD:
378 result=left+right;
379 break;
380 case SUB:
381 result=left-right;
382 break;
383 case MUL:
384 result=left*right;
385 break;
386 case DIV:
387 result=left/right;
388 break;
389 case SIN:
390 result=left.sin();
391 break;
392 case COS:
393 result=left.cos();
394 break;
395 case TAN:
396 result=left.tan();
397 break;
398 case ASIN:
399 result=left.asin();
400 break;
401 case ACOS:
402 result=left.acos();
403 break;
404 case ATAN:
405 result=left.atan();
406 break;
407 case SINH:
408 result=left.sinh();
409 break;
410 case COSH:
411 result=left.cosh();
412 break;
413 case TANH:
414 result=left.tanh();
415 break;
416 case ERF:
417 result=left.erf();
418 break;
419 case ASINH:
420 result=left.asinh();
421 break;
422 case ACOSH:
423 result=left.acosh();
424 break;
425 case ATANH:
426 result=left.atanh();
427 break;
428 case LOG10:
429 result=left.log10();
430 break;
431 case LOG:
432 result=left.log();
433 break;
434 case SIGN:
435 result=left.sign();
436 break;
437 case ABS:
438 result=left.abs();
439 break;
440 case NEG:
441 result=left.neg();
442 break;
443 case POS:
444 // it doesn't mean anything for delayed.
445 // it will just trigger a deep copy of the lazy object
446 throw DataException("Programmer error - POS not supported for lazy data.");
447 break;
448 case EXP:
449 result=left.exp();
450 break;
451 case SQRT:
452 result=left.sqrt();
453 break;
454 case RECIP:
455 result=left.oneOver();
456 break;
457 case GZ:
458 result=left.wherePositive();
459 break;
460 case LZ:
461 result=left.whereNegative();
462 break;
463 case GEZ:
464 result=left.whereNonNegative();
465 break;
466 case LEZ:
467 result=left.whereNonPositive();
468 break;
469 case SYM:
470 result=left.symmetric();
471 break;
472 case NSYM:
473 result=left.nonsymmetric();
474 break;
475 default:
476 throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
477 }
478 return result.borrowReadyPtr();
479 }
480
481 /*
482 \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
483 This method uses the original methods on the Data class to evaluate the expressions.
484 For this reason, it should not be used on DataExpanded instances. (To do so would defeat
485 the purpose of using DataLazy in the first place).
486 */
487 void
488 DataLazy::collapse()
489 {
490 if (m_op==IDENTITY)
491 {
492 return;
493 }
494 if (m_readytype=='E')
495 { // this is more an efficiency concern than anything else
496 throw DataException("Programmer Error - do not use collapse on Expanded data.");
497 }
498 m_id=collapseToReady();
499 m_op=IDENTITY;
500 }
501
502 /*
503 \brief Compute the value of the expression (unary operation) for the given sample.
504 \return Vector which stores the value of the subexpression for the given sample.
505 \param v A vector to store intermediate results.
506 \param offset Index in v to begin storing results.
507 \param sampleNo Sample number to evaluate.
508 \param roffset (output parameter) the offset in the return vector where the result begins.
509
510 The return value will be an existing vector so do not deallocate it.
511 If the result is stored in v it should be stored at the offset given.
512 Everything from offset to the end of v should be considered available for this method to use.
513 */
514 DataTypes::ValueType*
515 DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
516 {
517 // we assume that any collapsing has been done before we get here
518 // since we only have one argument we don't need to think about only
519 // processing single points.
520 if (m_readytype!='E')
521 {
522 throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
523 }
524 const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
525 const double* left=&((*vleft)[roffset]);
526 double* result=&(v[offset]);
527 roffset=offset;
528 switch (m_op)
529 {
530 case SIN:
531 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
532 break;
533 case COS:
534 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
535 break;
536 case TAN:
537 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
538 break;
539 case ASIN:
540 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
541 break;
542 case ACOS:
543 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
544 break;
545 case ATAN:
546 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
547 break;
548 case SINH:
549 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
550 break;
551 case COSH:
552 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
553 break;
554 case TANH:
555 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
556 break;
557 case ERF:
558 #ifdef _WIN32
559 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
560 #else
561 tensor_unary_operation(m_samplesize, left, result, ::erf);
562 break;
563 #endif
564 case ASINH:
565 #ifdef _WIN32
566 tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
567 #else
568 tensor_unary_operation(m_samplesize, left, result, ::asinh);
569 #endif
570 break;
571 case ACOSH:
572 #ifdef _WIN32
573 tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
574 #else
575 tensor_unary_operation(m_samplesize, left, result, ::acosh);
576 #endif
577 break;
578 case ATANH:
579 #ifdef _WIN32
580 tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
581 #else
582 tensor_unary_operation(m_samplesize, left, result, ::atanh);
583 #endif
584 break;
585 case LOG10:
586 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
587 break;
588 case LOG:
589 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
590 break;
591 case SIGN:
592 tensor_unary_operation(m_samplesize, left, result, escript::fsign);
593 break;
594 case ABS:
595 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
596 break;
597 case NEG:
598 tensor_unary_operation(m_samplesize, left, result, negate<double>());
599 break;
600 case POS:
601 // it doesn't mean anything for delayed.
602 // it will just trigger a deep copy of the lazy object
603 throw DataException("Programmer error - POS not supported for lazy data.");
604 break;
605 case EXP:
606 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
607 break;
608 case SQRT:
609 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
610 break;
611 case RECIP:
612 tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
613 break;
614 case GZ:
615 tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
616 break;
617 case LZ:
618 tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
619 break;
620 case GEZ:
621 tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
622 break;
623 case LEZ:
624 tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
625 break;
626
627 default:
628 throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
629 }
630 return &v;
631 }
632
633
634 /*
635 \brief Compute the value of the expression (unary operation) for the given sample.
636 \return Vector which stores the value of the subexpression for the given sample.
637 \param v A vector to store intermediate results.
638 \param offset Index in v to begin storing results.
639 \param sampleNo Sample number to evaluate.
640 \param roffset (output parameter) the offset in the return vector where the result begins.
641
642 The return value will be an existing vector so do not deallocate it.
643 If the result is stored in v it should be stored at the offset given.
644 Everything from offset to the end of v should be considered available for this method to use.
645 */
646 DataTypes::ValueType*
647 DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
648 {
649 // we assume that any collapsing has been done before we get here
650 // since we only have one argument we don't need to think about only
651 // processing single points.
652 if (m_readytype!='E')
653 {
654 throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
655 }
656 // since we can't write the result over the input, we need a result offset further along
657 size_t subroffset=roffset+m_samplesize;
658 const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);
659 roffset=offset;
660 switch (m_op)
661 {
662 case SYM:
663 DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
664 break;
665 case NSYM:
666 DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
667 break;
668 default:
669 throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
670 }
671 return &v;
672 }
673
674
675
676
677 #define PROC_OP(TYPE,X) \
678 for (int i=0;i<steps;++i,resultp+=resultStep) \
679 { \
680 tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
681 lroffset+=leftStep; \
682 rroffset+=rightStep; \
683 }
684
685 /*
686 \brief Compute the value of the expression (binary operation) for the given sample.
687 \return Vector which stores the value of the subexpression for the given sample.
688 \param v A vector to store intermediate results.
689 \param offset Index in v to begin storing results.
690 \param sampleNo Sample number to evaluate.
691 \param roffset (output parameter) the offset in the return vector where the result begins.
692
693 The return value will be an existing vector so do not deallocate it.
694 If the result is stored in v it should be stored at the offset given.
695 Everything from offset to the end of v should be considered available for this method to use.
696 */
697 // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
698 // have already been collapsed to IDENTITY. So we must have at least one expanded child.
699 // If both children are expanded, then we can process them in a single operation (we treat
700 // the whole sample as one big datapoint.
701 // If one of the children is not expanded, then we need to treat each point in the sample
702 // individually.
703 // There is an additional complication when scalar operations are considered.
704 // For example, 2+Vector.
705 // In this case each double within the point is treated individually
706 DataTypes::ValueType*
707 DataLazy::resolveBinary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
708 {
709 cout << "Resolve binary: " << toString() << endl;
710
711 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
712 // first work out which of the children are expanded
713 bool leftExp=(m_left->m_readytype=='E');
714 bool rightExp=(m_right->m_readytype=='E');
715 bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
716 int steps=(bigloops?1:getNumDPPSample());
717 size_t chunksize=(bigloops? m_samplesize : getNoValues()); // if bigloops, pretend the whole sample is a datapoint
718 if (m_left->getRank()!=m_right->getRank()) // need to deal with scalar * ? ops
719 {
720 EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");
721 steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());
722 chunksize=1; // for scalar
723 }
724 int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
725 int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
726 int resultStep=max(leftStep,rightStep); // only one (at most) should be !=0
727 // Get the values of sub-expressions
728 const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
729 const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
730 // the right child starts further along.
731 double* resultp=&(v[offset]); // results are stored at the vector offset we recieved
732 switch(m_op)
733 {
734 case ADD:
735 PROC_OP(NO_ARG,plus<double>());
736 break;
737 case SUB:
738 PROC_OP(NO_ARG,minus<double>());
739 break;
740 case MUL:
741 PROC_OP(NO_ARG,multiplies<double>());
742 break;
743 case DIV:
744 PROC_OP(NO_ARG,divides<double>());
745 break;
746 case POW:
747 PROC_OP(double (double,double),::pow);
748 break;
749 default:
750 throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
751 }
752 roffset=offset;
753 return &v;
754 }
755
756
757
758 /*
759 \brief Compute the value of the expression for the given sample.
760 \return Vector which stores the value of the subexpression for the given sample.
761 \param v A vector to store intermediate results.
762 \param offset Index in v to begin storing results.
763 \param sampleNo Sample number to evaluate.
764 \param roffset (output parameter) the offset in the return vector where the result begins.
765
766 The return value will be an existing vector so do not deallocate it.
767 */
768 // the vector and the offset are a place where the method could write its data if it wishes
769 // it is not obligated to do so. For example, if it has its own storage already, it can use that.
770 // Hence the return value to indicate where the data is actually stored.
771 // Regardless, the storage should be assumed to be used, even if it isn't.
772
773 // the roffset is the offset within the returned vector where the data begins
774 const DataTypes::ValueType*
775 DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
776 {
777 cout << "Resolve sample " << toString() << endl;
778 // collapse so we have a 'E' node or an IDENTITY for some other type
779 if (m_readytype!='E' && m_op!=IDENTITY)
780 {
781 collapse();
782 }
783 if (m_op==IDENTITY)
784 {
785 const ValueType& vec=m_id->getVector();
786 if (m_readytype=='C')
787 {
788 roffset=0;
789 return &(vec);
790 }
791 roffset=m_id->getPointOffset(sampleNo, 0);
792 return &(vec);
793 }
794 if (m_readytype!='E')
795 {
796 throw DataException("Programmer Error - Collapse did not produce an expanded node.");
797 }
798 switch (getOpgroup(m_op))
799 {
800 case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);
801 case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
802 default:
803 throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
804 }
805 }
806
807
808 // To simplify the memory management, all threads operate on one large vector, rather than one each.
809 // Each sample is evaluated independently and copied into the result DataExpanded.
810 DataReady_ptr
811 DataLazy::resolve()
812 {
813
814 cout << "Sample size=" << m_samplesize << endl;
815 cout << "Buffers=" << m_buffsRequired << endl;
816
817 if (m_readytype!='E') // if the whole sub-expression is Constant or Tagged, then evaluate it normally
818 {
819 collapse();
820 }
821 if (m_op==IDENTITY) // So a lazy expression of Constant or Tagged data will be returned here.
822 {
823 return m_id;
824 }
825 // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
826 size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough
827 // storage to evaluate its expression
828 int numthreads=1;
829 #ifdef _OPENMP
830 numthreads=getNumberOfThreads();
831 int threadnum=0;
832 #endif
833 ValueType v(numthreads*threadbuffersize);
834 cout << "Buffer created with size=" << v.size() << endl;
835 DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(), ValueType(getNoValues()));
836 ValueType& resvec=result->getVector();
837 DataReady_ptr resptr=DataReady_ptr(result);
838 int sample;
839 size_t outoffset; // offset in the output data
840 int totalsamples=getNumSamples();
841 const ValueType* res=0; // Vector storing the answer
842 size_t resoffset=0; // where in the vector to find the answer
843 #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)
844 for (sample=0;sample<totalsamples;++sample)
845 {
846 cout << "################################# " << sample << endl;
847 #ifdef _OPENMP
848 res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
849 #else
850 res=resolveSample(v,0,sample,resoffset); // res would normally be v, but not if its a single IDENTITY op.
851 #endif
852 cerr << "-------------------------------- " << endl;
853 outoffset=result->getPointOffset(sample,0);
854 cerr << "offset=" << outoffset << endl;
855 for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset) // copy values into the output vector
856 {
857 resvec[outoffset]=(*res)[resoffset];
858 }
859 cerr << "*********************************" << endl;
860 }
861 return resptr;
862 }
863
864 std::string
865 DataLazy::toString() const
866 {
867 ostringstream oss;
868 oss << "Lazy Data:";
869 intoString(oss);
870 return oss.str();
871 }
872
873
874 void
875 DataLazy::intoString(ostringstream& oss) const
876 {
877 switch (getOpgroup(m_op))
878 {
879 case G_IDENTITY:
880 if (m_id->isExpanded())
881 {
882 oss << "E";
883 }
884 else if (m_id->isTagged())
885 {
886 oss << "T";
887 }
888 else if (m_id->isConstant())
889 {
890 oss << "C";
891 }
892 else
893 {
894 oss << "?";
895 }
896 oss << '@' << m_id.get();
897 break;
898 case G_BINARY:
899 oss << '(';
900 m_left->intoString(oss);
901 oss << ' ' << opToString(m_op) << ' ';
902 m_right->intoString(oss);
903 oss << ')';
904 break;
905 case G_UNARY:
906 case G_NP1OUT:
907 oss << opToString(m_op) << '(';
908 m_left->intoString(oss);
909 oss << ')';
910 break;
911 default:
912 oss << "UNKNOWN";
913 }
914 }
915
916 DataAbstract*
917 DataLazy::deepCopy()
918 {
919 switch (getOpgroup(m_op))
920 {
921 case G_IDENTITY: return new DataLazy(m_id->deepCopy()->getPtr());
922 case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
923 case G_BINARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
924 default:
925 throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
926 }
927 }
928
929
930 DataTypes::ValueType::size_type
931 DataLazy::getLength() const
932 {
933 return m_length;
934 }
935
936
937 DataAbstract*
938 DataLazy::getSlice(const DataTypes::RegionType& region) const
939 {
940 throw DataException("getSlice - not implemented for Lazy objects.");
941 }
942
943
944 // To do this we need to rely on our child nodes
945 DataTypes::ValueType::size_type
946 DataLazy::getPointOffset(int sampleNo,
947 int dataPointNo)
948 {
949 if (m_op==IDENTITY)
950 {
951 return m_id->getPointOffset(sampleNo,dataPointNo);
952 }
953 if (m_readytype!='E')
954 {
955 collapse();
956 return m_id->getPointOffset(sampleNo,dataPointNo);
957 }
958 // at this point we do not have an identity node and the expression will be Expanded
959 // so we only need to know which child to ask
960 if (m_left->m_readytype=='E')
961 {
962 return m_left->getPointOffset(sampleNo,dataPointNo);
963 }
964 else
965 {
966 return m_right->getPointOffset(sampleNo,dataPointNo);
967 }
968 }
969
970 // To do this we need to rely on our child nodes
971 DataTypes::ValueType::size_type
972 DataLazy::getPointOffset(int sampleNo,
973 int dataPointNo) const
974 {
975 if (m_op==IDENTITY)
976 {
977 return m_id->getPointOffset(sampleNo,dataPointNo);
978 }
979 if (m_readytype=='E')
980 {
981 // at this point we do not have an identity node and the expression will be Expanded
982 // so we only need to know which child to ask
983 if (m_left->m_readytype=='E')
984 {
985 return m_left->getPointOffset(sampleNo,dataPointNo);
986 }
987 else
988 {
989 return m_right->getPointOffset(sampleNo,dataPointNo);
990 }
991 }
992 if (m_readytype=='C')
993 {
994 return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
995 }
996 throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
997 }
998
999 // It would seem that DataTagged will need to be treated differently since even after setting all tags
1000 // to zero, all the tags from all the DataTags would be in the result.
1001 // However since they all have the same value (0) whether they are there or not should not matter.
1002 // So I have decided that for all types this method will create a constant 0.
1003 // It can be promoted up as required.
1004 // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
1005 // but we can deal with that if it arrises.
1006 void
1007 DataLazy::setToZero()
1008 {
1009 DataTypes::ValueType v(getNoValues(),0);
1010 m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1011 m_op=IDENTITY;
1012 m_right.reset();
1013 m_left.reset();
1014 m_readytype='C';
1015 m_buffsRequired=1;
1016 }
1017
1018 } // end namespace

  ViewVC Help
Powered by ViewVC 1.1.26