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

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

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

revision 1889 by jfenwick, Thu Oct 16 05:57:09 2008 UTC revision 1901 by jfenwick, Wed Oct 22 02:44:34 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 98  resultShape(DataAbstract_ptr left, DataA Line 141  resultShape(DataAbstract_ptr left, DataA
141      return left->getShape();      return left->getShape();
142  }  }
143    
144    // determine the number of points in the result of "left op right"
145  size_t  size_t
146  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
147  {  {
# Line 110  resultLength(DataAbstract_ptr left, Data Line 154  resultLength(DataAbstract_ptr left, Data
154     }     }
155  }  }
156    
157    // determine the number of samples requires to evaluate an expression combining left and right
158  int  int
159  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
160  {  {
# Line 123  calcBuffs(const DataLazy_ptr& left, cons Line 168  calcBuffs(const DataLazy_ptr& left, cons
168     }     }
169  }  }
170    
171    
172  }   // end anonymous namespace  }   // end anonymous namespace
173    
174    
175    
176    // Return a string representing the operation
177  const std::string&  const std::string&
178  opToString(ES_optype op)  opToString(ES_optype op)
179  {  {
# Line 143  DataLazy::DataLazy(DataAbstract_ptr p) Line 191  DataLazy::DataLazy(DataAbstract_ptr p)
191  {  {
192     if (p->isLazy())     if (p->isLazy())
193     {     {
     // TODO: fix this.   We could make the new node a copy of p?  
194      // I don't want identity of Lazy.      // I don't want identity of Lazy.
195      // Question: Why would that be so bad?      // Question: Why would that be so bad?
196      // 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 163  DataLazy::DataLazy(DataAbstract_ptr p) Line 210  DataLazy::DataLazy(DataAbstract_ptr p)
210  cout << "(1)Lazy created with " << m_samplesize << endl;  cout << "(1)Lazy created with " << m_samplesize << endl;
211  }  }
212    
213    
214    
215    
216  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
217      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
218      m_op(op)      m_op(op)
# Line 188  DataLazy::DataLazy(DataAbstract_ptr left Line 238  DataLazy::DataLazy(DataAbstract_ptr left
238  }  }
239    
240    
 // 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;  
 // }  
   
241  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
242      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
243      m_op(op)      m_op(op)
# Line 212  DataLazy::DataLazy(DataAbstract_ptr left Line 246  DataLazy::DataLazy(DataAbstract_ptr left
246     {     {
247      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.");
248     }     }
249     if (left->isLazy())     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
250     {     {
251      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
252     }     }
# Line 243  DataLazy::DataLazy(DataAbstract_ptr left Line 277  DataLazy::DataLazy(DataAbstract_ptr left
277      m_readytype='C';      m_readytype='C';
278     }     }
279     m_length=resultLength(m_left,m_right,m_op);     m_length=resultLength(m_left,m_right,m_op);
280     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();    
281     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
282  cout << "(3)Lazy created with " << m_samplesize << endl;  cout << "(3)Lazy created with " << m_samplesize << endl;
283  }  }
# Line 261  DataLazy::getBuffsRequired() const Line 295  DataLazy::getBuffsRequired() const
295  }  }
296    
297    
298    /*
299      \brief Evaluates the expression using methods on Data.
300      This does the work for the collapse method.
301      For reasons of efficiency do not call this method on DataExpanded nodes.
302    */
303  DataReady_ptr  DataReady_ptr
304  DataLazy::collapseToReady()  DataLazy::collapseToReady()
305  {  {
# Line 380  DataLazy::collapseToReady() Line 419  DataLazy::collapseToReady()
419    return result.borrowReadyPtr();    return result.borrowReadyPtr();
420  }  }
421    
422    /*
423       \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
424       This method uses the original methods on the Data class to evaluate the expressions.
425       For this reason, it should not be used on DataExpanded instances. (To do so would defeat
426       the purpose of using DataLazy in the first place).
427    */
428  void  void
429  DataLazy::collapse()  DataLazy::collapse()
430  {  {
# Line 395  DataLazy::collapse() Line 440  DataLazy::collapse()
440    m_op=IDENTITY;    m_op=IDENTITY;
441  }  }
442    
443  const double*  /*
444  DataLazy::resolveUnary(ValueType& v,int sampleNo,  size_t offset) const    \brief Compute the value of the expression (binary operation) for the given sample.
445      \return Vector which stores the value of the subexpression for the given sample.
446      \param v A vector to store intermediate results.
447      \param offset Index in v to begin storing results.
448      \param sampleNo Sample number to evaluate.
449      \param roffset (output parameter) the offset in the return vector where the result begins.
450    
451      The return value will be an existing vector so do not deallocate it.
452      If the result is stored in v it should be stored at the offset given.
453      Everything from offset to the end of v should be considered available for this method to use.
454    */
455    DataTypes::ValueType*
456    DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
457  {  {
458      // we assume that any collapsing has been done before we get here      // we assume that any collapsing has been done before we get here
459      // since we only have one argument we don't need to think about only      // since we only have one argument we don't need to think about only
# Line 405  DataLazy::resolveUnary(ValueType& v,int Line 462  DataLazy::resolveUnary(ValueType& v,int
462    {    {
463      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");      throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
464    }    }
465    const double* left=m_left->resolveSample(v,sampleNo,offset);    const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
466      const double* left=&((*vleft)[roffset]);
467    double* result=&(v[offset]);    double* result=&(v[offset]);
468      roffset=offset;
469    switch (m_op)    switch (m_op)
470    {    {
471      case SIN:        case SIN:  
# Line 509  DataLazy::resolveUnary(ValueType& v,int Line 568  DataLazy::resolveUnary(ValueType& v,int
568      default:      default:
569      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
570    }    }
571    return result;    return &v;
572  }  }
573    
574    
575    
576    
577    
578  #define PROC_OP(X) \  #define PROC_OP(X) \
579      for (int i=0;i<steps;++i,resultp+=getNoValues()) \      for (int i=0;i<steps;++i,resultp+=getNoValues()) \
580      { \      { \
581  cout << "Step#" << i << " chunk=" << chunksize << endl; \         tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
582  cout << left[0] << left[1] << left[2] << endl; \         lroffset+=leftStep; \
583  cout << right[0] << right[1] << right[2] << endl; \         rroffset+=rightStep; \
        tensor_binary_operation(chunksize, left, right, resultp, X); \  
        left+=leftStep; \  
        right+=rightStep; \  
 cout << "Result=" << result << " " << result[0] << result[1] << result[2] << endl; \  
584      }      }
585    
586  const double*  /*
587  DataLazy::resolveBinary(ValueType& v,int sampleNo,  size_t offset) const    \brief Compute the value of the expression (binary operation) for the given sample.
588      \return Vector which stores the value of the subexpression for the given sample.
589      \param v A vector to store intermediate results.
590      \param offset Index in v to begin storing results.
591      \param sampleNo Sample number to evaluate.
592      \param roffset (output parameter) the offset in the return vector where the result begins.
593    
594      The return value will be an existing vector so do not deallocate it.
595      If the result is stored in v it should be stored at the offset given.
596      Everything from offset to the end of v should be considered available for this method to use.
597    */
598    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
599    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
600    // If both children are expanded, then we can process them in a single operation (we treat
601    // the whole sample as one big datapoint.
602    // If one of the children is not expanded, then we need to treat each point in the sample
603    // individually.
604    DataTypes::ValueType*
605    DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
606  {  {
     // again we assume that all collapsing has already been done  
     // so we have at least one expanded child.  
     // however, we could still have one of the children being not expanded.  
   
607  cout << "Resolve binary: " << toString() << endl;  cout << "Resolve binary: " << toString() << endl;
608    
609    const double* left=m_left->resolveSample(v,sampleNo,offset);    size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
610  cout << "Done Left " << left[0] << left[1] << left[2] << endl;      // first work out which of the children are expanded
   const double* right=m_right->resolveSample(v,sampleNo,offset);      
 cout << "Done Right"  << right[0] << right[1] << right[2] << endl;  
     // now we need to know which args are expanded  
611    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
612    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
613    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
614    int steps=(bigloops?1:getNumSamples());    int steps=(bigloops?1:getNumDPPSample());    
615    size_t chunksize=(bigloops? m_samplesize : getNoValues());    size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint
616    int leftStep=(rightExp? getNoValues() : 0);    int leftStep=((leftExp && !rightExp)? getNoValues() : 0);
617    int rightStep=(leftExp? getNoValues() : 0);    int rightStep=((rightExp && !leftExp)? getNoValues() : 0);
618  cout << "left=" << left << " right=" << right << endl;      // Get the values of sub-expressions
619    double* result=&(v[offset]);    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
620    double* resultp=result;    const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
621        // the right child starts further along.
622      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
623    switch(m_op)    switch(m_op)
624    {    {
625      case ADD:      case ADD:
626      PROC_OP(plus<double>())      PROC_OP(plus<double>());
627        break;
628        case SUB:
629        PROC_OP(minus<double>());
630        break;
631        case MUL:
632        PROC_OP(multiplies<double>());
633        break;
634        case DIV:
635        PROC_OP(divides<double>());
636      break;      break;
 // need to fill in the rest  
637      default:      default:
638      throw DataException("Programmer error - resolveBinay can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
639    }    }
640  cout << "About to return "  << result[0] << result[1] << result[2] << endl;;    roffset=offset;  
641    return result;    return &v;
642  }  }
643    
644    
645    
646    /*
647      \brief Compute the value of the expression for the given sample.
648      \return Vector which stores the value of the subexpression for the given sample.
649      \param v A vector to store intermediate results.
650      \param offset Index in v to begin storing results.
651      \param sampleNo Sample number to evaluate.
652      \param roffset (output parameter) the offset in the return vector where the result begins.
653    
654      The return value will be an existing vector so do not deallocate it.
655    */
656  // the vector and the offset are a place where the method could write its data if it wishes  // the vector and the offset are a place where the method could write its data if it wishes
657  // it is not obligated to do so. For example, if it has its own storage already, it can use that.  // it is not obligated to do so. For example, if it has its own storage already, it can use that.
658  // Hence the return value to indicate where the data is actually stored.  // Hence the return value to indicate where the data is actually stored.
659  // Regardless, the storage should be assumed to be used, even if it isn't.  // Regardless, the storage should be assumed to be used, even if it isn't.
660  const double*  
661  DataLazy::resolveSample(ValueType& v,int sampleNo,  size_t offset )  // the roffset is the offset within the returned vector where the data begins
662    const DataTypes::ValueType*
663    DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
664  {  {
665  cout << "Resolve sample " << toString() << endl;  cout << "Resolve sample " << toString() << endl;
666      // collapse so we have a 'E' node or an IDENTITY for some other type      // collapse so we have a 'E' node or an IDENTITY for some other type
# Line 575  cout << "Resolve sample " << toString() Line 668  cout << "Resolve sample " << toString()
668    {    {
669      collapse();      collapse();
670    }    }
 cout << "Post collapse check\n";  
671    if (m_op==IDENTITY)      if (m_op==IDENTITY)  
672    {    {
 cout << "In IDENTITY check\n";  
673      const ValueType& vec=m_id->getVector();      const ValueType& vec=m_id->getVector();
674      if (m_readytype=='C')      if (m_readytype=='C')
675      {      {
676      return &(vec[0]);      roffset=0;
677        return &(vec);
678      }      }
679      return &(vec[m_id->getPointOffset(sampleNo, 0)]);      roffset=m_id->getPointOffset(sampleNo, 0);
680        return &(vec);
681    }    }
682    if (m_readytype!='E')    if (m_readytype!='E')
683    {    {
684      throw DataException("Programmer Error - Collapse did not produce an expanded node.");      throw DataException("Programmer Error - Collapse did not produce an expanded node.");
685    }    }
 cout << "Calling sub resolvers\n";  
686    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
687    {    {
688    case G_UNARY: return resolveUnary(v,sampleNo,offset);    case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);
689    case G_BINARY: return resolveBinary(v,sampleNo,offset);    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
690    default:    default:
691      throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");      throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
692    }    }
693  }  }
694    
695    
696  // the vector and the offset are a place where the method could write its data if it wishes  // To simplify the memory management, all threads operate on one large vector, rather than one each.
697  // it is not obligated to do so. For example, if it has its own storage already, it can use that.  // Each sample is evaluated independently and copied into the result DataExpanded.
 // Hence the return value to indicate where the data is actually stored.  
 // Regardless, the storage should be assumed to be used, even if it isn't.  
 const double*  
 DataLazy::resolveSample2(ValueType& v,int sampleNo,  size_t offset )  
 {  
   if (m_readytype!='E')  
   {  
     throw DataException("Only supporting Expanded Data.");  
   }  
   if (m_op==IDENTITY)    
   {  
     const ValueType& vec=m_id->getVector();  
     return &(vec[m_id->getPointOffset(sampleNo, 0)]);  
   }  
   size_t rightoffset=offset+m_samplesize;  
   const double* left=m_left->resolveSample(v,sampleNo,offset);  
   const double* right=0;  
   if (getOpgroup(m_op)==G_BINARY)  
   {  
     right=m_right->resolveSample(v,sampleNo,rightoffset);  
   }  
   double* result=&(v[offset]);  
   {  
     switch(m_op)  
     {  
     case ADD:       // since these are pointwise ops, pretend each sample is one point  
     tensor_binary_operation(m_samplesize, left, right, result, plus<double>());  
     break;  
     case SUB:        
     tensor_binary_operation(m_samplesize, left, right, result, minus<double>());  
     break;  
     case MUL:        
     tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());  
     break;  
     case DIV:        
     tensor_binary_operation(m_samplesize, left, right, result, divides<double>());  
     break;  
 // unary ops  
     case SIN:  
     tensor_unary_operation(m_samplesize, left, result, ::sin);  
     break;  
     case COS:  
     tensor_unary_operation(m_samplesize, left, result, ::cos);  
     break;  
     case TAN:  
     tensor_unary_operation(m_samplesize, left, result, ::tan);  
     break;  
     case ASIN:  
     tensor_unary_operation(m_samplesize, left, result, ::asin);  
     break;  
     case ACOS:  
     tensor_unary_operation(m_samplesize, left, result, ::acos);  
     break;  
     case ATAN:  
     tensor_unary_operation(m_samplesize, left, result, ::atan);  
     break;  
     case SINH:  
     tensor_unary_operation(m_samplesize, left, result, ::sinh);  
     break;  
     case COSH:  
     tensor_unary_operation(m_samplesize, left, result, ::cosh);  
     break;  
     case TANH:  
     tensor_unary_operation(m_samplesize, left, result, ::tanh);  
     break;  
     case ERF:  
 #ifdef _WIN32  
     throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::erf);  
     break;  
 #endif  
    case ASINH:  
 #ifdef _WIN32  
     tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::asinh);  
 #endif    
     break;  
    case ACOSH:  
 #ifdef _WIN32  
     tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::acosh);  
 #endif    
     break;  
    case ATANH:  
 #ifdef _WIN32  
     tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);  
 #else  
     tensor_unary_operation(m_samplesize, left, result, ::atanh);  
 #endif    
     break;  
     case LOG10:  
     tensor_unary_operation(m_samplesize, left, result, ::log10);  
     break;  
     case LOG:  
     tensor_unary_operation(m_samplesize, left, result, ::log);  
     break;  
     case SIGN:  
     tensor_unary_operation(m_samplesize, left, result, escript::fsign);  
     break;  
     case ABS:  
     tensor_unary_operation(m_samplesize, left, result, ::fabs);  
     break;  
     case NEG:  
     tensor_unary_operation(m_samplesize, left, result, negate<double>());  
     break;  
     case POS:  
     // it doesn't mean anything for delayed.  
     // it will just trigger a deep copy of the lazy object  
     throw DataException("Programmer error - POS not supported for lazy data.");  
     break;  
     case EXP:  
     tensor_unary_operation(m_samplesize, left, result, ::exp);  
     break;  
     case SQRT:  
     tensor_unary_operation(m_samplesize, left, result, ::sqrt);  
     break;  
     case RECIP:  
     tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));  
     break;  
     case GZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));  
     break;  
     case LZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));  
     break;  
     case GEZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));  
     break;  
     case LEZ:  
     tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));  
     break;  
   
     default:  
     throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");  
     }  
   }  
   return result;  
 }  
   
