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

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

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

revision 1888 by jfenwick, Wed Oct 15 04:00:21 2008 UTC revision 1943 by jfenwick, Wed Oct 29 04:05:14 2008 UTC
# Line 27  Line 27 
27  #include "Data.h"  #include "Data.h"
28  #include "UnaryFuncs.h"     // for escript::fsign  #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;  using namespace boost;
77    
# Line 39  opToString(ES_optype op); Line 84  opToString(ES_optype op);
84  namespace  namespace
85  {  {
86    
   
   
87  enum ES_opgroup  enum ES_opgroup
88  {  {
89     G_UNKNOWN,     G_UNKNOWN,
90     G_IDENTITY,     G_IDENTITY,
91     G_BINARY,     G_BINARY,        // pointwise operations with two arguments
92     G_UNARY     G_UNARY      // pointwise operations with one argument
93  };  };
94    
95    
96    
97    
98  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","sin","cos","tan",  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
99                "sin","cos","tan",
100              "asin","acos","atan","sinh","cosh","tanh","erf",              "asin","acos","atan","sinh","cosh","tanh","erf",
101              "asinh","acosh","atanh",              "asinh","acosh","atanh",
102              "log10","log","sign","abs","neg","pos","exp","sqrt",              "log10","log","sign","abs","neg","pos","exp","sqrt",
103              "1/","where>0","where<0","where>=0","where<=0"};              "1/","where>0","where<0","where>=0","where<=0"};
104  int ES_opcount=32;  int ES_opcount=33;
105  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY,G_UNARY,G_UNARY,G_UNARY, //9  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
106              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 16              G_UNARY,G_UNARY,G_UNARY, //10
107              G_UNARY,G_UNARY,G_UNARY,                    // 19              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17
108              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 27              G_UNARY,G_UNARY,G_UNARY,                    // 20
109                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 28
110              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};
111  inline  inline
112  ES_opgroup  ES_opgroup
# Line 79  resultFS(DataAbstract_ptr left, DataAbst Line 124  resultFS(DataAbstract_ptr left, DataAbst
124      // 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
125      // programming error exception.      // programming error exception.
126    
127      FunctionSpace l=left->getFunctionSpace();
128      if (left->getFunctionSpace()!=right->getFunctionSpace())    FunctionSpace r=right->getFunctionSpace();
129      {    if (l!=r)
130          throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");    {
131      }      if (r.probeInterpolation(l))
132      return left->getFunctionSpace();      {
133        return l;
134        }
135        if (l.probeInterpolation(r))
136        {
137        return r;
138        }
139        throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
140      }
141      return l;
142  }  }
143    
144  // return the shape of the result of "left op right"  // return the shape of the result of "left op right"
# Line 93  resultShape(DataAbstract_ptr left, DataA Line 147  resultShape(DataAbstract_ptr left, DataA
147  {  {
148      if (left->getShape()!=right->getShape())      if (left->getShape()!=right->getShape())
149      {      {
150          throw DataException("Shapes not the same - shapes must match for lazy data.");        if (getOpgroup(op)!=G_BINARY)
151          {
152            throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
153          }
154          if (left->getRank()==0)   // we need to allow scalar * anything
155          {
156            return right->getShape();
157          }
158          if (right->getRank()==0)
159          {
160            return left->getShape();
161          }
162          throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
163      }      }
164      return left->getShape();      return left->getShape();
165  }  }
166    
167    // determine the number of points in the result of "left op right"
168  size_t  size_t
169  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
170  {  {
# Line 110  resultLength(DataAbstract_ptr left, Data Line 177  resultLength(DataAbstract_ptr left, Data
177     }     }
178  }  }
179    
180    // determine the number of samples requires to evaluate an expression combining left and right
181  int  int
182  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
183  {  {
# Line 124  calcBuffs(const DataLazy_ptr& left, cons Line 192  calcBuffs(const DataLazy_ptr& left, cons
192  }  }
193    
194    
   
