/[escript]/branches/clazy/escriptcore/src/DataLazy.cpp
ViewVC logotype

Diff of /branches/clazy/escriptcore/src/DataLazy.cpp

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

branches/schroedinger/escript/src/DataLazy.cpp revision 1884 by jfenwick, Tue Oct 14 04:54:59 2008 UTC branches/schroedinger_upto1946/escript/src/DataLazy.cpp revision 1950 by jfenwick, Thu Oct 30 00:59:34 2008 UTC
# Line 19  Line 19 
19  #ifdef PASO_MPI  #ifdef PASO_MPI
20  #include <mpi.h>  #include <mpi.h>
21  #endif  #endif
22    #ifdef _OPENMP
23    #include <omp.h>
24    #endif
25  #include "FunctionSpace.h"  #include "FunctionSpace.h"
26  #include "DataTypes.h"  #include "DataTypes.h"
27  #include "Data.h"  #include "Data.h"
28    #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 35  opToString(ES_optype op); Line 85  opToString(ES_optype op);
85  namespace  namespace
86  {  {
87    
88    enum ES_opgroup
89    {
90       G_UNKNOWN,
91       G_IDENTITY,
92       G_BINARY,        // pointwise operations with two arguments
93       G_UNARY      // pointwise operations with one argument
94    };
95    
96    
97    
98    
99    string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
100                "sin","cos","tan",
101                "asin","acos","atan","sinh","cosh","tanh","erf",
102                "asinh","acosh","atanh",
103                "log10","log","sign","abs","neg","pos","exp","sqrt",
104                "1/","where>0","where<0","where>=0","where<=0"};
105    int ES_opcount=33;
106    ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
107                G_UNARY,G_UNARY,G_UNARY, //10
108                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17
109                G_UNARY,G_UNARY,G_UNARY,                    // 20
110                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 28
111                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};
112    inline
113    ES_opgroup
114    getOpgroup(ES_optype op)
115    {
116      return opgroups[op];
117    }
118    
119  // return the FunctionSpace of the result of "left op right"  // return the FunctionSpace of the result of "left op right"
120  FunctionSpace  FunctionSpace
121  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
# Line 44  resultFS(DataAbstract_ptr left, DataAbst Line 125  resultFS(DataAbstract_ptr left, DataAbst
125      // 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
126      // programming error exception.      // programming error exception.
127    
128      FunctionSpace l=left->getFunctionSpace();
129      if (left->getFunctionSpace()!=right->getFunctionSpace())    FunctionSpace r=right->getFunctionSpace();
130      {    if (l!=r)
131          throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");    {
132      }      if (r.probeInterpolation(l))
133      return left->getFunctionSpace();      {
134        return l;
135        }
136        if (l.probeInterpolation(r))
137        {
138        return r;
139        }
140        throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
141      }
142      return l;
143  }  }
144    
145  // return the shape of the result of "left op right"  // return the shape of the result of "left op right"
# Line 58  resultShape(DataAbstract_ptr left, DataA Line 148  resultShape(DataAbstract_ptr left, DataA
148  {  {
149      if (left->getShape()!=right->getShape())      if (left->getShape()!=right->getShape())
150      {      {
151          throw DataException("Shapes not the same - shapes must match for lazy data.");        if (getOpgroup(op)!=G_BINARY)
152          {
153            throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
154          }
155          if (left->getRank()==0)   // we need to allow scalar * anything
156          {
157            return right->getShape();
158          }
159          if (right->getRank()==0)
160          {
161            return left->getShape();
162          }
163          throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
164      }      }
165      return left->getShape();      return left->getShape();
166  }  }
167    
168    // determine the number of points in the result of "left op right"
169  size_t  size_t
170  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
171  {  {
172     switch(op)     switch (getOpgroup(op))
173     {     {
174  //   case IDENTITY: return left->getLength();     case G_BINARY: return left->getLength();
175     case ADD:    // the length is preserved in these ops     case G_UNARY: return left->getLength();
    case SUB:  
    case MUL:  
    case DIV: return left->getLength();  
176     default:     default:
177      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
   
178     }     }
179  }  }
180    
181    // determine the number of samples requires to evaluate an expression combining left and right
182  int  int
183  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
184  {  {
185     switch(op)     switch(getOpgroup(op))
186     {     {
187     case IDENTITY: return 0;     case G_IDENTITY: return 1;
188     case ADD:    // the length is preserved in these ops     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
189     case SUB:     case G_UNARY: return max(left->getBuffsRequired(),1);
    case MUL:  
    case DIV: return max(left->getBuffsRequired(),right->getBuffsRequired());  
190     default:     default:
191      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
192     }     }
193  }  }
194    
 string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/"};  
 int ES_opcount=5;  
   
 // void performOp(ValueType& v, int startOffset, ES_optype op, int m_samplesize)  
 // {  
 //    switch(op)  
 //    {  
 //    case ADD:  DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),  
 //      startOffset+m_samplesize,plus<double>());  
 //        break;      
 //    case SUB:  DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),  
 //      startOffset+m_samplesize,minus<double>());  
 //        break;  
 //    case MUL:  DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),  
 //      startOffset+m_samplesize,multiplies<double>());  
 //        break;  
 //    case DIV:  DataMaths::binaryOp(v,getShape(),startOffset,v,getShape(),  
 //      startOffset+m_samplesize,divides<double>());  
 //        break;  
 //    default:  
 //  throw DataException("Programmer Error - attempt to performOp() for operator "+opToString(op)+".");  
 //    }  
 //  
 // }  
