/[escript]/trunk/escript/src/DataLazy.cpp
ViewVC logotype

Annotation of /trunk/escript/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2050 - (hide annotations)
Mon Nov 17 08:59:57 2008 UTC (10 years, 7 months ago) by phornby
File size: 30938 byte(s)
OK, I have to include this file in the experiment
to use the intelc math library on windows.


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 1950 #include "Utils.h"
30 jfenwick 1865
31 jfenwick 1899 /*
32     How does DataLazy work?
33     ~~~~~~~~~~~~~~~~~~~~~~~
34    
35     Each instance represents a single operation on one or two other DataLazy instances. These arguments are normally
36     denoted left and right.
37    
38     A special operation, IDENTITY, stores an instance of DataReady in the m_id member.
39     This means that all "internal" nodes in the structure are instances of DataLazy.
40    
41     Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
42     Note that IDENITY is not considered a unary operation.
43    
44     I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
45     It must however form a DAG (directed acyclic graph).
46     I will refer to individual DataLazy objects with the structure as nodes.
47    
48     Each node also stores:
49     - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
50     evaluated.
51     - m_length ~ how many values would be stored in the answer if the expression was evaluated.
52     - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to
53     evaluate the expression.
54     - m_samplesize ~ the number of doubles stored in a sample.
55    
56     When a new node is created, the above values are computed based on the values in the child nodes.
57     Eg: if left requires 4 samples and right requires 6 then left+right requires 7 samples.
58    
59     The resolve method, which produces a DataReady from a DataLazy, does the following:
60     1) Create a DataReady to hold the new result.
61     2) Allocate a vector (v) big enough to hold m_buffsrequired samples.
62     3) For each sample, call resolveSample with v, to get its values and copy them into the result object.
63    
64     (In the case of OMP, multiple samples are resolved in parallel so the vector needs to be larger.)
65    
66     resolveSample returns a Vector* and an offset within that vector where the result is stored.
67     Normally, this would be v, but for identity nodes their internal vector is returned instead.
68    
69     The convention that I use, is that the resolve methods should store their results starting at the offset they are passed.
70    
71     For expressions which evaluate to Constant or Tagged, there is a different evaluation method.
72     The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.
73 jfenwick 2037
74     To add a new operator you need to do the following (plus anything I might have forgotten):
75     1) Add to the ES_optype.
76     2) determine what opgroup your operation belongs to (X)
77     3) add a string for the op to the end of ES_opstrings
78     4) increase ES_opcount
79     5) add an entry (X) to opgroups
80     6) add an entry to the switch in collapseToReady
81     7) add an entry to resolveX
82 jfenwick 1899 */
83    
84    
85 jfenwick 1865 using namespace std;
86 jfenwick 1868 using namespace boost;
87 jfenwick 1865
88     namespace escript
89     {
90    
91     namespace
92     {
93    
94 jfenwick 1886 enum ES_opgroup
95     {
96     G_UNKNOWN,
97     G_IDENTITY,
98 jfenwick 1899 G_BINARY, // pointwise operations with two arguments
99 jfenwick 2037 G_UNARY, // pointwise operations with one argument
100     G_NP1OUT // non-pointwise op with one output
101 jfenwick 1886 };
102    
103    
104    
105    
106 jfenwick 1910 string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
107     "sin","cos","tan",
108 jfenwick 1886 "asin","acos","atan","sinh","cosh","tanh","erf",
109     "asinh","acosh","atanh",
110 jfenwick 1888 "log10","log","sign","abs","neg","pos","exp","sqrt",
111 jfenwick 2037 "1/","where>0","where<0","where>=0","where<=0",
112     "symmetric","nonsymmetric"};
113     int ES_opcount=35;
114 jfenwick 1910 ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
115     G_UNARY,G_UNARY,G_UNARY, //10
116     G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 17
117     G_UNARY,G_UNARY,G_UNARY, // 20
118     G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 28
119 jfenwick 2037 G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, // 33
120     G_NP1OUT,G_NP1OUT};
121 jfenwick 1886 inline
122     ES_opgroup
123     getOpgroup(ES_optype op)
124     {
125     return opgroups[op];
126     }
127    
128 jfenwick 1865 // return the FunctionSpace of the result of "left op right"
129     FunctionSpace
130     resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
131     {
132     // perhaps this should call interpolate and throw or something?
133     // maybe we need an interpolate node -
134     // that way, if interpolate is required in any other op we can just throw a
135     // programming error exception.
136 jfenwick 1879
137 jfenwick 1943 FunctionSpace l=left->getFunctionSpace();
138     FunctionSpace r=right->getFunctionSpace();
139     if (l!=r)
140     {
141     if (r.probeInterpolation(l))
142     {
143     return l;
144     }
145     if (l.probeInterpolation(r))
146     {
147     return r;
148     }
149     throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
150     }
151     return l;
152 jfenwick 1865 }
153    
154     // return the shape of the result of "left op right"
155     DataTypes::ShapeType
156     resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
157     {
158 jfenwick 1879 if (left->getShape()!=right->getShape())
159     {
160 jfenwick 2037 if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
161 jfenwick 1908 {
162     throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
163     }
164     if (left->getRank()==0) // we need to allow scalar * anything
165     {
166     return right->getShape();
167     }
168     if (right->getRank()==0)
169     {
170     return left->getShape();
171     }
172     throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
173 jfenwick 1879 }
174     return left->getShape();
175 jfenwick 1865 }
176    
177 jfenwick 1899 // determine the number of points in the result of "left op right"
178 jfenwick 1865 size_t
179     resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
180     {
181 jfenwick 1886 switch (getOpgroup(op))
182 jfenwick 1865 {
183 jfenwick 1886 case G_BINARY: return left->getLength();
184     case G_UNARY: return left->getLength();
185 jfenwick 2037 case G_NP1OUT: return left->getLength();
186 jfenwick 1865 default:
187     throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
188     }
189     }
190    
191 jfenwick 1899 // determine the number of samples requires to evaluate an expression combining left and right
192 jfenwick 2037 // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
193 jfenwick 1879 int
194     calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
195     {
196 jfenwick 1886 switch(getOpgroup(op))
197 jfenwick 1879 {
198 jfenwick 1886 case G_IDENTITY: return 1;
199     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
200     case G_UNARY: return max(left->getBuffsRequired(),1);
201 jfenwick 2037 case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
202 jfenwick 1879 default:
203     throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
204     }
205     }
206    
207 jfenwick 1899
208 jfenwick 1865 } // end anonymous namespace
209    
210    
211 jfenwick 1899
212     // Return a string representing the operation
213 jfenwick 1865 const std::string&
214     opToString(ES_optype op)
215     {
216     if (op<0 || op>=ES_opcount)
217     {
218     op=UNKNOWNOP;
219     }
220     return ES_opstrings[op];
221     }
222    
223    
224     DataLazy::DataLazy(DataAbstract_ptr p)
225     : parent(p->getFunctionSpace(),p->getShape()),
226     m_op(IDENTITY)
227     {
228 jfenwick 1879 if (p->isLazy())
229     {
230     // I don't want identity of Lazy.
231     // Question: Why would that be so bad?
232     // Answer: We assume that the child of ID is something we can call getVector on
233     throw DataException("Programmer error - attempt to create identity from a DataLazy.");
234     }
235     else
236     {
237     m_id=dynamic_pointer_cast<DataReady>(p);
238 jfenwick 1889 if(p->isConstant()) {m_readytype='C';}
239     else if(p->isExpanded()) {m_readytype='E';}
240     else if (p->isTagged()) {m_readytype='T';}
241     else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}
242 jfenwick 1879 }
243     m_length=p->getLength();
244 jfenwick 1886 m_buffsRequired=1;
245 jfenwick 1879 m_samplesize=getNumDPPSample()*getNoValues();
246     cout << "(1)Lazy created with " << m_samplesize << endl;
247 jfenwick 1865 }
248    
249 jfenwick 1899
250    
251 jfenwick 1901
252 jfenwick 1886 DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
253     : parent(left->getFunctionSpace(),left->getShape()),
254     m_op(op)
255     {
256 jfenwick 2037 if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
257 jfenwick 1886 {
258     throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
259     }
260     DataLazy_ptr lleft;
261     if (!left->isLazy())
262     {
263     lleft=DataLazy_ptr(new DataLazy(left));
264     }
265     else
266     {
267     lleft=dynamic_pointer_cast<DataLazy>(left);
268     }
269 jfenwick 1889 m_readytype=lleft->m_readytype;
270 jfenwick 1886 m_length=left->getLength();
271     m_left=lleft;
272 jfenwick 2037 m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
273 jfenwick 1886 m_samplesize=getNumDPPSample()*getNoValues();
274     }
275    
276    
277 jfenwick 1943 // In this constructor we need to consider interpolation
278 jfenwick 1879 DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
279     : parent(resultFS(left,right,op), resultShape(left,right,op)),
280     m_op(op)
281     {
282 jfenwick 2037 if ((getOpgroup(op)!=G_BINARY))
283 jfenwick 1886 {
284 jfenwick 2037 cout << opToString(op) << endl;
285 jfenwick 1889 throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
286 jfenwick 1886 }
287 jfenwick 1943
288     if (getFunctionSpace()!=left->getFunctionSpace()) // left needs to be interpolated
289     {
290     FunctionSpace fs=getFunctionSpace();
291     Data ltemp(left);
292     Data tmp(ltemp,fs);
293     left=tmp.borrowDataPtr();
294     }
295     if (getFunctionSpace()!=right->getFunctionSpace()) // left needs to be interpolated
296     {
297     Data tmp(Data(right),getFunctionSpace());
298     right=tmp.borrowDataPtr();
299     }
300     left->operandCheck(*right);
301    
302 jfenwick 1899 if (left->isLazy()) // the children need to be DataLazy. Wrap them in IDENTITY if required
303 jfenwick 1879 {
304     m_left=dynamic_pointer_cast<DataLazy>(left);
305     }
306     else
307     {
308     m_left=DataLazy_ptr(new DataLazy(left));
309     }
310     if (right->isLazy())
311     {
312     m_right=dynamic_pointer_cast<DataLazy>(right);
313     }
314     else
315     {
316     m_right=DataLazy_ptr(new DataLazy(right));
317     }
318 jfenwick 1889 char lt=m_left->m_readytype;
319     char rt=m_right->m_readytype;
320     if (lt=='E' || rt=='E')
321     {
322     m_readytype='E';
323     }
324     else if (lt=='T' || rt=='T')
325     {
326     m_readytype='T';
327     }
328     else
329     {
330     m_readytype='C';
331     }
332 jfenwick 1879 m_length=resultLength(m_left,m_right,m_op);
333 jfenwick 1899 m_samplesize=getNumDPPSample()*getNoValues();
334 jfenwick 1879 m_buffsRequired=calcBuffs(m_left, m_right,m_op);
335     cout << "(3)Lazy created with " << m_samplesize << endl;
336     }
337    
338    
339 jfenwick 1865 DataLazy::~DataLazy()
340     {
341     }
342    
343 jfenwick 1879
344     int
345     DataLazy::getBuffsRequired() const
346 jfenwick 1865 {
347 jfenwick 1879 return m_buffsRequired;
348     }
349    
350 jfenwick 1884
351 jfenwick 1899 /*
352     \brief Evaluates the expression using methods on Data.
353     This does the work for the collapse method.
354     For reasons of efficiency do not call this method on DataExpanded nodes.
355     */
356 jfenwick 1889 DataReady_ptr
357     DataLazy::collapseToReady()
358     {
359     if (m_readytype=='E')
360     { // this is more an efficiency concern than anything else
361     throw DataException("Programmer Error - do not use collapse on Expanded data.");
362     }
363     if (m_op==IDENTITY)
364     {
365     return m_id;
366     }
367     DataReady_ptr pleft=m_left->collapseToReady();
368     Data left(pleft);
369     Data right;
370     if (getOpgroup(m_op)==G_BINARY)
371     {
372     right=Data(m_right->collapseToReady());
373     }
374     Data result;
375     switch(m_op)
376     {
377     case ADD:
378     result=left+right;
379     break;
380     case SUB:
381     result=left-right;
382     break;
383     case MUL:
384     result=left*right;
385     break;
386     case DIV:
387     result=left/right;
388     break;
389     case SIN:
390     result=left.sin();
391     break;
392     case COS:
393     result=left.cos();
394     break;
395     case TAN:
396     result=left.tan();
397     break;
398     case ASIN:
399     result=left.asin();
400     break;
401     case ACOS:
402     result=left.acos();
403     break;
404     case ATAN:
405     result=left.atan();
406     break;
407     case SINH:
408     result=left.sinh();
409     break;
410     case COSH:
411     result=left.cosh();
412     break;
413     case TANH:
414     result=left.tanh();
415     break;
416     case ERF:
417     result=left.erf();
418     break;
419     case ASINH:
420     result=left.asinh();
421     break;
422     case ACOSH:
423     result=left.acosh();
424     break;
425     case ATANH:
426     result=left.atanh();
427     break;
428     case LOG10:
429     result=left.log10();
430     break;
431     case LOG:
432     result=left.log();
433     break;
434     case SIGN:
435     result=left.sign();
436     break;
437     case ABS:
438     result=left.abs();
439     break;
440     case NEG:
441     result=left.neg();
442     break;
443     case POS:
444     // it doesn't mean anything for delayed.
445     // it will just trigger a deep copy of the lazy object
446     throw DataException("Programmer error - POS not supported for lazy data.");
447     break;
448     case EXP:
449     result=left.exp();
450     break;
451     case SQRT:
452     result=left.sqrt();
453     break;
454     case RECIP:
455     result=left.oneOver();
456     break;
457     case GZ:
458     result=left.wherePositive();
459     break;
460     case LZ:
461     result=left.whereNegative();
462     break;
463     case GEZ:
464     result=left.whereNonNegative();
465     break;
466     case LEZ:
467     result=left.whereNonPositive();
468     break;
469 jfenwick 2037 case SYM:
470     result=left.symmetric();
471     break;
472     case NSYM:
473     result=left.nonsymmetric();
474     break;
475 jfenwick 1889 default:
476 jfenwick 2037 throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
477 jfenwick 1889 }
478     return result.borrowReadyPtr();
479     }
480    
481 jfenwick 1899 /*
482     \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
483     This method uses the original methods on the Data class to evaluate the expressions.
484     For this reason, it should not be used on DataExpanded instances. (To do so would defeat
485     the purpose of using DataLazy in the first place).
486     */
487 jfenwick 1889 void
488     DataLazy::collapse()
489     {
490     if (m_op==IDENTITY)
491     {
492     return;
493     }
494     if (m_readytype=='E')
495     { // this is more an efficiency concern than anything else
496     throw DataException("Programmer Error - do not use collapse on Expanded data.");
497     }
498     m_id=collapseToReady();
499     m_op=IDENTITY;
500     }
501    
502 jfenwick 1899 /*
503 jfenwick 2037 \brief Compute the value of the expression (unary operation) for the given sample.
504 jfenwick 1899 \return Vector which stores the value of the subexpression for the given sample.
505     \param v A vector to store intermediate results.
506     \param offset Index in v to begin storing results.
507     \param sampleNo Sample number to evaluate.
508     \param roffset (output parameter) the offset in the return vector where the result begins.
509    
510     The return value will be an existing vector so do not deallocate it.
511     If the result is stored in v it should be stored at the offset given.
512     Everything from offset to the end of v should be considered available for this method to use.
513     */
514 jfenwick 1898 DataTypes::ValueType*
515     DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
516 jfenwick 1889 {
517     // we assume that any collapsing has been done before we get here
518     // since we only have one argument we don't need to think about only
519     // processing single points.
520     if (m_readytype!='E')
521     {
522     throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
523     }
524 jfenwick 1898 const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
525     const double* left=&((*vleft)[roffset]);
526 jfenwick 1889 double* result=&(v[offset]);
527 jfenwick 1898 roffset=offset;
528 jfenwick 1889 switch (m_op)
529     {
530     case SIN:
531 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
532 jfenwick 1889 break;
533     case COS:
534 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
535 jfenwick 1889 break;
536     case TAN:
537 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
538 jfenwick 1889 break;
539     case ASIN:
540 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
541 jfenwick 1889 break;
542     case ACOS:
543 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
544 jfenwick 1889 break;
545     case ATAN:
546 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
547 jfenwick 1889 break;
548     case SINH:
549 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
550 jfenwick 1889 break;
551     case COSH:
552 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
553 jfenwick 1889 break;
554     case TANH:
555 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
556 jfenwick 1889 break;
557     case ERF:
558 phornby 2050 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
559 jfenwick 1889 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
560     #else
561     tensor_unary_operation(m_samplesize, left, result, ::erf);
562     break;
563     #endif
564     case ASINH:
565 phornby 2050 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
566 jfenwick 1889 tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
567     #else
568     tensor_unary_operation(m_samplesize, left, result, ::asinh);
569     #endif
570     break;
571     case ACOSH:
572 phornby 2050 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
573 jfenwick 1889 tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
574     #else
575     tensor_unary_operation(m_samplesize, left, result, ::acosh);
576     #endif
577     break;
578     case ATANH:
579 phornby 2050 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
580 jfenwick 1889 tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
581     #else
582     tensor_unary_operation(m_samplesize, left, result, ::atanh);
583     #endif
584     break;
585     case LOG10:
586 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
587 jfenwick 1889 break;
588     case LOG:
589 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
590 jfenwick 1889 break;
591     case SIGN:
592     tensor_unary_operation(m_samplesize, left, result, escript::fsign);
593     break;
594     case ABS:
595 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
596 jfenwick 1889 break;
597     case NEG:
598     tensor_unary_operation(m_samplesize, left, result, negate<double>());
599     break;
600     case POS:
601     // it doesn't mean anything for delayed.
602     // it will just trigger a deep copy of the lazy object
603     throw DataException("Programmer error - POS not supported for lazy data.");
604     break;
605     case EXP:
606 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
607 jfenwick 1889 break;
608     case SQRT:
609 jfenwick 1972 tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
610 jfenwick 1889 break;
611     case RECIP:
612     tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
613     break;
614     case GZ:
615     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
616     break;
617     case LZ:
618     tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
619     break;
620     case GEZ:
621     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
622     break;
623     case LEZ:
624     tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
625     break;
626    
627     default:
628     throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
629     }
630 jfenwick 1898 return &v;
631 jfenwick 1889 }
632    
633    
634 jfenwick 2037 /*
635     \brief Compute the value of the expression (unary operation) for the given sample.
636     \return Vector which stores the value of the subexpression for the given sample.
637     \param v A vector to store intermediate results.
638     \param offset Index in v to begin storing results.
639     \param sampleNo Sample number to evaluate.
640     \param roffset (output parameter) the offset in the return vector where the result begins.
641 jfenwick 1898
642 jfenwick 2037 The return value will be an existing vector so do not deallocate it.
643     If the result is stored in v it should be stored at the offset given.
644     Everything from offset to the end of v should be considered available for this method to use.
645     */
646     DataTypes::ValueType*
647     DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
648     {
649     // we assume that any collapsing has been done before we get here
650     // since we only have one argument we don't need to think about only
651     // processing single points.
652     if (m_readytype!='E')
653     {
654     throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
655     }
656     // since we can't write the result over the input, we need a result offset further along
657     size_t subroffset=roffset+m_samplesize;
658     const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);
659     roffset=offset;
660     switch (m_op)
661     {
662     case SYM:
663     DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
664     break;
665     case NSYM:
666     DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
667     break;
668     default:
669     throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
670     }
671     return &v;
672     }
673 jfenwick 1898
674 jfenwick 1899
675 jfenwick 2037
676    
677 phornby 1987 #define PROC_OP(TYPE,X) \
678 jfenwick 1908 for (int i=0;i<steps;++i,resultp+=resultStep) \
679 jfenwick 1889 { \
680 phornby 1994 tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
681 jfenwick 1899 lroffset+=leftStep; \
682     rroffset+=rightStep; \
683 jfenwick 1889 }
684    
685 jfenwick 1899 /*
686     \brief Compute the value of the expression (binary operation) for the given sample.
687     \return Vector which stores the value of the subexpression for the given sample.
688     \param v A vector to store intermediate results.
689     \param offset Index in v to begin storing results.
690     \param sampleNo Sample number to evaluate.
691     \param roffset (output parameter) the offset in the return vector where the result begins.
692    
693     The return value will be an existing vector so do not deallocate it.
694     If the result is stored in v it should be stored at the offset given.
695     Everything from offset to the end of v should be considered available for this method to use.
696     */
697     // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
698     // have already been collapsed to IDENTITY. So we must have at least one expanded child.
699     // If both children are expanded, then we can process them in a single operation (we treat
700     // the whole sample as one big datapoint.
701     // If one of the children is not expanded, then we need to treat each point in the sample
702     // individually.
703 jfenwick 1908 // There is an additional complication when scalar operations are considered.
704     // For example, 2+Vector.
705     // In this case each double within the point is treated individually
706 jfenwick 1898 DataTypes::ValueType*
707 jfenwick 1899 DataLazy::resolveBinary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
708 jfenwick 1889 {
709     cout << "Resolve binary: " << toString() << endl;
710    
711 jfenwick 1899 size_t lroffset=0, rroffset=0; // offsets in the left and right result vectors
712     // first work out which of the children are expanded
713 jfenwick 1889 bool leftExp=(m_left->m_readytype=='E');
714     bool rightExp=(m_right->m_readytype=='E');
715 jfenwick 1899 bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
716 jfenwick 1908 int steps=(bigloops?1:getNumDPPSample());
717 jfenwick 1899 size_t chunksize=(bigloops? m_samplesize : getNoValues()); // if bigloops, pretend the whole sample is a datapoint
718 jfenwick 1908 if (m_left->getRank()!=m_right->getRank()) // need to deal with scalar * ? ops
719     {
720     EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");
721     steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());
722     chunksize=1; // for scalar
723     }
724     int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
725     int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
726     int resultStep=max(leftStep,rightStep); // only one (at most) should be !=0
727 jfenwick 1899 // Get the values of sub-expressions
728 jfenwick 1898 const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
729 jfenwick 1899 const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
730     // the right child starts further along.
731     double* resultp=&(v[offset]); // results are stored at the vector offset we recieved
732 jfenwick 1889 switch(m_op)
733     {
734     case ADD:
735 phornby 2021 PROC_OP(NO_ARG,plus<double>());
736 jfenwick 1889 break;
737 jfenwick 1899 case SUB:
738 phornby 2021 PROC_OP(NO_ARG,minus<double>());
739 jfenwick 1899 break;
740     case MUL:
741 phornby 2021 PROC_OP(NO_ARG,multiplies<double>());
742 jfenwick 1899 break;
743     case DIV:
744 phornby 2021 PROC_OP(NO_ARG,divides<double>());
745 jfenwick 1899 break;
746 jfenwick 1910 case POW:
747 phornby 1994 PROC_OP(double (double,double),::pow);
748 jfenwick 1910 break;
749 jfenwick 1889 default:
750 jfenwick 1898 throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
751 jfenwick 1889 }
752 jfenwick 1899 roffset=offset;
753 jfenwick 1898 return &v;
754 jfenwick 1889 }
755    
756 jfenwick 1898
757    
758 jfenwick 1899 /*
759     \brief Compute the value of the expression for the given sample.
760     \return Vector which stores the value of the subexpression for the given sample.
761     \param v A vector to store intermediate results.
762     \param offset Index in v to begin storing results.
763     \param sampleNo Sample number to evaluate.
764     \param roffset (output parameter) the offset in the return vector where the result begins.
765 jfenwick 1898
766 jfenwick 1899 The return value will be an existing vector so do not deallocate it.
767     */
768 jfenwick 1884 // the vector and the offset are a place where the method could write its data if it wishes
769     // it is not obligated to do so. For example, if it has its own storage already, it can use that.
770 jfenwick 1885 // Hence the return value to indicate where the data is actually stored.
771     // Regardless, the storage should be assumed to be used, even if it isn't.
772 jfenwick 1898
773     // the roffset is the offset within the returned vector where the data begins
774     const DataTypes::ValueType*
775     DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
776 jfenwick 1879 {
777 jfenwick 1889 cout << "Resolve sample " << toString() << endl;
778     // collapse so we have a 'E' node or an IDENTITY for some other type
779     if (m_readytype!='E' && m_op!=IDENTITY)
780     {
781     collapse();
782     }
783 jfenwick 1885 if (m_op==IDENTITY)
784 jfenwick 1865 {
785 jfenwick 1879 const ValueType& vec=m_id->getVector();
786 jfenwick 1889 if (m_readytype=='C')
787     {
788 jfenwick 1898 roffset=0;
789     return &(vec);
790 jfenwick 1889 }
791 jfenwick 1898 roffset=m_id->getPointOffset(sampleNo, 0);
792     return &(vec);
793 jfenwick 1865 }
794 jfenwick 1889 if (m_readytype!='E')
795     {
796     throw DataException("Programmer Error - Collapse did not produce an expanded node.");
797     }
798     switch (getOpgroup(m_op))
799     {
800 jfenwick 1898 case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);
801     case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
802 jfenwick 1889 default:
803     throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
804     }
805     }
806    
807    
808 jfenwick 1899 // To simplify the memory management, all threads operate on one large vector, rather than one each.
809     // Each sample is evaluated independently and copied into the result DataExpanded.
810 jfenwick 1879 DataReady_ptr
811     DataLazy::resolve()
812     {
813    
814     cout << "Sample size=" << m_samplesize << endl;
815     cout << "Buffers=" << m_buffsRequired << endl;
816    
817 jfenwick 1899 if (m_readytype!='E') // if the whole sub-expression is Constant or Tagged, then evaluate it normally
818 jfenwick 1889 {
819     collapse();
820     }
821 jfenwick 1899 if (m_op==IDENTITY) // So a lazy expression of Constant or Tagged data will be returned here.
822 jfenwick 1889 {
823     return m_id;
824     }
825     // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
826 jfenwick 1899 size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough
827     // storage to evaluate its expression
828 jfenwick 1885 int numthreads=1;
829     #ifdef _OPENMP
830 jfenwick 1889 numthreads=getNumberOfThreads();
831 jfenwick 1885 int threadnum=0;
832     #endif
833 jfenwick 1886 ValueType v(numthreads*threadbuffersize);
834 jfenwick 1879 cout << "Buffer created with size=" << v.size() << endl;
835 jfenwick 1898 DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(), ValueType(getNoValues()));
836 jfenwick 1879 ValueType& resvec=result->getVector();
837     DataReady_ptr resptr=DataReady_ptr(result);
838     int sample;
839 jfenwick 1898 size_t outoffset; // offset in the output data
840 jfenwick 1885 int totalsamples=getNumSamples();
841 jfenwick 1899 const ValueType* res=0; // Vector storing the answer
842     size_t resoffset=0; // where in the vector to find the answer
843 jfenwick 1898 #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)
844 jfenwick 1885 for (sample=0;sample<totalsamples;++sample)
845 jfenwick 1879 {
846 jfenwick 1898 cout << "################################# " << sample << endl;
847 jfenwick 1885 #ifdef _OPENMP
848 jfenwick 1898 res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
849 jfenwick 1885 #else
850 jfenwick 1899 res=resolveSample(v,0,sample,resoffset); // res would normally be v, but not if its a single IDENTITY op.
851 jfenwick 1885 #endif
852 jfenwick 1898 cerr << "-------------------------------- " << endl;
853     outoffset=result->getPointOffset(sample,0);
854     cerr << "offset=" << outoffset << endl;
855     for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset) // copy values into the output vector
856 jfenwick 1879 {
857 jfenwick 1898 resvec[outoffset]=(*res)[resoffset];
858 jfenwick 1879 }
859 jfenwick 1898 cerr << "*********************************" << endl;
860 jfenwick 1879 }
861     return resptr;
862     }
863    
864 jfenwick 1865 std::string
865     DataLazy::toString() const
866     {
867 jfenwick 1886 ostringstream oss;
868     oss << "Lazy Data:";
869     intoString(oss);
870     return oss.str();
871 jfenwick 1865 }
872    
873 jfenwick 1899
874 jfenwick 1886 void
875     DataLazy::intoString(ostringstream& oss) const
876     {
877     switch (getOpgroup(m_op))
878     {
879 jfenwick 1889 case G_IDENTITY:
880     if (m_id->isExpanded())
881     {
882     oss << "E";
883     }
884     else if (m_id->isTagged())
885     {
886     oss << "T";
887     }
888     else if (m_id->isConstant())
889     {
890     oss << "C";
891     }
892     else
893     {
894     oss << "?";
895     }
896 jfenwick 1886 oss << '@' << m_id.get();
897     break;
898     case G_BINARY:
899     oss << '(';
900     m_left->intoString(oss);
901     oss << ' ' << opToString(m_op) << ' ';
902     m_right->intoString(oss);
903     oss << ')';
904     break;
905     case G_UNARY:
906 jfenwick 2037 case G_NP1OUT:
907 jfenwick 1886 oss << opToString(m_op) << '(';
908     m_left->intoString(oss);
909     oss << ')';
910     break;
911     default:
912     oss << "UNKNOWN";
913     }
914     }
915    
916 jfenwick 1865 DataAbstract*
917     DataLazy::deepCopy()
918     {
919 jfenwick 1935 switch (getOpgroup(m_op))
920 jfenwick 1865 {
921 jfenwick 1935 case G_IDENTITY: return new DataLazy(m_id->deepCopy()->getPtr());
922     case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
923     case G_BINARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
924     default:
925     throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
926 jfenwick 1865 }
927     }
928    
929    
930     DataTypes::ValueType::size_type
931     DataLazy::getLength() const
932     {
933 jfenwick 1879 return m_length;
934 jfenwick 1865 }
935    
936    
937     DataAbstract*
938     DataLazy::getSlice(const DataTypes::RegionType& region) const
939     {
940 jfenwick 1879 throw DataException("getSlice - not implemented for Lazy objects.");
941 jfenwick 1865 }
942    
943 jfenwick 1935
944     // To do this we need to rely on our child nodes
945 jfenwick 1865 DataTypes::ValueType::size_type
946     DataLazy::getPointOffset(int sampleNo,
947 jfenwick 1935 int dataPointNo)
948     {
949     if (m_op==IDENTITY)
950     {
951     return m_id->getPointOffset(sampleNo,dataPointNo);
952     }
953     if (m_readytype!='E')
954     {
955     collapse();
956     return m_id->getPointOffset(sampleNo,dataPointNo);
957     }
958     // at this point we do not have an identity node and the expression will be Expanded
959     // so we only need to know which child to ask
960     if (m_left->m_readytype=='E')
961     {
962     return m_left->getPointOffset(sampleNo,dataPointNo);
963     }
964     else
965     {
966     return m_right->getPointOffset(sampleNo,dataPointNo);
967     }
968     }
969    
970     // To do this we need to rely on our child nodes
971     DataTypes::ValueType::size_type
972     DataLazy::getPointOffset(int sampleNo,
973 jfenwick 1865 int dataPointNo) const
974     {
975 jfenwick 1935 if (m_op==IDENTITY)
976     {
977     return m_id->getPointOffset(sampleNo,dataPointNo);
978     }
979     if (m_readytype=='E')
980     {
981     // at this point we do not have an identity node and the expression will be Expanded
982     // so we only need to know which child to ask
983     if (m_left->m_readytype=='E')
984     {
985     return m_left->getPointOffset(sampleNo,dataPointNo);
986     }
987     else
988     {
989     return m_right->getPointOffset(sampleNo,dataPointNo);
990     }
991     }
992     if (m_readytype=='C')
993     {
994     return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
995     }
996     throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
997 jfenwick 1865 }
998    
999 jfenwick 1901 // It would seem that DataTagged will need to be treated differently since even after setting all tags
1000     // to zero, all the tags from all the DataTags would be in the result.
1001     // However since they all have the same value (0) whether they are there or not should not matter.
1002     // So I have decided that for all types this method will create a constant 0.
1003     // It can be promoted up as required.
1004     // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
1005     // but we can deal with that if it arrises.
1006     void
1007     DataLazy::setToZero()
1008     {
1009     DataTypes::ValueType v(getNoValues(),0);
1010     m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1011     m_op=IDENTITY;
1012     m_right.reset();
1013     m_left.reset();
1014     m_readytype='C';
1015     m_buffsRequired=1;
1016     }
1017    
1018 jfenwick 1865 } // end namespace

  ViewVC Help
Powered by ViewVC 1.1.26