/[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 1879 by jfenwick, Tue Oct 14 03:54:42 2008 UTC trunk/escript/src/DataLazy.cpp revision 2082 by caltinay, Fri Nov 21 01:46:05 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_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to
52        evaluate the expression.
53    - m_samplesize ~ the number of doubles stored in a sample.
54    
55    When a new node is created, the above values are computed based on the values in the child nodes.
56    Eg: if left requires 4 samples and right requires 6 then left+right requires 7 samples.
57    
58    The resolve method, which produces a DataReady from a DataLazy, does the following:
59    1) Create a DataReady to hold the new result.
60    2) Allocate a vector (v) big enough to hold m_buffsrequired samples.
61    3) For each sample, call resolveSample with v, to get its values and copy them into the result object.
62    
63    (In the case of OMP, multiple samples are resolved in parallel so the vector needs to be larger.)
64    
65    resolveSample returns a Vector* and an offset within that vector where the result is stored.
66    Normally, this would be v, but for identity nodes their internal vector is returned instead.
67    
68    The convention that I use, is that the resolve methods should store their results starting at the offset they are passed.
69    
70    For expressions which evaluate to Constant or Tagged, there is a different evaluation method.
71    The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.
72    
73    To add a new operator you need to do the following (plus anything I might have forgotten):
74    1) Add to the ES_optype.
75    2) determine what opgroup your operation belongs to (X)
76    3) add a string for the op to the end of ES_opstrings
77    4) increase ES_opcount
78    5) add an entry (X) to opgroups
79    6) add an entry to the switch in collapseToReady
80    7) add an entry to resolveX
81    */
82    
83    
84  using namespace std;  using namespace std;
85  using namespace boost;  using namespace boost;
# Line 29  using namespace boost; Line 87  using namespace boost;
87  namespace escript  namespace escript
88  {  {
89    
 const std::string&  
 opToString(ES_optype op);  
   
90  namespace  namespace
91  {  {
92    
93    enum ES_opgroup
94    {
95       G_UNKNOWN,
96       G_IDENTITY,
97       G_BINARY,        // pointwise operations with two arguments
98       G_UNARY,     // pointwise operations with one argument
99       G_NP1OUT,        // non-pointwise op with one output
100       G_TENSORPROD     // general tensor product
101    };
102    
103    
104    
105    
106    string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
107                "sin","cos","tan",
108                "asin","acos","atan","sinh","cosh","tanh","erf",
109                "asinh","acosh","atanh",
110                "log10","log","sign","abs","neg","pos","exp","sqrt",
111                "1/","where>0","where<0","where>=0","where<=0",
112                "symmetric","nonsymmetric",
113                "prod"};
114    int ES_opcount=36;
115    ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
116                G_UNARY,G_UNARY,G_UNARY, //10
117                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17
118                G_UNARY,G_UNARY,G_UNARY,                    // 20
119                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 28
120                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,            // 33
121                G_NP1OUT,G_NP1OUT,
122                G_TENSORPROD};
123    inline
124    ES_opgroup
125    getOpgroup(ES_optype op)
126    {
127      return opgroups[op];
128    }
129    
130  // return the FunctionSpace of the result of "left op right"  // return the FunctionSpace of the result of "left op right"
131  FunctionSpace  FunctionSpace
132  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 136  resultFS(DataAbstract_ptr left, DataAbst
136      // 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
137      // programming error exception.      // programming error exception.
138    
139      FunctionSpace l=left->getFunctionSpace();
140      if (left->getFunctionSpace()!=right->getFunctionSpace())    FunctionSpace r=right->getFunctionSpace();
141      {    if (l!=r)
142          throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");    {
143      }      if (r.probeInterpolation(l))
144      return left->getFunctionSpace();      {
145        return l;
146        }
147        if (l.probeInterpolation(r))
148        {
149        return r;
150        }
151        throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
152      }
153      return l;
154  }  }
155    
156  // return the shape of the result of "left op right"  // return the shape of the result of "left op right"
157    // the shapes resulting from tensor product are more complex to compute so are worked out elsewhere
158  DataTypes::ShapeType  DataTypes::ShapeType
159  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
160  {  {
161      if (left->getShape()!=right->getShape())      if (left->getShape()!=right->getShape())
162      {      {
163          throw DataException("Shapes not the same - shapes must match for lazy data.");        if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
164          {
165            throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
166          }
167          if (left->getRank()==0)   // we need to allow scalar * anything
168          {
169            return right->getShape();
170          }
171          if (right->getRank()==0)
172          {
173            return left->getShape();
174          }
175          throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
176      }      }
177      return left->getShape();      return left->getShape();
178  }  }
179    
180  size_t  // determine the output shape for the general tensor product operation
181  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  // the additional parameters return information required later for the product
182    // the majority of this code is copy pasted from C_General_Tensor_Product
183    DataTypes::ShapeType
184    GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
185  {  {
186     switch(op)      
187     {    // Get rank and shape of inputs
188  //   case IDENTITY: return left->getLength();    int rank0 = left->getRank();
189     case ADD:    // the length is preserved in these ops    int rank1 = right->getRank();
190     case SUB:    const DataTypes::ShapeType& shape0 = left->getShape();
191     case MUL:    const DataTypes::ShapeType& shape1 = right->getShape();
192     case DIV: return left->getLength();  
193     default:    // Prepare for the loops of the product and verify compatibility of shapes
194      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");    int start0=0, start1=0;
195      if (transpose == 0)       {}
196      else if (transpose == 1)  { start0 = axis_offset; }
197      else if (transpose == 2)  { start1 = rank1-axis_offset; }
198      else              { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
199    
200    
201      // Adjust the shapes for transpose
202      DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather
203      DataTypes::ShapeType tmpShape1(rank1);    // than using push_back
204      for (int i=0; i<rank0; i++)   { tmpShape0[i]=shape0[(i+start0)%rank0]; }
205      for (int i=0; i<rank1; i++)   { tmpShape1[i]=shape1[(i+start1)%rank1]; }
206    
207      // Prepare for the loops of the product
208      SL=1, SM=1, SR=1;
209      for (int i=0; i<rank0-axis_offset; i++)   {
210        SL *= tmpShape0[i];
211      }
212      for (int i=rank0-axis_offset; i<rank0; i++)   {
213        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
214          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
215        }
216        SM *= tmpShape0[i];
217      }
218      for (int i=axis_offset; i<rank1; i++)     {
219        SR *= tmpShape1[i];
220      }
221    
222     }    // Define the shape of the output (rank of shape is the sum of the loop ranges below)
223      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);  
224      {         // block to limit the scope of out_index
225         int out_index=0;
226         for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
227         for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
228      }
229      return shape2;
230  }  }
231    
232    
233    // determine the number of points in the result of "left op right"
234    // note that determining the resultLength for G_TENSORPROD is more complex and will not be processed here
235    // size_t
236    // resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
237    // {
238    //    switch (getOpgroup(op))
239    //    {
240    //    case G_BINARY: return left->getLength();
241    //    case G_UNARY: return left->getLength();
242    //    case G_NP1OUT: return left->getLength();
243    //    default:
244    //  throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
245    //    }
246    // }
247    
248    // determine the number of samples requires to evaluate an expression combining left and right
249    // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
250    // The same goes for G_TENSORPROD
251  int  int
252  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
253  {  {
254     switch(op)     switch(getOpgroup(op))
255     {     {
256     case IDENTITY: return 0;     case G_IDENTITY: return 1;
257     case ADD:    // the length is preserved in these ops     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
258     case SUB:     case G_UNARY: return max(left->getBuffsRequired(),1);
259     case MUL:     case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
260     case DIV: return max(left->getBuffsRequired(),right->getBuffsRequired());     case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
261     default:     default:
262      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
263     }     }
264  }  }
265    
 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)+".");  
 //    }  
 //  
 // }  
