/[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 1898 by jfenwick, Mon Oct 20 01:20:18 2008 UTC revision 1899 by jfenwick, Mon Oct 20 05:13:24 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  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
216      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
217      m_op(op)      m_op(op)
# Line 188  DataLazy::DataLazy(DataAbstract_ptr left Line 237  DataLazy::DataLazy(DataAbstract_ptr left
237  }  }
238    
239    
 // 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;  
 // }  
   
240  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
241      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
242      m_op(op)      m_op(op)
# Line 212  DataLazy::DataLazy(DataAbstract_ptr left Line 245  DataLazy::DataLazy(DataAbstract_ptr left
245     {     {
246      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.");
247     }     }
248     if (left->isLazy())     if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
249     {     {
250      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
251     }     }
# Line 243  DataLazy::DataLazy(DataAbstract_ptr left Line 276  DataLazy::DataLazy(DataAbstract_ptr left
276      m_readytype='C';      m_readytype='C';
277     }     }
278     m_length=resultLength(m_left,m_right,m_op);     m_length=resultLength(m_left,m_right,m_op);
279     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();    
280     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
281  cout << "(3)Lazy created with " << m_samplesize << endl;  cout << "(3)Lazy created with " << m_samplesize << endl;
282  }  }
# Line 261  DataLazy::getBuffsRequired() const Line 294  DataLazy::getBuffsRequired() const
294  }  }
295    
296    
297    /*
298      \brief Evaluates the expression using methods on Data.
299      This does the work for the collapse method.
300      For reasons of efficiency do not call this method on DataExpanded nodes.
301    */
302  DataReady_ptr  DataReady_ptr
303  DataLazy::collapseToReady()  DataLazy::collapseToReady()
304  {  {
# Line 380  DataLazy::collapseToReady() Line 418  DataLazy::collapseToReady()
418    return result.borrowReadyPtr();    return result.borrowReadyPtr();
419  }  }
420    
421    /*
422       \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
423       This method uses the original methods on the Data class to evaluate the expressions.
424       For this reason, it should not be used on DataExpanded instances. (To do so would defeat
425       the purpose of using DataLazy in the first place).
426    */
427  void  void
428  DataLazy::collapse()  DataLazy::collapse()
429  {  {
# Line 395  DataLazy::collapse() Line 439  DataLazy::collapse()
439    m_op=IDENTITY;    m_op=IDENTITY;
440  }  }
441    
442    /*
443      \brief Compute the value of the expression (binary operation) for the given sample.
444      \return Vector which stores the value of the subexpression for the given sample.
445      \param v A vector to store intermediate results.
446      \param offset Index in v to begin storing results.
447      \param sampleNo Sample number to evaluate.
448      \param roffset (output parameter) the offset in the return vector where the result begins.
449    
450      The return value will be an existing vector so do not deallocate it.
451      If the result is stored in v it should be stored at the offset given.
452      Everything from offset to the end of v should be considered available for this method to use.
453    */
454  DataTypes::ValueType*  DataTypes::ValueType*
455  DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
456  {  {
# Line 516  DataLazy::resolveUnary(ValueType& v, siz Line 572  DataLazy::resolveUnary(ValueType& v, siz
572    
573    
574    
575  // const double*  
 // DataLazy::resolveUnary(ValueType& v,int sampleNo,  size_t offset) const  
 // {  
 //  // we assume that any collapsing has been done before we get here  
 //  // since we only have one argument we don't need to think about only  
 //  // processing single points.  
 //   if (m_readytype!='E')  
 //   {  
 //     throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");  
 //   }  
 //   const double* left=m_left->resolveSample(v,sampleNo,offset);  
 //   double* result=&(v[offset]);  
 //   switch (m_op)  
 //   {  
 //     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 - resolveUnary can not resolve operator "+opToString(m_op)+".");  
 //   }  
 //   return result;  
 // }  
576    
577  #define PROC_OP(X) \  #define PROC_OP(X) \
578      for (int i=0;i<steps;++i,resultp+=getNoValues()) \      for (int i=0;i<steps;++i,resultp+=getNoValues()) \
579      { \      { \
580  cout << "Step#" << i << " chunk=" << chunksize << endl; \         tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
581  cout << left[0] << left[1] << left[2] << endl; \         lroffset+=leftStep; \
582  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; \  
583      }      }
584    
585    /*
586      \brief Compute the value of the expression (binary operation) for the given sample.
587      \return Vector which stores the value of the subexpression for the given sample.
588      \param v A vector to store intermediate results.
589      \param offset Index in v to begin storing results.
590      \param sampleNo Sample number to evaluate.
591      \param roffset (output parameter) the offset in the return vector where the result begins.
592    
593      The return value will be an existing vector so do not deallocate it.
594      If the result is stored in v it should be stored at the offset given.
595      Everything from offset to the end of v should be considered available for this method to use.
596    */
597    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
598    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
599    // If both children are expanded, then we can process them in a single operation (we treat
600    // the whole sample as one big datapoint.
601    // If one of the children is not expanded, then we need to treat each point in the sample
602    // individually.
603  DataTypes::ValueType*  DataTypes::ValueType*
604  DataLazy::resolveBinary(ValueType& v,  size_t offset ,int sampleNo, size_t& roffset) const  DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
605  {  {
     // 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.  
   
606  cout << "Resolve binary: " << toString() << endl;  cout << "Resolve binary: " << toString() << endl;
607    
608    size_t lroffset=0, rroffset=0;    size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
609        // first work out which of the children are expanded
610    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
611    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
612    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
613    int steps=(bigloops?1:getNumDPPSample());    int steps=(bigloops?1:getNumDPPSample());    
614    size_t chunksize=(bigloops? m_samplesize : getNoValues());    size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint
615    int leftStep=((leftExp && !rightExp)? getNoValues() : 0);    int leftStep=((leftExp && !rightExp)? getNoValues() : 0);
616    int rightStep=((rightExp && !leftExp)? getNoValues() : 0);    int rightStep=((rightExp && !leftExp)? getNoValues() : 0);
617        // Get the values of sub-expressions
618    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
619    const ValueType* right=m_right->resolveSample(v,offset,sampleNo,rroffset);        const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
620      // now we need to know which args are expanded      // the right child starts further along.
621  cout << "left=" << left << " right=" << right << endl;    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
 cout << "(Length) l=" << left->size() << " r=" << right->size() << " res=" << v.size() << endl;  
   double* resultp=&(v[offset]);  
622    switch(m_op)    switch(m_op)
623    {    {
624      case ADD:      case ADD:
625      for (int i=0;i<steps;++i,resultp+=getNoValues())      PROC_OP(plus<double>());
626      {      break;
627  cerr << "Step#" << i << " chunk=" << chunksize << endl;      case SUB:
628  cerr << left << "[" << lroffset << "] " << right << "[" << rroffset << "]" << endl;      PROC_OP(minus<double>());
629         tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, plus<double>());      break;
630         lroffset+=leftStep;      case MUL:
631         rroffset+=rightStep;      PROC_OP(multiplies<double>());
632  cerr << "left=" << lroffset << " right=" << rroffset << endl;      break;
633      }      case DIV:
634        PROC_OP(divides<double>());
635      break;      break;
 // need to fill in the rest  
636      default:      default:
637      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
638    }    }
639    roffset=offset;    roffset=offset;  
640    return &v;    return &v;
641  }  }
642    
643    
644    
645  // #define PROC_OP(X) \  /*
646  //  for (int i=0;i<steps;++i,resultp+=getNoValues()) \    \brief Compute the value of the expression for the given sample.
647  //  { \    \return Vector which stores the value of the subexpression for the given sample.
648  // cout << "Step#" << i << " chunk=" << chunksize << endl; \    \param v A vector to store intermediate results.
649  // cout << left[0] << left[1] << left[2] << endl; \    \param offset Index in v to begin storing results.
650  // cout << right[0] << right[1] << right[2] << endl; \    \param sampleNo Sample number to evaluate.
651  //     tensor_binary_operation(chunksize, left, right, resultp, X); \    \param roffset (output parameter) the offset in the return vector where the result begins.
 //     left+=leftStep; \  
 //     right+=rightStep; \  
 // cout << "Result=" << result << " " << result[0] << result[1] << result[2] << endl; \  
 //  }  
 //  
 // const double*  
 // DataLazy::resolveBinary(ValueType& v,int sampleNo,  size_t offset) const  
 // {  
 //  // 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.  
 //  
 // cout << "Resolve binary: " << toString() << endl;  
 //  
 //   const double* left=m_left->resolveSample(v,sampleNo,offset);  
 // // cout << "Done Left " << /*left[0] << left[1] << left[2] << */endl;  
 //   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  
 //   bool leftExp=(m_left->m_readytype=='E');  
 //   bool rightExp=(m_right->m_readytype=='E');  
 //   bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp));  // is processing in single step  
 //   int steps=(bigloops?1:getNumSamples());  
 //   size_t chunksize=(bigloops? m_samplesize : getNoValues());  
 //   int leftStep=((leftExp && !rightExp)? getNoValues() : 0);  
 //   int rightStep=((rightExp && !leftExp)? getNoValues() : 0);  
 // cout << "left=" << left << " right=" << right << endl;  
 //   double* result=&(v[offset]);  
 //   double* resultp=result;  
 //   switch(m_op)  
 //   {  
 //     case ADD:  
 //  for (int i=0;i<steps;++i,resultp+=getNoValues())  
 //  {  
 // cout << "Step#" << i << " chunk=" << chunksize << endl; \  
 // // cout << left[0] << left[1] << left[2] << endl;  
 // // cout << right[0] << right[1] << right[2] << endl;  
 //     tensor_binary_operation(chunksize, left, right, resultp, plus<double>());  
 // cout << "left=" << left << " right=" << right << " resp=" << resultp << endl;  
 //     left+=leftStep;  
 //     right+=rightStep;  
 // cout << "left=" << left << " right=" << right << endl;  
 // // cout << "Result=" << result << " " << result[0] << result[1] << result[2] << endl;  
 //  }  
 //  break;  
 // // need to fill in the rest  
 //     default:  
 //  throw DataException("Programmer error - resolveBinay can not resolve operator "+opToString(m_op)+".");  
 //   }  
 // // cout << "About to return "  << result[0] << result[1] << result[2] << endl;;  
 //   return result;  
 // }  
   
 // // the vector and the offset are a place where the method could write its data if it wishes  
 // // it is not obligated to do so. For example, if it has its own storage already, it can use that.  
 // // 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::resolveSample(ValueType& v,int sampleNo,  size_t offset )  
 // {  
 // cout << "Resolve sample " << toString() << endl;  
 //  // collapse so we have a 'E' node or an IDENTITY for some other type  
 //   if (m_readytype!='E' && m_op!=IDENTITY)  
 //   {  
 //  collapse();  
 //   }  
 //   if (m_op==IDENTITY)      
 //   {  
 //     const ValueType& vec=m_id->getVector();  
 //     if (m_readytype=='C')  
 //     {  
 //  return &(vec[0]);  
 //     }  
 //     return &(vec[m_id->getPointOffset(sampleNo, 0)]);  
 //   }  
 //   if (m_readytype!='E')  
 //   {  
 //     throw DataException("Programmer Error - Collapse did not produce an expanded node.");  
 //   }  
 //   switch (getOpgroup(m_op))  
 //   {  
 //   case G_UNARY: return resolveUnary(v,sampleNo,offset);  
 //   case G_BINARY: return resolveBinary(v,sampleNo,offset);  
 //   default:  
 //     throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");  
 //   }  
 // }  
   
   
652    
653      The return value will be an existing vector so do not deallocate it.
654    */
655  // 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
656  // 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.
657  // Hence the return value to indicate where the data is actually stored.  // Hence the return value to indicate where the data is actually stored.
# Line 830  cout << "Resolve sample " << toString() Line 692  cout << "Resolve sample " << toString()
692  }  }
693    
694    
695    // To simplify the memory management, all threads operate on one large vector, rather than one each.
696    // Each sample is evaluated independently and copied into the result DataExpanded.
 // This version uses double* trying again with vectors  
 // DataReady_ptr  
 // DataLazy::resolve()  
 // {  
 //  
 // cout << "Sample size=" << m_samplesize << endl;  
 // cout << "Buffers=" << m_buffsRequired << endl;  
 //  
 //   if (m_readytype!='E')  
 //   {  
 //     collapse();  
 //   }  
 //   if (m_op==IDENTITY)  
 //   {  
 //     return m_id;  
 //   }  
 //      // from this point on we must have m_op!=IDENTITY and m_readytype=='E'  
 //   size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);  
 //   int numthreads=1;  
 // #ifdef _OPENMP  
 //   numthreads=getNumberOfThreads();  
 //   int threadnum=0;  
 // #endif  
 //   ValueType v(numthreads*threadbuffersize);    
 // cout << "Buffer created with size=" << v.size() << endl;  
 //   DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));  
 //   ValueType& resvec=result->getVector();  
 //   DataReady_ptr resptr=DataReady_ptr(result);  
 //   int sample;  
 //   int resoffset;  
 //   int totalsamples=getNumSamples();  
 //   const double* res=0;  
 //   #pragma omp parallel for private(sample,resoffset,threadnum,res) schedule(static)  
 //   for (sample=0;sample<totalsamples;++sample)  
 //   {  
 // cout << "################################# " << sample << endl;  
 // #ifdef _OPENMP  
 //     res=resolveSample(v,sample,threadbuffersize*omp_get_thread_num());  
 // #else  
 //     res=resolveSample(v,sample,0);   // this would normally be v, but not if its a single IDENTITY op.  
 // #endif  
 // cerr << "-------------------------------- " << endl;  
 //     resoffset=result->getPointOffset(sample,0);  
 // cerr << "offset=" << resoffset << endl;  
 //     for (unsigned int i=0;i<m_samplesize;++i,++resoffset)    // copy values into the output vector  
 //     {  
 //  resvec[resoffset]=res[i];  
 //     }  
 // cerr << "*********************************" << endl;  
 //   }  
 //   return resptr;  
 // }  
   
   
