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

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

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

branches/schroedinger/escript/src/DataLazy.cpp revision 1898 by jfenwick, Mon Oct 20 01:20:18 2008 UTC trunk/escript/src/DataLazy.cpp revision 2153 by jfenwick, Fri Dec 12 00:18:18 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    // #define LAZYDEBUG(X) X;
32    #define LAZYDEBUG(X)
33    
34    /*
35    How does DataLazy work?
36    ~~~~~~~~~~~~~~~~~~~~~~~
37    
38    Each instance represents a single operation on one or two other DataLazy instances. These arguments are normally
39    denoted left and right.
40    
41    A special operation, IDENTITY, stores an instance of DataReady in the m_id member.
42    This means that all "internal" nodes in the structure are instances of DataLazy.
43    
44    Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
45    Note that IDENITY is not considered a unary operation.
46    
47    I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
48    It must however form a DAG (directed acyclic graph).
49    I will refer to individual DataLazy objects with the structure as nodes.
50    
51    Each node also stores:
52    - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
53        evaluated.
54    - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to
55        evaluate the expression.
56    - m_samplesize ~ the number of doubles stored in a sample.
57    
58    When a new node is created, the above values are computed based on the values in the child nodes.
59    Eg: if left requires 4 samples and right requires 6 then left+right requires 7 samples.
60    
61    The resolve method, which produces a DataReady from a DataLazy, does the following:
62    1) Create a DataReady to hold the new result.
63    2) Allocate a vector (v) big enough to hold m_buffsrequired samples.
64    3) For each sample, call resolveSample with v, to get its values and copy them into the result object.
65    
66    (In the case of OMP, multiple samples are resolved in parallel so the vector needs to be larger.)
67    
68    resolveSample returns a Vector* and an offset within that vector where the result is stored.
69    Normally, this would be v, but for identity nodes their internal vector is returned instead.
70    
71    The convention that I use, is that the resolve methods should store their results starting at the offset they are passed.
72    
73    For expressions which evaluate to Constant or Tagged, there is a different evaluation method.
74    The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.
75    
76    To add a new operator you need to do the following (plus anything I might have forgotten - adding a new group for example):
77    1) Add to the ES_optype.
78    2) determine what opgroup your operation belongs to (X)
79    3) add a string for the op to the end of ES_opstrings
80    4) increase ES_opcount
81    5) add an entry (X) to opgroups
82    6) add an entry to the switch in collapseToReady
83    7) add an entry to resolveX
84    */
85    
86    
87  using namespace std;  using namespace std;
88  using namespace boost;  using namespace boost;
# Line 33  using namespace boost; Line 90  using namespace boost;
90  namespace escript  namespace escript
91  {  {
92    
 const std::string&  
 opToString(ES_optype op);  
   
93  namespace  namespace
94  {  {
95    
   
   
96  enum ES_opgroup  enum ES_opgroup
97  {  {
98     G_UNKNOWN,     G_UNKNOWN,
99     G_IDENTITY,     G_IDENTITY,
100     G_BINARY,     G_BINARY,        // pointwise operations with two arguments
101     G_UNARY     G_UNARY,     // pointwise operations with one argument
102       G_UNARY_P,       // pointwise operations with one argument, requiring a parameter
103       G_NP1OUT,        // non-pointwise op with one output
104       G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
105       G_TENSORPROD     // general tensor product
106  };  };
107    
108    
109    
110    
111  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","sin","cos","tan",  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
112                "sin","cos","tan",
113              "asin","acos","atan","sinh","cosh","tanh","erf",              "asin","acos","atan","sinh","cosh","tanh","erf",
114              "asinh","acosh","atanh",              "asinh","acosh","atanh",
115              "log10","log","sign","abs","neg","pos","exp","sqrt",              "log10","log","sign","abs","neg","pos","exp","sqrt",
116              "1/","where>0","where<0","where>=0","where<=0"};              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
117  int ES_opcount=32;              "symmetric","nonsymmetric",
118  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY,G_UNARY,G_UNARY,G_UNARY, //9              "prod",
119              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 16              "transpose", "trace"};
120              G_UNARY,G_UNARY,G_UNARY,                    // 19  int ES_opcount=40;
121              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,        // 27  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
122              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};              G_UNARY,G_UNARY,G_UNARY, //10
123                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17
124                G_UNARY,G_UNARY,G_UNARY,                    // 20
125                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 28
126                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P,      // 35
127                G_NP1OUT,G_NP1OUT,
128                G_TENSORPROD,
129                G_NP1OUT_P, G_NP1OUT_P};
130  inline  inline
131  ES_opgroup  ES_opgroup
132  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 79  resultFS(DataAbstract_ptr left, DataAbst Line 143  resultFS(DataAbstract_ptr left, DataAbst
143      // 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
144      // programming error exception.      // programming error exception.
145    
146      FunctionSpace l=left->getFunctionSpace();
147      if (left->getFunctionSpace()!=right->getFunctionSpace())    FunctionSpace r=right->getFunctionSpace();
148      {    if (l!=r)
149          throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");    {
150      }      if (r.probeInterpolation(l))
151      return left->getFunctionSpace();      {
152        return l;
153        }
154        if (l.probeInterpolation(r))
155        {
156        return r;
157        }
158        throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
159      }
160      return l;
161  }  }
162    
163  // return the shape of the result of "left op right"  // return the shape of the result of "left op right"
164    // the shapes resulting from tensor product are more complex to compute so are worked out elsewhere
165  DataTypes::ShapeType  DataTypes::ShapeType
166  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
167  {  {
168      if (left->getShape()!=right->getShape())      if (left->getShape()!=right->getShape())
169      {      {
170          throw DataException("Shapes not the same - shapes must match for lazy data.");        if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
171          {
172            throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
173          }
174          if (left->getRank()==0)   // we need to allow scalar * anything
175          {
176            return right->getShape();
177          }
178          if (right->getRank()==0)
179          {
180            return left->getShape();
181          }
182          throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
183      }      }
184      return left->getShape();      return left->getShape();
185  }  }
186    
187  size_t  // return the shape for "op left"
188  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  
189    DataTypes::ShapeType
190    resultShape(DataAbstract_ptr left, ES_optype op)
191  {  {
192     switch (getOpgroup(op))      switch(op)
193     {      {
194     case G_BINARY: return left->getLength();          case TRANS:
195     case G_UNARY: return left->getLength();          return left->getShape();
196     default:      break;
197      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");      case TRACE:
198     }          return DataTypes::scalarShape;
199        break;
200            default:
201        throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
202        }
203    }
204    
205    // determine the output shape for the general tensor product operation
206    // the additional parameters return information required later for the product
207    // the majority of this code is copy pasted from C_General_Tensor_Product
208    DataTypes::ShapeType
209    GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
210    {
211        
212      // Get rank and shape of inputs
213      int rank0 = left->getRank();
214      int rank1 = right->getRank();
215      const DataTypes::ShapeType& shape0 = left->getShape();
216      const DataTypes::ShapeType& shape1 = right->getShape();
217    
218      // Prepare for the loops of the product and verify compatibility of shapes
219      int start0=0, start1=0;
220      if (transpose == 0)       {}
221      else if (transpose == 1)  { start0 = axis_offset; }
222      else if (transpose == 2)  { start1 = rank1-axis_offset; }
223      else              { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
224    
225      if (rank0<axis_offset)
226      {
227        throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
228      }
229    
230      // Adjust the shapes for transpose
231      DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather
232      DataTypes::ShapeType tmpShape1(rank1);    // than using push_back
233      for (int i=0; i<rank0; i++)   { tmpShape0[i]=shape0[(i+start0)%rank0]; }
234      for (int i=0; i<rank1; i++)   { tmpShape1[i]=shape1[(i+start1)%rank1]; }
235    
236      // Prepare for the loops of the product
237      SL=1, SM=1, SR=1;
238      for (int i=0; i<rank0-axis_offset; i++)   {
239        SL *= tmpShape0[i];
240      }
241      for (int i=rank0-axis_offset; i<rank0; i++)   {
242        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
243          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
244        }
245        SM *= tmpShape0[i];
246      }
247      for (int i=axis_offset; i<rank1; i++)     {
248        SR *= tmpShape1[i];
249      }
250    
251      // Define the shape of the output (rank of shape is the sum of the loop ranges below)
252      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);  
253      {         // block to limit the scope of out_index
254         int out_index=0;
255         for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
256         for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
257      }
258    
259      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
260      {
261         ostringstream os;
262         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
263         throw DataException(os.str());
264      }
265    
266      return shape2;
267  }  }
268    
269    
270    // determine the number of points in the result of "left op right"
271    // note that determining the resultLength for G_TENSORPROD is more complex and will not be processed here
272    // size_t
273    // resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
274    // {
275    //    switch (getOpgroup(op))
276    //    {
277    //    case G_BINARY: return left->getLength();
278    //    case G_UNARY: return left->getLength();
279    //    case G_NP1OUT: return left->getLength();
280    //    default:
281    //  throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");
282    //    }
283    // }
284    
285    // determine the number of samples requires to evaluate an expression combining left and right
286    // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
287    // The same goes for G_TENSORPROD
288  int  int
289  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
290  {  {
# Line 117  calcBuffs(const DataLazy_ptr& left, cons Line 292  calcBuffs(const DataLazy_ptr& left, cons
292     {     {
293     case G_IDENTITY: return 1;     case G_IDENTITY: return 1;
294     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);
295     case G_UNARY: return max(left->getBuffsRequired(),1);     case G_UNARY:
296       case G_UNARY_P: return max(left->getBuffsRequired(),1);
297       case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
298       case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);
299       case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
300     default:     default:
301      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
302     }     }
303  }  }
304    
305    
306  }   // end anonymous namespace  }   // end anonymous namespace
307    
308    
309    
310    // Return a string representing the operation
311  const std::string&  const std::string&
312  opToString(ES_optype op)  opToString(ES_optype op)
313  {  {
# Line 139  opToString(ES_optype op) Line 321  opToString(ES_optype op)
321    
322  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
323      : parent(p->getFunctionSpace(),p->getShape()),      : parent(p->getFunctionSpace(),p->getShape()),
324      m_op(IDENTITY)      m_op(IDENTITY),
325        m_axis_offset(0),
326        m_transpose(0),
327        m_SL(0), m_SM(0), m_SR(0)
328  {  {
329     if (p->isLazy())     if (p->isLazy())
330     {     {
     // TODO: fix this.   We could make the new node a copy of p?  
331      // I don't want identity of Lazy.      // I don't want identity of Lazy.
332      // Question: Why would that be so bad?      // Question: Why would that be so bad?
333      // 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 157  DataLazy::DataLazy(DataAbstract_ptr p) Line 341  DataLazy::DataLazy(DataAbstract_ptr p)
341      else if (p->isTagged()) {m_readytype='T';}      else if (p->isTagged()) {m_readytype='T';}
342      else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}      else {throw DataException("Unknown DataReady instance in DataLazy constructor.");}
343     }     }
    m_length=p->getLength();  
344     m_buffsRequired=1;     m_buffsRequired=1;
345     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
346  cout << "(1)Lazy created with " << m_samplesize << endl;     m_maxsamplesize=m_samplesize;
347    LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
348  }  }
349    
350    
351    
352    
353  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
354      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
355      m_op(op)      m_op(op),
356        m_axis_offset(0),
357        m_transpose(0),
358        m_SL(0), m_SM(0), m_SR(0)
359  {  {
360     if (getOpgroup(op)!=G_UNARY)     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
361     {     {
362      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.");
363     }     }
364    
365     DataLazy_ptr lleft;     DataLazy_ptr lleft;
366     if (!left->isLazy())     if (!left->isLazy())
367     {     {
# Line 181  DataLazy::DataLazy(DataAbstract_ptr left Line 372  DataLazy::DataLazy(DataAbstract_ptr left
372      lleft=dynamic_pointer_cast<DataLazy>(left);      lleft=dynamic_pointer_cast<DataLazy>(left);
373     }     }
374     m_readytype=lleft->m_readytype;     m_readytype=lleft->m_readytype;
    m_length=left->getLength();  
375     m_left=lleft;     m_left=lleft;
376     m_buffsRequired=1;     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
377     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
378       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
379  }  }
380    
381    
382  // DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)  // In this constructor we need to consider interpolation
 //  : parent(resultFS(left,right,op), resultShape(left,right,op)),  
 //  m_left(left),  
 //  m_right(right),  
 //  m_op(op)  
 // {  
 //    if (getOpgroup(op)!=G_BINARY)  
 //    {  
 //  throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");  
 //    }  
 //    m_length=resultLength(m_left,m_right,m_op);  
 //    m_samplesize=getNumDPPSample()*getNoValues();  
 //    m_buffsRequired=calcBuffs(m_left, m_right, m_op);  
 // cout << "(2)Lazy created with " << m_samplesize << endl;  
 // }  
   
383  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
384      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
385      m_op(op)      m_op(op),
386        m_SL(0), m_SM(0), m_SR(0)
387  {  {
388     if (getOpgroup(op)!=G_BINARY)     if ((getOpgroup(op)!=G_BINARY))
389     {     {
390      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.");
391     }     }
392     if (left->isLazy())  
393       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
394       {
395        FunctionSpace fs=getFunctionSpace();
396        Data ltemp(left);
397        Data tmp(ltemp,fs);
398        left=tmp.borrowDataPtr();
399       }
400       if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
401       {
402        Data tmp(Data(right),getFunctionSpace());
403        right=tmp.borrowDataPtr();
404       }
405       left->operandCheck(*right);
406    
407       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
408       {
409        m_left=dynamic_pointer_cast<DataLazy>(left);
410       }
411       else
412       {
413        m_left=DataLazy_ptr(new DataLazy(left));
414       }
415       if (right->isLazy())
416       {
417        m_right=dynamic_pointer_cast<DataLazy>(right);
418       }
419       else
420       {
421        m_right=DataLazy_ptr(new DataLazy(right));
422       }
423       char lt=m_left->m_readytype;
424       char rt=m_right->m_readytype;
425       if (lt=='E' || rt=='E')
426       {
427        m_readytype='E';
428       }
429       else if (lt=='T' || rt=='T')
430       {
431        m_readytype='T';
432       }
433       else
434       {
435        m_readytype='C';
436       }
437       m_samplesize=getNumDPPSample()*getNoValues();
438       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
439       m_buffsRequired=calcBuffs(m_left, m_right,m_op);
440    LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
441    }
442    
443    DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
444        : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
445        m_op(op),
446        m_axis_offset(axis_offset),
447        m_transpose(transpose)
448    {
449       if ((getOpgroup(op)!=G_TENSORPROD))
450       {
451        throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
452       }
453       if ((transpose>2) || (transpose<0))
454       {
455        throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
456       }
457       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
458       {
459        FunctionSpace fs=getFunctionSpace();
460        Data ltemp(left);
461        Data tmp(ltemp,fs);
462        left=tmp.borrowDataPtr();
463       }
464       if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
465       {
466        Data tmp(Data(right),getFunctionSpace());
467        right=tmp.borrowDataPtr();
468       }
469       left->operandCheck(*right);
470    
471       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
472     {     {
473      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
474     }     }
# Line 242  DataLazy::DataLazy(DataAbstract_ptr left Line 498  DataLazy::DataLazy(DataAbstract_ptr left
498     {     {
499      m_readytype='C';      m_readytype='C';
500     }     }
    m_length=resultLength(m_left,m_right,m_op);  
501     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
502       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
503     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
504  cout << "(3)Lazy created with " << m_samplesize << endl;  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
505  }  }
506    
507    
508    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
509        : parent(left->getFunctionSpace(), resultShape(left,op)),
510        m_op(op),
511        m_axis_offset(axis_offset),
512        m_transpose(0),
513        m_tol(0)
514    {
515       if ((getOpgroup(op)!=G_NP1OUT_P))
516       {
517        throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
518       }
519       DataLazy_ptr lleft;
520       if (!left->isLazy())
521       {
522        lleft=DataLazy_ptr(new DataLazy(left));
523       }
524       else
525       {
526        lleft=dynamic_pointer_cast<DataLazy>(left);
527       }
528       m_readytype=lleft->m_readytype;
529       m_left=lleft;
530       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
531       m_samplesize=getNumDPPSample()*getNoValues();
532       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
533    LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
534    }
535    
536    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
537        : parent(left->getFunctionSpace(), left->getShape()),
538        m_op(op),
539        m_axis_offset(0),
540        m_transpose(0),
541        m_tol(tol)
542    {
543       if ((getOpgroup(op)!=G_UNARY_P))
544       {
545        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
546       }
547       DataLazy_ptr lleft;
548       if (!left->isLazy())
549       {
550        lleft=DataLazy_ptr(new DataLazy(left));
551       }
552       else
553       {
554        lleft=dynamic_pointer_cast<DataLazy>(left);
555       }
556       m_readytype=lleft->m_readytype;
557       m_left=lleft;
558       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
559       m_samplesize=getNumDPPSample()*getNoValues();
560       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
561    LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
562    }
563    
564  DataLazy::~DataLazy()  DataLazy::~DataLazy()
565  {  {
566  }  }
# Line 261  DataLazy::getBuffsRequired() const Line 573  DataLazy::getBuffsRequired() const
573  }  }
574    
575    
576    size_t
577    DataLazy::getMaxSampleSize() const
578    {
579        return m_maxsamplesize;
580    }
581    
582    /*
583      \brief Evaluates the expression using methods on Data.
584      This does the work for the collapse method.
585      For reasons of efficiency do not call this method on DataExpanded nodes.
586    */
587  DataReady_ptr  DataReady_ptr
588  DataLazy::collapseToReady()  DataLazy::collapseToReady()
589  {  {
# Line 275  DataLazy::collapseToReady() Line 598  DataLazy::collapseToReady()
598    DataReady_ptr pleft=m_left->collapseToReady();    DataReady_ptr pleft=m_left->collapseToReady();
599    Data left(pleft);    Data left(pleft);
600    Data right;    Data right;
601    if (getOpgroup(m_op)==G_BINARY)    if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
602    {    {
603      right=Data(m_right->collapseToReady());      right=Data(m_right->collapseToReady());
604    }    }
# Line 374  DataLazy::collapseToReady() Line 697  DataLazy::collapseToReady()
697      case LEZ:      case LEZ:
698      result=left.whereNonPositive();      result=left.whereNonPositive();
699      break;      break;
700        case NEZ:
701        result=left.whereNonZero(m_tol);
702        break;
703        case EZ:
704        result=left.whereZero(m_tol);
705        break;
706        case SYM:
707        result=left.symmetric();
708        break;
709        case NSYM:
710        result=left.nonsymmetric();
711        break;
712        case PROD:
713        result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
714        break;
715        case TRANS:
716        result=left.transpose(m_axis_offset);
717        break;
718        case TRACE:
719        result=left.trace(m_axis_offset);
720        break;
721      default:      default:
722      throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
723    }    }
724    return result.borrowReadyPtr();    return result.borrowReadyPtr();
725  }  }
726    
727    /*
728       \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
729       This method uses the original methods on the Data class to evaluate the expressions.
730       For this reason, it should not be used on DataExpanded instances. (To do so would defeat
731       the purpose of using DataLazy in the first place).
732    */
733  void  void
734  DataLazy::collapse()  DataLazy::collapse()
735  {  {
# Line 395  DataLazy::collapse() Line 745  DataLazy::collapse()
745    m_op=IDENTITY;    m_op=IDENTITY;
746  }  }
747    
748    /*
749      \brief Compute the value of the expression (unary operation) for the given sample.
750      \return Vector which stores the value of the subexpression for the given sample.
751      \param v A vector to store intermediate results.
752      \param offset Index in v to begin storing results.
753      \param sampleNo Sample number to evaluate.
754      \param roffset (output parameter) the offset in the return vector where the result begins.
755    
756      The return value will be an existing vector so do not deallocate it.
757      If the result is stored in v it should be stored at the offset given.
758      Everything from offset to the end of v should be considered available for this method to use.
759    */
760  DataTypes::ValueType*  DataTypes::ValueType*
761  DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const  DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
762  {  {
# Line 412  DataLazy::resolveUnary(ValueType& v, siz Line 774  DataLazy::resolveUnary(ValueType& v, siz
774    switch (m_op)    switch (m_op)
775    {    {
776      case SIN:        case SIN:  
777      tensor_unary_operation(m_samplesize, left, result, ::sin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
778      break;      break;
779      case COS:      case COS:
780      tensor_unary_operation(m_samplesize, left, result, ::cos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
781      break;      break;
782      case TAN:      case TAN:
783      tensor_unary_operation(m_samplesize, left, result, ::tan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
784      break;      break;
785      case ASIN:      case ASIN:
786      tensor_unary_operation(m_samplesize, left, result, ::asin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
787      break;      break;
788      case ACOS:      case ACOS:
789      tensor_unary_operation(m_samplesize, left, result, ::acos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
790      break;      break;
791      case ATAN:      case ATAN:
792      tensor_unary_operation(m_samplesize, left, result, ::atan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
793      break;      break;
794      case SINH:      case SINH:
795      tensor_unary_operation(m_samplesize, left, result, ::sinh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
796      break;      break;
797      case COSH:      case COSH:
798      tensor_unary_operation(m_samplesize, left, result, ::cosh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
799      break;      break;
800      case TANH:      case TANH:
801      tensor_unary_operation(m_samplesize, left, result, ::tanh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
802      break;      break;
803      case ERF:      case ERF:
804  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
805      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
806  #else  #else
807      tensor_unary_operation(m_samplesize, left, result, ::erf);      tensor_unary_operation(m_samplesize, left, result, ::erf);
808      break;      break;
809  #endif  #endif
810     case ASINH:     case ASINH:
811  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
812      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
813  #else  #else
814      tensor_unary_operation(m_samplesize, left, result, ::asinh);      tensor_unary_operation(m_samplesize, left, result, ::asinh);
815  #endif    #endif  
816      break;      break;
817     case ACOSH:     case ACOSH:
818  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
819      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
820  #else  #else
821      tensor_unary_operation(m_samplesize, left, result, ::acosh);      tensor_unary_operation(m_samplesize, left, result, ::acosh);
822  #endif    #endif  
823      break;      break;
824     case ATANH:     case ATANH:
825  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
826      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
827  #else  #else
828      tensor_unary_operation(m_samplesize, left, result, ::atanh);      tensor_unary_operation(m_samplesize, left, result, ::atanh);
829  #endif    #endif  
830      break;      break;
831      case LOG10:      case LOG10:
832      tensor_unary_operation(m_samplesize, left, result, ::log10);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
833      break;      break;
834      case LOG:      case LOG:
835      tensor_unary_operation(m_samplesize, left, result, ::log);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
836      break;      break;
837      case SIGN:      case SIGN:
838      tensor_unary_operation(m_samplesize, left, result, escript::fsign);      tensor_unary_operation(m_samplesize, left, result, escript::fsign);
839      break;      break;
840      case ABS:      case ABS:
841      tensor_unary_operation(m_samplesize, left, result, ::fabs);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
842      break;      break;
843      case NEG:      case NEG:
844      tensor_unary_operation(m_samplesize, left, result, negate<double>());      tensor_unary_operation(m_samplesize, left, result, negate<double>());
# Line 487  DataLazy::resolveUnary(ValueType& v, siz Line 849  DataLazy::resolveUnary(ValueType& v, siz
849      throw DataException("Programmer error - POS not supported for lazy data.");      throw DataException("Programmer error - POS not supported for lazy data.");
850      break;      break;
851      case EXP:      case EXP:
852      tensor_unary_operation(m_samplesize, left, result, ::exp);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
853      break;      break;
854      case SQRT:      case SQRT:
855      tensor_unary_operation(m_samplesize, left, result, ::sqrt);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
856      break;      break;
857      case RECIP:      case RECIP:
858      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
# Line 507  DataLazy::resolveUnary(ValueType& v, siz Line 869  DataLazy::resolveUnary(ValueType& v, siz
869      case LEZ:      case LEZ:
870      tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));      tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
871      break;      break;
872    // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
873        case NEZ:
874        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
875        break;
876        case EZ:
877        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
878        break;
879    
880      default:      default:
881      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
# Line 516  DataLazy::resolveUnary(ValueType& v, siz Line 885  DataLazy::resolveUnary(ValueType& v, siz
885    
886    
887    
 // const double*  
 // DataLazy::resolveUnary(ValueType& v,int sampleNo,  size_t offset) const  
 // {  
 //  // we assume that any collapsing has been done before we get here  
 //  // since we only have one argument we don't need to think about only  
 //  // processing single points.  
 //   if (m_readytype!='E')  
 //   {  
 //     throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");  
 //   }  
 //   const double* left=m_left->resolveSample(v,sampleNo,offset);  
 //   double* result=&(v[offset]);  
 //   switch (m_op)  
 //   {  
 //     case SIN:      
 //  tensor_unary_operation(m_samplesize, left, result, ::sin);  
 //  break;  
 //     case COS:  
 //  tensor_unary_operation(m_samplesize, left, result, ::cos);  
 //  break;  
 //     case TAN:  
 //  tensor_unary_operation(m_samplesize, left, result, ::tan);  
 //  break;  
 //     case ASIN:  
 //  tensor_unary_operation(m_samplesize, left, result, ::asin);  
 //  break;  
 //     case ACOS:  
 //  tensor_unary_operation(m_samplesize, left, result, ::acos);  
 //  break;  
 //     case ATAN:  
 //  tensor_unary_operation(m_samplesize, left, result, ::atan);  
 //  break;  
 //     case SINH:  
 //  tensor_unary_operation(m_samplesize, left, result, ::sinh);  
 //  break;  
 //     case COSH:  
 //  tensor_unary_operation(m_samplesize, left, result, ::cosh);  
 //  break;  
 //     case TANH:  
 //  tensor_unary_operation(m_samplesize, left, result, ::tanh);  
 //  break;  
 //     case ERF:  
 // #ifdef _WIN32  
 //  throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");  
 // #else  
 //  tensor_unary_operation(m_samplesize, left, result, ::erf);  
 //  break;  
 // #endif  
 //    case ASINH:  
 // #ifdef _WIN32  
 //  tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);  
 // #else  
 //  tensor_unary_operation(m_samplesize, left, result, ::asinh);  
 // #endif    
 //  break;  
 //    case ACOSH:  
 // #ifdef _WIN32  
 //  tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);  
 // #else  
 //  tensor_unary_operation(m_samplesize, left, result, ::acosh);  
 // #endif    
 //  break;  
 //    case ATANH:  
 // #ifdef _WIN32  
 //  tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);  
 // #else  
 //  tensor_unary_operation(m_samplesize, left, result, ::atanh);  
 // #endif    
 //  break;  
 //     case LOG10:  
 //  tensor_unary_operation(m_samplesize, left, result, ::log10);  
 //  break;  
 //     case LOG:  
 //  tensor_unary_operation(m_samplesize, left, result, ::log);  
 //  break;  
 //     case SIGN:  
 //  tensor_unary_operation(m_samplesize, left, result, escript::fsign);  
 //  break;  
 //     case ABS:  
 //  tensor_unary_operation(m_samplesize, left, result, ::fabs);  
 //  break;  
 //     case NEG:  
 //  tensor_unary_operation(m_samplesize, left, result, negate<double>());  
 //  break;  
 //     case POS:  
 //  // it doesn't mean anything for delayed.  
 //  // it will just trigger a deep copy of the lazy object  
 //  throw DataException("Programmer error - POS not supported for lazy data.");  
 //  break;  
 //     case EXP:  
 //  tensor_unary_operation(m_samplesize, left, result, ::exp);  
 //  break;  
 //     case SQRT:  
 //  tensor_unary_operation(m_samplesize, left, result, ::sqrt);  
 //  break;  
 //     case RECIP:  
 //  tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));  
 //  break;  
 //     case GZ:  
 //  tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));  
 //  break;  
 //     case LZ:  
 //  tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));  
 //  break;  
 //     case GEZ:  
 //  tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));  
 //  break;  
 //     case LEZ:  
 //  tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));  
 //  break;  
 //  
 //     default:  
 //  throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");  
 //   }  
 //   return result;  
 // }  
888    
 #define PROC_OP(X) \  
     for (int i=0;i<steps;++i,resultp+=getNoValues()) \  
     { \  
 cout << "Step#" << i << " chunk=" << chunksize << endl; \  
 cout << left[0] << left[1] << left[2] << endl; \  
 cout << right[0] << right[1] << right[2] << endl; \  
        tensor_binary_operation(chunksize, left, right, resultp, X); \  
        left+=leftStep; \  
        right+=rightStep; \  
 cout << "Result=" << result << " " << result[0] << result[1] << result[2] << endl; \  
     }  
889    
890    
891    /*
892      \brief Compute the value of the expression (unary operation) for the given sample.
893      \return Vector which stores the value of the subexpression for the given sample.
894      \param v A vector to store intermediate results.
895      \param offset Index in v to begin storing results.
896      \param sampleNo Sample number to evaluate.
897      \param roffset (output parameter) the offset in the return vector where the result begins.
898    
899      The return value will be an existing vector so do not deallocate it.
900      If the result is stored in v it should be stored at the offset given.
901      Everything from offset to the end of v should be considered available for this method to use.
902    */
903  DataTypes::ValueType*  DataTypes::ValueType*
904  DataLazy::resolveBinary(ValueType& v,  size_t offset ,int sampleNo, size_t& roffset) const  DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
905  {  {
906      // again we assume that all collapsing has already been done      // we assume that any collapsing has been done before we get here
907      // so we have at least one expanded child.      // since we only have one argument we don't need to think about only
908      // however, we could still have one of the children being not expanded.      // processing single points.
909      if (m_readytype!='E')
910      {
911        throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
912      }
913        // since we can't write the result over the input, we need a result offset further along
914      size_t subroffset=roffset+m_samplesize;
915      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);
916      roffset=offset;
917      switch (m_op)
918      {
919        case SYM:
920        DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
921        break;
922        case NSYM:
923        DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
924        break;
925        default:
926        throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
927      }
928      return &v;
929    }
930    
931    /*
932      \brief Compute the value of the expression (unary operation) for the given sample.
933      \return Vector which stores the value of the subexpression for the given sample.
934      \param v A vector to store intermediate results.
935      \param offset Index in v to begin storing results.
936      \param sampleNo Sample number to evaluate.
937      \param roffset (output parameter) the offset in the return vector where the result begins.
938    
939      The return value will be an existing vector so do not deallocate it.
940      If the result is stored in v it should be stored at the offset given.
941      Everything from offset to the end of v should be considered available for this method to use.
942    */
943    DataTypes::ValueType*
944    DataLazy::resolveNP1OUT_P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
945    {
946        // we assume that any collapsing has been done before we get here
947        // since we only have one argument we don't need to think about only
948        // processing single points.
949      if (m_readytype!='E')
950      {
951        throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");
952      }
953        // since we can't write the result over the input, we need a result offset further along
954      size_t subroffset=roffset+m_samplesize;
955      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,subroffset);
956      roffset=offset;
957      switch (m_op)
958      {
959        case TRACE:
960             DataMaths::trace(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
961        break;
962        case TRANS:
963             DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
964        break;
965        default:
966        throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
967      }
968      return &v;
969    }
970    
 cout << "Resolve binary: " << toString() << endl;  
971    
972    size_t lroffset=0, rroffset=0;  #define PROC_OP(TYPE,X)                               \
973        for (int j=0;j<onumsteps;++j)\
974        {\
975          for (int i=0;i<numsteps;++i,resultp+=resultStep) \
976          { \
977    LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
978             tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
979             lroffset+=leftstep; \
980             rroffset+=rightstep; \
981          }\
982          lroffset+=oleftstep;\
983          rroffset+=orightstep;\
984        }
985    
986    /*
987      \brief Compute the value of the expression (binary operation) for the given sample.
988      \return Vector which stores the value of the subexpression for the given sample.
989      \param v A vector to store intermediate results.
990      \param offset Index in v to begin storing results.
991      \param sampleNo Sample number to evaluate.
992      \param roffset (output parameter) the offset in the return vector where the result begins.
993    
994      The return value will be an existing vector so do not deallocate it.
995      If the result is stored in v it should be stored at the offset given.
996      Everything from offset to the end of v should be considered available for this method to use.
997    */
998    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
999    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1000    // If both children are expanded, then we can process them in a single operation (we treat
1001    // the whole sample as one big datapoint.
1002    // If one of the children is not expanded, then we need to treat each point in the sample
1003    // individually.
1004    // There is an additional complication when scalar operations are considered.
1005    // For example, 2+Vector.
1006    // In this case each double within the point is treated individually
1007    DataTypes::ValueType*
1008    DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
1009    {
1010    LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1011    
1012      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
1013        // first work out which of the children are expanded
1014    bool leftExp=(m_left->m_readytype=='E');    bool leftExp=(m_left->m_readytype=='E');
1015    bool rightExp=(m_right->m_readytype=='E');    bool rightExp=(m_right->m_readytype=='E');
1016    bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp)); // is processing in single step    if (!leftExp && !rightExp)
1017    int steps=(bigloops?1:getNumDPPSample());    {
1018    size_t chunksize=(bigloops? m_samplesize : getNoValues());      throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1019    int leftStep=((leftExp && !rightExp)? getNoValues() : 0);    }
1020    int rightStep=((rightExp && !leftExp)? getNoValues() : 0);    bool leftScalar=(m_left->getRank()==0);
1021      bool rightScalar=(m_right->getRank()==0);
1022      if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1023      {
1024        throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1025      }
1026      size_t leftsize=m_left->getNoValues();
1027      size_t rightsize=m_right->getNoValues();
1028      size_t chunksize=1;           // how many doubles will be processed in one go
1029      int leftstep=0;       // how far should the left offset advance after each step
1030      int rightstep=0;
1031      int numsteps=0;       // total number of steps for the inner loop
1032      int oleftstep=0;  // the o variables refer to the outer loop
1033      int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1034      int onumsteps=1;
1035      
1036      bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1037      bool RES=(rightExp && rightScalar);
1038      bool LS=(!leftExp && leftScalar); // left is a single scalar
1039      bool RS=(!rightExp && rightScalar);
1040      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1041      bool RN=(!rightExp && !rightScalar);
1042      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1043      bool REN=(rightExp && !rightScalar);
1044    
1045      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1046      {
1047        chunksize=m_left->getNumDPPSample()*leftsize;
1048        leftstep=0;
1049        rightstep=0;
1050        numsteps=1;
1051      }
1052      else if (LES || RES)
1053      {
1054        chunksize=1;
1055        if (LES)        // left is an expanded scalar
1056        {
1057            if (RS)
1058            {
1059               leftstep=1;
1060               rightstep=0;
1061               numsteps=m_left->getNumDPPSample();
1062            }
1063            else        // RN or REN
1064            {
1065               leftstep=0;
1066               oleftstep=1;
1067               rightstep=1;
1068               orightstep=(RN?-rightsize:0);
1069               numsteps=rightsize;
1070               onumsteps=m_left->getNumDPPSample();
1071            }
1072        }
1073        else        // right is an expanded scalar
1074        {
1075            if (LS)
1076            {
1077               rightstep=1;
1078               leftstep=0;
1079               numsteps=m_right->getNumDPPSample();
1080            }
1081            else
1082            {
1083               rightstep=0;
1084               orightstep=1;
1085               leftstep=1;
1086               oleftstep=(LN?-leftsize:0);
1087               numsteps=leftsize;
1088               onumsteps=m_right->getNumDPPSample();
1089            }
1090        }
1091      }
1092      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1093      {
1094        if (LEN)    // and Right will be a single value
1095        {
1096            chunksize=rightsize;
1097            leftstep=rightsize;
1098            rightstep=0;
1099            numsteps=m_left->getNumDPPSample();
1100            if (RS)
1101            {
1102               numsteps*=leftsize;
1103            }
1104        }
1105        else    // REN
1106        {
1107            chunksize=leftsize;
1108            rightstep=leftsize;
1109            leftstep=0;
1110            numsteps=m_right->getNumDPPSample();
1111            if (LS)
1112            {
1113               numsteps*=rightsize;
1114            }
1115        }
1116      }
1117    
1118      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
1119        // Get the values of sub-expressions
1120    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);    const ValueType* left=m_left->resolveSample(v,offset,sampleNo,lroffset);
1121    const ValueType* right=m_right->resolveSample(v,offset,sampleNo,rroffset);        const ValueType* right=m_right->resolveSample(v,offset+m_left->getMaxSampleSize(),sampleNo,rroffset); // Note
1122      // now we need to know which args are expanded      // the right child starts further along.
1123  cout << "left=" << left << " right=" << right << endl;  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1124  cout << "(Length) l=" << left->size() << " r=" << right->size() << " res=" << v.size() << endl;  LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1125    double* resultp=&(v[offset]);  LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1126    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1127    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1128    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1129    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1130      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
1131    switch(m_op)    switch(m_op)
1132    {    {
1133      case ADD:      case ADD:
1134      for (int i=0;i<steps;++i,resultp+=getNoValues())          PROC_OP(NO_ARG,plus<double>());
1135      {      break;
1136  cerr << "Step#" << i << " chunk=" << chunksize << endl;      case SUB:
1137  cerr << left << "[" << lroffset << "] " << right << "[" << rroffset << "]" << endl;      PROC_OP(NO_ARG,minus<double>());
1138         tensor_binary_operation(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, plus<double>());      break;
1139         lroffset+=leftStep;      case MUL:
1140         rroffset+=rightStep;      PROC_OP(NO_ARG,multiplies<double>());
1141  cerr << "left=" << lroffset << " right=" << rroffset << endl;      break;
1142      }      case DIV:
1143        PROC_OP(NO_ARG,divides<double>());
1144        break;
1145        case POW:
1146           PROC_OP(double (double,double),::pow);
1147      break;      break;
 // need to fill in the rest  
1148      default:      default:
1149      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1150    }    }
1151    roffset=offset;    roffset=offset;  
1152    return &v;    return &v;
1153  }  }
1154    
1155    
1156    
1157  // #define PROC_OP(X) \  /*
1158  //  for (int i=0;i<steps;++i,resultp+=getNoValues()) \    \brief Compute the value of the expression (tensor product) for the given sample.
1159  //  { \    \return Vector which stores the value of the subexpression for the given sample.
1160  // cout << "Step#" << i << " chunk=" << chunksize << endl; \    \param v A vector to store intermediate results.
1161  // cout << left[0] << left[1] << left[2] << endl; \    \param offset Index in v to begin storing results.
1162  // cout << right[0] << right[1] << right[2] << endl; \    \param sampleNo Sample number to evaluate.
1163  //     tensor_binary_operation(chunksize, left, right, resultp, X); \    \param roffset (output parameter) the offset in the return vector where the result begins.
1164  //     left+=leftStep; \  
1165  //     right+=rightStep; \    The return value will be an existing vector so do not deallocate it.
1166  // cout << "Result=" << result << " " << result[0] << result[1] << result[2] << endl; \    If the result is stored in v it should be stored at the offset given.
1167  //  }    Everything from offset to the end of v should be considered available for this method to use.
1168  //  */
1169  // const double*  // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1170  // DataLazy::resolveBinary(ValueType& v,int sampleNo,  size_t offset) const  // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1171  // {  // unlike the other resolve helpers, we must treat these datapoints separately.
1172  //  // again we assume that all collapsing has already been done  DataTypes::ValueType*
1173  //  // so we have at least one expanded child.  DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
1174  //  // however, we could still have one of the children being not expanded.  {
1175  //  LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
1176  // cout << "Resolve binary: " << toString() << endl;  
1177  //    size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
1178  //   const double* left=m_left->resolveSample(v,sampleNo,offset);      // first work out which of the children are expanded
1179  // // cout << "Done Left " << /*left[0] << left[1] << left[2] << */endl;    bool leftExp=(m_left->m_readytype=='E');
1180  //   const double* right=m_right->resolveSample(v,sampleNo,offset);    bool rightExp=(m_right->m_readytype=='E');
1181  // // cout << "Done Right"  << /*right[0] << right[1] << right[2] <<*/ endl;    int steps=getNumDPPSample();
1182  //      // now we need to know which args are expanded    int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
1183  //   bool leftExp=(m_left->m_readytype=='E');    int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);
1184  //   bool rightExp=(m_right->m_readytype=='E');    int resultStep=max(leftStep,rightStep);   // only one (at most) should be !=0
1185  //   bool bigloops=((leftExp && rightExp) || (!leftExp && !rightExp));  // is processing in single step      // Get the values of sub-expressions (leave a gap of one sample for the result).
1186  //   int steps=(bigloops?1:getNumSamples());    int gap=offset+m_left->getMaxSampleSize();    // actually only needs to be m_left->m_samplesize
1187  //   size_t chunksize=(bigloops? m_samplesize : getNoValues());    const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);
1188  //   int leftStep=((leftExp && !rightExp)? getNoValues() : 0);    gap+=m_right->getMaxSampleSize();
1189  //   int rightStep=((rightExp && !leftExp)? getNoValues() : 0);    const ValueType* right=m_right->resolveSample(v,gap,sampleNo,rroffset);
1190  // cout << "left=" << left << " right=" << right << endl;  LAZYDEBUG(cout << "Post sub calls: " << toString() << endl;)
1191  //   double* result=&(v[offset]);  LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1192  //   double* resultp=result;  LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1193  //   switch(m_op)  LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1194  //   {  LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1195  //     case ADD:  LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1196  //  for (int i=0;i<steps;++i,resultp+=getNoValues())    double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
1197  //  {    switch(m_op)
1198  // cout << "Step#" << i << " chunk=" << chunksize << endl; \    {
1199  // // cout << left[0] << left[1] << left[2] << endl;      case PROD:
1200  // // cout << right[0] << right[1] << right[2] << endl;      for (int i=0;i<steps;++i,resultp+=resultStep)
1201  //     tensor_binary_operation(chunksize, left, right, resultp, plus<double>());      {
1202  // cout << "left=" << left << " right=" << right << " resp=" << resultp << endl;  LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)
1203  //     left+=leftStep;  LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)
1204  //     right+=rightStep;  LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)
1205  // cout << "left=" << left << " right=" << right << endl;            const double *ptr_0 = &((*left)[lroffset]);
1206  // // cout << "Result=" << result << " " << result[0] << result[1] << result[2] << endl;            const double *ptr_1 = &((*right)[rroffset]);
1207  //  }  LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1208  //  break;  LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1209  // // need to fill in the rest            matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
1210  //     default:        lroffset+=leftStep;
1211  //  throw DataException("Programmer error - resolveBinay can not resolve operator "+opToString(m_op)+".");        rroffset+=rightStep;
1212  //   }      }
1213  // // cout << "About to return "  << result[0] << result[1] << result[2] << endl;;      break;
1214  //   return result;      default:
1215  // }      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1216      }
1217      roffset=offset;
1218      return &v;
1219    }
1220    
 // // the vector and the offset are a place where the method could write its data if it wishes  
 // // it is not obligated to do so. For example, if it has its own storage already, it can use that.  
 // // Hence the return value to indicate where the data is actually stored.  
 // // Regardless, the storage should be assumed to be used, even if it isn't.  
 // const double*  
 // DataLazy::resolveSample(ValueType& v,int sampleNo,  size_t offset )  
 // {  
 // cout << "Resolve sample " << toString() << endl;  
 //  // collapse so we have a 'E' node or an IDENTITY for some other type  
 //   if (m_readytype!='E' && m_op!=IDENTITY)  
 //   {  
 //  collapse();  
 //   }  
 //   if (m_op==IDENTITY)      
 //   {  
 //     const ValueType& vec=m_id->getVector();  
 //     if (m_readytype=='C')  
 //     {  
 //  return &(vec[0]);  
 //     }  
 //     return &(vec[m_id->getPointOffset(sampleNo, 0)]);  
 //   }  
 //   if (m_readytype!='E')  
 //   {  
 //     throw DataException("Programmer Error - Collapse did not produce an expanded node.");  
 //   }  
 //   switch (getOpgroup(m_op))  
 //   {  
 //   case G_UNARY: return resolveUnary(v,sampleNo,offset);  
 //   case G_BINARY: return resolveBinary(v,sampleNo,offset);  
 //   default:  
 //     throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");  
 //   }  
 // }  
1221    
1222    
1223    /*
1224      \brief Compute the value of the expression for the given sample.
1225      \return Vector which stores the value of the subexpression for the given sample.
1226      \param v A vector to store intermediate results.
1227      \param offset Index in v to begin storing results.
1228      \param sampleNo Sample number to evaluate.
1229      \param roffset (output parameter) the offset in the return vector where the result begins.
1230    
1231      The return value will be an existing vector so do not deallocate it.
1232    */
1233  // the vector and the offset are a place where the method could write its data if it wishes  // the vector and the offset are a place where the method could write its data if it wishes
1234  // it is not obligated to do so. For example, if it has its own storage already, it can use that.  // it is not obligated to do so. For example, if it has its own storage already, it can use that.
1235  // Hence the return value to indicate where the data is actually stored.  // Hence the return value to indicate where the data is actually stored.
# Line 799  cerr << "left=" << lroffset << " right=" Line 1239  cerr << "left=" << lroffset << " right="
1239  const DataTypes::ValueType*  const DataTypes::ValueType*
1240  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)  DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
1241  {  {
1242  cout << "Resolve sample " << toString() << endl;  LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1243      // collapse so we have a 'E' node or an IDENTITY for some other type      // collapse so we have a 'E' node or an IDENTITY for some other type
1244    if (m_readytype!='E' && m_op!=IDENTITY)    if (m_readytype!='E' && m_op!=IDENTITY)
1245    {    {
# Line 811  cout << "Resolve sample " << toString() Line 1251  cout << "Resolve sample " << toString()
1251      if (m_readytype=='C')      if (m_readytype=='C')
1252      {      {
1253      roffset=0;      roffset=0;
1254    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1255      return &(vec);      return &(vec);
1256      }      }
1257      roffset=m_id->getPointOffset(sampleNo, 0);      roffset=m_id->getPointOffset(sampleNo, 0);
1258    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1259      return &(vec);      return &(vec);
1260    }    }
1261    if (m_readytype!='E')    if (m_readytype!='E')
# Line 822  cout << "Resolve sample " << toString() Line 1264  cout << "Resolve sample " << toString()
1264    }    }
1265    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1266    {    {
1267    case G_UNARY: return resolveUnary(v, offset,sampleNo,roffset);    case G_UNARY:
1268      case G_UNARY_P: return resolveUnary(v, offset,sampleNo,roffset);
1269    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);    case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
1270      case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
1271      case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);
1272      case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
1273    default:    default:
1274      throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");      throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
1275    }    }
 }  
   
1276    
1277    }
1278    
1279    
1280  // This version uses double* trying again with vectors  // To simplify the memory management, all threads operate on one large vector, rather than one each.
1281  // DataReady_ptr  // Each sample is evaluated independently and copied into the result DataExpanded.
 // DataLazy::resolve()  
 // {  
 //  
 // cout << "Sample size=" << m_samplesize << endl;  
 // cout << "Buffers=" << m_buffsRequired << endl;  
 //  
 //   if (m_readytype!='E')  
 //   {  
 //     collapse();  
 //   }  
 //   if (m_op==IDENTITY)  
 //   {  
 //     return m_id;  
 //   }  
 //      // from this point on we must have m_op!=IDENTITY and m_readytype=='E'  
 //   size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);  
 //   int numthreads=1;  
 // #ifdef _OPENMP  
 //   numthreads=getNumberOfThreads();  
 //   int threadnum=0;  
 // #endif  
 //   ValueType v(numthreads*threadbuffersize);    
 // cout << "Buffer created with size=" << v.size() << endl;  
 //   DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));  
 //   ValueType& resvec=result->getVector();  
 //   DataReady_ptr resptr=DataReady_ptr(result);  
 //   int sample;  
 //   int resoffset;  
 //   int totalsamples=getNumSamples();  
 //   const double* res=0;  
 //   #pragma omp parallel for private(sample,resoffset,threadnum,res) schedule(static)  
 //   for (sample=0;sample<totalsamples;++sample)  
 //   {  
 // cout << "################################# " << sample << endl;  
 // #ifdef _OPENMP  
 //     res=resolveSample(v,sample,threadbuffersize*omp_get_thread_num());  
 // #else  
 //     res=resolveSample(v,sample,0);   // this would normally be v, but not if its a single IDENTITY op.  
 // #endif  
 // cerr << "-------------------------------- " << endl;  
 //     resoffset=result->getPointOffset(sample,0);  
 // cerr << "offset=" << resoffset << endl;  
 //     for (unsigned int i=0;i<m_samplesize;++i,++resoffset)    // copy values into the output vector  
 //     {  
 //  resvec[resoffset]=res[i];  
 //     }  
 // cerr << "*********************************" << endl;  
 //   }  
 //   return resptr;  
 // }  
   
   
1282  DataReady_ptr  DataReady_ptr
1283  DataLazy::resolve()  DataLazy::resolve()
1284  {  {
1285    
1286  cout << "Sample size=" << m_samplesize << endl;  LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)
1287  cout << "Buffers=" << m_buffsRequired << endl;  LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)
1288    
1289    if (m_readytype!='E')    if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally
1290    {    {
1291      collapse();      collapse();
1292    }    }
1293    if (m_op==IDENTITY)    if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
1294    {    {
1295      return m_id;      return m_id;
1296    }    }
1297      // from this point on we must have m_op!=IDENTITY and m_readytype=='E'      // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
1298    size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);    size_t threadbuffersize=m_maxsamplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough
1299        // storage to evaluate its expression
1300    int numthreads=1;    int numthreads=1;
1301  #ifdef _OPENMP  #ifdef _OPENMP
1302    numthreads=getNumberOfThreads();    numthreads=getNumberOfThreads();
   int threadnum=0;  