195  }   // end anonymous namespace  }   // end anonymous namespace
196    
197    
198    
199    // Return a string representing the operation
200  const std::string&  const std::string&
201  opToString(ES_optype op)  opToString(ES_optype op)
202  {  {
# Line 145  DataLazy::DataLazy(DataAbstract_ptr p) Line 214  DataLazy::DataLazy(DataAbstract_ptr p)
214  {  {
215     if (p->isLazy())     if (p->isLazy())
216     {     {
     // TODO: fix this.   We could make the new node a copy of p?  
217      // I don't want identity of Lazy.      // I don't want identity of Lazy.
218      // Question: Why would that be so bad?      // Question: Why would that be so bad?
219      // Answer: We assume that the child of ID is something we can call getVector on      // Answer: We assume that the child of ID is something we can call getVector on
# Line 154  DataLazy::DataLazy(DataAbstract_ptr p) Line 222  DataLazy::DataLazy(DataAbstract_ptr p)
222     else     else
223     {     {
224      m_id=dynamic_pointer_cast<DataReady>(p);      m_id=dynamic_pointer_cast<DataReady>(p);
225        if(p->isConstant()) {m_readytype='C';}
226        else if(p->isExpanded()) {m_readytype='E';}
227        else if (p->isTagged()) {m_readytype='T';}
228        else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}
229     }     }
230     m_length=p->getLength();     m_length=p->getLength();
231     m_buffsRequired=1;     m_buffsRequired=1;
# Line 161  DataLazy::DataLazy(DataAbstract_ptr p) Line 233  DataLazy::DataLazy(DataAbstract_ptr p)
233  cout << "(1)Lazy created with " << m_samplesize << endl;  cout << "(1)Lazy created with " << m_samplesize << endl;
234  }  }
235    
236    
237    
238    
239  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
240      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
241      m_op(op)      m_op(op)
# Line 178  DataLazy::DataLazy(DataAbstract_ptr left Line 253  DataLazy::DataLazy(DataAbstract_ptr left
253     {     {
254      lleft=dynamic_pointer_cast<DataLazy>(left);      lleft=dynamic_pointer_cast<DataLazy>(left);
255     }     }
256       m_readytype=lleft->m_readytype;
257     m_length=left->getLength();     m_length=left->getLength();
258     m_left=lleft;     m_left=lleft;
259     m_buffsRequired=1;     m_buffsRequired=1;
# Line 185  DataLazy::DataLazy(DataAbstract_ptr left Line 261  DataLazy::DataLazy(DataAbstract_ptr left
261  }  }
262    
263    
264  DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)  // In this constructor we need to consider interpolation
265    DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
266      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
     m_left(left),  
     m_right(right),  
