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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1889 - (hide annotations)
Thu Oct 16 05:57:09 2008 UTC (10 years, 6 months ago) by jfenwick
Original Path: branches/schroedinger/escript/src/DataLazy.cpp
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 jfenwick 1865
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 jfenwick 1885 #ifdef _OPENMP
23     #include <omp.h>
24     #endif
25 jfenwick 1865 #include "FunctionSpace.h"
26     #include "DataTypes.h"
27 jfenwick 1868 #include "Data.h"
28 jfenwick 1886 #include "UnaryFuncs.h" // for escript::fsign
29 jfenwick 1865
30     using namespace std;
31 jfenwick 1868 using namespace boost;
32 jfenwick 1865
33     namespace escript
34     {
35    
36     const std::string&
37     opToString(ES_optype op);
38    
39     namespace
40     {
41    
42 jfenwick 1886
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 jfenwick 1888 "log10","log","sign","abs","neg","pos","exp","sqrt",
59     "1/","where>0","where<0","where>=0","where<=0"};
60     int ES_opcount=32;
61 jfenwick 1886 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 jfenwick 1888 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 jfenwick 1886 inline
67     ES_opgroup
68     getOpgroup(ES_optype op)
69     {
70     return opgroups[op];
71     }
72    
73 jfenwick 1865 // 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 jfenwick 1879
82    
83     if (left->getFunctionSpace()!=right->getFunctionSpace())
84     {
85     throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");
86     }
87 jfenwick 1865 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 jfenwick 1879 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 jfenwick 1865 }
100    
101     size_t
102     resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
103     {
104 jfenwick 1886 switch (getOpgroup(op))
105 jfenwick 1865 {
106 jfenwick 1886 case G_BINARY: return left->getLength();
107     case G_UNARY: return left->getLength();
108 jfenwick 1865 default:
109     throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
110     }
111     }
112    
113 jfenwick 1879 int
114     calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
115     {
116 jfenwick 1886 switch(getOpgroup(op))
117 jfenwick 1879 {
118 jfenwick 1886 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 jfenwick 1879 default:
122     throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
123     }
124     }
125    
126 jfenwick 1865 } // 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 jfenwick 1879 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 jfenwick 1889 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 jfenwick 1879 }
160     m_length=p->getLength();
161 jfenwick 1886 m_buffsRequired=1;
162 jfenwick 1879 m_samplesize=getNumDPPSample()*getNoValues();
163     cout << "(1)Lazy created with " << m_samplesize << endl;
164 jfenwick 1865 }
165    
166 jfenwick 1886 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 jfenwick 1889 m_readytype=lleft->m_readytype;
184 jfenwick 1886 m_length=left->getLength();
185     m_left=lleft;
186     m_buffsRequired=1;
187     m_samplesize=getNumDPPSample()*getNoValues();
188     }
189    
190    
191 jfenwick 1889 // 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 jfenwick 1865
207 jfenwick 1879 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 jfenwick 1886 if (getOpgroup(op)!=G_BINARY)
212     {
213 jfenwick 1889 throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
214 jfenwick 1886 }
215 jfenwick 1879 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 jfenwick 1889 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 jfenwick 1879 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 jfenwick 1865 DataLazy::~DataLazy()
253     {
254     }
255    
256 jfenwick 1879
257     int
258     DataLazy::getBuffsRequired() const
259 jfenwick 1865 {
260 jfenwick 1879 return m_buffsRequired;
261     }
262    
263 jfenwick 1884
264 jfenwick 1889 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 jfenwick 1884 // 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 jfenwick 1885 // 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 jfenwick 1884 const double*
570 jfenwick 1889 DataLazy::resolveSample(ValueType& v,int sampleNo, size_t offset )
571 jfenwick 1879 {
572 jfenwick 1889 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 jfenwick 1885 if (m_op==IDENTITY)
580 jfenwick 1865 {
581 jfenwick 1889 cout << "In IDENTITY check\n";
582 jfenwick 1879 const ValueType& vec=m_id->getVector();
583 jfenwick 1889 if (m_readytype=='C')
584     {
585     return &(vec[0]);
586     }
587 jfenwick 1884 return &(vec[m_id->getPointOffset(sampleNo, 0)]);
588 jfenwick 1865 }
589 jfenwick 1889 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 jfenwick 1879 size_t rightoffset=offset+m_samplesize;
621 jfenwick 1884 const double* left=m_left->resolveSample(v,sampleNo,offset);
622 jfenwick 1886 const double* right=0;
623     if (getOpgroup(m_op)==G_BINARY)
624     {
625     right=m_right->resolveSample(v,sampleNo,rightoffset);
626     }
627 jfenwick 1884 double* result=&(v[offset]);
628 jfenwick 1865 {
629 jfenwick 1879 switch(m_op)
630     {
631     case ADD: // since these are pointwise ops, pretend each sample is one point
632 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, plus<double>());
633 jfenwick 1879 break;
634     case SUB:
635 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, minus<double>());
636 jfenwick 1879 break;
637     case MUL:
638 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());
639 jfenwick 1879 break;
640     case DIV:
641 jfenwick 1884 tensor_binary_operation(m_samplesize, left, right, result, divides<double>());
642 jfenwick 1879 break;
643 jfenwick 1886 // 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 jfenwick 1888 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 jfenwick 1868 default:
742     throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
743 jfenwick 1879 }
744 jfenwick 1868 }
745 jfenwick 1884 return result;
746 jfenwick 1865 }
747    
748 jfenwick 1879 DataReady_ptr
749     DataLazy::resolve()
750     {
751    
752     cout << "Sample size=" << m_samplesize << endl;
753     cout << "Buffers=" << m_buffsRequired << endl;
754    
755 jfenwick 1889 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 jfenwick 1885 size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);
765     int numthreads=1;
766     #ifdef _OPENMP
767 jfenwick 1889 numthreads=getNumberOfThreads();
768 jfenwick 1885 int threadnum=0;
769     #endif
770 jfenwick 1886 ValueType v(numthreads*threadbuffersize);
771 jfenwick 1879 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 jfenwick 1885 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 jfenwick 1879 {
782 jfenwick 1885 #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 jfenwick 1889 for (unsigned int i=0;i<m_samplesize;++i,++resoffset) // copy values into the output vector
789 jfenwick 1879 {
790 jfenwick 1885 resvec[resoffset]=res[i];
791 jfenwick 1879 }
792     }
793     return resptr;
794     }
795    
796 jfenwick 1865 std::string
797     DataLazy::toString() const
798     {
799 jfenwick 1886 ostringstream oss;
800     oss << "Lazy Data:";
801     intoString(oss);
802     return oss.str();
803 jfenwick 1865 }
804    
805 jfenwick 1886 void
806     DataLazy::intoString(ostringstream& oss) const
807     {
808     switch (getOpgroup(m_op))
809     {
810 jfenwick 1889 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 jfenwick 1886 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 jfenwick 1865 // 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 jfenwick 1879 return m_length;
863 jfenwick 1865 }
864    
865    
866     DataAbstract*
867     DataLazy::getSlice(const DataTypes::RegionType& region) const
868     {
869 jfenwick 1879 throw DataException("getSlice - not implemented for Lazy objects.");
870 jfenwick 1865 }
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