195    
196  }   // end anonymous namespace  }   // end anonymous namespace
197    
198    
199    
200    // Return a string representing the operation
201  const std::string&  const std::string&
202  opToString(ES_optype op)  opToString(ES_optype op)
203  {  {
# Line 139  DataLazy::DataLazy(DataAbstract_ptr p) Line 215  DataLazy::DataLazy(DataAbstract_ptr p)
215  {  {
216     if (p->isLazy())     if (p->isLazy())
217     {     {
     // TODO: fix this.   We could make the new node a copy of p?  
218      // I don't want identity of Lazy.      // I don't want identity of Lazy.
219      // Question: Why would that be so bad?      // Question: Why would that be so bad?
220      // 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 148  DataLazy::DataLazy(DataAbstract_ptr p) Line 223  DataLazy::DataLazy(DataAbstract_ptr p)
223     else     else
224     {     {
225      m_id=dynamic_pointer_cast<DataReady>(p);      m_id=dynamic_pointer_cast<DataReady>(p);
226        if(p->isConstant()) {m_readytype='C';}
227        else if(p->isExpanded()) {m_readytype='E';}
228        else if (p->isTagged()) {m_readytype='T';}
229        else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}
230     }     }
231     m_length=p->getLength();     m_length=p->getLength();
232     m_buffsRequired=0;     m_buffsRequired=1;
233     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
234  cout << "(1)Lazy created with " << m_samplesize << endl;  cout << "(1)Lazy created with " << m_samplesize << endl;
235  }  }
236    
237  DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)  
238      : parent(resultFS(left,right,op), resultShape(left,right,op)),  
239      m_left(left),  
240      m_right(right),  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
241        : parent(left->getFunctionSpace(),left->getShape()),
242      m_op(op)      m_op(op)
243  {  {
244     m_length=resultLength(m_left,m_right,m_op);     if (getOpgroup(op)!=G_UNARY)
245       {
246        throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
247       }
248       DataLazy_ptr lleft;
249       if (!left->isLazy())
250       {
251        lleft=DataLazy_ptr(new DataLazy(left));
252       }
253       else
254       {
255        lleft=dynamic_pointer_cast<DataLazy>(left);
256       }
257       m_readytype=lleft->m_readytype;
258       m_length=left->getLength();
259       m_left=lleft;
260       m_buffsRequired=1;
261     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
    m_buffsRequired=calcBuffs(m_left, m_right, m_op);  
 cout << "(2)Lazy created with " << m_samplesize << endl;  
262  }  }
263    
264    
265    // In this constructor we need to consider interpolation
266  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
267      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
268      m_op(op)      m_op(op)
269  {  {
270     if (left->isLazy())     if (getOpgroup(op)!=G_BINARY)
271       {
272        throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
273       }
274    
275       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
276       {
277        FunctionSpace fs=getFunctionSpace();
278        Data ltemp(left);
279        Data tmp(ltemp,fs);
280        left=tmp.borrowDataPtr();
281       }
282       if (getFunctionSpace()!=right->getFunctionSpace())   // left needs to be interpolated
283       {
284        Data tmp(Data(right),getFunctionSpace());
285        right=tmp.borrowDataPtr();
286       }
287       left->operandCheck(*right);
288    
289       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
290     {     {
291      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
292     }     }
# Line 187  DataLazy::DataLazy(DataAbstract_ptr left Line 302  DataLazy::DataLazy(DataAbstract_ptr left
302     {     {
303      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
304     }     }
305       char lt=m_left->m_readytype;
306       char rt=m_right->m_readytype;
307       if (lt=='E' || rt=='E')
308       {
309        m_readytype='E';
310       }
311       else if (lt=='T' || rt=='T')
312       {
313        m_readytype='T';
314       }
315       else
316       {
317        m_readytype='C';
318       }
319     m_length=resultLength(m_left,m_right,m_op);     m_length=resultLength(m_left,m_right,m_op);
320     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();    
321     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
322  cout << "(3)Lazy created with " << m_samplesize << endl;  cout << "(3)Lazy created with " << m_samplesize << endl;
323  }  }
# Line 207  DataLazy::getBuffsRequired() const Line 335  DataLazy::getBuffsRequired() const
335  }  }
336    
337    
338  // the vector and the offset are a place where the method could write its data if it wishes  /*
339  // it is not obligated to do so. For example, if it has its own storage already, it can use that.    \brief Evaluates the expression using methods on Data.
340  // hence the return value to indicate where the data is actually stored.    This does the work for the collapse method.
341  // regardless, the storage should be assumed to be used, even if it isn't.    For reasons of efficiency do not call this method on DataExpanded nodes.
342  const double*  */
343  DataLazy::resolveSample(ValueType& v,int sampleNo,  size_t offset ) const  DataReady_ptr
344    DataLazy::collapseToReady()
345  {  {
346    if (m_op==IDENTITY)   // copy the contents into the vector    if (m_readytype=='E')
347      { // this is more an efficiency concern than anything else
348        throw DataException("Programmer Error - do not use collapse on Expanded data.");
349      }
350      if (m_op==IDENTITY)
351    {    {
352  cout << "Begin ID" << endl;      return m_id;
353  cout << "dpps=" << getNumDPPSample() << " novals=" << getNoValues() << endl;    }
354      const ValueType& vec=m_id->getVector();    DataReady_ptr pleft=m_left->collapseToReady();
355  //     size_t srcOffset=m_id->getPointOffset(sampleNo, 0);    Data left(pleft);
356  // cout << "v.size()=" << v.size() << " vec=" << vec.size() << endl;    Data right;
357  //     for (size_t i=0;i<m_samplesize;++i,++srcOffset,++offset)    if (getOpgroup(m_op)==G_BINARY)
 //     {  
 // cout << "Trying offset=" << offset << " srcOffset=" << srcOffset << endl;  
 //  v[offset]=vec[srcOffset];    
 //     }  
 cout << "End ID - returning offset " << m_id->getPointOffset(sampleNo, 0) << " of vector@" << &vec<<endl;  
     return &(vec[m_id->getPointOffset(sampleNo, 0)]);  
 //     return;  
   }  
 cout << "Begin op";  
   size_t rightoffset=offset+m_samplesize;  
   const double* left=m_left->resolveSample(v,sampleNo,offset);  
   const double* right=m_right->resolveSample(v,sampleNo,rightoffset);  
   double* result=&(v[offset]);  
 cout << "left=" << left << " right=" << right << " result=" << result << endl;  
 //  for (int i=0;i<getNumDPPSample();++i)  
358    {    {
359      switch(m_op)      right=Data(m_right->collapseToReady());
360      {    }
361      case ADD:       // since these are pointwise ops, pretend each sample is one point    Data result;
362      tensor_binary_operation(m_samplesize, left, right, result, plus<double>());    switch(m_op)
363      {
364        case ADD:
365        result=left+right;
366      break;      break;
367      case SUB:            case SUB:      
368      tensor_binary_operation(m_samplesize, left, right, result, minus<double>());      result=left-right;
369      break;      break;
370      case MUL:            case MUL:      
371      tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());      result=left*right;
372      break;      break;
373      case DIV:            case DIV:      
374      tensor_binary_operation(m_samplesize, left, right, result, divides<double>());      result=left/right;
375        break;
376        case SIN:
377        result=left.sin();  
378        break;
379        case COS:
380        result=left.cos();
381        break;
382        case TAN:
383        result=left.tan();
384        break;
385        case ASIN:
386        result=left.asin();
387        break;
388        case ACOS:
389        result=left.acos();
390        break;
391        case ATAN:
392        result=left.atan();
393        break;
394        case SINH:
395        result=left.sinh();
396        break;
397        case COSH:
398        result=left.cosh();
399        break;
400        case TANH:
401        result=left.tanh();
402        break;
403        case ERF:
404        result=left.erf();
405        break;
406       case ASINH:
407        result=left.asinh();
408        break;
409       case ACOSH:
410        result=left.acosh();
411        break;
412       case ATANH:
413        result=left.atanh();
414        break;
415        case LOG10:
416        result=left.log10();
417        break;
418        case LOG:
419        result=left.log();
420        break;
421        case SIGN:
422        result=left.sign();
423        break;
424        case ABS:
425        result=left.abs();
426        break;
427        case NEG:
428        result=left.neg();
429        break;
430        case POS:
431        // it doesn't mean anything for delayed.
432        // it will just trigger a deep copy of the lazy object
433        throw DataException("Programmer error - POS not supported for lazy data.");
434        break;
435        case EXP:
436        result=left.exp();
437        break;
438        case SQRT:
439        result=left.sqrt();
440        break;
441        case RECIP:
442        result=left.oneOver();
443        break;
444        case GZ:
445        result=left.wherePositive();
446        break;
447        case LZ:
448        result=left.whereNegative();
449        break;
450        case GEZ:
451        result=left.whereNonNegative();
452        break;
453        case LEZ:
454        result=left.whereNonPositive();
455      break;      break;
456      default:      default:
457      throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");
458      }
459      return result.borrowReadyPtr();
460    }
461    
462    /*
463       \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
464       This method uses the original methods on the Data class to evaluate the expressions.
465       For this reason, it should not be used on DataExpanded instances. (To do so would defeat
466       the purpose of using DataLazy in the first place).
467    */
468    void
469    DataLazy::collapse()
470    {
471      if (m_op==IDENTITY)
472      {
473        return;
474      }
475      if (m_readytype=='E')
476      { // this is more an efficiency concern than anything else
477        throw DataException("Programmer Error - do not use collapse on Expanded data.");
478      }
479      m_id=collapseToReady();
480      m_op=IDENTITY;
481    }
482    
483    /*
484      \brief Compute the value of the expression (binary operation) for the given sample.
485      \return Vector which stores the value of the subexpression for the given sample.
486      \param v A vector to store intermediate results.
487      \param offset Index in v to begin storing results.
488      \param sampleNo Sample number to evaluate.
489      \param roffset (output parameter) the offset in the return vector where the result begins.
490    
491      The return value will be an existing vector so do not deallocate it.
492      If the result is stored in v it should be stored at the offset given.
493      Everything from offset to the end of v should be considered available for this method to use.
494    */
495    DataTypes::ValueType*
496    DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
497    {
498        // we assume that any collapsing has been done before we get here
499        // since we only have one argument we don't need to think about only
500        // processing single points.
501      if (m_readytype!='E')
502      {
503        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
504      }
505      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
506      const double* left=&((*vleft)[roffset]);
507      double* result=&(v[offset]);
508      roffset=offset;
509      switch (m_op)
510      {
511        case SIN:  
512        tensor_unary_operation(m_samplesize, left, result, ::sin);
513        break;
514        case COS:
515        tensor_unary_operation(m_samplesize, left, result, ::cos);
516        break;
517        case TAN:
518        tensor_unary_operation(m_samplesize, left, result, ::tan);
519        break;
520        case ASIN:
521        tensor_unary_operation(m_samplesize, left, result, ::asin);
522        break;
523        case ACOS:
524        tensor_unary_operation(m_samplesize, left, result, ::acos);
525        break;
526        case ATAN:
527        tensor_unary_operation(m_samplesize, left, result, ::atan);
528        break;
529        case SINH:
530        tensor_unary_operation(m_samplesize, left, result, ::sinh);
531        break;
532        case COSH:
533        tensor_unary_operation(m_samplesize, left, result, ::cosh);
534        break;
535        case TANH:
536        tensor_unary_operation(m_samplesize, left, result, ::tanh);
537        break;
538        case ERF:
539    #ifdef _WIN32
540        throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
541    #else
542        tensor_unary_operation(m_samplesize, left, result, ::erf);
543        break;
544    #endif
545       case ASINH:
546    #ifdef _WIN32
547        tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
548    #else
549        tensor_unary_operation(m_samplesize, left, result, ::asinh);
550    #endif  
551        break;
552       case ACOSH:
553    #ifdef _WIN32
554        tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
555    #else
556        tensor_unary_operation(m_samplesize, left, result, ::acosh);
557    #endif  
558        break;
559       case ATANH:
560    #ifdef _WIN32
561        tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
562    #else
563        tensor_unary_operation(m_samplesize, left, result, ::atanh);
564    #endif  
565        break;
566        case LOG10:
567        tensor_unary_operation(m_samplesize, left, result, ::log10);
568        break;
569        case LOG:
570        tensor_unary_operation(m_samplesize, left, result, ::log);
571        break;
572        case SIGN:
573        tensor_unary_operation(m_samplesize, left, result, escript::fsign);
574        break;
575        case ABS:
576        tensor_unary_operation(m_samplesize, left, result, ::fabs);
577        break;
578        case NEG:
579        tensor_unary_operation(m_samplesize, left, result, negate<double>());
580        break;
581        case POS:
582        // it doesn't mean anything for delayed.
583        // it will just trigger a deep copy of the lazy object
584        throw DataException("Programmer error - POS not supported for lazy data.");
585        break;
586        case EXP:
587        tensor_unary_operation(m_samplesize, left, result, ::exp);
588        break;
589        case SQRT:
590        tensor_unary_operation(m_samplesize, left, result, ::sqrt);
591        break;
592        case RECIP:
593        tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
594        break;
595        case GZ:
596        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
597        break;
598        case LZ:
599        tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
600        break;
601        case GEZ:
602        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
603        break;
604        case LEZ:
605        tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
606        break;
607    
608        default:
609        throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
610      }
611      return &v;
612    }
613    
614    
615    
616    
617    
618    #define PROC_OP(X) \
619        for (int i=0;i<steps;++i,resultp+=resultStep) \
620        { \
621           tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
622           lroffset+=leftStep; \
623           rroffset+=rightStep; \
624        }
625    
626    /*
627      \brief Compute the value of the expression (binary operation) for the given sample.
628      \return Vector which stores the value of the subexpression for the given sample.
629      \param v A vector to store intermediate results.
630      \param offset Index in v to begin storing results.
631      \param sampleNo Sample number to evaluate.
632      \param roffset (output parameter) the offset in the return vector where the result begins.
633    
634      The return value will be an existing vector so do not deallocate it.
635      If the result is stored in v it should be stored at the offset given.
636      Everything from offset to the end of v should be considered available for this method to use.
637    */
638    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
639    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
640    // If both children are expanded, then we can process them in a single operation (we treat
641    // the whole sample as one big datapoint.
642    // If one of the children is not expanded, then we need to treat each point in the sample
643    // individually.
644    // There is an additional complication when scalar operations are considered.
645    // For example, 2+Vector.
646    // In this case each double within the point is treated individually
647    DataTypes::ValueType*
648    DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
649    {
650    cout << "Resolve binary: " << toString() << endl;
651    
652      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
653        // first work out which of the children are expanded
654      bool leftExp=(m_left->m_readytype=='E');
655      bool rightExp=(m_right->m_readytype=='E');
656      bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
657      int steps=(bigloops?1:getNumDPPSample());
658      size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint
659      if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops
660      {
661        EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");
662        steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());
663        chunksize=1;    // for scalar
664      }    
665      int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
666      int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
667      int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0
668        // Get the values of sub-expressions
669      const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
670      const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
671        // the right child starts further along.
672      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
673      switch(m_op)
674      {
675        case ADD:
676        PROC_OP(plus<double>());
677        break;
678        case SUB:
679        PROC_OP(minus<double>());
680        break;
681        case MUL:
682        PROC_OP(multiplies<double>());
683        break;
684        case DIV:
685        PROC_OP(divides<double>());
686        break;
687        case POW:
688        PROC_OP(::pow);
689        break;
690        default:
691        throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
692      }
693      roffset=offset;  
694      return &v;
695    }
696    
697    
698    
699    /*
700      \brief Compute the value of the expression for the given sample.
701      \return Vector which stores the value of the subexpression for the given sample.
702      \param v A vector to store intermediate results.
703      \param offset Index in v to begin storing results.
704      \param sampleNo Sample number to evaluate.
705      \param roffset (output parameter) the offset in the return vector where the result begins.
706    
707      The return value will be an existing vector so do not deallocate it.
708    */
709    // the vector and the offset are a place where the method could write its data if it wishes
710    // it is not obligated to do so. For example, if it has its own storage already, it can use that.
711    // Hence the return value to indicate where the data is actually stored.
712    // Regardless, the storage should be assumed to be used, even if it isn't.
713    
714    // the roffset is the offset within the returned vector where the data begins
715    const DataTypes::ValueType*
716    DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
717    {
718    cout << "Resolve sample " << toString() << endl;
719        // collapse so we have a 'E' node or an IDENTITY for some other type
720      if (m_readytype!='E' && m_op!=IDENTITY)
721      {
722        collapse();
723      }
724      if (m_op==IDENTITY)  
725      {
726        const ValueType& vec=m_id->getVector();
727        if (m_readytype=='C')
728        {
729        roffset=0;
730        return &(vec);
731      }      }
732        roffset=m_id->getPointOffset(sampleNo, 0);
733        return &(vec);
734      }
735      if (m_readytype!='E')
736      {
737        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
738      }
739      switch (getOpgroup(m_op))
740      {
741      case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);
742      case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
743      default:
744        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
745    }    }
   return result;  
