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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2021 - (show annotations)
Mon Nov 10 14:07:49 2008 UTC (11 years ago) by phornby
File size: 28145 byte(s)
Use the NO_ARG trick on the PROC_OP macro expansions.


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

  ViewVC Help
Powered by ViewVC 1.1.26