/[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 1888 by jfenwick, Wed Oct 15 04:00:21 2008 UTC trunk/escript/src/DataLazy.cpp revision 2066 by jfenwick, Thu Nov 20 05:31:33 2008 UTC
# Line 26  Line 26 
26  #include "DataTypes.h"  #include "DataTypes.h"
27  #include "Data.h"  #include "Data.h"
28  #include "UnaryFuncs.h"     // for escript::fsign  #include "UnaryFuncs.h"     // for escript::fsign
29    #include "Utils.h"
30    
31    /*
32    How does DataLazy work?
33    ~~~~~~~~~~~~~~~~~~~~~~~
34    
35    Each instance represents a single operation on one or two other DataLazy instances. These arguments are normally
36    denoted left and right.
37    
38    A special operation, IDENTITY, stores an instance of DataReady in the m_id member.
39    This means that all "internal" nodes in the structure are instances of DataLazy.
40    
41    Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
42    Note that IDENITY is not considered a unary operation.
43    
44    I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
45    It must however form a DAG (directed acyclic graph).
46    I will refer to individual DataLazy objects with the structure as nodes.
47    
48    Each node also stores:
49    - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
50        evaluated.
51    - m_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 33  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  enum ES_opgroup
94  {  {
95     G_UNKNOWN,     G_UNKNOWN,
96     G_IDENTITY,     G_IDENTITY,
97     G_BINARY,     G_BINARY,        // pointwise operations with two arguments
98     G_UNARY     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","+","-","*","/","sin","cos","tan",  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
107                "sin","cos","tan",
108              "asin","acos","atan","sinh","cosh","tanh","erf",              "asin","acos","atan","sinh","cosh","tanh","erf",
109              "asinh","acosh","atanh",              "asinh","acosh","atanh",
110              "log10","log","sign","abs","neg","pos","exp","sqrt",              "log10","log","sign","abs","neg","pos","exp","sqrt",
111              "1/","where>0","where<0","where>=0","where<=0"};              "1/","where>0","where<0","where>=0","where<=0",
112  int ES_opcount=32;              "symmetric","nonsymmetric",
113  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY,G_UNARY,G_UNARY,G_UNARY, //9              "prod"};
114              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 16  int ES_opcount=36;
115              G_UNARY,G_UNARY,G_UNARY,                    // 19  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
116              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 27              G_UNARY,G_UNARY,G_UNARY, //10
117              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};              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  inline
124  ES_opgroup  ES_opgroup
125  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 79  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 (getOpgroup(op))      
187     {    // Get rank and shape of inputs
188     case G_BINARY: return left->getLength();    int rank0 = left->getRank();
189     case G_UNARY: return left->getLength();    int rank1 = right->getRank();
190     default:    const DataTypes::ShapeType& shape0 = left->getShape();
191      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");    const DataTypes::ShapeType& shape1 = right->getShape();
192     }  
193      // Prepare for the loops of the product and verify compatibility of shapes
194      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  {  {
# Line 118  calcBuffs(const DataLazy_ptr& left, cons Line 256  calcBuffs(const DataLazy_ptr& left, cons
256     case G_IDENTITY: return 1;     case G_IDENTITY: return 1;
257     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
258     case G_UNARY: return max(left->getBuffsRequired(),1);     case G_UNARY: return max(left->getBuffsRequired(),1);
259       case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
260       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    
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 141  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 154  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     }     }
    m_length=p->getLength();  
305     m_buffsRequired=1;     m_buffsRequired=1;
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    
312    
313    
314  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
315      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
316      m_op(op)      m_op(op),
317        m_axis_offset(0),
318        m_transpose(0),
319        m_SL(0), m_SM(0), m_SR(0)
320  {  {
321     if (getOpgroup(op)!=G_UNARY)     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
322     {     {
323      throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");      throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
324     }     }
325    
326     DataLazy_ptr lleft;     DataLazy_ptr lleft;
327     if (!left->isLazy())     if (!left->isLazy())
328     {     {
# Line 178  DataLazy::DataLazy(DataAbstract_ptr left Line 332  DataLazy::DataLazy(DataAbstract_ptr left
332     {     {
333      lleft=dynamic_pointer_cast<DataLazy>(left);      lleft=dynamic_pointer_cast<DataLazy>(left);
334     }     }
335     m_length=left->getLength();     m_readytype=lleft->m_readytype;
336     m_left=lleft;     m_left=lleft;
337     m_buffsRequired=1;     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_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
340  }  }
341    
342    
343  DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)  // In this constructor we need to consider interpolation
344    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_left(left),      m_op(op),
347      m_right(right),      m_SL(0), m_SM(0), m_SR(0)
     m_op(op)  
348  {  {
349     if (getOpgroup(op)!=G_BINARY)     if ((getOpgroup(op)!=G_BINARY))
350     {     {
351      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");      throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
352     }     }
353     m_length=resultLength(m_left,m_right,m_op);  
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);
371       }
372       else
373       {
374        m_left=DataLazy_ptr(new DataLazy(left));
375       }
376       if (right->isLazy())
377       {
378        m_right=dynamic_pointer_cast<DataLazy>(right);
379       }
380       else
381       {
382        m_right=DataLazy_ptr(new DataLazy(right));
383       }
384       char lt=m_left->m_readytype;
385       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_buffsRequired=calcBuffs(m_left, m_right, m_op);     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
400  cout << "(2)Lazy created with " << m_samplesize << endl;     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
401    cout << "(3)Lazy created with " << m_samplesize << endl;
402  }  }
403    
404  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
405      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
406      m_op(op)      m_op(op),
407        m_axis_offset(axis_offset),
408        m_transpose(transpose)
409  {  {
410     if (getOpgroup(op)!=G_BINARY)     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      throw DataException("Programmer error - constructor DataLazy(left, op) will only process BINARY operations.");      FunctionSpace fs=getFunctionSpace();
421        Data ltemp(left);
422        Data tmp(ltemp,fs);
423        left=tmp.borrowDataPtr();
424     }     }
425     if (left->isLazy())     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);      m_left=dynamic_pointer_cast<DataLazy>(left);
435     }     }
# Line 225  DataLazy::DataLazy(DataAbstract_ptr left Line 445  DataLazy::DataLazy(DataAbstract_ptr left
445     {     {
446      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
447     }     }
448       char lt=m_left->m_readytype;
449     m_length=resultLength(m_left,m_right,m_op);     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();     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);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
465  cout << "(3)Lazy created with " << m_samplesize << endl;  cout << "(4)Lazy created with " << m_samplesize << endl;
466  }  }
467    
468    
# Line 245  DataLazy::getBuffsRequired() const Line 478  DataLazy::getBuffsRequired() const
478  }  }
479    
480    
481  // the vector and the offset are a place where the method could write its data if it wishes  size_t
482  // it is not obligated to do so. For example, if it has its own storage already, it can use that.  DataLazy::getMaxSampleSize() const
 // Hence the return value to indicate where the data is actually stored.  
 // Regardless, the storage should be assumed to be used, even if it isn't.  
 const double*  
 DataLazy::resolveSample(ValueType& v,int sampleNo,  size_t offset ) const  
483  {  {
484    if (m_op==IDENTITY)        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_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      const ValueType& vec=m_id->getVector();      return m_id;
     return &(vec[m_id->getPointOffset(sampleNo, 0)]);  
502    }    }
503    size_t rightoffset=offset+m_samplesize;    DataReady_ptr pleft=m_left->collapseToReady();
504    const double* left=m_left->resolveSample(v,sampleNo,offset);    Data left(pleft);
505    const double* right=0;    Data right;
506    if (getOpgroup(m_op)==G_BINARY)    if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
507    {    {
508      right=m_right->resolveSample(v,sampleNo,rightoffset);      right=Data(m_right->collapseToReady());
509    }    }
510    double* result=&(v[offset]);    Data result;
511      switch(m_op)
512    {    {
513      switch(m_op)      case ADD:
514      {      result=left+right;
     case ADD:       // since these are pointwise ops, pretend each sample is one point  
     tensor_binary_operation(m_samplesize, left, right, result, plus<double>());  
515      break;      break;
516      case SUB:            case SUB:      
517      tensor_binary_operation(m_samplesize, left, right, result, minus<double>());      result=left-right;
518      break;      break;
519      case MUL:            case MUL:      
520      tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());      result=left*right;
521      break;      break;
522      case DIV:            case DIV:      
523      tensor_binary_operation(m_samplesize, left, right, result, divides<double>());      result=left/right;
524      break;      break;
 // unary ops  
525      case SIN:      case SIN:
526      tensor_unary_operation(m_samplesize, left, result, ::sin);      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;      break;
672      case COS:      case COS:
673      tensor_unary_operation(m_samplesize, left, result, ::cos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
674      break;      break;
675      case TAN:      case TAN:
676      tensor_unary_operation(m_samplesize, left, result, ::tan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
677      break;      break;
678      case ASIN:      case ASIN:
679      tensor_unary_operation(m_samplesize, left, result, ::asin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
680      break;      break;
681      case ACOS:      case ACOS:
682      tensor_unary_operation(m_samplesize, left, result, ::acos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
683      break;      break;
684      case ATAN:      case ATAN:
685      tensor_unary_operation(m_samplesize, left, result, ::atan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
686      break;      break;
687      case SINH:      case SINH:
688      tensor_unary_operation(m_samplesize, left, result, ::sinh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
689      break;      break;
690      case COSH:      case COSH:
691      tensor_unary_operation(m_samplesize, left, result, ::cosh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
692      break;      break;
693      case TANH:      case TANH:
694      tensor_unary_operation(m_samplesize, left, result, ::tanh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
695      break;      break;
696      case ERF:      case ERF:
697  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
698      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
699  #else  #else
700      tensor_unary_operation(m_samplesize, left, result, ::erf);      tensor_unary_operation(m_samplesize, left, result, ::erf);
701      break;      break;
702  #endif  #endif
703     case ASINH:     case ASINH:
704  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
705      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
706  #else  #else
707      tensor_unary_operation(m_samplesize, left, result, ::asinh);      tensor_unary_operation(m_samplesize, left, result, ::asinh);
708  #endif    #endif  
709      break;      break;
710     case ACOSH:     case ACOSH:
711  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
712      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
713  #else  #else
714      tensor_unary_operation(m_samplesize, left, result, ::acosh);      tensor_unary_operation(m_samplesize, left, result, ::acosh);
715  #endif    #endif  
716      break;      break;
717     case ATANH:     case ATANH:
718  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
719      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
720  #else  #else
721      tensor_unary_operation(m_samplesize, left, result, ::atanh);      tensor_unary_operation(m_samplesize, left, result, ::atanh);
722  #endif    #endif  
723      break;      break;
724      case LOG10:      case LOG10:
725      tensor_unary_operation(m_samplesize, left, result, ::log10);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
726      break;      break;
727      case LOG:      case LOG:
728      tensor_unary_operation(m_samplesize, left, result, ::log);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
729      break;      break;
730      case SIGN:      case SIGN:
731      tensor_unary_operation(m_samplesize, left, result, escript::fsign);      tensor_unary_operation(m_samplesize, left, result, escript::fsign);
732      break;      break;
733      case ABS:      case ABS:
734      tensor_unary_operation(m_samplesize, left, result, ::fabs);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
735      break;      break;
736      case NEG:      case NEG:
737      tensor_unary_operation(m_samplesize, left, result, negate<double>());      tensor_unary_operation(m_samplesize, left, result, negate<double>());
# Line 357  DataLazy::resolveSample(ValueType& v,int Line 742  DataLazy::resolveSample(ValueType& v,int
742      throw DataException("Programmer error - POS not supported for lazy data.");      throw DataException("Programmer error - POS not supported for lazy data.");
743      break;      break;
744      case EXP:      case EXP:
745      tensor_unary_operation(m_samplesize, left, result, ::exp);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
746      break;      break;
747      case SQRT:      case SQRT:
748      tensor_unary_operation(m_samplesize, left, result, ::sqrt);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
749      break;      break;
750      case RECIP:      case RECIP:
751      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
# Line 379  DataLazy::resolveSample(ValueType& v,int Line 764  DataLazy::resolveSample(ValueType& v,int
764      break;      break;
765    
766      default:      default:
767      throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");      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;
940        default:
941        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    }    }
   return result;  
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    size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);    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;    int numthreads=1;
1022  #ifdef _OPENMP  #ifdef _OPENMP
1023    numthreads=omp_get_max_threads();    numthreads=getNumberOfThreads();
1024    int threadnum=0;    int threadnum=0;
1025  #endif  #endif
1026    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
1027  cout << "Buffer created with size=" << v.size() << endl;  cout << "Buffer created with size=" << v.size() << endl;
1028    ValueType dummy(getNoValues());    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
   DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),dummy);  
1029    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVector();
1030    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1031    int sample;    int sample;
1032    int resoffset;    size_t outoffset;     // offset in the output data
1033    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1034    #pragma omp parallel for private(sample,resoffset,threadnum) schedule(static)    const ValueType* res=0;   // Vector storing the answer
1035      size_t resoffset=0;       // where in the vector to find the answer
1036      #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)
1037    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
1038    {    {
1039    cout << "################################# " << sample << endl;
1040  #ifdef _OPENMP  #ifdef _OPENMP
1041      const double* res=resolveSample(v,sample,threadbuffersize*omp_get_thread_num());      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
1042  #else  #else
1043      const double* res=resolveSample(v,sample,0);   // this would normally be v, but not if its a single IDENTITY op.      res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.
1044  #endif  #endif
1045      resoffset=result->getPointOffset(sample,0);  cerr << "-------------------------------- " << endl;
1046      for (int i=0;i<m_samplesize;++i,++resoffset)    // copy values into the output vector      outoffset=result->getPointOffset(sample,0);
1047    cerr << "offset=" << outoffset << endl;
1048        for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector
1049      {      {
1050      resvec[resoffset]=res[i];      resvec[outoffset]=(*res)[resoffset];
1051      }      }
1052    cerr << "*********************************" << endl;
1053    }    }
1054    return resptr;    return resptr;
1055  }  }
# Line 435  DataLazy::toString() const Line 1063  DataLazy::toString() const
1063    return oss.str();    return oss.str();
1064  }  }
1065    
1066    
1067  void  void
1068  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
1069  {  {
1070    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1071    {    {
1072    case G_IDENTITY:    case G_IDENTITY:
1073        if (m_id->isExpanded())
1074        {
1075           oss << "E";
1076        }
1077        else if (m_id->isTagged())
1078        {
1079          oss << "T";
1080        }
1081        else if (m_id->isConstant())
1082        {
1083          oss << "C";
1084        }
1085        else
1086        {
1087          oss << "?";
1088        }
1089      oss << '@' << m_id.get();      oss << '@' << m_id.get();
1090      break;      break;
1091    case G_BINARY:    case G_BINARY:
# Line 451  DataLazy::intoString(ostringstream& oss) Line 1096  DataLazy::intoString(ostringstream& oss)
1096      oss << ')';      oss << ')';
1097      break;      break;
1098    case G_UNARY:    case G_UNARY:
1099      case G_NP1OUT:
1100      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
1101      m_left->intoString(oss);      m_left->intoString(oss);
1102      oss << ')';      oss << ')';
1103      break;      break;
1104      case G_TENSORPROD:
1105        oss << opToString(m_op) << '(';
1106        m_left->intoString(oss);
1107        oss << ", ";
1108        m_right->intoString(oss);
1109        oss << ')';
1110        break;
1111    default:    default:
1112      oss << "UNKNOWN";      oss << "UNKNOWN";
1113    }    }
1114  }  }
1115    
 // 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.  
1116  DataAbstract*  DataAbstract*
1117  DataLazy::deepCopy()  DataLazy::deepCopy()
1118  {  {
1119    if (m_op==IDENTITY)    switch (getOpgroup(m_op))
1120    {    {
1121      return new DataLazy(m_left);    // we don't need to copy the child here    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1122      case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1123      case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1124      case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1125      case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1126      default:
1127        throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1128    }    }
   return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);  
1129  }  }
1130    
1131    
1132    // There is no single, natural interpretation of getLength on DataLazy.
1133    // Instances of DataReady can look at the size of their vectors.
1134    // For lazy though, it could be the size the data would be if it were resolved;
1135    // or it could be some function of the lengths of the DataReady instances which
1136    // form part of the expression.
1137    // Rather than have people making assumptions, I have disabled the method.
1138  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1139  DataLazy::getLength() const  DataLazy::getLength() const
1140  {  {
1141    return m_length;    throw DataException("getLength() does not make sense for lazy data.");
1142  }  }
1143    
1144    
# Line 486  DataLazy::getSlice(const DataTypes::Regi Line 1148  DataLazy::getSlice(const DataTypes::Regi
1148    throw DataException("getSlice - not implemented for Lazy objects.");    throw DataException("getSlice - not implemented for Lazy objects.");
1149  }  }
1150    
1151    
1152    // To do this we need to rely on our child nodes
1153    DataTypes::ValueType::size_type
1154    DataLazy::getPointOffset(int sampleNo,
1155                     int dataPointNo)
1156    {
1157      if (m_op==IDENTITY)
1158      {
1159        return m_id->getPointOffset(sampleNo,dataPointNo);
1160      }
1161      if (m_readytype!='E')
1162      {
1163        collapse();
1164        return m_id->getPointOffset(sampleNo,dataPointNo);
1165      }
1166      // at this point we do not have an identity node and the expression will be Expanded
1167      // so we only need to know which child to ask
1168      if (m_left->m_readytype=='E')
1169      {
1170        return m_left->getPointOffset(sampleNo,dataPointNo);
1171      }
1172      else
1173      {
1174        return m_right->getPointOffset(sampleNo,dataPointNo);
1175      }
1176    }
1177    
1178    // To do this we need to rely on our child nodes
1179  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1180  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
1181                   int dataPointNo) const                   int dataPointNo) const
1182  {  {
1183    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");    if (m_op==IDENTITY)
1184      {
1185        return m_id->getPointOffset(sampleNo,dataPointNo);
1186      }
1187      if (m_readytype=='E')
1188      {
1189        // at this point we do not have an identity node and the expression will be Expanded
1190        // so we only need to know which child to ask
1191        if (m_left->m_readytype=='E')
1192        {
1193        return m_left->getPointOffset(sampleNo,dataPointNo);
1194        }
1195        else
1196        {
1197        return m_right->getPointOffset(sampleNo,dataPointNo);
1198        }
1199      }
1200      if (m_readytype=='C')
1201      {
1202        return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
1203      }
1204      throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1205    }
1206    
1207    // It would seem that DataTagged will need to be treated differently since even after setting all tags
1208    // to zero, all the tags from all the DataTags would be in the result.
1209    // However since they all have the same value (0) whether they are there or not should not matter.
1210    // So I have decided that for all types this method will create a constant 0.
1211    // It can be promoted up as required.
1212    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
1213    // but we can deal with that if it arrises.
1214    void
1215    DataLazy::setToZero()
1216    {
1217      DataTypes::ValueType v(getNoValues(),0);
1218      m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1219      m_op=IDENTITY;
1220      m_right.reset();  
1221      m_left.reset();
1222      m_readytype='C';
1223      m_buffsRequired=1;
1224  }  }
1225    
1226  }   // end namespace  }   // end namespace

Legend:
Removed from v.1888  
changed lines
  Added in v.2066

  ViewVC Help
Powered by ViewVC 1.1.26