746  }  }
747    
748    
749    // To simplify the memory management, all threads operate on one large vector, rather than one each.
750    // Each sample is evaluated independently and copied into the result DataExpanded.
751  DataReady_ptr  DataReady_ptr
752  DataLazy::resolve()  DataLazy::resolve()
753  {  {
   // This is broken!     We need to have a buffer per thread!  
   // so the allocation of v needs to move inside the loop somehow  
754    
755  cout << "Sample size=" << m_samplesize << endl;  cout << "Sample size=" << m_samplesize << endl;
756  cout << "Buffers=" << m_buffsRequired << endl;  cout << "Buffers=" << m_buffsRequired << endl;
757    
758    ValueType v(m_samplesize*(max(1,m_buffsRequired)+1)); // the +1 comes from the fact that I want to have a safe    if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally
759                              // space for the RHS of ops to write to even if they don't    {
760                              // need it.      collapse();
761      }
762      if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
763      {
764        return m_id;
765      }
766        // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
767      size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired));    // Each thread needs to have enough
768        // storage to evaluate its expression
769      int numthreads=1;
770    #ifdef _OPENMP
771      numthreads=getNumberOfThreads();
772      int threadnum=0;
773    #endif
774      ValueType v(numthreads*threadbuffersize);
775  cout << "Buffer created with size=" << v.size() << endl;  cout << "Buffer created with size=" << v.size() << endl;
776    ValueType dummy(getNoValues());    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
   DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),dummy);  