698  DataReady_ptr  DataReady_ptr
699  DataLazy::resolve()  DataLazy::resolve()
700  {  {
# Line 752  DataLazy::resolve() Line 702  DataLazy::resolve()
702  cout << "Sample size=" << m_samplesize << endl;  cout << "Sample size=" << m_samplesize << endl;
703  cout << "Buffers=" << m_buffsRequired << endl;  cout << "Buffers=" << m_buffsRequired << endl;
704    
705    if (m_readytype!='E')    if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally
706    {    {
707      collapse();      collapse();
708    }    }
709    if (m_op==IDENTITY)    if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
710    {    {
711      return m_id;      return m_id;
712    }    }
713      // from this point on we must have m_op!=IDENTITY and m_readytype=='E'      // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
714    size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);    size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired));    // Each thread needs to have enough
715        // storage to evaluate its expression
716    int numthreads=1;    int numthreads=1;
717  #ifdef _OPENMP  #ifdef _OPENMP
718    numthreads=getNumberOfThreads();    numthreads=getNumberOfThreads();
# Line 769  cout << "Buffers=" << m_buffsRequired << Line 720  cout << "Buffers=" << m_buffsRequired <<
720  #endif  #endif
721    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
722  cout << "Buffer created with size=" << v.size() << endl;  cout << "Buffer created with size=" << v.size() << endl;
723    ValueType dummy(getNoValues());    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
   DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),dummy);  
