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

Diff of /branches/schroedinger_upto1946/escript/src/DataLazy.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.1865  
changed lines
  Added in v.1901

  ViewVC Help
Powered by ViewVC 1.1.26