777    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVector();
778    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
779    int sample;    int sample;
780    #pragma omp parallel for private(sample) schedule(static)    size_t outoffset;     // offset in the output data
781    for (sample=0;sample<getNumSamples();++sample)    int totalsamples=getNumSamples();
782      const ValueType* res=0;   // Vector storing the answer
783      size_t resoffset=0;       // where in the vector to find the answer
784      #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)
785      for (sample=0;sample<totalsamples;++sample)
786    {    {
787  cout << "Processing sample#" << sample << endl;  cout << "################################# " << sample << endl;
788      resolveSample(v,sample,0);  #ifdef _OPENMP
789  cout << "Copying#" << sample << endl;      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
790      for (int i=0;i<m_samplesize;++i)    // copy values into the output vector  #else
791        res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.
792    #endif
793    cerr << "-------------------------------- " << endl;
794        outoffset=result->getPointOffset(sample,0);
795    cerr << "offset=" << outoffset << endl;
796        for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector
797      {      {
798      resvec[i]=v[i];      resvec[outoffset]=(*res)[resoffset];
799      }      }
800    cerr << "*********************************" << endl;
801    }    }
802    return resptr;    return resptr;
803  }  }
# Line 294  cout << "Copying#" << sample << endl; Line 805  cout << "Copying#" << sample << endl;
805  std::string  std::string
806  DataLazy::toString() const  DataLazy::toString() const
807  {  {
808    return "Lazy evaluation object. No details available.";    ostringstream oss;
809      oss << "Lazy Data:";
810      intoString(oss);
811      return oss.str();
812    }
813    
814    
815    void
816    DataLazy::intoString(ostringstream& oss) const
817    {
818      switch (getOpgroup(m_op))
819      {
820      case G_IDENTITY:
821        if (m_id->isExpanded())
822        {
823           oss << "E";
824        }
825        else if (m_id->isTagged())
826        {
827          oss << "T";
828        }
829        else if (m_id->isConstant())
830        {
831          oss << "C";
832        }
833        else
834        {
835          oss << "?";
836        }
837        oss << '@' << m_id.get();
838        break;
839      case G_BINARY:
840        oss << '(';
841        m_left->intoString(oss);
842        oss << ' ' << opToString(m_op) << ' ';
843        m_right->intoString(oss);
844        oss << ')';
845        break;
846      case G_UNARY:
847        oss << opToString(m_op) << '(';
848        m_left->intoString(oss);
849        oss << ')';
850        break;
851      default:
852        oss << "UNKNOWN";
853      }
854  }  }
855    
 // 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.  
