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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1947 - (show annotations)
Wed Oct 29 23:19:45 2008 UTC (10 years, 4 months ago) by jfenwick
File size: 27785 byte(s)
This does not actually have the changes in it yet.

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

  ViewVC Help
Powered by ViewVC 1.1.26