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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1910 - (show annotations)
Thu Oct 23 03:05:28 2008 UTC (10 years, 6 months ago) by jfenwick
File size: 25840 byte(s)
Branch commit.
Support for ** added.

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

  ViewVC Help
Powered by ViewVC 1.1.26