/[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 1910 by jfenwick, Thu Oct 23 03:05:28 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=32;
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 93  resultShape(DataAbstract_ptr left, DataA Line 138  resultShape(DataAbstract_ptr left, DataA
138  {  {
139      if (left->getShape()!=right->getShape())      if (left->getShape()!=right->getShape())
140      {      {
141          throw DataException("Shapes not the same - shapes must match for lazy data.");        if (getOpgroup(op)!=G_BINARY)
142          {
143            throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
144          }
145          if (left->getRank()==0)   // we need to allow scalar * anything
146          {
147            return right->getShape();
148          }
149          if (right->getRank()==0)
150          {
151            return left->getShape();
152          }
153          throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
154      }      }
155      return left->getShape();      return left->getShape();
156  }  }
157    
158    // determine the number of points in the result of "left op right"
159  size_t  size_t
160  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
161  {  {
# Line 110  resultLength(DataAbstract_ptr left, Data Line 168  resultLength(DataAbstract_ptr left, Data
168     }     }
169  }  }
170    
171    // determine the number of samples requires to evaluate an expression combining left and right
172  int  int
173  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
174  {  {
# Line 124  calcBuffs(const DataLazy_ptr& left, cons Line 183  calcBuffs(const DataLazy_ptr& left, cons
183  }  }
184    
185    
   
186  }   // end anonymous namespace  }   // end anonymous namespace
187    
188    
189    
190    // Return a string representing the operation
191  const std::string&  const std::string&
192  opToString(ES_optype op)  opToString(ES_optype op)
193  {  {
# Line 145  DataLazy::DataLazy(DataAbstract_ptr p) Line 205  DataLazy::DataLazy(DataAbstract_ptr p)
205  {  {
206     if (p->isLazy())     if (p->isLazy())
207     {     {
     // TODO: fix this.   We could make the new node a copy of p?  
208      // I don't want identity of Lazy.      // I don't want identity of Lazy.
209      // Question: Why would that be so bad?      // Question: Why would that be so bad?
210      // 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 213  DataLazy::DataLazy(DataAbstract_ptr p)
213     else     else
214     {     {
215      m_id=dynamic_pointer_cast<DataReady>(p);      m_id=dynamic_pointer_cast<DataReady>(p);
216        if(p->isConstant()) {m_readytype='C';}
217        else if(p->isExpanded()) {m_readytype='E';}
218        else if (p->isTagged()) {m_readytype='T';}
219        else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}
220     }     }
221     m_length=p->getLength();     m_length=p->getLength();
222     m_buffsRequired=1;     m_buffsRequired=1;
# Line 161  DataLazy::DataLazy(DataAbstract_ptr p) Line 224  DataLazy::DataLazy(DataAbstract_ptr p)
224  cout << "(1)Lazy created with " << m_samplesize << endl;  cout << "(1)Lazy created with " << m_samplesize << endl;
225  }  }
226    
227    
228    
229    
230  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
231      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
232      m_op(op)      m_op(op)
# Line 178  DataLazy::DataLazy(DataAbstract_ptr left Line 244  DataLazy::DataLazy(DataAbstract_ptr left
244     {     {
245      lleft=dynamic_pointer_cast<DataLazy>(left);      lleft=dynamic_pointer_cast<DataLazy>(left);
246     }     }
247       m_readytype=lleft->m_readytype;
248     m_length=left->getLength();     m_length=left->getLength();
249     m_left=lleft;     m_left=lleft;
250     m_buffsRequired=1;     m_buffsRequired=1;
# Line 185  DataLazy::DataLazy(DataAbstract_ptr left Line 252  DataLazy::DataLazy(DataAbstract_ptr left
252  }  }
253    
254    
 DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)  
     : parent(resultFS(left,right,op), resultShape(left,right,op)),  
     m_left(left),  
     m_right(right),  
     m_op(op)  
 {  
    if (getOpgroup(op)!=G_BINARY)  
    {  
     throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");  
    }  
    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;  
 }  
   