266    
267  }   // end anonymous namespace  }   // end anonymous namespace
268    
269    
270    
271    // Return a string representing the operation
272  const std::string&  const std::string&
273  opToString(ES_optype op)  opToString(ES_optype op)
274  {  {
# Line 135  opToString(ES_optype op) Line 282  opToString(ES_optype op)
282    
283  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
284      : parent(p->getFunctionSpace(),p->getShape()),      : parent(p->getFunctionSpace(),p->getShape()),
285      m_op(IDENTITY)      m_op(IDENTITY),
286        m_axis_offset(0),
287        m_transpose(0),
288        m_SL(0), m_SM(0), m_SR(0)
289  {  {
290     if (p->isLazy())     if (p->isLazy())
291     {     {
     // TODO: fix this.   We could make the new node a copy of p?  
292      // I don't want identity of Lazy.      // I don't want identity of Lazy.
293      // Question: Why would that be so bad?      // Question: Why would that be so bad?
294      // 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 297  DataLazy::DataLazy(DataAbstract_ptr p)
297     else     else
298     {     {
299      m_id=dynamic_pointer_cast<DataReady>(p);      m_id=dynamic_pointer_cast<DataReady>(p);
300        if(p->isConstant()) {m_readytype='C';}
301        else if(p->isExpanded()) {m_readytype='E';}
302        else if (p->isTagged()) {m_readytype='T';}
303        else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}
304     }     }
305     m_length=p->getLength();     m_buffsRequired=1;
    m_buffsRequired=0;  
306     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
307       m_maxsamplesize=m_samplesize;
308  cout << "(1)Lazy created with " << m_samplesize << endl;  cout << "(1)Lazy created with " << m_samplesize << endl;
309  }  }
310    
311  DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)  
312      : parent(resultFS(left,right,op), resultShape(left,right,op)),  
313      m_left(left),  
314      m_right(right),  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
315      m_op(op)      : parent(left->getFunctionSpace(),left->getShape()),
316        m_op(op),
317        m_axis_offset(0),
318        m_transpose(0),
319        m_SL(0), m_SM(0), m_SR(0)
320  {  {
321     m_length=resultLength(m_left,m_right,m_op);     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
322       {
323        throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
324       }
325    
326       DataLazy_ptr lleft;
327       if (!left->isLazy())
328       {
329        lleft=DataLazy_ptr(new DataLazy(left));
330       }
331       else
332       {
333        lleft=dynamic_pointer_cast<DataLazy>(left);
334       }
335       m_readytype=lleft->m_readytype;
336       m_left=lleft;
337       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
338     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
339     m_buffsRequired=calcBuffs(m_left, m_right, m_op);     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
 cout << "(2)Lazy created with " << m_samplesize << endl;  
340  }  }
341    
342    
343    // In this constructor we need to consider interpolation
344  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
345      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
346      m_op(op)      m_op(op),
347        m_SL(0), m_SM(0), m_SR(0)
348  {  {
349     if (left->isLazy())     if ((getOpgroup(op)!=G_BINARY))
350       {
351        throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
352       }
353    
354       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
355       {
356        FunctionSpace fs=getFunctionSpace();
357        Data ltemp(left);
358        Data tmp(ltemp,fs);
359        left=tmp.borrowDataPtr();
360       }
361       if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
362       {
363        Data tmp(Data(right),getFunctionSpace());
364        right=tmp.borrowDataPtr();
365       }
366       left->operandCheck(*right);
367    
368       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
369     {     {
370      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
371     }     }
# Line 187  DataLazy::DataLazy(DataAbstract_ptr left Line 381  DataLazy::DataLazy(DataAbstract_ptr left
381     {     {
382      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
383     }     }
384       char lt=m_left->m_readytype;
385     m_length=resultLength(m_left,m_right,m_op);     char rt=m_right->m_readytype;
386       if (lt=='E' || rt=='E')
387       {
388        m_readytype='E';
389       }
390       else if (lt=='T' || rt=='T')
391       {
392        m_readytype='T';
393       }
394       else
395       {
396        m_readytype='C';
397       }
398     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
399       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
400     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
401  cout << "(3)Lazy created with " << m_samplesize << endl;  cout << "(3)Lazy created with " << m_samplesize << endl;
402  }  }
403    
404    DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
405        : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
406        m_op(op),
407        m_axis_offset(axis_offset),
408        m_transpose(transpose)
409    {
410       if ((getOpgroup(op)!=G_TENSORPROD))
411       {
412        throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
413       }
414       if ((transpose>2) || (transpose<0))
415       {
416        throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
417       }
418       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
419       {
420        FunctionSpace fs=getFunctionSpace();
421        Data ltemp(left);
422        Data tmp(ltemp,fs);
423        left=tmp.borrowDataPtr();
424       }
425       if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
426       {
427        Data tmp(Data(right),getFunctionSpace());
428        right=tmp.borrowDataPtr();
429       }
430       left->operandCheck(*right);
431    
432       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
433       {
434        m_left=dynamic_pointer_cast<DataLazy>(left);
435       }
436       else
437       {
438        m_left=DataLazy_ptr(new DataLazy(left));
439       }
440       if (right->isLazy())
441       {
442        m_right=dynamic_pointer_cast<DataLazy>(right);
443       }
444       else
445       {
446        m_right=DataLazy_ptr(new DataLazy(right));
447       }
448       char lt=m_left->m_readytype;
449       char rt=m_right->m_readytype;
450       if (lt=='E' || rt=='E')
451       {
452        m_readytype='E';
453       }
454       else if (lt=='T' || rt=='T')
455       {
456        m_readytype='T';
457       }
458       else
459       {
460        m_readytype='C';
461       }
462       m_samplesize=getNumDPPSample()*getNoValues();
463       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
464       m_buffsRequired=calcBuffs(m_left, m_right,m_op);
465    cout << "(4)Lazy created with " << m_samplesize << endl;
466    }
467    
468    
469  DataLazy::~DataLazy()  DataLazy::~DataLazy()
470  {  {
# Line 206  DataLazy::getBuffsRequired() const Line 477  DataLazy::getBuffsRequired() const
477      return m_buffsRequired;      return m_buffsRequired;
478  }  }
479    
480  void  
481  DataLazy::resolveSample(ValueType& v,int sampleNo,  size_t offset ) const  size_t
482    DataLazy::getMaxSampleSize() const
483    {
484        return m_maxsamplesize;
485    }
486    
487    /*
488      \brief Evaluates the expression using methods on Data.
489      This does the work for the collapse method.
490      For reasons of efficiency do not call this method on DataExpanded nodes.
491    */
492    DataReady_ptr
493    DataLazy::collapseToReady()
494  {  {
495    if (m_op==IDENTITY)   // copy the contents into the vector    if (m_readytype=='E')
496      { // this is more an efficiency concern than anything else
497        throw DataException("Programmer Error - do not use collapse on Expanded data.");
498      }
499      if (m_op==IDENTITY)
500    {    {
501  cout << "Begin ID" << endl;      return m_id;
 cout << "dpps=" << getNumDPPSample() << " novals=" << getNoValues() << endl;  
     const ValueType& vec=m_id->getVector();  
     size_t srcOffset=m_id->getPointOffset(sampleNo, 0);  
 cout << "v.size()=" << v.size() << " vec=" << vec.size() << endl;  
     for (size_t i=0;i<m_samplesize;++i,++srcOffset,++offset)  
     {  
 cout << "Trying offset=" << offset << " srcOffset=" << srcOffset << endl;  
     v[offset]=vec[srcOffset];    
     }  
 cout << "End ID" << endl;  
     return;  
502    }    }
503    size_t rightoffset=offset+m_samplesize;    DataReady_ptr pleft=m_left->collapseToReady();
504    m_left->resolveSample(v,sampleNo,offset);    Data left(pleft);
505    m_right->resolveSample(v,sampleNo,rightoffset);    Data right;
506  //  for (int i=0;i<getNumDPPSample();++i)    if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
507    {    {
508      switch(m_op)      right=Data(m_right->collapseToReady());
509      {    }
510      case ADD:       // since these are pointwise ops, pretend each sample is one point    Data result;
511      tensor_binary_operation(m_samplesize,&(v[offset]),&(v[rightoffset]),&(v[offset]),plus<double>());    switch(m_op)
512      {
513        case ADD:
514        result=left+right;
515      break;      break;
516      case SUB:            case SUB:      
517      tensor_binary_operation(m_samplesize,&(v[offset]),&(v[rightoffset]),&(v[offset]),minus<double>());      result=left-right;
518      break;      break;
519      case MUL:            case MUL:      
520      tensor_binary_operation(m_samplesize,&(v[offset]),&(v[rightoffset]),&(v[offset]),multiplies<double>());      result=left*right;
521      break;      break;
522      case DIV:            case DIV:      
523      tensor_binary_operation(m_samplesize,&(v[offset]),&(v[rightoffset]),&(v[offset]),divides<double>());      result=left/right;
524        break;
525        case SIN:
526        result=left.sin();  
527        break;
528        case COS:
529        result=left.cos();
530        break;
531        case TAN:
532        result=left.tan();
533        break;
534        case ASIN:
535        result=left.asin();
536        break;
537        case ACOS:
538        result=left.acos();
539        break;
540        case ATAN:
541        result=left.atan();
542        break;
543        case SINH:
544        result=left.sinh();
545        break;
546        case COSH:
547        result=left.cosh();
548        break;
549        case TANH:
550        result=left.tanh();
551        break;
552        case ERF:
553        result=left.erf();
554        break;
555       case ASINH:
556        result=left.asinh();
557        break;
558       case ACOSH:
559        result=left.acosh();
560        break;
561       case ATANH:
562        result=left.atanh();
563        break;
564        case LOG10:
565        result=left.log10();
566        break;
567        case LOG:
568        result=left.log();
569        break;
570        case SIGN:
571        result=left.sign();
572        break;
573        case ABS:
574        result=left.abs();
575        break;
576        case NEG:
577        result=left.neg();
578        break;
579        case POS:
580        // it doesn't mean anything for delayed.
581        // it will just trigger a deep copy of the lazy object
582        throw DataException("Programmer error - POS not supported for lazy data.");
583        break;
584        case EXP:
585        result=left.exp();
586        break;
587        case SQRT:
588        result=left.sqrt();
589        break;
590        case RECIP:
591        result=left.oneOver();
592        break;
593        case GZ:
594        result=left.wherePositive();
595        break;
596        case LZ:
597        result=left.whereNegative();
598        break;
599        case GEZ:
600        result=left.whereNonNegative();
601        break;
602        case LEZ:
603        result=left.whereNonPositive();
604        break;
605        case SYM:
606        result=left.symmetric();
607        break;
608        case NSYM:
609        result=left.nonsymmetric();
610        break;
611        case PROD:
612        result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
613        break;
614        default:
615        throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
616      }
617      return result.borrowReadyPtr();
618    }
619    
620    /*
621       \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
622       This method uses the original methods on the Data class to evaluate the expressions.
623       For this reason, it should not be used on DataExpanded instances. (To do so would defeat
624       the purpose of using DataLazy in the first place).
625    */
626    void
627    DataLazy::collapse()
628    {
629      if (m_op==IDENTITY)
630      {
631        return;
632      }
633      if (m_readytype=='E')
634      { // this is more an efficiency concern than anything else
635        throw DataException("Programmer Error - do not use collapse on Expanded data.");
636      }
637      m_id=collapseToReady();
638      m_op=IDENTITY;
639    }
640    
641    /*
642      \brief Compute the value of the expression (unary operation) for the given sample.
643      \return Vector which stores the value of the subexpression for the given sample.
644      \param v A vector to store intermediate results.
645      \param offset Index in v to begin storing results.
646      \param sampleNo Sample number to evaluate.
647      \param roffset (output parameter) the offset in the return vector where the result begins.
648    
649      The return value will be an existing vector so do not deallocate it.
650      If the result is stored in v it should be stored at the offset given.
651      Everything from offset to the end of v should be considered available for this method to use.
652    */
653    DataTypes::ValueType*
654    DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
655    {
656        // we assume that any collapsing has been done before we get here
657        // since we only have one argument we don't need to think about only
658        // processing single points.
659      if (m_readytype!='E')
660      {
661        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
662      }
663      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
664      const double* left=&((*vleft)[roffset]);
665      double* result=&(v[offset]);
666      roffset=offset;
667      switch (m_op)
668      {
669        case SIN:  
670        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
671        break;
672        case COS:
673        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
674        break;
675        case TAN:
676        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
677        break;
678        case ASIN:
679        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
680        break;
681        case ACOS:
682        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
683        break;
684        case ATAN:
685        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
686        break;
687        case SINH:
688        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
689        break;
690        case COSH:
691        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
692        break;
693        case TANH:
694        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
695        break;
696        case ERF:
697    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
698        throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
699    #else
700        tensor_unary_operation(m_samplesize, left, result, ::erf);
701        break;
702    #endif
703       case ASINH:
704    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
705        tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
706    #else
707        tensor_unary_operation(m_samplesize, left, result, ::asinh);
708    #endif  
709        break;
710       case ACOSH:
711    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
712        tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
713    #else
714        tensor_unary_operation(m_samplesize, left, result, ::acosh);
715    #endif  
716        break;
717       case ATANH:
718    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
719        tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
720    #else
721        tensor_unary_operation(m_samplesize, left, result, ::atanh);
722    #endif  
723        break;
724        case LOG10:
725        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
726        break;
727        case LOG:
728        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
729        break;
730        case SIGN:
731        tensor_unary_operation(m_samplesize, left, result, escript::fsign);
732        break;
733        case ABS:
734        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
735        break;
736        case NEG:
737        tensor_unary_operation(m_samplesize, left, result, negate<double>());
738        break;
739        case POS:
740        // it doesn't mean anything for delayed.
741        // it will just trigger a deep copy of the lazy object
742        throw DataException("Programmer error - POS not supported for lazy data.");
743        break;
744        case EXP:
745        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
746        break;
747        case SQRT:
748        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
749        break;
750        case RECIP:
751        tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
752        break;
753        case GZ:
754        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
755        break;
756        case LZ:
757        tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
758        break;
759        case GEZ:
760        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
761        break;
762        case LEZ:
763        tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
764        break;
765    
766        default:
767        throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
768      }
769      return &v;
770    }
771    
772    
773    /*
774      \brief Compute the value of the expression (unary operation) for the given sample.
775      \return Vector which stores the value of the subexpression for the given sample.
776      \param v A vector to store intermediate results.
777      \param offset Index in v to begin storing results.
778      \param sampleNo Sample number to evaluate.
779      \param roffset (output parameter) the offset in the return vector where the result begins.
780    
781      The return value will be an existing vector so do not deallocate it.
782      If the result is stored in v it should be stored at the offset given.
783      Everything from offset to the end of v should be considered available for this method to use.
784    */
785    DataTypes::ValueType*
786    DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
787    {
788        // we assume that any collapsing has been done before we get here
789        // since we only have one argument we don't need to think about only
790        // processing single points.
791      if (m_readytype!='E')
792      {
793        throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
794      }
795        // since we can't write the result over the input, we need a result offset further along
796      size_t subroffset=roffset+m_samplesize;
797      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);
798      roffset=offset;
799      switch (m_op)
800      {
801        case SYM:
802        DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
803        break;
804        case NSYM:
805        DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
806        break;
807        default:
808        throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
809      }
810      return &v;
811    }
812    
813    
814    
815    
816    #define PROC_OP(TYPE,X)                               \
817        for (int i=0;i<steps;++i,resultp+=resultStep) \
818        { \
819           tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
820           lroffset+=leftStep; \
821           rroffset+=rightStep; \
822        }
823    
824    /*
825      \brief Compute the value of the expression (binary operation) for the given sample.
826      \return Vector which stores the value of the subexpression for the given sample.
827      \param v A vector to store intermediate results.
828      \param offset Index in v to begin storing results.
829      \param sampleNo Sample number to evaluate.
830      \param roffset (output parameter) the offset in the return vector where the result begins.
831    
832      The return value will be an existing vector so do not deallocate it.
833      If the result is stored in v it should be stored at the offset given.
834      Everything from offset to the end of v should be considered available for this method to use.
835    */
836    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
837    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
838    // If both children are expanded, then we can process them in a single operation (we treat
839    // the whole sample as one big datapoint.
840    // If one of the children is not expanded, then we need to treat each point in the sample
841    // individually.
842    // There is an additional complication when scalar operations are considered.
843    // For example, 2+Vector.
844    // In this case each double within the point is treated individually
845    DataTypes::ValueType*
846    DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
847    {
848    cout << "Resolve binary: " << toString() << endl;
849    
850      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
851        // first work out which of the children are expanded
852      bool leftExp=(m_left->m_readytype=='E');
853      bool rightExp=(m_right->m_readytype=='E');
854      bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step?
855      int steps=(bigloops?1:getNumDPPSample());
856      size_t chunksize=(bigloops? m_samplesize : getNoValues());    // if bigloops, pretend the whole sample is a datapoint
857      if (m_left->getRank()!=m_right->getRank())    // need to deal with scalar * ? ops
858      {
859        EsysAssert((m_left->getRank()==0) || (m_right->getRank()==0), "Error - Ranks must match unless one is 0.");
860        steps=getNumDPPSample()*max(m_left->getNoValues(),m_right->getNoValues());
861        chunksize=1;    // for scalar
862      }    
863      int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
864      int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
865      int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0
866        // Get the values of sub-expressions
867      const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
868      const ValueType* right=m_right->resolveSample(v,offset+m_samplesize,sampleNo,rroffset); // Note
869        // the right child starts further along.
870      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
871      switch(m_op)
872      {
873        case ADD:
874            PROC_OP(NO_ARG,plus<double>());
875        break;
876        case SUB:
877        PROC_OP(NO_ARG,minus<double>());
878        break;
879        case MUL:
880        PROC_OP(NO_ARG,multiplies<double>());
881        break;
882        case DIV:
883        PROC_OP(NO_ARG,divides<double>());
884        break;
885        case POW:
886           PROC_OP(double (double,double),::pow);
887        break;
888        default:
889        throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
890      }
891      roffset=offset;  
892      return &v;
893    }
894    
895    
896    /*
897      \brief Compute the value of the expression (tensor product) for the given sample.
898      \return Vector which stores the value of the subexpression for the given sample.
899      \param v A vector to store intermediate results.
900      \param offset Index in v to begin storing results.
901      \param sampleNo Sample number to evaluate.
902      \param roffset (output parameter) the offset in the return vector where the result begins.
903    
904      The return value will be an existing vector so do not deallocate it.
905      If the result is stored in v it should be stored at the offset given.
906      Everything from offset to the end of v should be considered available for this method to use.
907    */
908    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
909    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
910    // unlike the other resolve helpers, we must treat these datapoints separately.
911    DataTypes::ValueType*
912    DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
913    {
914    cout << "Resolve TensorProduct: " << toString() << endl;
915    
916      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
917        // first work out which of the children are expanded
918      bool leftExp=(m_left->m_readytype=='E');
919      bool rightExp=(m_right->m_readytype=='E');
920      int steps=getNumDPPSample();
921      int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
922      int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
923      int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0
924        // Get the values of sub-expressions (leave a gap of one sample for the result).
925      const ValueType* left=m_left->resolveSample(v,offset+m_samplesize,sampleNo,lroffset);
926      const ValueType* right=m_right->resolveSample(v,offset+2*m_samplesize,sampleNo,rroffset);
927      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
928      switch(m_op)
929      {
930        case PROD:
931        for (int i=0;i<steps;++i,resultp+=resultStep)
932        {
933              const double *ptr_0 = &((*left)[lroffset]);
934              const double *ptr_1 = &((*right)[rroffset]);
935              matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
936          lroffset+=leftStep;
937          rroffset+=rightStep;
938        }
939      break;      break;
940      default:      default:
941      throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
942      }
943      roffset=offset;
944      return &v;
945    }
946    
947    
948    
949    /*
950      \brief Compute the value of the expression for the given sample.
951      \return Vector which stores the value of the subexpression for the given sample.
952      \param v A vector to store intermediate results.
953      \param offset Index in v to begin storing results.
954      \param sampleNo Sample number to evaluate.
955      \param roffset (output parameter) the offset in the return vector where the result begins.
956    
957      The return value will be an existing vector so do not deallocate it.
958    */
959    // the vector and the offset are a place where the method could write its data if it wishes
960    // it is not obligated to do so. For example, if it has its own storage already, it can use that.
961    // Hence the return value to indicate where the data is actually stored.
962    // Regardless, the storage should be assumed to be used, even if it isn't.
963    
964    // the roffset is the offset within the returned vector where the data begins
965    const DataTypes::ValueType*
966    DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
967    {
968    cout << "Resolve sample " << toString() << endl;
969        // collapse so we have a 'E' node or an IDENTITY for some other type
970      if (m_readytype!='E' && m_op!=IDENTITY)
971      {
972        collapse();
973      }
974      if (m_op==IDENTITY)  
975      {
976        const ValueType& vec=m_id->getVector();
977        if (m_readytype=='C')
978        {
979        roffset=0;
980        return &(vec);
981      }      }
982        roffset=m_id->getPointOffset(sampleNo, 0);
983        return &(vec);
984      }
985      if (m_readytype!='E')
986      {
987        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
988      }
989      switch (getOpgroup(m_op))
990      {
991      case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);
992      case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
993      case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
994      case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
995      default:
996        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
997    }    }
998  }  }
999    
1000    
1001    // To simplify the memory management, all threads operate on one large vector, rather than one each.
1002    // Each sample is evaluated independently and copied into the result DataExpanded.
1003  DataReady_ptr  DataReady_ptr
1004  DataLazy::resolve()  DataLazy::resolve()
1005  {  {
   // This is broken!     We need to have a buffer per thread!  
   // so the allocation of v needs to move inside the loop somehow  
1006    
1007  cout << "Sample size=" << m_samplesize << endl;  cout << "Sample size=" << m_samplesize << endl;
1008  cout << "Buffers=" << m_buffsRequired << endl;  cout << "Buffers=" << m_buffsRequired << endl;
1009    
1010    ValueType v(m_samplesize*max(1,m_buffsRequired));    if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally
1011      {
1012        collapse();
1013      }
1014      if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
1015      {
1016        return m_id;
1017      }
1018        // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
1019      size_t threadbuffersize=m_maxsamplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough
1020        // storage to evaluate its expression
1021      int numthreads=1;
1022    #ifdef _OPENMP
1023      numthreads=getNumberOfThreads();
1024    #endif
1025      ValueType v(numthreads*threadbuffersize);
1026  cout << "Buffer created with size=" << v.size() << endl;  cout << "Buffer created with size=" << v.size() << endl;
1027    ValueType dummy(getNoValues());    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
   DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),dummy);  
