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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1972 - (show annotations)
Thu Nov 6 01:49:39 2008 UTC (10 years, 6 months ago) by jfenwick
Original Path: branches/schroedinger_upto1946/escript/src/DataLazy.cpp
File size: 28084 byte(s)
Branch commit.
Fixed to help windows determine template params.

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

  ViewVC Help
Powered by ViewVC 1.1.26