267      m_op(op)      m_op(op)
268  {  {
269     if (getOpgroup(op)!=G_BINARY)     if (getOpgroup(op)!=G_BINARY)
270     {     {
271      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
272     }     }
    m_length=resultLength(m_left,m_right,m_op);  
    m_samplesize=getNumDPPSample()*getNoValues();  
    m_buffsRequired=calcBuffs(m_left, m_right, m_op);  
 cout << "(2)Lazy created with " << m_samplesize << endl;  
 }  
273    
274  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)     if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
     : parent(resultFS(left,right,op), resultShape(left,right,op)),  
     m_op(op)  
 {  
    if (getOpgroup(op)!=G_BINARY)  
275     {     {
276      throw DataException("Programmer error - constructor DataLazy(left, op) will only process BINARY operations.");      FunctionSpace fs=getFunctionSpace();
277        Data ltemp(left);
278        Data tmp(ltemp,fs);
279        left=tmp.borrowDataPtr();
280     }     }
281     if (left->isLazy())     if (getFunctionSpace()!=right->getFunctionSpace())   // left needs to be interpolated
282       {
283        Data tmp(Data(right),getFunctionSpace());
284        right=tmp.borrowDataPtr();
285       }
286       left->operandCheck(*right);
287    
288       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
289     {     {
290      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
291     }     }
# Line 225  DataLazy::DataLazy(DataAbstract_ptr left Line 301  DataLazy::DataLazy(DataAbstract_ptr left
301     {     {
302      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
303     }     }
304       char lt=m_left->m_readytype;
305       char rt=m_right->m_readytype;
306       if (lt=='E' || rt=='E')
307       {
308        m_readytype='E';
309       }
310       else if (lt=='T' || rt=='T')
311       {
312        m_readytype='T';
313       }
314       else
315       {
316        m_readytype='C';
317       }
318     m_length=resultLength(m_left,m_right,m_op);     m_length=resultLength(m_left,m_right,m_op);
319     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();    
320     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
321  cout << "(3)Lazy created with " << m_samplesize << endl;  cout << "(3)Lazy created with " << m_samplesize << endl;
322  }  }
# Line 245  DataLazy::getBuffsRequired() const Line 334  DataLazy::getBuffsRequired() const
334  }  }
335    
336    
337  // the vector and the offset are a place where the method could write its data if it wishes  /*
338  // it is not obligated to do so. For example, if it has its own storage already, it can use that.    \brief Evaluates the expression using methods on Data.
339  // Hence the return value to indicate where the data is actually stored.    This does the work for the collapse method.
340  // Regardless, the storage should be assumed to be used, even if it isn't.    For reasons of efficiency do not call this method on DataExpanded nodes.
341  const double*  */
342  DataLazy::resolveSample(ValueType& v,int sampleNo,  size_t offset ) const  DataReady_ptr
343    DataLazy::collapseToReady()
344  {  {
345    if (m_op==IDENTITY)      if (m_readytype=='E')
346      { // this is more an efficiency concern than anything else
347        throw DataException("Programmer Error - do not use collapse on Expanded data.");
348      }
349      if (m_op==IDENTITY)
350    {    {
351      const ValueType& vec=m_id->getVector();      return m_id;
     return &(vec[m_id->getPointOffset(sampleNo, 0)]);  
352    }    }
353    size_t rightoffset=offset+m_samplesize;    DataReady_ptr pleft=m_left->collapseToReady();
354    const double* left=m_left->resolveSample(v,sampleNo,offset);    Data left(pleft);
355    const double* right=0;    Data right;
356    if (getOpgroup(m_op)==G_BINARY)    if (getOpgroup(m_op)==G_BINARY)
357    {    {
358      right=m_right->resolveSample(v,sampleNo,rightoffset);      right=Data(m_right->collapseToReady());
359    }    }
360    double* result=&(v[offset]);    Data result;
361      switch(m_op)
362    {    {
363      switch(m_op)      case ADD:
364      {      result=left+right;
     case ADD:       // since these are pointwise ops, pretend each sample is one point  
     tensor_binary_operation(m_samplesize, left, right, result, plus<double>());  
365      break;      break;
366      case SUB:            case SUB:      
367      tensor_binary_operation(m_samplesize, left, right, result, minus<double>());      result=left-right;
368      break;      break;
369      case MUL:            case MUL:      
370      tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());      result=left*right;
371      break;      break;
372      case DIV:            case DIV:      
373      tensor_binary_operation(m_samplesize, left, right, result, divides<double>());      result=left/right;
374      break;      break;
 // unary ops  
375      case SIN:      case SIN:
376        result=left.sin();  
377        break;
378        case COS:
379        result=left.cos();
380        break;
381        case TAN:
382        result=left.tan();
383        break;
384        case ASIN:
385        result=left.asin();
386        break;
387        case ACOS:
388        result=left.acos();
389        break;
390        case ATAN:
391        result=left.atan();
392        break;
393        case SINH:
394        result=left.sinh();
395        break;
396        case COSH:
397        result=left.cosh();
398        break;
399        case TANH:
400        result=left.tanh();
401        break;
402        case ERF:
403        result=left.erf();
404        break;
405       case ASINH:
406        result=left.asinh();
407        break;
408       case ACOSH:
409        result=left.acosh();
410        break;
411       case ATANH:
412        result=left.atanh();
413        break;
414        case LOG10:
415        result=left.log10();
416        break;
417        case LOG:
418        result=left.log();
419        break;
420        case SIGN:
421        result=left.sign();
422        break;
423        case ABS:
424        result=left.abs();
425        break;
426        case NEG:
427        result=left.neg();
428        break;
429        case POS:
430        // it doesn't mean anything for delayed.
431        // it will just trigger a deep copy of the lazy object
432        throw DataException("Programmer error - POS not supported for lazy data.");
433        break;
434        case EXP:
435        result=left.exp();
436        break;
437        case SQRT:
438        result=left.sqrt();
439        break;
440        case RECIP:
441        result=left.oneOver();
442        break;
443        case GZ:
444        result=left.wherePositive();
445        break;
446        case LZ:
447        result=left.whereNegative();
448        break;
449        case GEZ:
450        result=left.whereNonNegative();
451        break;
452        case LEZ:
453        result=left.whereNonPositive();
454        break;
455        default:
456        throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
457      }
458      return result.borrowReadyPtr();
459    }
460    
461    /*
462       \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
463       This method uses the original methods on the Data class to evaluate the expressions.
464       For this reason, it should not be used on DataExpanded instances. (To do so would defeat
465       the purpose of using DataLazy in the first place).
466    */
467    void
468    DataLazy::collapse()
469    {
470      if (m_op==IDENTITY)
471      {
472        return;
473      }
474      if (m_readytype=='E')
475      { // this is more an efficiency concern than anything else
476        throw DataException("Programmer Error - do not use collapse on Expanded data.");
477      }
478      m_id=collapseToReady();
479      m_op=IDENTITY;
480    }
481    
482    /*
483      \brief Compute the value of the expression (binary operation) for the given sample.
484      \return Vector which stores the value of the subexpression for the given sample.
485      \param v A vector to store intermediate results.
486      \param offset Index in v to begin storing results.
487      \param sampleNo Sample number to evaluate.
488      \param roffset (output parameter) the offset in the return vector where the result begins.
489    
490      The return value will be an existing vector so do not deallocate it.
491      If the result is stored in v it should be stored at the offset given.
492      Everything from offset to the end of v should be considered available for this method to use.
493    */
494    DataTypes::ValueType*
495    DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
496    {
497        // we assume that any collapsing has been done before we get here
498        // since we only have one argument we don't need to think about only
499        // processing single points.
500      if (m_readytype!='E')
501      {
502        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
503      }
504      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
505      const double* left=&((*vleft)[roffset]);
506      double* result=&(v[offset]);
507      roffset=offset;
508      switch (m_op)
509      {
510        case SIN:  
511      tensor_unary_operation(m_samplesize, left, result, ::sin);      tensor_unary_operation(m_samplesize, left, result, ::sin);
512      break;      break;
513      case COS:      case COS:
# Line 379  DataLazy::resolveSample(ValueType& v,int Line 605  DataLazy::resolveSample(ValueType& v,int
605      break;      break;
606    
607      default:      default:
608      throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
609      }
610      return &v;
611    }
612    
613    
614    
615    
616    
617    #define PROC_OP(X) \
618        for (int i=0;i<steps;++i,resultp+=resultStep) \
619        { \
620           tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
621           lroffset+=leftStep; \
622           rroffset+=rightStep; \
623        }
624    
625    /*
626      \brief Compute the value of the expression (binary operation) for the given sample.
627      \return Vector which stores the value of the subexpression for the given sample.
628      \param v A vector to store intermediate results.
629      \param offset Index in v to begin storing results.
630      \param sampleNo Sample number to evaluate.
631      \param roffset (output parameter) the offset in the return vector where the result begins.
632    
633      The return value will be an existing vector so do not deallocate it.
634      If the result is stored in v it should be stored at the offset given.
635      Everything from offset to the end of v should be considered available for this method to use.
636    */
637    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
638    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
639    // If both children are expanded, then we can process them in a single operation (we treat
640    // the whole sample as one big datapoint.
641    // If one of the children is not expanded, then we need to treat each point in the sample
642    // individually.
643    // There is an additional complication when scalar operations are considered.
644    // For example, 2+Vector.
645    // In this case each double within the point is treated individually
646    DataTypes::ValueType*
647    DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
648    {
649    cout << "Resolve binary: " << toString() << endl;
650    
651      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
652        // first work out which of the children are expanded
653      bool leftExp=(m_left->m_readytype=='E');
654      bool rightExp=(m_right->m_readytype=='E');
655      bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
656      int steps=(bigloops?1:getNumDPPSample());
657      size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint
658      if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops
659      {
660        EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");
661        steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());
662        chunksize=1;    // for scalar
663      }    
664      int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
665      int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
666      int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0
667        // Get the values of sub-expressions
668      const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
669      const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
670        // the right child starts further along.
671      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
672      switch(m_op)
673      {
674        case ADD:
675        PROC_OP(plus<double>());
676        break;
677        case SUB:
678        PROC_OP(minus<double>());
679        break;
680        case MUL:
681        PROC_OP(multiplies<double>());
682        break;
683        case DIV:
684        PROC_OP(divides<double>());
685        break;
686        case POW:
687        PROC_OP(::pow);
688        break;
689        default:
690        throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
691      }
692      roffset=offset;  
693      return &v;
694    }
695    
696    
697    
698    /*
699      \brief Compute the value of the expression for the given sample.
700      \return Vector which stores the value of the subexpression for the given sample.
701      \param v A vector to store intermediate results.
702      \param offset Index in v to begin storing results.
703      \param sampleNo Sample number to evaluate.
704      \param roffset (output parameter) the offset in the return vector where the result begins.
705    
706      The return value will be an existing vector so do not deallocate it.
707    */
708    // the vector and the offset are a place where the method could write its data if it wishes
709    // it is not obligated to do so. For example, if it has its own storage already, it can use that.
710    // Hence the return value to indicate where the data is actually stored.
711    // Regardless, the storage should be assumed to be used, even if it isn't.
712    
713    // the roffset is the offset within the returned vector where the data begins
714    const DataTypes::ValueType*
715    DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
716    {
717    cout << "Resolve sample " << toString() << endl;
718        // collapse so we have a 'E' node or an IDENTITY for some other type
719      if (m_readytype!='E' && m_op!=IDENTITY)
720      {
721        collapse();
722      }
723      if (m_op==IDENTITY)  
724      {
725        const ValueType& vec=m_id->getVector();
726        if (m_readytype=='C')
727        {
728        roffset=0;
729        return &(vec);
730      }      }
731        roffset=m_id->getPointOffset(sampleNo, 0);
732        return &(vec);
733      }
734      if (m_readytype!='E')
735      {
736        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
737      }
738      switch (getOpgroup(m_op))
739      {
740      case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);
741      case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
742      default:
743        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
744    }    }
   return result;  