724    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVector();
725    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
726    int sample;    int sample;
727    int resoffset;    size_t outoffset;     // offset in the output data
728    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
729    #pragma omp parallel for private(sample,resoffset,threadnum) schedule(static)    const ValueType* res=0;   // Vector storing the answer
730      size_t resoffset=0;       // where in the vector to find the answer
731      #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)
732    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
733    {    {
734    cout << "################################# " << sample << endl;
735  #ifdef _OPENMP  #ifdef _OPENMP
736      const double* res=resolveSample(v,sample,threadbuffersize*omp_get_thread_num());      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
737  #else  #else
738      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.
739  #endif  #endif
740      resoffset=result->getPointOffset(sample,0);  cerr << "-------------------------------- " << endl;
741      for (unsigned int i=0;i<m_samplesize;++i,++resoffset)   // copy values into the output vector      outoffset=result->getPointOffset(sample,0);
742    cerr << "offset=" << outoffset << endl;
743        for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector
744      {      {
745      resvec[resoffset]=res[i];      resvec[outoffset]=(*res)[resoffset];
746      }      }
747    cerr << "*********************************" << endl;
748    }    }
749    return resptr;    return resptr;
750  }  }
# Line 802  DataLazy::toString() const Line 758  DataLazy::toString() const
758    return oss.str();    return oss.str();
759  }  }
760    
761    
762  void  void
763  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
764  {  {
# Line 876  DataLazy::getPointOffset(int sampleNo, Line 833  DataLazy::getPointOffset(int sampleNo,
833    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");
834  }  }
835    
836    // It would seem that DataTagged will need to be treated differently since even after setting all tags
837    // to zero, all the tags from all the DataTags would be in the result.
838    // However since they all have the same value (0) whether they are there or not should not matter.
839    // So I have decided that for all types this method will create a constant 0.
840    // It can be promoted up as required.
841    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
842    // but we can deal with that if it arrises.
843    void
844    DataLazy::setToZero()
845    {
846      DataTypes::ValueType v(getNoValues(),0);
847      m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
848      m_op=IDENTITY;
849      m_right.reset();  
850      m_left.reset();
851      m_readytype='C';
852      m_buffsRequired=1;
853    }
854    
855  }   // end namespace  }   // end namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26