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

  ViewVC Help
Powered by ViewVC 1.1.26