255  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
256      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
257      m_op(op)      m_op(op)
258  {  {
259     if (getOpgroup(op)!=G_BINARY)     if (getOpgroup(op)!=G_BINARY)
260     {     {
261      throw DataException("Programmer error - constructor DataLazy(left, op) will only process BINARY operations.");      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
262     }     }
263     if (left->isLazy())     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
264     {     {
265      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
266     }     }
# Line 225  DataLazy::DataLazy(DataAbstract_ptr left Line 276  DataLazy::DataLazy(DataAbstract_ptr left
276     {     {
277      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
278     }     }
279       char lt=m_left->m_readytype;
280       char rt=m_right->m_readytype;
281       if (lt=='E' || rt=='E')
282       {
283        m_readytype='E';
284       }
285       else if (lt=='T' || rt=='T')
286       {
287        m_readytype='T';
288       }
289       else
290       {
291        m_readytype='C';
292       }
293     m_length=resultLength(m_left,m_right,m_op);     m_length=resultLength(m_left,m_right,m_op);
294     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();    
295     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
296  cout << "(3)Lazy created with " << m_samplesize << endl;  cout << "(3)Lazy created with " << m_samplesize << endl;
297  }  }
# Line 245  DataLazy::getBuffsRequired() const Line 309  DataLazy::getBuffsRequired() const
309  }  }
310    
311    
312  // the vector and the offset are a place where the method could write its data if it wishes  /*
313  // 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.
314  // Hence the return value to indicate where the data is actually stored.    This does the work for the collapse method.
315  // 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.
316  const double*  */
317  DataLazy::resolveSample(ValueType& v,int sampleNo,  size_t offset ) const  DataReady_ptr
318    DataLazy::collapseToReady()
319  {  {
320    if (m_op==IDENTITY)      if (m_readytype=='E')
321      { // this is more an efficiency concern than anything else
322        throw DataException("Programmer Error - do not use collapse on Expanded data.");
323      }
324      if (m_op==IDENTITY)
325    {    {
326      const ValueType& vec=m_id->getVector();      return m_id;
     return &(vec[m_id->getPointOffset(sampleNo, 0)]);  
327    }    }
328    size_t rightoffset=offset+m_samplesize;    DataReady_ptr pleft=m_left->collapseToReady();
329    const double* left=m_left->resolveSample(v,sampleNo,offset);    Data left(pleft);
330    const double* right=0;    Data right;
331    if (getOpgroup(m_op)==G_BINARY)    if (getOpgroup(m_op)==G_BINARY)
332    {    {
333      right=m_right->resolveSample(v,sampleNo,rightoffset);      right=Data(m_right->collapseToReady());
334    }    }
335    double* result=&(v[offset]);    Data result;
336      switch(m_op)
337    {    {
338      switch(m_op)      case ADD:
339      {      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>());  
340      break;      break;
341      case SUB:            case SUB:      
342      tensor_binary_operation(m_samplesize, left, right, result, minus<double>());      result=left-right;
343      break;      break;
344      case MUL:            case MUL:      
345      tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());      result=left*right;
346      break;      break;
347      case DIV:            case DIV:      
348      tensor_binary_operation(m_samplesize, left, right, result, divides<double>());      result=left/right;
349      break;      break;
 // unary ops  
350      case SIN:      case SIN:
351        result=left.sin();  
352        break;
353        case COS:
354        result=left.cos();
355        break;
356        case TAN:
357        result=left.tan();
358        break;
359        case ASIN:
360        result=left.asin();
361        break;
362        case ACOS:
363        result=left.acos();
364        break;
365        case ATAN:
366        result=left.atan();
367        break;
368        case SINH:
369        result=left.sinh();
370        break;
371        case COSH:
372        result=left.cosh();
373        break;
374        case TANH:
375        result=left.tanh();
376        break;
377        case ERF:
378        result=left.erf();
379        break;
380       case ASINH:
381        result=left.asinh();
382        break;
383       case ACOSH:
384        result=left.acosh();
385        break;
386       case ATANH:
387        result=left.atanh();
388        break;
389        case LOG10:
390        result=left.log10();
391        break;
392        case LOG:
393        result=left.log();
394        break;
395        case SIGN:
396        result=left.sign();
397        break;
398        case ABS:
399        result=left.abs();
400        break;
401        case NEG:
402        result=left.neg();
403        break;
404        case POS:
405        // it doesn't mean anything for delayed.
406        // it will just trigger a deep copy of the lazy object
407        throw DataException("Programmer error - POS not supported for lazy data.");
408        break;
409        case EXP:
410        result=left.exp();
411        break;
412        case SQRT:
413        result=left.sqrt();
414        break;
415        case RECIP:
416        result=left.oneOver();
417        break;
418        case GZ:
419        result=left.wherePositive();
420        break;
421        case LZ:
422        result=left.whereNegative();
423        break;
424        case GEZ:
425        result=left.whereNonNegative();
426        break;
427        case LEZ:
428        result=left.whereNonPositive();
429        break;
430        default:
431        throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
432      }
433      return result.borrowReadyPtr();
434    }
435    
436    /*
437       \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
438       This method uses the original methods on the Data class to evaluate the expressions.
439       For this reason, it should not be used on DataExpanded instances. (To do so would defeat
440       the purpose of using DataLazy in the first place).
441    */
442    void
443    DataLazy::collapse()
444    {
445      if (m_op==IDENTITY)
446      {
447        return;
448      }
449      if (m_readytype=='E')
450      { // this is more an efficiency concern than anything else
451        throw DataException("Programmer Error - do not use collapse on Expanded data.");
452      }
453      m_id=collapseToReady();
454      m_op=IDENTITY;
455    }
456    
457    /*
458      \brief Compute the value of the expression (binary operation) for the given sample.
459      \return Vector which stores the value of the subexpression for the given sample.
460      \param v A vector to store intermediate results.
461      \param offset Index in v to begin storing results.
462      \param sampleNo Sample number to evaluate.
463      \param roffset (output parameter) the offset in the return vector where the result begins.
464    
465      The return value will be an existing vector so do not deallocate it.
466      If the result is stored in v it should be stored at the offset given.
467      Everything from offset to the end of v should be considered available for this method to use.
468    */
469    DataTypes::ValueType*
470    DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
471    {
472        // we assume that any collapsing has been done before we get here
473        // since we only have one argument we don't need to think about only
474        // processing single points.
475      if (m_readytype!='E')
476      {
477        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
478      }
479      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
480      const double* left=&((*vleft)[roffset]);
481      double* result=&(v[offset]);
482      roffset=offset;
483      switch (m_op)
484      {
485        case SIN:  
486      tensor_unary_operation(m_samplesize, left, result, ::sin);      tensor_unary_operation(m_samplesize, left, result, ::sin);
487      break;      break;
488      case COS:      case COS:
# Line 379  DataLazy::resolveSample(ValueType& v,int Line 580  DataLazy::resolveSample(ValueType& v,int
580      break;      break;
581    
582      default:      default:
583      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)+".");
584      }
585      return &v;
586    }
587    
588    
589    
590    
591    
592    #define PROC_OP(X) \
593        for (int i=0;i<steps;++i,resultp+=resultStep) \
594        { \
595           tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
596           lroffset+=leftStep; \
597           rroffset+=rightStep; \
598        }
599    
600    /*
601      \brief Compute the value of the expression (binary operation) for the given sample.
602      \return Vector which stores the value of the subexpression for the given sample.
603      \param v A vector to store intermediate results.
604      \param offset Index in v to begin storing results.
605      \param sampleNo Sample number to evaluate.
606      \param roffset (output parameter) the offset in the return vector where the result begins.
607    
608      The return value will be an existing vector so do not deallocate it.
609      If the result is stored in v it should be stored at the offset given.
610      Everything from offset to the end of v should be considered available for this method to use.
611    */
612    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
613    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
614    // If both children are expanded, then we can process them in a single operation (we treat
615    // the whole sample as one big datapoint.
616    // If one of the children is not expanded, then we need to treat each point in the sample
617    // individually.
618    // There is an additional complication when scalar operations are considered.
619    // For example, 2+Vector.
620    // In this case each double within the point is treated individually
621    DataTypes::ValueType*
622    DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
623    {
624    cout << "Resolve binary: " << toString() << endl;
625    
626      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
627        // first work out which of the children are expanded
628      bool leftExp=(m_left->m_readytype=='E');
629      bool rightExp=(m_right->m_readytype=='E');
630      bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
631      int steps=(bigloops?1:getNumDPPSample());
632      size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint
633      if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops
634      {
635        EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");
636        steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());
637        chunksize=1;    // for scalar
638      }    
639      int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
640      int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
641      int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0
642        // Get the values of sub-expressions
643      const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
644      const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
645        // the right child starts further along.
646      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
647      switch(m_op)
648      {
649        case ADD:
650        PROC_OP(plus<double>());
651        break;
652        case SUB:
653        PROC_OP(minus<double>());
654        break;
655        case MUL:
656        PROC_OP(multiplies<double>());
657        break;
658        case DIV:
659        PROC_OP(divides<double>());
660        break;
661        case POW:
662        PROC_OP(::pow);
663        break;
664        default:
665        throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
666      }
667      roffset=offset;  
668      return &v;
669    }
670    
671    
672    
673    /*
674      \brief Compute the value of the expression for the given sample.
675      \return Vector which stores the value of the subexpression for the given sample.
676      \param v A vector to store intermediate results.
677      \param offset Index in v to begin storing results.
678      \param sampleNo Sample number to evaluate.
679      \param roffset (output parameter) the offset in the return vector where the result begins.
680    
681      The return value will be an existing vector so do not deallocate it.
682    */
683    // the vector and the offset are a place where the method could write its data if it wishes
684    // it is not obligated to do so. For example, if it has its own storage already, it can use that.
685    // Hence the return value to indicate where the data is actually stored.
686    // Regardless, the storage should be assumed to be used, even if it isn't.
687    
688    // the roffset is the offset within the returned vector where the data begins
689    const DataTypes::ValueType*
690    DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
691    {
692    cout << "Resolve sample " << toString() << endl;
693        // collapse so we have a 'E' node or an IDENTITY for some other type
694      if (m_readytype!='E' && m_op!=IDENTITY)
695      {
696        collapse();
697      }
698      if (m_op==IDENTITY)  
699      {
700        const ValueType& vec=m_id->getVector();
701        if (m_readytype=='C')
702        {
703        roffset=0;
704        return &(vec);
705      }      }
706        roffset=m_id->getPointOffset(sampleNo, 0);
707        return &(vec);
708      }
709      if (m_readytype!='E')
710      {
711        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
712      }
713      switch (getOpgroup(m_op))
714      {
715      case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);
716      case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
717      default:
718        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
719    }    }
   return result;  