745  }  }
746    
747    
748    // To simplify the memory management, all threads operate on one large vector, rather than one each.
749    // Each sample is evaluated independently and copied into the result DataExpanded.
750  DataReady_ptr  DataReady_ptr
751  DataLazy::resolve()  DataLazy::resolve()
752  {  {
   // This is broken!     We need to have a buffer per thread!  
   // so the allocation of v needs to move inside the loop somehow  
753    
754  cout << "Sample size=" << m_samplesize << endl;  cout << "Sample size=" << m_samplesize << endl;
755  cout << "Buffers=" << m_buffsRequired << endl;  cout << "Buffers=" << m_buffsRequired << endl;
756    
757    size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);    if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally
758      {
759        collapse();
760      }
761      if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
762      {
763        return m_id;
764      }
765        // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
766      size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired));    // Each thread needs to have enough
767        // storage to evaluate its expression
768    int numthreads=1;    int numthreads=1;
769  #ifdef _OPENMP  #ifdef _OPENMP
770    numthreads=omp_get_max_threads();    numthreads=getNumberOfThreads();
771    int threadnum=0;    int threadnum=0;
772  #endif  #endif
773    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
774  cout << "Buffer created with size=" << v.size() << endl;  cout << "Buffer created with size=" << v.size() << endl;
775    ValueType dummy(getNoValues());    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
   DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),dummy);  