1303  #endif  #endif
1304    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
1305  cout << "Buffer created with size=" << v.size() << endl;  LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)
1306    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1307    ValueType& resvec=result->getVector();    ValueType& resvec=result->getVector();
1308    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1309    int sample;    int sample;
1310    size_t outoffset;     // offset in the output data    size_t outoffset;     // offset in the output data
1311    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1312    const ValueType* res=0;    const ValueType* res=0;   // Vector storing the answer
1313    size_t resoffset=0;    size_t resoffset=0;       // where in the vector to find the answer
1314    #pragma omp parallel for private(sample,resoffset,outoffset,threadnum,res) schedule(static)  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1315      #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
1316    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
1317    {    {
1318  cout << "################################# " << sample << endl;  LAZYDEBUG(cout << "################################# " << sample << endl;)
1319  #ifdef _OPENMP  #ifdef _OPENMP
1320      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
1321  #else  #else
1322      res=resolveSample(v,0,sample,resoffset);   // this would normally be v, but not if its a single IDENTITY op.      res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.
1323  #endif  #endif
1324  cerr << "-------------------------------- " << endl;  LAZYDEBUG(cerr << "-------------------------------- " << endl;)
1325      outoffset=result->getPointOffset(sample,0);      outoffset=result->getPointOffset(sample,0);
1326  cerr << "offset=" << outoffset << endl;  LAZYDEBUG(cerr << "offset=" << outoffset << endl;)
1327      for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector      for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector
1328      {      {
1329      resvec[outoffset]=(*res)[resoffset];      resvec[outoffset]=(*res)[resoffset];
1330      }      }
1331  cerr << "*********************************" << endl;  LAZYDEBUG(cerr << "*********************************" << endl;)
1332    }    }
1333    return resptr;    return resptr;
1334  }  }
# Line 948  DataLazy::toString() const Line 1342  DataLazy::toString() const
1342    return oss.str();    return oss.str();
1343  }  }
1344    
1345    
1346  void  void
1347  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
1348  {  {
# Line 980  DataLazy::intoString(ostringstream& oss) Line 1375  DataLazy::intoString(ostringstream& oss)
1375      oss << ')';      oss << ')';
1376      break;      break;
1377    case G_UNARY:    case G_UNARY:
1378      case G_UNARY_P:
1379      case G_NP1OUT:
1380      case G_NP1OUT_P:
1381      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
1382      m_left->intoString(oss);      m_left->intoString(oss);
1383      oss << ')';      oss << ')';
1384      break;      break;
1385      case G_TENSORPROD:
1386        oss << opToString(m_op) << '(';
1387        m_left->intoString(oss);
1388        oss << ", ";
1389        m_right->intoString(oss);
1390        oss << ')';
1391        break;
1392    default:    default:
1393      oss << "UNKNOWN";      oss << "UNKNOWN";
1394    }    }
1395  }  }
1396    
 // 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.  