720  }  }
721    
722    
723    // To simplify the memory management, all threads operate on one large vector, rather than one each.
724    // Each sample is evaluated independently and copied into the result DataExpanded.
725  DataReady_ptr  DataReady_ptr
726  DataLazy::resolve()  DataLazy::resolve()
727  {  {
   // This is broken!     We need to have a buffer per thread!  
   // so the allocation of v needs to move inside the loop somehow  
728    
729  cout << "Sample size=" << m_samplesize << endl;  cout << "Sample size=" << m_samplesize << endl;
730  cout << "Buffers=" << m_buffsRequired << endl;  cout << "Buffers=" << m_buffsRequired << endl;
731    
732    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
733      {
734        collapse();
735      }
736      if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
737      {
738        return m_id;
739      }
740        // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
741      size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired));    // Each thread needs to have enough
742        // storage to evaluate its expression
743    int numthreads=1;    int numthreads=1;
744  #ifdef _OPENMP  #ifdef _OPENMP
745    numthreads=omp_get_max_threads();    numthreads=getNumberOfThreads();
746    int threadnum=0;    int threadnum=0;
747  #endif  #endif
748    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
749  cout << "Buffer created with size=" << v.size() << endl;  cout << "Buffer created with size=" << v.size() << endl;
750    ValueType dummy(getNoValues());    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
   DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),dummy);  