856  DataAbstract*  DataAbstract*
857  DataLazy::deepCopy()  DataLazy::deepCopy()
858  {  {
859    if (m_op==IDENTITY)    switch (getOpgroup(m_op))
860    {    {
861      return new DataLazy(m_left);    // we don't need to copy the child here    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
862      case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
863      case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
864      default:
865        throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
866    }    }
   return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);  
867  }  }
868    
869    
# Line 323  DataLazy::getSlice(const DataTypes::Regi Line 880  DataLazy::getSlice(const DataTypes::Regi
880    throw DataException("getSlice - not implemented for Lazy objects.");    throw DataException("getSlice - not implemented for Lazy objects.");
881  }  }
882    
883    
884    // To do this we need to rely on our child nodes
885    DataTypes::ValueType::size_type
886    DataLazy::getPointOffset(int sampleNo,
887                     int dataPointNo)
888    {
889      if (m_op==IDENTITY)
890      {
891        return m_id->getPointOffset(sampleNo,dataPointNo);
892      }
893      if (m_readytype!='E')
894      {
895        collapse();
896        return m_id->getPointOffset(sampleNo,dataPointNo);
897      }
898      // at this point we do not have an identity node and the expression will be Expanded
899      // so we only need to know which child to ask
900      if (m_left->m_readytype=='E')
901      {
902        return m_left->getPointOffset(sampleNo,dataPointNo);
903      }
904      else
905      {
906        return m_right->getPointOffset(sampleNo,dataPointNo);
907      }
908    }
909    
910    // To do this we need to rely on our child nodes
911  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
912  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
913                   int dataPointNo) const                   int dataPointNo) const
914  {  {
915    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");    if (m_op==IDENTITY)
916      {
917        return m_id->getPointOffset(sampleNo,dataPointNo);
918      }
919      if (m_readytype=='E')
920      {
921        // at this point we do not have an identity node and the expression will be Expanded
922        // so we only need to know which child to ask
923        if (m_left->m_readytype=='E')
924        {
925        return m_left->getPointOffset(sampleNo,dataPointNo);
926        }
927        else
928        {
929        return m_right->getPointOffset(sampleNo,dataPointNo);
930        }
931      }
932      if (m_readytype=='C')
933      {
934        return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
935      }
936      throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
937    }
938    
939    // It would seem that DataTagged will need to be treated differently since even after setting all tags
940    // to zero, all the tags from all the DataTags would be in the result.
941    // However since they all have the same value (0) whether they are there or not should not matter.
942    // So I have decided that for all types this method will create a constant 0.
943    // It can be promoted up as required.
944    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
945    // but we can deal with that if it arrises.
946    void
947    DataLazy::setToZero()
948    {
949      DataTypes::ValueType v(getNoValues(),0);
950      m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
951      m_op=IDENTITY;
952      m_right.reset();  
953      m_left.reset();
954      m_readytype='C';
955      m_buffsRequired=1;
956  }  }
957    
958  }   // end namespace  }   // end namespace

Legend:
Removed from v.1884  
changed lines
  Added in v.1950

  ViewVC Help
Powered by ViewVC 1.1.26