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

  ViewVC Help
Powered by ViewVC 1.1.26