697  DataReady_ptr  DataReady_ptr
698  DataLazy::resolve()  DataLazy::resolve()
699  {  {
# Line 893  DataLazy::resolve() Line 701  DataLazy::resolve()
701  cout << "Sample size=" << m_samplesize << endl;  cout << "Sample size=" << m_samplesize << endl;
702  cout << "Buffers=" << m_buffsRequired << endl;  cout << "Buffers=" << m_buffsRequired << endl;
703    
704    if (m_readytype!='E')    if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally
705    {    {
706      collapse();      collapse();
707    }    }
708    if (m_op==IDENTITY)    if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
709    {    {
710      return m_id;      return m_id;
711    }    }
712      // 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'
713    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
714        // storage to evaluate its expression
715    int numthreads=1;    int numthreads=1;
716  #ifdef _OPENMP  #ifdef _OPENMP
717    numthreads=getNumberOfThreads();    numthreads=getNumberOfThreads();
# Line 916  cout << "Buffer created with size=" << v Line 725  cout << "Buffer created with size=" << v
725    int sample;    int sample;
726    size_t outoffset;     // offset in the output data    size_t outoffset;     // offset in the output data
727    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
728    const ValueType* res=0;    const ValueType* res=0;   // Vector storing the answer
729    size_t resoffset=0;    size_t resoffset=0;       // where in the vector to find the answer
730    #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)    #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)
731    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
732    {    {
# Line 925  cout << "############################### Line 734  cout << "###############################
734  #ifdef _OPENMP  #ifdef _OPENMP
735      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
736  #else  #else
737      res=resolveSample(v,0,sample,resoffset);   // 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.
738  #endif  #endif
739  cerr << "-------------------------------- " << endl;  cerr << "-------------------------------- " << endl;
740      outoffset=result->getPointOffset(sample,0);      outoffset=result->getPointOffset(sample,0);
# Line 948  DataLazy::toString() const Line 757  DataLazy::toString() const
757    return oss.str();    return oss.str();
758  }  }
759    
760    
761  void  void
762  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
763  {  {

Legend:
Removed from v.1898  
changed lines
  Added in v.1899

  ViewVC Help
Powered by ViewVC 1.1.26