1028    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVector();
1029    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1030    int sample;    int sample;
1031    #pragma omp parallel for private(sample) schedule(static)    size_t outoffset;     // offset in the output data
1032    for (sample=0;sample<getNumSamples();++sample)    int totalsamples=getNumSamples();
1033      const ValueType* res=0;   // Vector storing the answer
1034      size_t resoffset=0;       // where in the vector to find the answer
1035      #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
1036      for (sample=0;sample<totalsamples;++sample)
1037    {    {
1038      resolveSample(v,sample,0);  cout << "################################# " << sample << endl;
1039      for (int i=0;i<m_samplesize;++i)    // copy values into the output vector  #ifdef _OPENMP
1040        res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
1041    #else
1042        res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.
1043    #endif
1044    cerr << "-------------------------------- " << endl;
1045        outoffset=result->getPointOffset(sample,0);
1046    cerr << "offset=" << outoffset << endl;
1047        for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector
1048      {      {
1049      resvec[i]=v[i];      resvec[outoffset]=(*res)[resoffset];
1050      }      }
1051    cerr << "*********************************" << endl;
1052    }    }
1053    return resptr;    return resptr;
1054  }  }
# Line 280  cout << "Buffer created with size=" << v Line 1056  cout << "Buffer created with size=" << v
1056  std::string  std::string
1057  DataLazy::toString() const  DataLazy::toString() const
1058  {  {
1059    return "Lazy evaluation object. No details available.";    ostringstream oss;
1060      oss << "Lazy Data:";
1061      intoString(oss);
1062      return oss.str();
1063    }
1064    
1065    
1066    void
1067    DataLazy::intoString(ostringstream& oss) const
1068    {
1069      switch (getOpgroup(m_op))
1070      {
1071      case G_IDENTITY:
1072        if (m_id->isExpanded())
1073        {
1074           oss << "E";
1075        }
1076        else if (m_id->isTagged())
1077        {
1078          oss << "T";
1079        }
1080        else if (m_id->isConstant())
1081        {
1082          oss << "C";
1083        }
1084        else
1085        {
1086          oss << "?";
1087        }
1088        oss << '@' << m_id.get();
1089        break;
1090      case G_BINARY:
1091        oss << '(';
1092        m_left->intoString(oss);
1093        oss << ' ' << opToString(m_op) << ' ';
1094        m_right->intoString(oss);
1095        oss << ')';
1096        break;
1097      case G_UNARY:
1098      case G_NP1OUT:
1099        oss << opToString(m_op) << '(';
1100        m_left->intoString(oss);
1101        oss << ')';
1102        break;
1103      case G_TENSORPROD:
1104        oss << opToString(m_op) << '(';
1105        m_left->intoString(oss);
1106        oss << ", ";
1107        m_right->intoString(oss);
1108        oss << ')';
1109        break;
1110      default:
1111        oss << "UNKNOWN";
1112      }
1113  }  }
1114    
 // 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.  