751    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVector();
752    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
753    int sample;    int sample;
754    int resoffset;    size_t outoffset;     // offset in the output data
755    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
756    #pragma omp parallel for private(sample,resoffset,threadnum) schedule(static)    const ValueType* res=0;   // Vector storing the answer
757      size_t resoffset=0;       // where in the vector to find the answer
758      #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)
759    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
760    {    {
761    cout << "################################# " << sample << endl;
762  #ifdef _OPENMP  #ifdef _OPENMP
763      const double* res=resolveSample(v,sample,threadbuffersize*omp_get_thread_num());      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
764  #else  #else
765      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.
766  #endif  #endif
767      resoffset=result->getPointOffset(sample,0);  cerr << "-------------------------------- " << endl;
768      for (int i=0;i<m_samplesize;++i,++resoffset)    // copy values into the output vector      outoffset=result->getPointOffset(sample,0);
769    cerr << "offset=" << outoffset << endl;
770        for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector
771      {      {
772      resvec[resoffset]=res[i];      resvec[outoffset]=(*res)[resoffset];
773      }      }
774    cerr << "*********************************" << endl;
775    }    }
776    return resptr;    return resptr;
777  }  }
# Line 435  DataLazy::toString() const Line 785  DataLazy::toString() const
785    return oss.str();    return oss.str();
786  }  }
787    
788    
789  void  void
790  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
791  {  {
792    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
793    {    {
794    case G_IDENTITY:    case G_IDENTITY:
795        if (m_id->isExpanded())
796        {
797           oss << "E";
798        }
799        else if (m_id->isTagged())
800        {
801          oss << "T";
802        }
803        else if (m_id->isConstant())
804        {
805          oss << "C";
806        }
807        else
808        {
809          oss << "?";
810        }
811      oss << '@' << m_id.get();      oss << '@' << m_id.get();
812      break;      break;
813    case G_BINARY:    case G_BINARY:
# Line 493  DataLazy::getPointOffset(int sampleNo, Line 860  DataLazy::getPointOffset(int sampleNo,
860    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");
861  }  }
862    
863    // It would seem that DataTagged will need to be treated differently since even after setting all tags
864    // to zero, all the tags from all the DataTags would be in the result.
865    // However since they all have the same value (0) whether they are there or not should not matter.
866    // So I have decided that for all types this method will create a constant 0.
867    // It can be promoted up as required.
868    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
869    // but we can deal with that if it arrises.
870    void
871    DataLazy::setToZero()
872    {
873      DataTypes::ValueType v(getNoValues(),0);
874      m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
875      m_op=IDENTITY;
876      m_right.reset();  
877      m_left.reset();
878      m_readytype='C';
879      m_buffsRequired=1;
880    }
881    
882  }   // end namespace  }   // end namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26