1397  DataAbstract*  DataAbstract*
1398  DataLazy::deepCopy()  DataLazy::deepCopy()
1399  {  {
1400    if (m_op==IDENTITY)    switch (getOpgroup(m_op))
1401    {    {
1402      return new DataLazy(m_left);    // we don't need to copy the child here    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1403      case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1404      case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1405      case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1406      case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1407      default:
1408        throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1409    }    }
   return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);  
1410  }  }
1411    
1412    
1413    // There is no single, natural interpretation of getLength on DataLazy.
1414    // Instances of DataReady can look at the size of their vectors.
1415    // For lazy though, it could be the size the data would be if it were resolved;
1416    // or it could be some function of the lengths of the DataReady instances which
1417    // form part of the expression.
1418    // Rather than have people making assumptions, I have disabled the method.
1419  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1420  DataLazy::getLength() const  DataLazy::getLength() const
1421  {  {
1422    return m_length;    throw DataException("getLength() does not make sense for lazy data.");
1423  }  }
1424    
1425    
# Line 1015  DataLazy::getSlice(const DataTypes::Regi Line 1429  DataLazy::getSlice(const DataTypes::Regi
1429    throw DataException("getSlice - not implemented for Lazy objects.");    throw DataException("getSlice - not implemented for Lazy objects.");
1430  }  }
1431    
1432    
1433    // To do this we need to rely on our child nodes
1434    DataTypes::ValueType::size_type
1435    DataLazy::getPointOffset(int sampleNo,
1436                     int dataPointNo)
1437    {
1438      if (m_op==IDENTITY)
1439      {
1440        return m_id->getPointOffset(sampleNo,dataPointNo);
1441      }
1442      if (m_readytype!='E')
1443      {
1444        collapse();
1445        return m_id->getPointOffset(sampleNo,dataPointNo);
1446      }
1447      // at this point we do not have an identity node and the expression will be Expanded
1448      // so we only need to know which child to ask
1449      if (m_left->m_readytype=='E')
1450      {
1451        return m_left->getPointOffset(sampleNo,dataPointNo);
1452      }
1453      else
1454      {
1455        return m_right->getPointOffset(sampleNo,dataPointNo);
1456      }
1457    }
1458    
1459    // To do this we need to rely on our child nodes
1460  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1461  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
1462                   int dataPointNo) const                   int dataPointNo) const
1463  {  {
1464    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");    if (m_op==IDENTITY)
1465      {
1466        return m_id->getPointOffset(sampleNo,dataPointNo);
1467      }
1468      if (m_readytype=='E')
1469      {
1470        // at this point we do not have an identity node and the expression will be Expanded
1471        // so we only need to know which child to ask
1472        if (m_left->m_readytype=='E')
1473        {
1474        return m_left->getPointOffset(sampleNo,dataPointNo);
1475        }
1476        else
1477        {
1478        return m_right->getPointOffset(sampleNo,dataPointNo);
1479        }
1480      }
1481      if (m_readytype=='C')
1482      {
1483        return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
1484      }
1485      throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1486    }
1487    
1488    // It would seem that DataTagged will need to be treated differently since even after setting all tags
1489    // to zero, all the tags from all the DataTags would be in the result.
1490    // However since they all have the same value (0) whether they are there or not should not matter.
1491    // So I have decided that for all types this method will create a constant 0.
1492    // It can be promoted up as required.
1493    // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
1494    // but we can deal with that if it arrises.
1495    void
1496    DataLazy::setToZero()
1497    {
1498      DataTypes::ValueType v(getNoValues(),0);
1499      m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1500      m_op=IDENTITY;
1501      m_right.reset();  
1502      m_left.reset();
1503      m_readytype='C';
1504      m_buffsRequired=1;
1505  }  }
1506    
1507  }   // end namespace  }   // end namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26