776    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVector();
777    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
778    int sample;    int sample;
779    int resoffset;    size_t outoffset;     // offset in the output data
780    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
781    #pragma omp parallel for private(sample,resoffset,threadnum) schedule(static)    const ValueType* res=0;   // Vector storing the answer
782      size_t resoffset=0;       // where in the vector to find the answer
783      #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)
784    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
785    {    {
786    cout << "################################# " << sample << endl;
787  #ifdef _OPENMP  #ifdef _OPENMP
788      const double* res=resolveSample(v,sample,threadbuffersize*omp_get_thread_num());      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
789  #else  #else
790      const double* res=resolveSample(v,sample,0);   // this would normally be v, but not if its a single IDENTITY op.      res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.
791  #endif  #endif
792      resoffset=result->getPointOffset(sample,0);  cerr << "-------------------------------- " << endl;
793      for (int i=0;i<m_samplesize;++i,++resoffset)    // copy values into the output vector      outoffset=result->getPointOffset(sample,0);
794    cerr << "offset=" << outoffset << endl;
795        for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector
796      {      {
797      resvec[resoffset]=res[i];      resvec[outoffset]=(*res)[resoffset];
798      }      }
799    cerr << "*********************************" << endl;
800    }    }
801    return resptr;    return resptr;
802  }  }
# Line 435  DataLazy::toString() const Line 810  DataLazy::toString() const
810    return oss.str();    return oss.str();
811  }  }
812    
813    
814  void  void
815  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
816  {  {
817    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
818    {    {
819    case G_IDENTITY:    case G_IDENTITY:
820        if (m_id->isExpanded())
821        {
822           oss << "E";
823        }
824        else if (m_id->isTagged())
825        {
826          oss << "T";
827        }
828        else if (m_id->isConstant())
829        {
830          oss << "C";
831        }
832        else
833        {
834          oss << "?";
835        }
836      oss << '@' << m_id.get();      oss << '@' << m_id.get();
837      break;      break;
838    case G_BINARY:    case G_BINARY:
# Line 460  DataLazy::intoString(ostringstream& oss) Line 852  DataLazy::intoString(ostringstream& oss)
852    }    }
853  }  }
854    
 // Note that in this case, deepCopy does not make copies of the leaves.  
 // Hopefully copy on write (or whatever we end up using) will take care of this.  
