/[escript]/trunk/escript/src/DataLazy.cpp
ViewVC logotype

Diff of /trunk/escript/src/DataLazy.cpp

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

branches/schroedinger/escript/src/DataLazy.cpp revision 1898 by jfenwick, Mon Oct 20 01:20:18 2008 UTC branches/schroedinger_upto1946/escript/src/DataLazy.cpp revision 1993 by phornby, Fri Nov 7 04:52:15 2008 UTC
# Line 26  Line 26 
26  #include "DataTypes.h"  #include "DataTypes.h"
27  #include "Data.h"  #include "Data.h"
28  #include "UnaryFuncs.h"     // for escript::fsign  #include "UnaryFuncs.h"     // for escript::fsign
29    #include "Utils.h"
30    
31    /*
32    How does DataLazy work?
33    ~~~~~~~~~~~~~~~~~~~~~~~
34    
35    Each instance represents a single operation on one or two other DataLazy instances. These arguments are normally
36    denoted left and right.
37    
38    A special operation, IDENTITY, stores an instance of DataReady in the m_id member.
39    This means that all "internal" nodes in the structure are instances of DataLazy.
40    
41    Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
42    Note that IDENITY is not considered a unary operation.
43    
44    I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
45    It must however form a DAG (directed acyclic graph).
46    I will refer to individual DataLazy objects with the structure as nodes.
47    
48    Each node also stores:
49    - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
50        evaluated.
51    - m_length ~ how many values would be stored in the answer if the expression was evaluated.
52    - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to
53        evaluate the expression.
54    - m_samplesize ~ the number of doubles stored in a sample.
55    
56    When a new node is created, the above values are computed based on the values in the child nodes.
57    Eg: if left requires 4 samples and right requires 6 then left+right requires 7 samples.
58    
59    The resolve method, which produces a DataReady from a DataLazy, does the following:
60    1) Create a DataReady to hold the new result.
61    2) Allocate a vector (v) big enough to hold m_buffsrequired samples.
62    3) For each sample, call resolveSample with v, to get its values and copy them into the result object.
63    
64    (In the case of OMP, multiple samples are resolved in parallel so the vector needs to be larger.)
65    
66    resolveSample returns a Vector* and an offset within that vector where the result is stored.
67    Normally, this would be v, but for identity nodes their internal vector is returned instead.
68    
69    The convention that I use, is that the resolve methods should store their results starting at the offset they are passed.
70    
71    For expressions which evaluate to Constant or Tagged, there is a different evaluation method.
72    The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.
73    */
74    
75    
76  using namespace std;  using namespace std;
77  using namespace boost;  using namespace boost;
# Line 33  using namespace boost; Line 79  using namespace boost;
79  namespace escript  namespace escript
80  {  {
81    
 const std::string&  
 opToString(ES_optype op);  
   
82  namespace  namespace
83  {  {
84    
   
   
85  enum ES_opgroup  enum ES_opgroup
86  {  {
87     G_UNKNOWN,     G_UNKNOWN,
88     G_IDENTITY,     G_IDENTITY,
89     G_BINARY,     G_BINARY,        // pointwise operations with two arguments
90     G_UNARY     G_UNARY      // pointwise operations with one argument
91  };  };
92    
93    
94    
95    
96  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","sin","cos","tan",  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
97                "sin","cos","tan",
98              "asin","acos","atan","sinh","cosh","tanh","erf",              "asin","acos","atan","sinh","cosh","tanh","erf",
99              "asinh","acosh","atanh",              "asinh","acosh","atanh",
100              "log10","log","sign","abs","neg","pos","exp","sqrt",              "log10","log","sign","abs","neg","pos","exp","sqrt",
101              "1/","where>0","where<0","where>=0","where<=0"};              "1/","where>0","where<0","where>=0","where<=0"};
102  int ES_opcount=32;  int ES_opcount=33;
103  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY,G_UNARY,G_UNARY,G_UNARY, //9  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
104              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 16              G_UNARY,G_UNARY,G_UNARY, //10
105              G_UNARY,G_UNARY,G_UNARY,                    // 19              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17
106              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 27              G_UNARY,G_UNARY,G_UNARY,                    // 20
107                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 28
108              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};
109  inline  inline
110  ES_opgroup  ES_opgroup
# Line 79  resultFS(DataAbstract_ptr left, DataAbst Line 122  resultFS(DataAbstract_ptr left, DataAbst
122      // that way, if interpolate is required in any other op we can just throw a      // that way, if interpolate is required in any other op we can just throw a
123      // programming error exception.      // programming error exception.
124    
125      FunctionSpace l=left->getFunctionSpace();
126      if (left->getFunctionSpace()!=right->getFunctionSpace())    FunctionSpace r=right->getFunctionSpace();
127      {    if (l!=r)
128          throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");    {
129      }      if (r.probeInterpolation(l))
130      return left->getFunctionSpace();      {
131        return l;
132        }
133        if (l.probeInterpolation(r))
134        {
135        return r;
136        }
137        throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
138      }
139      return l;
140  }  }
141    
142  // return the shape of the result of "left op right"  // return the shape of the result of "left op right"
# Line 93  resultShape(DataAbstract_ptr left, DataA Line 145  resultShape(DataAbstract_ptr left, DataA
145  {  {
146      if (left->getShape()!=right->getShape())      if (left->getShape()!=right->getShape())
147      {      {
148          throw DataException("Shapes not the same - shapes must match for lazy data.");        if (getOpgroup(op)!=G_BINARY)
149          {
150            throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
151          }
152          if (left->getRank()==0)   // we need to allow scalar * anything
153          {
154            return right->getShape();
155          }
156          if (right->getRank()==0)
157          {
158            return left->getShape();
159          }
160          throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
161      }      }
162      return left->getShape();      return left->getShape();
163  }  }
164    
165    // determine the number of points in the result of "left op right"
166  size_t  size_t
167  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
168  {  {
# Line 110  resultLength(DataAbstract_ptr left, Data Line 175  resultLength(DataAbstract_ptr left, Data
175     }     }
176  }  }
177    
178    // determine the number of samples requires to evaluate an expression combining left and right
179  int  int
180  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
181  {  {
# Line 123  calcBuffs(const DataLazy_ptr& left, cons Line 189  calcBuffs(const DataLazy_ptr& left, cons
189     }     }
190  }  }
191    
192    
193  }   // end anonymous namespace  }   // end anonymous namespace
194    
195    
196    
197    // Return a string representing the operation
198  const std::string&  const std::string&
199  opToString(ES_optype op)  opToString(ES_optype op)
200  {  {
# Line 143  DataLazy::DataLazy(DataAbstract_ptr p) Line 212  DataLazy::DataLazy(DataAbstract_ptr p)
212  {  {
213     if (p->isLazy())     if (p->isLazy())
214     {     {
     // TODO: fix this.   We could make the new node a copy of p?  
215      // I don't want identity of Lazy.      // I don't want identity of Lazy.
216      // Question: Why would that be so bad?      // Question: Why would that be so bad?
217      // 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 231  DataLazy::DataLazy(DataAbstract_ptr p)
231  cout << "(1)Lazy created with " << m_samplesize << endl;  cout << "(1)Lazy created with " << m_samplesize << endl;
232  }  }
233    
234    
235    
236    
237  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
238      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
239      m_op(op)      m_op(op)
# Line 188  DataLazy::DataLazy(DataAbstract_ptr left Line 259  DataLazy::DataLazy(DataAbstract_ptr left
259  }  }
260    
261    
262  // DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)  // In this constructor we need to consider interpolation
 //  : 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;  
 // }  
   
263  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
264      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
265      m_op(op)      m_op(op)
# Line 212  DataLazy::DataLazy(DataAbstract_ptr left Line 268  DataLazy::DataLazy(DataAbstract_ptr left
268     {     {
269      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.");
270     }     }
271     if (left->isLazy())  
272       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
273       {
274        FunctionSpace fs=getFunctionSpace();
275        Data ltemp(left);
276        Data tmp(ltemp,fs);
277        left=tmp.borrowDataPtr();
278       }
279       if (getFunctionSpace()!=right->getFunctionSpace())   // left needs to be interpolated
280       {
281        Data tmp(Data(right),getFunctionSpace());
282        right=tmp.borrowDataPtr();
283       }
284       left->operandCheck(*right);
285    
286       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
287     {     {
288      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
289     }     }
# Line 243  DataLazy::DataLazy(DataAbstract_ptr left Line 314  DataLazy::DataLazy(DataAbstract_ptr left
314      m_readytype='C';      m_readytype='C';
315     }     }
316     m_length=resultLength(m_left,m_right,m_op);     m_length=resultLength(m_left,m_right,m_op);
317     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();    
318     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
319  cout << "(3)Lazy created with " << m_samplesize << endl;  cout << "(3)Lazy created with " << m_samplesize << endl;
320  }  }
# Line 261  DataLazy::getBuffsRequired() const Line 332  DataLazy::getBuffsRequired() const
332  }  }
333    
334    
335    /*
336      \brief Evaluates the expression using methods on Data.
337      This does the work for the collapse method.
338      For reasons of efficiency do not call this method on DataExpanded nodes.
339    */
340  DataReady_ptr  DataReady_ptr
341  DataLazy::collapseToReady()  DataLazy::collapseToReady()
342  {  {
# Line 380  DataLazy::collapseToReady() Line 456  DataLazy::collapseToReady()
456    return result.borrowReadyPtr();    return result.borrowReadyPtr();
457  }  }
458    
459    /*
460       \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
461       This method uses the original methods on the Data class to evaluate the expressions.
462       For this reason, it should not be used on DataExpanded instances. (To do so would defeat
463       the purpose of using DataLazy in the first place).
464    */
465  void  void
466  DataLazy::collapse()  DataLazy::collapse()
467  {  {
# Line 395  DataLazy::collapse() Line 477  DataLazy::collapse()
477    m_op=IDENTITY;    m_op=IDENTITY;
478  }  }
479    
480    /*
481      \brief Compute the value of the expression (binary operation) for the given sample.
482      \return Vector which stores the value of the subexpression for the given sample.
483      \param v A vector to store intermediate results.
484      \param offset Index in v to begin storing results.
485      \param sampleNo Sample number to evaluate.
486      \param roffset (output parameter) the offset in the return vector where the result begins.
487    
488      The return value will be an existing vector so do not deallocate it.
489      If the result is stored in v it should be stored at the offset given.
490      Everything from offset to the end of v should be considered available for this method to use.
491    */
492  DataTypes::ValueType*  DataTypes::ValueType*
493  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
494  {  {
# Line 412  DataLazy::resolveUnary(ValueType& v, siz Line 506  DataLazy::resolveUnary(ValueType& v, siz
506    switch (m_op)    switch (m_op)
507    {    {
508      case SIN:        case SIN:  
509      tensor_unary_operation(m_samplesize, left, result, ::sin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
510      break;      break;
511      case COS:      case COS:
512      tensor_unary_operation(m_samplesize, left, result, ::cos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
513      break;      break;
514      case TAN:      case TAN:
515      tensor_unary_operation(m_samplesize, left, result, ::tan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
516      break;      break;
517      case ASIN:      case ASIN:
518      tensor_unary_operation(m_samplesize, left, result, ::asin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
519      break;      break;
520      case ACOS:      case ACOS:
521      tensor_unary_operation(m_samplesize, left, result, ::acos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
522      break;      break;
523      case ATAN:      case ATAN:
524      tensor_unary_operation(m_samplesize, left, result, ::atan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
525      break;      break;
526      case SINH:      case SINH:
527      tensor_unary_operation(m_samplesize, left, result, ::sinh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
528      break;      break;
529      case COSH:      case COSH:
530      tensor_unary_operation(m_samplesize, left, result, ::cosh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
531      break;      break;
532      case TANH:      case TANH:
533      tensor_unary_operation(m_samplesize, left, result, ::tanh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
534      break;      break;
535      case ERF:      case ERF:
536  #ifdef _WIN32  #ifdef _WIN32
# Line 467  DataLazy::resolveUnary(ValueType& v, siz Line 561  DataLazy::resolveUnary(ValueType& v, siz
561  #endif    #endif  
562      break;      break;
563      case LOG10:      case LOG10:
564      tensor_unary_operation(m_samplesize, left, result, ::log10);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
565      break;      break;
566      case LOG:      case LOG:
567      tensor_unary_operation(m_samplesize, left, result, ::log);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
568      break;      break;
569      case SIGN:      case SIGN:
570      tensor_unary_operation(m_samplesize, left, result, escript::fsign);      tensor_unary_operation(m_samplesize, left, result, escript::fsign);
571      break;      break;
572      case ABS:      case ABS:
573      tensor_unary_operation(m_samplesize, left, result, ::fabs);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
574      break;      break;
575      case NEG:      case NEG:
576      tensor_unary_operation(m_samplesize, left, result, negate<double>());      tensor_unary_operation(m_samplesize, left, result, negate<double>());
# Line 487  DataLazy::resolveUnary(ValueType& v, siz Line 581  DataLazy::resolveUnary(ValueType& v, siz
581      throw DataException("Programmer error - POS not supported for lazy data.");      throw DataException("Programmer error - POS not supported for lazy data.");
582      break;      break;
583      case EXP:      case EXP:
584      tensor_unary_operation(m_samplesize, left, result, ::exp);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
585      break;      break;
586      case SQRT:      case SQRT:
587      tensor_unary_operation(m_samplesize, left, result, ::sqrt);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
588      break;      break;
589      case RECIP:      case RECIP:
590      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
# Line 516  DataLazy::resolveUnary(ValueType& v, siz Line 610  DataLazy::resolveUnary(ValueType& v, siz
610    
611    
612    
 // 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;  
 // }  
613    
614  #define PROC_OP(X) \  
615      for (int i=0;i<steps;++i,resultp+=getNoValues()) \  #define PROC_OP(TYPE,X)                               \
616        for (int i=0;i<steps;++i,resultp+=resultStep) \
617      { \      { \
618  cout << "Step#" << i << " chunk=" << chunksize << endl; \         tensor_binary_operation##TYPE(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
619  cout << left[0] << left[1] << left[2] << endl; \         lroffset+=leftStep; \
620  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; \  
621      }      }
622    
623    /*
624      \brief Compute the value of the expression (binary operation) for the given sample.
625      \return Vector which stores the value of the subexpression for the given sample.
626      \param v A vector to store intermediate results.
627      \param offset Index in v to begin storing results.
628      \param sampleNo Sample number to evaluate.
629      \param roffset (output parameter) the offset in the return vector where the result begins.
630    
631      The return value will be an existing vector so do not deallocate it.
632      If the result is stored in v it should be stored at the offset given.
633      Everything from offset to the end of v should be considered available for this method to use.
634    */
635    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
636    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
637    // If both children are expanded, then we can process them in a single operation (we treat
638    // the whole sample as one big datapoint.
639    // If one of the children is not expanded, then we need to treat each point in the sample
640    // individually.
641    // There is an additional complication when scalar operations are considered.
642    // For example, 2+Vector.
643    // In this case each double within the point is treated individually
644  DataTypes::ValueType*  DataTypes::ValueType*
645  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
646  {  {
     // 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.  
   
647  cout << "Resolve binary: " << toString() << endl;  cout << "Resolve binary: " << toString() << endl;
648    
649    size_t lroffset=0, rroffset=0;    size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
650        // first work out which of the children are expanded
651    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
652    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
653    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
654    int steps=(bigloops?1:getNumDPPSample());    int steps=(bigloops?1:getNumDPPSample());
655    size_t chunksize=(bigloops? m_samplesize : getNoValues());    size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint
656    int leftStep=((leftExp && !rightExp)? getNoValues() : 0);    if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops
657    int rightStep=((rightExp && !leftExp)? getNoValues() : 0);    {
658        EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");
659        steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());
660        chunksize=1;    // for scalar
661      }    
662      int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
663      int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
664      int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0
665        // Get the values of sub-expressions
666    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
667    const ValueType* right=m_right->resolveSample(v,offset,sampleNo,rroffset);        const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
668      // now we need to know which args are expanded      // the right child starts further along.
669  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]);  
670    switch(m_op)    switch(m_op)
671    {    {
672      case ADD:      case ADD:
673      for (int i=0;i<steps;++i,resultp+=getNoValues())          PROC_OP(/**/,plus<double>());
674      {      break;
675  cerr << "Step#" << i << " chunk=" << chunksize << endl;      case SUB:
676  cerr << left << "[" << lroffset << "] " << right << "[" << rroffset << "]" << endl;      PROC_OP(/**/,minus<double>());
677         tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, plus<double>());      break;
678         lroffset+=leftStep;      case MUL:
679         rroffset+=rightStep;      PROC_OP(/**/,multiplies<double>());
680  cerr << "left=" << lroffset << " right=" << rroffset << endl;      break;
681      }      case DIV:
682        PROC_OP(/**/,divides<double>());
683        break;
684        case POW:
685           PROC_OP(<double (double,double)>,::pow);
686      break;      break;
 // need to fill in the rest  
687      default:      default:
688      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
689    }    }
690    roffset=offset;    roffset=offset;  
691    return &v;    return &v;
692  }  }
693    
694    
695    
696  // #define PROC_OP(X) \  /*
697  //  for (int i=0;i<steps;++i,resultp+=getNoValues()) \    \brief Compute the value of the expression for the given sample.
698  //  { \    \return Vector which stores the value of the subexpression for the given sample.
699  // cout << "Step#" << i << " chunk=" << chunksize << endl; \    \param v A vector to store intermediate results.
700  // cout << left[0] << left[1] << left[2] << endl; \    \param offset Index in v to begin storing results.
701  // cout << right[0] << right[1] << right[2] << endl; \    \param sampleNo Sample number to evaluate.
702  //     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)+".");  
 //   }  
 // }  
   
   
703    
704      The return value will be an existing vector so do not deallocate it.
705    */
706  // 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
707  // 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.
708  // 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 743  cout << "Resolve sample " << toString()
743  }  }
744    
745    
746    // To simplify the memory management, all threads operate on one large vector, rather than one each.
747    // 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;  
 // }  
   
   
748  DataReady_ptr  DataReady_ptr
749  DataLazy::resolve()  DataLazy::resolve()
750  {  {
# Line 893  DataLazy::resolve() Line 752  DataLazy::resolve()
752  cout << "Sample size=" << m_samplesize << endl;  cout << "Sample size=" << m_samplesize << endl;
753  cout << "Buffers=" << m_buffsRequired << endl;  cout << "Buffers=" << m_buffsRequired << endl;
754    
755    if (m_readytype!='E')    if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally
756    {    {
757      collapse();      collapse();
758    }    }
759    if (m_op==IDENTITY)    if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
760    {    {
761      return m_id;      return m_id;
762    }    }
763      // 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'
764    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
765        // storage to evaluate its expression
766    int numthreads=1;    int numthreads=1;
767  #ifdef _OPENMP  #ifdef _OPENMP
768    numthreads=getNumberOfThreads();    numthreads=getNumberOfThreads();
# Line 916  cout << "Buffer created with size=" << v Line 776  cout << "Buffer created with size=" << v
776    int sample;    int sample;
777    size_t outoffset;     // offset in the output data    size_t outoffset;     // offset in the output data
778    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
779    const ValueType* res=0;    const ValueType* res=0;   // Vector storing the answer
780    size_t resoffset=0;    size_t resoffset=0;       // where in the vector to find the answer
781    #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)    #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)
782    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
783    {    {
# Line 925  cout << "############################### Line 785  cout << "###############################
785  #ifdef _OPENMP  #ifdef _OPENMP
786      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
787  #else  #else
788      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.
789  #endif  #endif
790  cerr << "-------------------------------- " << endl;  cerr << "-------------------------------- " << endl;
791      outoffset=result->getPointOffset(sample,0);      outoffset=result->getPointOffset(sample,0);
# Line 948  DataLazy::toString() const Line 808  DataLazy::toString() const
808    return oss.str();    return oss.str();
809  }  }
810    
811    
812  void  void
813  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
814  {  {
# Line 989  DataLazy::intoString(ostringstream& oss) Line 850  DataLazy::intoString(ostringstream& oss)
850    }    }
851  }  }
852    
 // Note that in this case, deepCopy does not make copies of the leaves.  
 // Hopefully copy on write (or whatever we end up using) will take care of this.  
853  DataAbstract*  DataAbstract*
854  DataLazy::deepCopy()  DataLazy::deepCopy()
855  {  {
856    if (m_op==IDENTITY)    switch (getOpgroup(m_op))
857    {    {
858      return new DataLazy(m_left);    // we don't need to copy the child here    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
859      case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
860      case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
861      default:
862        throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
863    }    }
   return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);  
864  }  }
865    
866    
# Line 1015  DataLazy::getSlice(const DataTypes::Regi Line 877  DataLazy::getSlice(const DataTypes::Regi
877    throw DataException("getSlice - not implemented for Lazy objects.");    throw DataException("getSlice - not implemented for Lazy objects.");
878  }  }
879    
880    
881    // To do this we need to rely on our child nodes
882    DataTypes::ValueType::size_type
883    DataLazy::getPointOffset(int sampleNo,
884                     int dataPointNo)
885    {
886      if (m_op==IDENTITY)
887      {
888        return m_id->getPointOffset(sampleNo,dataPointNo);
889      }
890      if (m_readytype!='E')
891      {
892        collapse();
893        return m_id->getPointOffset(sampleNo,dataPointNo);
894      }
895      // at this point we do not have an identity node and the expression will be Expanded
896      // so we only need to know which child to ask
897      if (m_left->m_readytype=='E')
898      {
899        return m_left->getPointOffset(sampleNo,dataPointNo);
900      }
901      else
902      {
903        return m_right->getPointOffset(sampleNo,dataPointNo);
904      }
905    }
906    
907    // To do this we need to rely on our child nodes
908  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
909  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
910                   int dataPointNo) const                   int dataPointNo) const
911  {  {
912    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");    if (m_op==IDENTITY)
913      {
914        return m_id->getPointOffset(sampleNo,dataPointNo);
915      }
916      if (m_readytype=='E')
917      {
918        // at this point we do not have an identity node and the expression will be Expanded
919        // so we only need to know which child to ask
920        if (m_left->m_readytype=='E')
921        {
922        return m_left->getPointOffset(sampleNo,dataPointNo);
923        }
924        else
925        {
926        return m_right->getPointOffset(sampleNo,dataPointNo);
927        }
928      }
929      if (m_readytype=='C')
930      {
931        return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
932      }
933      throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
934    }
935    
936    // It would seem that DataTagged will need to be treated differently since even after setting all tags
937    // to zero, all the tags from all the DataTags would be in the result.
938    // However since they all have the same value (0) whether they are there or not should not matter.
939    // So I have decided that for all types this method will create a constant 0.
940    // It can be promoted up as required.
941    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
942    // but we can deal with that if it arrises.
943    void
944    DataLazy::setToZero()
945    {
946      DataTypes::ValueType v(getNoValues(),0);
947      m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
948      m_op=IDENTITY;
949      m_right.reset();  
950      m_left.reset();
951      m_readytype='C';
952      m_buffsRequired=1;
953  }  }
954    
955  }   // end namespace  }   // end namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26