/[escript]/branches/clazy/escriptcore/src/DataLazy.cpp
ViewVC logotype

Contents of /branches/clazy/escriptcore/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1908 - (show annotations)
Thu Oct 23 01:35:31 2008 UTC (10 years, 6 months ago) by jfenwick
Original Path: branches/schroedinger/escript/src/DataLazy.cpp
File size: 25778 byte(s)
Branch commit.
Fixed to support scalar op non-scalar operations.

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

  ViewVC Help
Powered by ViewVC 1.1.26