1115  DataAbstract*  DataAbstract*
1116  DataLazy::deepCopy()  DataLazy::deepCopy()
1117  {  {
1118    if (m_op==IDENTITY)    switch (getOpgroup(m_op))
1119    {    {
1120      return new DataLazy(m_left);    // we don't need to copy the child here    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1121      case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1122      case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1123      case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1124      case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1125      default:
1126        throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1127    }    }
   return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);  
1128  }  }
1129    
1130    
1131    // There is no single, natural interpretation of getLength on DataLazy.
1132    // Instances of DataReady can look at the size of their vectors.
1133    // For lazy though, it could be the size the data would be if it were resolved;
1134    // or it could be some function of the lengths of the DataReady instances which
1135    // form part of the expression.
1136    // Rather than have people making assumptions, I have disabled the method.
1137  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1138  DataLazy::getLength() const  DataLazy::getLength() const
1139  {  {
1140    return m_length;    throw DataException("getLength() does not make sense for lazy data.");
1141  }  }
1142    
1143    
# Line 309  DataLazy::getSlice(const DataTypes::Regi Line 1147  DataLazy::getSlice(const DataTypes::Regi
1147    throw DataException("getSlice - not implemented for Lazy objects.");    throw DataException("getSlice - not implemented for Lazy objects.");
1148  }  }
1149    
1150    
1151    // To do this we need to rely on our child nodes
1152  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1153  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
1154                   int dataPointNo) const                   int dataPointNo)
1155  {  {
1156    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");    if (m_op==IDENTITY)
1157      {
1158        return m_id->getPointOffset(sampleNo,dataPointNo);
1159      }
1160      if (m_readytype!='E')
1161      {
1162        collapse();
1163        return m_id->getPointOffset(sampleNo,dataPointNo);
1164      }
1165      // at this point we do not have an identity node and the expression will be Expanded
1166      // so we only need to know which child to ask
1167      if (m_left->m_readytype=='E')
1168      {
1169        return m_left->getPointOffset(sampleNo,dataPointNo);
1170      }
1171      else
1172      {
1173        return m_right->getPointOffset(sampleNo,dataPointNo);
1174      }
1175  }  }
1176    
1177  // // The startOffset is where to write results in the output vector v  // To do this we need to rely on our child nodes
1178  // void  DataTypes::ValueType::size_type
1179  // DataLazy::processSample(ValueType& v, int sampleNo, size_t startOffset)  DataLazy::getPointOffset(int sampleNo,
1180  // {                   int dataPointNo) const
1181  //     m_left.processSample(v,sampleNo,startOffset);  {
1182  //     m_right.processSample(v,sampleNo,startOffset+getSampleSize());    if (m_op==IDENTITY)
1183  //     int i;    {
1184  //     #pragma omp parallel for private(i) schedule(static)      return m_id->getPointOffset(sampleNo,dataPointNo);
1185  //     for (i=0;i<getSampleSize();++i)    }
1186  //     {    if (m_readytype=='E')
1187  //      performOp(v,startOffset+i*m_pointsize,ES_optype,m_samplesize);    {
1188  //     }      // at this point we do not have an identity node and the expression will be Expanded
1189  // }      // so we only need to know which child to ask
1190        if (m_left->m_readytype=='E')
1191        {
1192        return m_left->getPointOffset(sampleNo,dataPointNo);
1193        }
1194        else
1195        {
1196        return m_right->getPointOffset(sampleNo,dataPointNo);
1197        }
1198      }
1199      if (m_readytype=='C')
1200      {
1201        return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
1202      }
1203      throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1204    }
1205    
1206    // It would seem that DataTagged will need to be treated differently since even after setting all tags
1207    // to zero, all the tags from all the DataTags would be in the result.
1208    // However since they all have the same value (0) whether they are there or not should not matter.
1209    // So I have decided that for all types this method will create a constant 0.
1210    // It can be promoted up as required.
1211    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
1212    // but we can deal with that if it arrises.
1213    void
1214    DataLazy::setToZero()
1215    {
1216      DataTypes::ValueType v(getNoValues(),0);
1217      m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1218      m_op=IDENTITY;
1219      m_right.reset();  
1220      m_left.reset();
1221      m_readytype='C';
1222      m_buffsRequired=1;
1223    }
1224    
1225  }   // end namespace  }   // end namespace

Legend:
Removed from v.1879  
changed lines
  Added in v.2082

  ViewVC Help
Powered by ViewVC 1.1.26