855  DataAbstract*  DataAbstract*
856  DataLazy::deepCopy()  DataLazy::deepCopy()
857  {  {
858    if (m_op==IDENTITY)    switch (getOpgroup(m_op))
859    {    {
860      return new DataLazy(m_left);    // we don't need to copy the child here    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
861      case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
862      case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
863      default:
864        throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
865    }    }
   return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);  
866  }  }
867    
868    
# Line 486  DataLazy::getSlice(const DataTypes::Regi Line 879  DataLazy::getSlice(const DataTypes::Regi
879    throw DataException("getSlice - not implemented for Lazy objects.");    throw DataException("getSlice - not implemented for Lazy objects.");
880  }  }
881    
882    
883    // To do this we need to rely on our child nodes
884    DataTypes::ValueType::size_type
885    DataLazy::getPointOffset(int sampleNo,
886                     int dataPointNo)
887    {
888      if (m_op==IDENTITY)
889      {
890        return m_id->getPointOffset(sampleNo,dataPointNo);
891      }
892      if (m_readytype!='E')
893      {
894        collapse();
895        return m_id->getPointOffset(sampleNo,dataPointNo);
896      }
897      // at this point we do not have an identity node and the expression will be Expanded
898      // so we only need to know which child to ask
899      if (m_left->m_readytype=='E')
900      {
901        return m_left->getPointOffset(sampleNo,dataPointNo);
902      }
903      else
904      {
905        return m_right->getPointOffset(sampleNo,dataPointNo);
906      }
907    }
908    
909    // To do this we need to rely on our child nodes
910  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
911  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
912                   int dataPointNo) const                   int dataPointNo) const
913  {  {
914    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");    if (m_op==IDENTITY)
915      {
916        return m_id->getPointOffset(sampleNo,dataPointNo);
917      }
918      if (m_readytype=='E')
919      {
920        // at this point we do not have an identity node and the expression will be Expanded
921        // so we only need to know which child to ask
922        if (m_left->m_readytype=='E')
923        {
924        return m_left->getPointOffset(sampleNo,dataPointNo);
925        }
926        else
927        {
928        return m_right->getPointOffset(sampleNo,dataPointNo);
929        }
930      }
931      if (m_readytype=='C')
932      {
933        return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
934      }
935      throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
936    }
937    
938    // It would seem that DataTagged will need to be treated differently since even after setting all tags
939    // to zero, all the tags from all the DataTags would be in the result.
940    // However since they all have the same value (0) whether they are there or not should not matter.
941    // So I have decided that for all types this method will create a constant 0.
942    // It can be promoted up as required.
943    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
944    // but we can deal with that if it arrises.
945    void
946    DataLazy::setToZero()
947    {
948      DataTypes::ValueType v(getNoValues(),0);
949      m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
950      m_op=IDENTITY;
951      m_right.reset();  
952      m_left.reset();
953      m_readytype='C';
954      m_buffsRequired=1;
955  }  }
956    
957  }   // end namespace  }   // end namespace

Legend:
Removed from v.1888  
changed lines
  Added in v.1943

  ViewVC Help
Powered by ViewVC 1.1.26