/[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 1993 by phornby, Fri Nov 7 04:52:15 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 29  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
86    {
87       G_UNKNOWN,
88       G_IDENTITY,
89       G_BINARY,        // pointwise operations with two arguments
90       G_UNARY      // pointwise operations with one argument
91    };
92    
93    
94    
95    
96    string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
97                "sin","cos","tan",
98                "asin","acos","atan","sinh","cosh","tanh","erf",
99                "asinh","acosh","atanh",
100                "log10","log","sign","abs","neg","pos","exp","sqrt",
101                "1/","where>0","where<0","where>=0","where<=0"};
102    int ES_opcount=33;
103    ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
104                G_UNARY,G_UNARY,G_UNARY, //10
105                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17
106                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};
109    inline
110    ES_opgroup
111    getOpgroup(ES_optype op)
112    {
113      return opgroups[op];
114    }
115    
116  // return the FunctionSpace of the result of "left op right"  // return the FunctionSpace of the result of "left op right"
117  FunctionSpace  FunctionSpace
118  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 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 58  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  {  {
169     switch(op)     switch (getOpgroup(op))
170     {     {
171  //   case IDENTITY: return left->getLength();     case G_BINARY: return left->getLength();
172     case ADD:    // the length is preserved in these ops     case G_UNARY: return left->getLength();
    case SUB:  
    case MUL:  
    case DIV: return left->getLength();  
173     default:     default:
174      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
   
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  {  {
182     switch(op)     switch(getOpgroup(op))
183     {     {
184     case IDENTITY: return 0;     case G_IDENTITY: return 1;
185     case ADD:    // the length is preserved in these ops     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
186     case SUB:     case G_UNARY: return max(left->getBuffsRequired(),1);
    case MUL:  
    case DIV: return max(left->getBuffsRequired(),right->getBuffsRequired());  
187     default:     default:
188      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
189     }     }
190  }  }
191    
 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)+".");  
 //    }  
 //  
 // }  
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 139  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 148  DataLazy::DataLazy(DataAbstract_ptr p) Line 220  DataLazy::DataLazy(DataAbstract_ptr p)
220     else     else
221     {     {
222      m_id=dynamic_pointer_cast<DataReady>(p);      m_id=dynamic_pointer_cast<DataReady>(p);
223        if(p->isConstant()) {m_readytype='C';}
224        else if(p->isExpanded()) {m_readytype='E';}
225        else if (p->isTagged()) {m_readytype='T';}
226        else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}
227     }     }
228     m_length=p->getLength();     m_length=p->getLength();
229     m_buffsRequired=0;     m_buffsRequired=1;
230     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
231  cout << "(1)Lazy created with " << m_samplesize << endl;  cout << "(1)Lazy created with " << m_samplesize << endl;
232  }  }
233    
234  DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)  
235      : parent(resultFS(left,right,op), resultShape(left,right,op)),  
236      m_left(left),  
237      m_right(right),  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
238        : parent(left->getFunctionSpace(),left->getShape()),
239      m_op(op)      m_op(op)
240  {  {
241     m_length=resultLength(m_left,m_right,m_op);     if (getOpgroup(op)!=G_UNARY)
242       {
243        throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
244       }
245       DataLazy_ptr lleft;
246       if (!left->isLazy())
247       {
248        lleft=DataLazy_ptr(new DataLazy(left));
249       }
250       else
251       {
252        lleft=dynamic_pointer_cast<DataLazy>(left);
253       }
254       m_readytype=lleft->m_readytype;
255       m_length=left->getLength();
256       m_left=lleft;
257       m_buffsRequired=1;
258     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;  
259  }  }
260    
261    
262    // In this constructor we need to consider interpolation
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)
266  {  {
267     if (left->isLazy())     if (getOpgroup(op)!=G_BINARY)
268       {
269        throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
270       }
271    
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 187  DataLazy::DataLazy(DataAbstract_ptr left Line 299  DataLazy::DataLazy(DataAbstract_ptr left
299     {     {
300      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
301     }     }
302       char lt=m_left->m_readytype;
303       char rt=m_right->m_readytype;
304       if (lt=='E' || rt=='E')
305       {
306        m_readytype='E';
307       }
308       else if (lt=='T' || rt=='T')
309       {
310        m_readytype='T';
311       }
312       else
313       {
314        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 207  DataLazy::getBuffsRequired() const Line 332  DataLazy::getBuffsRequired() const
332  }  }
333    
334    
335  // the vector and the offset are a place where the method could write its data if it wishes  /*
336  // 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.
337  // hence the return value to indicate where the data is actually stored.    This does the work for the collapse method.
338  // 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.
339  const double*  */
340  DataLazy::resolveSample(ValueType& v,int sampleNo,  size_t offset ) const  DataReady_ptr
341    DataLazy::collapseToReady()
342  {  {
343    if (m_op==IDENTITY)   // copy the contents into the vector    if (m_readytype=='E')
344      { // this is more an efficiency concern than anything else
345        throw DataException("Programmer Error - do not use collapse on Expanded data.");
346      }
347      if (m_op==IDENTITY)
348    {    {
349  cout << "Begin ID" << endl;      return m_id;
350  cout << "dpps=" << getNumDPPSample() << " novals=" << getNoValues() << endl;    }
351      const ValueType& vec=m_id->getVector();    DataReady_ptr pleft=m_left->collapseToReady();
352  //     size_t srcOffset=m_id->getPointOffset(sampleNo, 0);    Data left(pleft);
353  // cout << "v.size()=" << v.size() << " vec=" << vec.size() << endl;    Data right;
354  //     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)  
355    {    {
356      switch(m_op)      right=Data(m_right->collapseToReady());
357      {    }
358      case ADD:       // since these are pointwise ops, pretend each sample is one point    Data result;
359      tensor_binary_operation(m_samplesize, left, right, result, plus<double>());    switch(m_op)
360      {
361        case ADD:
362        result=left+right;
363      break;      break;
364      case SUB:            case SUB:      
365      tensor_binary_operation(m_samplesize, left, right, result, minus<double>());      result=left-right;
366      break;      break;
367      case MUL:            case MUL:      
368      tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());      result=left*right;
369      break;      break;
370      case DIV:            case DIV:      
371      tensor_binary_operation(m_samplesize, left, right, result, divides<double>());      result=left/right;
372        break;
373        case SIN:
374        result=left.sin();  
375        break;
376        case COS:
377        result=left.cos();
378        break;
379        case TAN:
380        result=left.tan();
381        break;
382        case ASIN:
383        result=left.asin();
384        break;
385        case ACOS:
386        result=left.acos();
387        break;
388        case ATAN:
389        result=left.atan();
390        break;
391        case SINH:
392        result=left.sinh();
393        break;
394        case COSH:
395        result=left.cosh();
396        break;
397        case TANH:
398        result=left.tanh();
399        break;
400        case ERF:
401        result=left.erf();
402        break;
403       case ASINH:
404        result=left.asinh();
405        break;
406       case ACOSH:
407        result=left.acosh();
408        break;
409       case ATANH:
410        result=left.atanh();
411        break;
412        case LOG10:
413        result=left.log10();
414        break;
415        case LOG:
416        result=left.log();
417        break;
418        case SIGN:
419        result=left.sign();
420        break;
421        case ABS:
422        result=left.abs();
423        break;
424        case NEG:
425        result=left.neg();
426        break;
427        case POS:
428        // it doesn't mean anything for delayed.
429        // it will just trigger a deep copy of the lazy object
430        throw DataException("Programmer error - POS not supported for lazy data.");
431        break;
432        case EXP:
433        result=left.exp();
434        break;
435        case SQRT:
436        result=left.sqrt();
437        break;
438        case RECIP:
439        result=left.oneOver();
440        break;
441        case GZ:
442        result=left.wherePositive();
443        break;
444        case LZ:
445        result=left.whereNegative();
446        break;
447        case GEZ:
448        result=left.whereNonNegative();
449        break;
450        case LEZ:
451        result=left.whereNonPositive();
452      break;      break;
453      default:      default:
454      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)+".");
455      }
456      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
466    DataLazy::collapse()
467    {
468      if (m_op==IDENTITY)
469      {
470        return;
471      }
472      if (m_readytype=='E')
473      { // this is more an efficiency concern than anything else
474        throw DataException("Programmer Error - do not use collapse on Expanded data.");
475      }
476      m_id=collapseToReady();
477      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*
493    DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
494    {
495        // we assume that any collapsing has been done before we get here
496        // since we only have one argument we don't need to think about only
497        // processing single points.
498      if (m_readytype!='E')
499      {
500        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
501      }
502      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
503      const double* left=&((*vleft)[roffset]);
504      double* result=&(v[offset]);
505      roffset=offset;
506      switch (m_op)
507      {
508        case SIN:  
509        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
510        break;
511        case COS:
512        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
513        break;
514        case TAN:
515        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
516        break;
517        case ASIN:
518        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
519        break;
520        case ACOS:
521        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
522        break;
523        case ATAN:
524        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
525        break;
526        case SINH:
527        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
528        break;
529        case COSH:
530        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
531        break;
532        case TANH:
533        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
534        break;
535        case ERF:
536    #ifdef _WIN32
537        throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
538    #else
539        tensor_unary_operation(m_samplesize, left, result, ::erf);
540        break;
541    #endif
542       case ASINH:
543    #ifdef _WIN32
544        tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
545    #else
546        tensor_unary_operation(m_samplesize, left, result, ::asinh);
547    #endif  
548        break;
549       case ACOSH:
550    #ifdef _WIN32
551        tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
552    #else
553        tensor_unary_operation(m_samplesize, left, result, ::acosh);
554    #endif  
555        break;
556       case ATANH:
557    #ifdef _WIN32
558        tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
559    #else
560        tensor_unary_operation(m_samplesize, left, result, ::atanh);
561    #endif  
562        break;
563        case LOG10:
564        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
565        break;
566        case LOG:
567        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
568        break;
569        case SIGN:
570        tensor_unary_operation(m_samplesize, left, result, escript::fsign);
571        break;
572        case ABS:
573        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
574        break;
575        case NEG:
576        tensor_unary_operation(m_samplesize, left, result, negate<double>());
577        break;
578        case POS:
579        // it doesn't mean anything for delayed.
580        // it will just trigger a deep copy of the lazy object
581        throw DataException("Programmer error - POS not supported for lazy data.");
582        break;
583        case EXP:
584        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
585        break;
586        case SQRT:
587        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
588        break;
589        case RECIP:
590        tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
591        break;
592        case GZ:
593        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
594        break;
595        case LZ:
596        tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
597        break;
598        case GEZ:
599        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
600        break;
601        case LEZ:
602        tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
603        break;
604    
605        default:
606        throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
607      }
608      return &v;
609    }
610    
611    
612    
613    
614    
615    #define PROC_OP(TYPE,X)                               \
616        for (int i=0;i<steps;++i,resultp+=resultStep) \
617        { \
618           tensor_binary_operation##TYPE(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
619           lroffset+=leftStep; \
620           rroffset+=rightStep; \
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*
645    DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
646    {
647    cout << "Resolve binary: " << toString() << endl;
648    
649      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');
652      bool rightExp=(m_right->m_readytype=='E');
653      bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
654      int steps=(bigloops?1:getNumDPPSample());
655      size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint
656      if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops
657      {
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);
667      const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
668        // the right child starts further along.
669      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
670      switch(m_op)
671      {
672        case ADD:
673            PROC_OP(/**/,plus<double>());
674        break;
675        case SUB:
676        PROC_OP(/**/,minus<double>());
677        break;
678        case MUL:
679        PROC_OP(/**/,multiplies<double>());
680        break;
681        case DIV:
682        PROC_OP(/**/,divides<double>());
683        break;
684        case POW:
685           PROC_OP(<double (double,double)>,::pow);
686        break;
687        default:
688        throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
689      }
690      roffset=offset;  
691      return &v;
692    }
693    
694    
695    
696    /*
697      \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      \param v A vector to store intermediate results.
700      \param offset Index in v to begin storing results.
701      \param sampleNo Sample number to evaluate.
702      \param roffset (output parameter) the offset in the return vector where the result begins.
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
707    // 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.
709    // Regardless, the storage should be assumed to be used, even if it isn't.
710    
711    // the roffset is the offset within the returned vector where the data begins
712    const DataTypes::ValueType*
713    DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
714    {
715    cout << "Resolve sample " << toString() << endl;
716        // collapse so we have a 'E' node or an IDENTITY for some other type
717      if (m_readytype!='E' && m_op!=IDENTITY)
718      {
719        collapse();
720      }
721      if (m_op==IDENTITY)  
722      {
723        const ValueType& vec=m_id->getVector();
724        if (m_readytype=='C')
725        {
726        roffset=0;
727        return &(vec);
728      }      }
729        roffset=m_id->getPointOffset(sampleNo, 0);
730        return &(vec);
731      }
732      if (m_readytype!='E')
733      {
734        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
735      }
736      switch (getOpgroup(m_op))
737      {
738      case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);
739      case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
740      default:
741        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
742    }    }
   return result;  
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.
748  DataReady_ptr  DataReady_ptr
749  DataLazy::resolve()  DataLazy::resolve()
750  {  {
   // This is broken!     We need to have a buffer per thread!  
   // so the allocation of v needs to move inside the loop somehow  
751    
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    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
756                              // space for the RHS of ops to write to even if they don't    {
757                              // need it.      collapse();
758      }
759      if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
760      {
761        return m_id;
762      }
763        // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
764      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;
767    #ifdef _OPENMP
768      numthreads=getNumberOfThreads();
769      int threadnum=0;
770    #endif
771      ValueType v(numthreads*threadbuffersize);
772  cout << "Buffer created with size=" << v.size() << endl;  cout << "Buffer created with size=" << v.size() << endl;
773    ValueType dummy(getNoValues());    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
   DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),dummy);  
774    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVector();
775    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
776    int sample;    int sample;
777    #pragma omp parallel for private(sample) schedule(static)    size_t outoffset;     // offset in the output data
778    for (sample=0;sample<getNumSamples();++sample)    int totalsamples=getNumSamples();
779      const ValueType* res=0;   // Vector storing the answer
780      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)
782      for (sample=0;sample<totalsamples;++sample)
783    {    {
784  cout << "Processing sample#" << sample << endl;  cout << "################################# " << sample << endl;
785      resolveSample(v,sample,0);  #ifdef _OPENMP
786  cout << "Copying#" << sample << endl;      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
787      for (int i=0;i<m_samplesize;++i)    // copy values into the output vector  #else
788        res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.
789    #endif
790    cerr << "-------------------------------- " << endl;
791        outoffset=result->getPointOffset(sample,0);
792    cerr << "offset=" << outoffset << endl;
793        for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector
794      {      {
795      resvec[i]=v[i];      resvec[outoffset]=(*res)[resoffset];
796      }      }
797    cerr << "*********************************" << endl;
798    }    }
799    return resptr;    return resptr;
800  }  }
# Line 294  cout << "Copying#" << sample << endl; Line 802  cout << "Copying#" << sample << endl;
802  std::string  std::string
803  DataLazy::toString() const  DataLazy::toString() const
804  {  {
805    return "Lazy evaluation object. No details available.";    ostringstream oss;
806      oss << "Lazy Data:";
807      intoString(oss);
808      return oss.str();
809    }
810    
811    
812    void
813    DataLazy::intoString(ostringstream& oss) const
814    {
815      switch (getOpgroup(m_op))
816      {
817      case G_IDENTITY:
818        if (m_id->isExpanded())
819        {
820           oss << "E";
821        }
822        else if (m_id->isTagged())
823        {
824          oss << "T";
825        }
826        else if (m_id->isConstant())
827        {
828          oss << "C";
829        }
830        else
831        {
832          oss << "?";
833        }
834        oss << '@' << m_id.get();
835        break;
836      case G_BINARY:
837        oss << '(';
838        m_left->intoString(oss);
839        oss << ' ' << opToString(m_op) << ' ';
840        m_right->intoString(oss);
841        oss << ')';
842        break;
843      case G_UNARY:
844        oss << opToString(m_op) << '(';
845        m_left->intoString(oss);
846        oss << ')';
847        break;
848      default:
849        oss << "UNKNOWN";
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 323  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.1884  
changed lines
  Added in v.1993

  ViewVC Help
Powered by ViewVC 1.1.26