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

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

  ViewVC Help
Powered by ViewVC 1.1.26