/[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 2271 by jfenwick, Mon Feb 16 05:08:29 2009 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    #include <iomanip>      // for some fancy formatting in debug
32    
33    // #define LAZYDEBUG(X) if (privdebug){X;}
34    #define LAZYDEBUG(X)
35    namespace
36    {
37    bool privdebug=false;
38    
39    #define ENABLEDEBUG privdebug=true;
40    #define DISABLEDEBUG privdebug=false;
41    }
42    
43    // #define SIZELIMIT
44    // #define SIZELIMIT if ((m_height>7) || (m_children>127)) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl; resolveToIdentity();}
45    #define SIZELIMIT if ((m_height>7) || (m_children>127)) {resolveToIdentity();}
46    
47    /*
48    How does DataLazy work?
49    ~~~~~~~~~~~~~~~~~~~~~~~
50    
51    Each instance represents a single operation on one or two other DataLazy instances. These arguments are normally
52    denoted left and right.
53    
54    A special operation, IDENTITY, stores an instance of DataReady in the m_id member.
55    This means that all "internal" nodes in the structure are instances of DataLazy.
56    
57    Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
58    Note that IDENITY is not considered a unary operation.
59    
60    I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
61    It must however form a DAG (directed acyclic graph).
62    I will refer to individual DataLazy objects with the structure as nodes.
63    
64    Each node also stores:
65    - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
66        evaluated.
67    - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to
68        evaluate the expression.
69    - m_samplesize ~ the number of doubles stored in a sample.
70    
71    When a new node is created, the above values are computed based on the values in the child nodes.
72    Eg: if left requires 4 samples and right requires 6 then left+right requires 7 samples.
73    
74    The resolve method, which produces a DataReady from a DataLazy, does the following:
75    1) Create a DataReady to hold the new result.
76    2) Allocate a vector (v) big enough to hold m_buffsrequired samples.
77    3) For each sample, call resolveSample with v, to get its values and copy them into the result object.
78    
79    (In the case of OMP, multiple samples are resolved in parallel so the vector needs to be larger.)
80    
81    resolveSample returns a Vector* and an offset within that vector where the result is stored.
82    Normally, this would be v, but for identity nodes their internal vector is returned instead.
83    
84    The convention that I use, is that the resolve methods should store their results starting at the offset they are passed.
85    
86    For expressions which evaluate to Constant or Tagged, there is a different evaluation method.
87    The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.
88    
89    To add a new operator you need to do the following (plus anything I might have forgotten - adding a new group for example):
90    1) Add to the ES_optype.
91    2) determine what opgroup your operation belongs to (X)
92    3) add a string for the op to the end of ES_opstrings
93    4) increase ES_opcount
94    5) add an entry (X) to opgroups
95    6) add an entry to the switch in collapseToReady
96    7) add an entry to resolveX
97    */
98    
99    
100  using namespace std;  using namespace std;
101  using namespace boost;  using namespace boost;
# Line 33  using namespace boost; Line 103  using namespace boost;
103  namespace escript  namespace escript
104  {  {
105    
 const std::string&  
 opToString(ES_optype op);  
   
106  namespace  namespace
107  {  {
108    
   
   
109  enum ES_opgroup  enum ES_opgroup
110  {  {
111     G_UNKNOWN,     G_UNKNOWN,
112     G_IDENTITY,     G_IDENTITY,
113     G_BINARY,     G_BINARY,        // pointwise operations with two arguments
114     G_UNARY     G_UNARY,     // pointwise operations with one argument
115       G_UNARY_P,       // pointwise operations with one argument, requiring a parameter
116       G_NP1OUT,        // non-pointwise op with one output
117       G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
118       G_TENSORPROD     // general tensor product
119  };  };
120    
121    
122    
123    
124  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","sin","cos","tan",  string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
125                "sin","cos","tan",
126              "asin","acos","atan","sinh","cosh","tanh","erf",              "asin","acos","atan","sinh","cosh","tanh","erf",
127              "asinh","acosh","atanh",              "asinh","acosh","atanh",
128              "log10","log","sign","abs","neg","pos","exp","sqrt",              "log10","log","sign","abs","neg","pos","exp","sqrt",
129              "1/","where>0","where<0","where>=0","where<=0"};              "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
130  int ES_opcount=32;              "symmetric","nonsymmetric",
131  ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY,G_UNARY,G_UNARY,G_UNARY, //9              "prod",
132              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 16              "transpose", "trace"};
133              G_UNARY,G_UNARY,G_UNARY,                    // 19  int ES_opcount=40;
134              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,
135              G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY};              G_UNARY,G_UNARY,G_UNARY, //10
136                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17
137                G_UNARY,G_UNARY,G_UNARY,                    // 20
138                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 28
139                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P,      // 35
140                G_NP1OUT,G_NP1OUT,
141                G_TENSORPROD,
142                G_NP1OUT_P, G_NP1OUT_P};
143  inline  inline
144  ES_opgroup  ES_opgroup
145  getOpgroup(ES_optype op)  getOpgroup(ES_optype op)
# Line 79  resultFS(DataAbstract_ptr left, DataAbst Line 156  resultFS(DataAbstract_ptr left, DataAbst
156      // 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
157      // programming error exception.      // programming error exception.
158    
159      FunctionSpace l=left->getFunctionSpace();
160      if (left->getFunctionSpace()!=right->getFunctionSpace())    FunctionSpace r=right->getFunctionSpace();
161      {    if (l!=r)
162          throw DataException("FunctionSpaces not equal - interpolation not supported on lazy data.");    {
163      }      if (r.probeInterpolation(l))
164      return left->getFunctionSpace();      {
165        return l;
166        }
167        if (l.probeInterpolation(r))
168        {
169        return r;
170        }
171        throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
172      }
173      return l;
174  }  }
175    
176  // return the shape of the result of "left op right"  // return the shape of the result of "left op right"
177    // the shapes resulting from tensor product are more complex to compute so are worked out elsewhere
178  DataTypes::ShapeType  DataTypes::ShapeType
179  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
180  {  {
181      if (left->getShape()!=right->getShape())      if (left->getShape()!=right->getShape())
182      {      {
183          throw DataException("Shapes not the same - shapes must match for lazy data.");        if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
184          {
185            throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
186          }
187          if (left->getRank()==0)   // we need to allow scalar * anything
188          {
189            return right->getShape();
190          }
191          if (right->getRank()==0)
192          {
193            return left->getShape();
194          }
195          throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
196      }      }
197      return left->getShape();      return left->getShape();
198  }  }
199    
200  size_t  // return the shape for "op left"
201  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  
202    DataTypes::ShapeType
203    resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
204  {  {
205     switch (getOpgroup(op))      switch(op)
206     {      {
207     case G_BINARY: return left->getLength();          case TRANS:
208     case G_UNARY: return left->getLength();         {            // for the scoping of variables
209     default:          const DataTypes::ShapeType& s=left->getShape();
210      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");          DataTypes::ShapeType sh;
211     }          int rank=left->getRank();
212            if (axis_offset<0 || axis_offset>rank)
213            {
214                   throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
215                }
216                for (int i=0; i<rank; i++)
217            {
218               int index = (axis_offset+i)%rank;
219                   sh.push_back(s[index]); // Append to new shape
220                }
221            return sh;
222           }
223        break;
224        case TRACE:
225           {
226            int rank=left->getRank();
227            if (rank<2)
228            {
229               throw DataException("Trace can only be computed for objects with rank 2 or greater.");
230            }
231            if ((axis_offset>rank-2) || (axis_offset<0))
232            {
233               throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
234            }
235            if (rank==2)
236            {
237               return DataTypes::scalarShape;
238            }
239            else if (rank==3)
240            {
241               DataTypes::ShapeType sh;
242                   if (axis_offset==0)
243               {
244                    sh.push_back(left->getShape()[2]);
245                   }
246                   else     // offset==1
247               {
248                sh.push_back(left->getShape()[0]);
249                   }
250               return sh;
251            }
252            else if (rank==4)
253            {
254               DataTypes::ShapeType sh;
255               const DataTypes::ShapeType& s=left->getShape();
256                   if (axis_offset==0)
257               {
258                    sh.push_back(s[2]);
259                    sh.push_back(s[3]);
260                   }
261                   else if (axis_offset==1)
262               {
263                    sh.push_back(s[0]);
264                    sh.push_back(s[3]);
265                   }
266               else     // offset==2
267               {
268                sh.push_back(s[0]);
269                sh.push_back(s[1]);
270               }
271               return sh;
272            }
273            else        // unknown rank
274            {
275               throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
276            }
277           }
278        break;
279            default:
280        throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
281        }
282    }
283    
284    // determine the output shape for the general tensor product operation
285    // the additional parameters return information required later for the product
286    // the majority of this code is copy pasted from C_General_Tensor_Product
287    DataTypes::ShapeType
288    GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
289    {
290        
291      // Get rank and shape of inputs
292      int rank0 = left->getRank();
293      int rank1 = right->getRank();
294      const DataTypes::ShapeType& shape0 = left->getShape();
295      const DataTypes::ShapeType& shape1 = right->getShape();
296    
297      // Prepare for the loops of the product and verify compatibility of shapes
298      int start0=0, start1=0;
299      if (transpose == 0)       {}
300      else if (transpose == 1)  { start0 = axis_offset; }
301      else if (transpose == 2)  { start1 = rank1-axis_offset; }
302      else              { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
303    
304      if (rank0<axis_offset)
305      {
306        throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
307      }
308    
309      // Adjust the shapes for transpose
310      DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather
311      DataTypes::ShapeType tmpShape1(rank1);    // than using push_back
312      for (int i=0; i<rank0; i++)   { tmpShape0[i]=shape0[(i+start0)%rank0]; }
313      for (int i=0; i<rank1; i++)   { tmpShape1[i]=shape1[(i+start1)%rank1]; }
314    
315      // Prepare for the loops of the product
316      SL=1, SM=1, SR=1;
317      for (int i=0; i<rank0-axis_offset; i++)   {
318        SL *= tmpShape0[i];
319      }
320      for (int i=rank0-axis_offset; i<rank0; i++)   {
321        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
322          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
323        }
324        SM *= tmpShape0[i];
325      }
326      for (int i=axis_offset; i<rank1; i++)     {
327        SR *= tmpShape1[i];
328      }
329    
330      // Define the shape of the output (rank of shape is the sum of the loop ranges below)
331      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);  
332      {         // block to limit the scope of out_index
333         int out_index=0;
334         for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
335         for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
336      }
337    
338      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
339      {
340         ostringstream os;
341         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
342         throw DataException(os.str());
343      }
344    
345      return shape2;
346  }  }
347    
348    // determine the number of samples requires to evaluate an expression combining left and right
349    // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
350    // The same goes for G_TENSORPROD
351    // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefts.
352    // This would be true were it not for the possibility that the LHS could be a scalar which needs to be examined
353    // multiple times
354  int  int
355  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)  calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
356  {  {
357     switch(getOpgroup(op))     switch(getOpgroup(op))
358     {     {
359     case G_IDENTITY: return 1;     case G_IDENTITY: return 1;
360     case G_BINARY: return max(left->getBuffsRequired(),right->getBuffsRequired()+1);     case G_BINARY: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
361     case G_UNARY: return max(left->getBuffsRequired(),1);     case G_UNARY:
362       case G_UNARY_P: return max(left->getBuffsRequired(),1);
363       case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
364       case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);
365       case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
366     default:     default:
367      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
368     }     }
369  }  }
370    
371    
   
372  }   // end anonymous namespace  }   // end anonymous namespace
373    
374    
375    
376    // Return a string representing the operation
377  const std::string&  const std::string&
378  opToString(ES_optype op)  opToString(ES_optype op)
379  {  {
# Line 140  opToString(ES_optype op) Line 386  opToString(ES_optype op)
386    
387    
388  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
389      : parent(p->getFunctionSpace(),p->getShape()),      : parent(p->getFunctionSpace(),p->getShape())
     m_op(IDENTITY)  
390  {  {
391     if (p->isLazy())     if (p->isLazy())
392     {     {
     // TODO: fix this.   We could make the new node a copy of p?  
393      // I don't want identity of Lazy.      // I don't want identity of Lazy.
394      // Question: Why would that be so bad?      // Question: Why would that be so bad?
395      // 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 153  DataLazy::DataLazy(DataAbstract_ptr p) Line 397  DataLazy::DataLazy(DataAbstract_ptr p)
397     }     }
398     else     else
399     {     {
400      m_id=dynamic_pointer_cast<DataReady>(p);      p->makeLazyShared();
401        DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
402        makeIdentity(dr);
403    LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
404     }     }
405     m_length=p->getLength();  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
    m_buffsRequired=1;  
    m_samplesize=getNumDPPSample()*getNoValues();  
 cout << "(1)Lazy created with " << m_samplesize << endl;  
406  }  }
407    
408    
409    
410    
411  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
412      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
413      m_op(op)      m_op(op),
414        m_axis_offset(0),
415        m_transpose(0),
416        m_SL(0), m_SM(0), m_SR(0)
417  {  {
418     if (getOpgroup(op)!=G_UNARY)     if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
419     {     {
420      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.");
421     }     }
422    
423     DataLazy_ptr lleft;     DataLazy_ptr lleft;
424     if (!left->isLazy())     if (!left->isLazy())
425     {     {
# Line 178  DataLazy::DataLazy(DataAbstract_ptr left Line 429  DataLazy::DataLazy(DataAbstract_ptr left
429     {     {
430      lleft=dynamic_pointer_cast<DataLazy>(left);      lleft=dynamic_pointer_cast<DataLazy>(left);
431     }     }
432     m_length=left->getLength();     m_readytype=lleft->m_readytype;
433     m_left=lleft;     m_left=lleft;
434     m_buffsRequired=1;     m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
435     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
436       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
437       m_children=m_left->m_children+1;
438       m_height=m_left->m_height+1;
439       SIZELIMIT
440  }  }
441    
442    
443  DataLazy::DataLazy(DataLazy_ptr left, DataLazy_ptr right, ES_optype op)  // In this constructor we need to consider interpolation
444    DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
445      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
446      m_left(left),      m_op(op),
447      m_right(right),      m_SL(0), m_SM(0), m_SR(0)
     m_op(op)  
448  {  {
449     if (getOpgroup(op)!=G_BINARY)  LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
450       if ((getOpgroup(op)!=G_BINARY))
451     {     {
452      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.");
453     }     }
454     m_length=resultLength(m_left,m_right,m_op);  
455       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
456       {
457        FunctionSpace fs=getFunctionSpace();
458        Data ltemp(left);
459        Data tmp(ltemp,fs);
460        left=tmp.borrowDataPtr();
461       }
462       if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
463       {
464        Data tmp(Data(right),getFunctionSpace());
465        right=tmp.borrowDataPtr();
466    LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
467       }
468       left->operandCheck(*right);
469    
470       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
471       {
472        m_left=dynamic_pointer_cast<DataLazy>(left);
473    LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
474       }
475       else
476       {
477        m_left=DataLazy_ptr(new DataLazy(left));
478    LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
479       }
480       if (right->isLazy())
481       {
482        m_right=dynamic_pointer_cast<DataLazy>(right);
483    LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
484       }
485       else
486       {
487        m_right=DataLazy_ptr(new DataLazy(right));
488    LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
489       }
490       char lt=m_left->m_readytype;
491       char rt=m_right->m_readytype;
492       if (lt=='E' || rt=='E')
493       {
494        m_readytype='E';
495       }
496       else if (lt=='T' || rt=='T')
497       {
498        m_readytype='T';
499       }
500       else
501       {
502        m_readytype='C';
503       }
504     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
505     m_buffsRequired=calcBuffs(m_left, m_right, m_op);     m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
506  cout << "(2)Lazy created with " << m_samplesize << endl;     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
507       m_children=m_left->m_children+m_right->m_children+2;
508       m_height=max(m_left->m_height,m_right->m_height)+1;
509       SIZELIMIT
510    LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
511  }  }
512    
513  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)
514      : 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)),
515      m_op(op)      m_op(op),
516        m_axis_offset(axis_offset),
517        m_transpose(transpose)
518  {  {
519     if (getOpgroup(op)!=G_BINARY)     if ((getOpgroup(op)!=G_TENSORPROD))
520       {
521        throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
522       }
523       if ((transpose>2) || (transpose<0))
524       {
525        throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
526       }
527       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
528     {     {
529      throw DataException("Programmer error - constructor DataLazy(left, op) will only process BINARY operations.");      FunctionSpace fs=getFunctionSpace();
530        Data ltemp(left);
531        Data tmp(ltemp,fs);
532        left=tmp.borrowDataPtr();
533     }     }
534     if (left->isLazy())     if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
535       {
536        Data tmp(Data(right),getFunctionSpace());
537        right=tmp.borrowDataPtr();
538       }
539    //    left->operandCheck(*right);
540    
541       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
542     {     {
543      m_left=dynamic_pointer_cast<DataLazy>(left);      m_left=dynamic_pointer_cast<DataLazy>(left);
544     }     }
# Line 225  DataLazy::DataLazy(DataAbstract_ptr left Line 554  DataLazy::DataLazy(DataAbstract_ptr left
554     {     {
555      m_right=DataLazy_ptr(new DataLazy(right));      m_right=DataLazy_ptr(new DataLazy(right));
556     }     }
557       char lt=m_left->m_readytype;
558     m_length=resultLength(m_left,m_right,m_op);     char rt=m_right->m_readytype;
559       if (lt=='E' || rt=='E')
560       {
561        m_readytype='E';
562       }
563       else if (lt=='T' || rt=='T')
564       {
565        m_readytype='T';
566       }
567       else
568       {
569        m_readytype='C';
570       }
571     m_samplesize=getNumDPPSample()*getNoValues();     m_samplesize=getNumDPPSample()*getNoValues();
572       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
573     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
574  cout << "(3)Lazy created with " << m_samplesize << endl;     m_children=m_left->m_children+m_right->m_children+2;
575       m_height=max(m_left->m_height,m_right->m_height)+1;
576       SIZELIMIT
577    LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
578    }
579    
580    
581    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
582        : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
583        m_op(op),
584        m_axis_offset(axis_offset),
585        m_transpose(0),
586        m_tol(0)
587    {
588       if ((getOpgroup(op)!=G_NP1OUT_P))
589       {
590        throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
591       }
592       DataLazy_ptr lleft;
593       if (!left->isLazy())
594       {
595        lleft=DataLazy_ptr(new DataLazy(left));
596       }
597       else
598       {
599        lleft=dynamic_pointer_cast<DataLazy>(left);
600       }
601       m_readytype=lleft->m_readytype;
602       m_left=lleft;
603       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
604       m_samplesize=getNumDPPSample()*getNoValues();
605       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
606       m_children=m_left->m_children+1;
607       m_height=m_left->m_height+1;
608       SIZELIMIT
609    LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
610  }  }
611    
612    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
613        : parent(left->getFunctionSpace(), left->getShape()),
614        m_op(op),
615        m_axis_offset(0),
616        m_transpose(0),
617        m_tol(tol)
618    {
619       if ((getOpgroup(op)!=G_UNARY_P))
620       {
621        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
622       }
623       DataLazy_ptr lleft;
624       if (!left->isLazy())
625       {
626        lleft=DataLazy_ptr(new DataLazy(left));
627       }
628       else
629       {
630        lleft=dynamic_pointer_cast<DataLazy>(left);
631       }
632       m_readytype=lleft->m_readytype;
633       m_left=lleft;
634       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
635       m_samplesize=getNumDPPSample()*getNoValues();
636       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
637       m_children=m_left->m_children+1;
638       m_height=m_left->m_height+1;
639       SIZELIMIT
640    LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
641    }
642    
643  DataLazy::~DataLazy()  DataLazy::~DataLazy()
644  {  {
# Line 245  DataLazy::getBuffsRequired() const Line 652  DataLazy::getBuffsRequired() const
652  }  }
653    
654    
655  // the vector and the offset are a place where the method could write its data if it wishes  size_t
656  // 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  
657  {  {
658    if (m_op==IDENTITY)        return m_maxsamplesize;
659    }
660    
661    
662    
663    size_t
664    DataLazy::getSampleBufferSize() const
665    {
666        return m_maxsamplesize*(max(1,m_buffsRequired));
667    }
668    
669    /*
670      \brief Evaluates the expression using methods on Data.
671      This does the work for the collapse method.
672      For reasons of efficiency do not call this method on DataExpanded nodes.
673    */
674    DataReady_ptr
675    DataLazy::collapseToReady()
676    {
677      if (m_readytype=='E')
678      { // this is more an efficiency concern than anything else
679        throw DataException("Programmer Error - do not use collapse on Expanded data.");
680      }
681      if (m_op==IDENTITY)
682    {    {
683      const ValueType& vec=m_id->getVector();      return m_id;
     return &(vec[m_id->getPointOffset(sampleNo, 0)]);  
684    }    }
685    size_t rightoffset=offset+m_samplesize;    DataReady_ptr pleft=m_left->collapseToReady();
686    const double* left=m_left->resolveSample(v,sampleNo,offset);    Data left(pleft);
687    const double* right=0;    Data right;
688    if (getOpgroup(m_op)==G_BINARY)    if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
689    {    {
690      right=m_right->resolveSample(v,sampleNo,rightoffset);      right=Data(m_right->collapseToReady());
691    }    }
692    double* result=&(v[offset]);    Data result;
693      switch(m_op)
694    {    {
695      switch(m_op)      case ADD:
696      {      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>());  
697      break;      break;
698      case SUB:            case SUB:      
699      tensor_binary_operation(m_samplesize, left, right, result, minus<double>());      result=left-right;
700      break;      break;
701      case MUL:            case MUL:      
702      tensor_binary_operation(m_samplesize, left, right, result, multiplies<double>());      result=left*right;
703      break;      break;
704      case DIV:            case DIV:      
705      tensor_binary_operation(m_samplesize, left, right, result, divides<double>());      result=left/right;
706      break;      break;
 // unary ops  
707      case SIN:      case SIN:
708      tensor_unary_operation(m_samplesize, left, result, ::sin);      result=left.sin();  
709        break;
710        case COS:
711        result=left.cos();
712        break;
713        case TAN:
714        result=left.tan();
715        break;
716        case ASIN:
717        result=left.asin();
718        break;
719        case ACOS:
720        result=left.acos();
721        break;
722        case ATAN:
723        result=left.atan();
724        break;
725        case SINH:
726        result=left.sinh();
727        break;
728        case COSH:
729        result=left.cosh();
730        break;
731        case TANH:
732        result=left.tanh();
733        break;
734        case ERF:
735        result=left.erf();
736        break;
737       case ASINH:
738        result=left.asinh();
739        break;
740       case ACOSH:
741        result=left.acosh();
742        break;
743       case ATANH:
744        result=left.atanh();
745        break;
746        case LOG10:
747        result=left.log10();
748        break;
749        case LOG:
750        result=left.log();
751        break;
752        case SIGN:
753        result=left.sign();
754        break;
755        case ABS:
756        result=left.abs();
757        break;
758        case NEG:
759        result=left.neg();
760        break;
761        case POS:
762        // it doesn't mean anything for delayed.
763        // it will just trigger a deep copy of the lazy object
764        throw DataException("Programmer error - POS not supported for lazy data.");
765        break;
766        case EXP:
767        result=left.exp();
768        break;
769        case SQRT:
770        result=left.sqrt();
771        break;
772        case RECIP:
773        result=left.oneOver();
774        break;
775        case GZ:
776        result=left.wherePositive();
777        break;
778        case LZ:
779        result=left.whereNegative();
780        break;
781        case GEZ:
782        result=left.whereNonNegative();
783        break;
784        case LEZ:
785        result=left.whereNonPositive();
786        break;
787        case NEZ:
788        result=left.whereNonZero(m_tol);
789        break;
790        case EZ:
791        result=left.whereZero(m_tol);
792        break;
793        case SYM:
794        result=left.symmetric();
795        break;
796        case NSYM:
797        result=left.nonsymmetric();
798        break;
799        case PROD:
800        result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
801        break;
802        case TRANS:
803        result=left.transpose(m_axis_offset);
804        break;
805        case TRACE:
806        result=left.trace(m_axis_offset);
807        break;
808        default:
809        throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
810      }
811      return result.borrowReadyPtr();
812    }
813    
814    /*
815       \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
816       This method uses the original methods on the Data class to evaluate the expressions.
817       For this reason, it should not be used on DataExpanded instances. (To do so would defeat
818       the purpose of using DataLazy in the first place).
819    */
820    void
821    DataLazy::collapse()
822    {
823      if (m_op==IDENTITY)
824      {
825        return;
826      }
827      if (m_readytype=='E')
828      { // this is more an efficiency concern than anything else
829        throw DataException("Programmer Error - do not use collapse on Expanded data.");
830      }
831      m_id=collapseToReady();
832      m_op=IDENTITY;
833    }
834    
835    /*
836      \brief Compute the value of the expression (unary operation) for the given sample.
837      \return Vector which stores the value of the subexpression for the given sample.
838      \param v A vector to store intermediate results.
839      \param offset Index in v to begin storing results.
840      \param sampleNo Sample number to evaluate.
841      \param roffset (output parameter) the offset in the return vector where the result begins.
842    
843      The return value will be an existing vector so do not deallocate it.
844      If the result is stored in v it should be stored at the offset given.
845      Everything from offset to the end of v should be considered available for this method to use.
846    */
847    DataTypes::ValueType*
848    DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
849    {
850        // we assume that any collapsing has been done before we get here
851        // since we only have one argument we don't need to think about only
852        // processing single points.
853      if (m_readytype!='E')
854      {
855        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
856      }
857      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
858      const double* left=&((*vleft)[roffset]);
859      double* result=&(v[offset]);
860      roffset=offset;
861      switch (m_op)
862      {
863        case SIN:  
864        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
865      break;      break;
866      case COS:      case COS:
867      tensor_unary_operation(m_samplesize, left, result, ::cos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
868      break;      break;
869      case TAN:      case TAN:
870      tensor_unary_operation(m_samplesize, left, result, ::tan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
871      break;      break;
872      case ASIN:      case ASIN:
873      tensor_unary_operation(m_samplesize, left, result, ::asin);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
874      break;      break;
875      case ACOS:      case ACOS:
876      tensor_unary_operation(m_samplesize, left, result, ::acos);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
877      break;      break;
878      case ATAN:      case ATAN:
879      tensor_unary_operation(m_samplesize, left, result, ::atan);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
880      break;      break;
881      case SINH:      case SINH:
882      tensor_unary_operation(m_samplesize, left, result, ::sinh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
883      break;      break;
884      case COSH:      case COSH:
885      tensor_unary_operation(m_samplesize, left, result, ::cosh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
886      break;      break;
887      case TANH:      case TANH:
888      tensor_unary_operation(m_samplesize, left, result, ::tanh);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
889      break;      break;
890      case ERF:      case ERF:
891  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
892      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
893  #else  #else
894      tensor_unary_operation(m_samplesize, left, result, ::erf);      tensor_unary_operation(m_samplesize, left, result, ::erf);
895      break;      break;
896  #endif  #endif
897     case ASINH:     case ASINH:
898  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
899      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
900  #else  #else
901      tensor_unary_operation(m_samplesize, left, result, ::asinh);      tensor_unary_operation(m_samplesize, left, result, ::asinh);
902  #endif    #endif  
903      break;      break;
904     case ACOSH:     case ACOSH:
905  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
906      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
907  #else  #else
908      tensor_unary_operation(m_samplesize, left, result, ::acosh);      tensor_unary_operation(m_samplesize, left, result, ::acosh);
909  #endif    #endif  
910      break;      break;
911     case ATANH:     case ATANH:
912  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
913      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);      tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
914  #else  #else
915      tensor_unary_operation(m_samplesize, left, result, ::atanh);      tensor_unary_operation(m_samplesize, left, result, ::atanh);
916  #endif    #endif  
917      break;      break;
918      case LOG10:      case LOG10:
919      tensor_unary_operation(m_samplesize, left, result, ::log10);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
920      break;      break;
921      case LOG:      case LOG:
922      tensor_unary_operation(m_samplesize, left, result, ::log);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
923      break;      break;
924      case SIGN:      case SIGN:
925      tensor_unary_operation(m_samplesize, left, result, escript::fsign);      tensor_unary_operation(m_samplesize, left, result, escript::fsign);
926      break;      break;
927      case ABS:      case ABS:
928      tensor_unary_operation(m_samplesize, left, result, ::fabs);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
929      break;      break;
930      case NEG:      case NEG:
931      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 936  DataLazy::resolveSample(ValueType& v,int
936      throw DataException("Programmer error - POS not supported for lazy data.");      throw DataException("Programmer error - POS not supported for lazy data.");
937      break;      break;
938      case EXP:      case EXP:
939      tensor_unary_operation(m_samplesize, left, result, ::exp);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
940      break;      break;
941      case SQRT:      case SQRT:
942      tensor_unary_operation(m_samplesize, left, result, ::sqrt);      tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
943      break;      break;
944      case RECIP:      case RECIP:
945      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));      tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
# Line 377  DataLazy::resolveSample(ValueType& v,int Line 956  DataLazy::resolveSample(ValueType& v,int
956      case LEZ:      case LEZ:
957      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));
958      break;      break;
959    // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
960        case NEZ:
961        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
962        break;
963        case EZ:
964        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
965        break;
966    
967        default:
968        throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
969      }
970      return &v;
971    }
972    
973    
974    
975    
976    
977    
978    /*
979      \brief Compute the value of the expression (unary operation) for the given sample.
980      \return Vector which stores the value of the subexpression for the given sample.
981      \param v A vector to store intermediate results.
982      \param offset Index in v to begin storing results.
983      \param sampleNo Sample number to evaluate.
984      \param roffset (output parameter) the offset in the return vector where the result begins.
985    
986      The return value will be an existing vector so do not deallocate it.
987      If the result is stored in v it should be stored at the offset given.
988      Everything from offset to the end of v should be considered available for this method to use.
989    */
990    DataTypes::ValueType*
991    DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
992    {
993        // we assume that any collapsing has been done before we get here
994        // since we only have one argument we don't need to think about only
995        // processing single points.
996      if (m_readytype!='E')
997      {
998        throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
999      }
1000        // since we can't write the result over the input, we need a result offset further along
1001      size_t subroffset=roffset+m_samplesize;
1002    LAZYDEBUG(cerr << "subroffset=" << subroffset << endl;)
1003      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1004      roffset=offset;
1005      size_t loop=0;
1006      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1007      size_t step=getNoValues();
1008      switch (m_op)
1009      {
1010        case SYM:
1011        for (loop=0;loop<numsteps;++loop)
1012        {
1013            DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
1014            subroffset+=step;
1015            offset+=step;
1016        }
1017        break;
1018        case NSYM:
1019        for (loop=0;loop<numsteps;++loop)
1020        {
1021            DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
1022            subroffset+=step;
1023            offset+=step;
1024        }
1025        break;
1026        default:
1027        throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1028      }
1029      return &v;
1030    }
1031    
1032    /*
1033      \brief Compute the value of the expression (unary operation) for the given sample.
1034      \return Vector which stores the value of the subexpression for the given sample.
1035      \param v A vector to store intermediate results.
1036      \param offset Index in v to begin storing results.
1037      \param sampleNo Sample number to evaluate.
1038      \param roffset (output parameter) the offset in the return vector where the result begins.
1039    
1040      The return value will be an existing vector so do not deallocate it.
1041      If the result is stored in v it should be stored at the offset given.
1042      Everything from offset to the end of v should be considered available for this method to use.
1043    */
1044    DataTypes::ValueType*
1045    DataLazy::resolveNP1OUT_P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1046    {
1047        // we assume that any collapsing has been done before we get here
1048        // since we only have one argument we don't need to think about only
1049        // processing single points.
1050      if (m_readytype!='E')
1051      {
1052        throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");
1053      }
1054        // since we can't write the result over the input, we need a result offset further along
1055      size_t subroffset;
1056      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1057    LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)
1058    LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)
1059      roffset=offset;
1060      size_t loop=0;
1061      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1062      size_t outstep=getNoValues();
1063      size_t instep=m_left->getNoValues();
1064    LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)
1065      switch (m_op)
1066      {
1067        case TRACE:
1068        for (loop=0;loop<numsteps;++loop)
1069        {
1070    size_t zz=sampleNo*getNumDPPSample()+loop;
1071    if (zz==5800)
1072    {
1073    LAZYDEBUG(cerr << "point=" <<  zz<< endl;)
1074    LAZYDEBUG(cerr << "Input to  trace=" << DataTypes::pointToString(*vleft,m_left->getShape(),subroffset,"") << endl;)
1075    LAZYDEBUG(cerr << "Offset for point=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << " vs ";)
1076    LAZYDEBUG(cerr << subroffset << endl;)
1077    LAZYDEBUG(cerr << "output=" << offset << endl;)
1078    }
1079                DataMaths::trace(*vleft,m_left->getShape(),subroffset, v ,getShape(),offset,m_axis_offset);
1080    if (zz==5800)
1081    {
1082    LAZYDEBUG(cerr << "Result of trace=" << DataTypes::pointToString(v,getShape(),offset,"") << endl;)
1083    }
1084            subroffset+=instep;
1085            offset+=outstep;
1086        }
1087        break;
1088        case TRANS:
1089        for (loop=0;loop<numsteps;++loop)
1090        {
1091                DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
1092            subroffset+=instep;
1093            offset+=outstep;
1094        }
1095        break;
1096        default:
1097        throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1098      }
1099      return &v;
1100    }
1101    
1102    
1103    #define PROC_OP(TYPE,X)                               \
1104        for (int j=0;j<onumsteps;++j)\
1105        {\
1106          for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1107          { \
1108    LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1109    LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1110             tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1111    LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \
1112             lroffset+=leftstep; \
1113             rroffset+=rightstep; \
1114          }\
1115          lroffset+=oleftstep;\
1116          rroffset+=orightstep;\
1117        }
1118    
1119    /*
1120      \brief Compute the value of the expression (binary operation) for the given sample.
1121      \return Vector which stores the value of the subexpression for the given sample.
1122      \param v A vector to store intermediate results.
1123      \param offset Index in v to begin storing results.
1124      \param sampleNo Sample number to evaluate.
1125      \param roffset (output parameter) the offset in the return vector where the result begins.
1126    
1127      The return value will be an existing vector so do not deallocate it.
1128      If the result is stored in v it should be stored at the offset given.
1129      Everything from offset to the end of v should be considered available for this method to use.
1130    */
1131    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1132    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1133    // If both children are expanded, then we can process them in a single operation (we treat
1134    // the whole sample as one big datapoint.
1135    // If one of the children is not expanded, then we need to treat each point in the sample
1136    // individually.
1137    // There is an additional complication when scalar operations are considered.
1138    // For example, 2+Vector.
1139    // In this case each double within the point is treated individually
1140    DataTypes::ValueType*
1141    DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
1142    {
1143    LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1144    
1145      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
1146        // first work out which of the children are expanded
1147      bool leftExp=(m_left->m_readytype=='E');
1148      bool rightExp=(m_right->m_readytype=='E');
1149      if (!leftExp && !rightExp)
1150      {
1151        throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1152      }
1153      bool leftScalar=(m_left->getRank()==0);
1154      bool rightScalar=(m_right->getRank()==0);
1155      if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1156      {
1157        throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1158      }
1159      size_t leftsize=m_left->getNoValues();
1160      size_t rightsize=m_right->getNoValues();
1161      size_t chunksize=1;           // how many doubles will be processed in one go
1162      int leftstep=0;       // how far should the left offset advance after each step
1163      int rightstep=0;
1164      int numsteps=0;       // total number of steps for the inner loop
1165      int oleftstep=0;  // the o variables refer to the outer loop
1166      int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1167      int onumsteps=1;
1168      
1169      bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1170      bool RES=(rightExp && rightScalar);
1171      bool LS=(!leftExp && leftScalar); // left is a single scalar
1172      bool RS=(!rightExp && rightScalar);
1173      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1174      bool RN=(!rightExp && !rightScalar);
1175      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1176      bool REN=(rightExp && !rightScalar);
1177    
1178      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1179      {
1180        chunksize=m_left->getNumDPPSample()*leftsize;
1181        leftstep=0;
1182        rightstep=0;
1183        numsteps=1;
1184      }
1185      else if (LES || RES)
1186      {
1187        chunksize=1;
1188        if (LES)        // left is an expanded scalar
1189        {
1190            if (RS)
1191            {
1192               leftstep=1;
1193               rightstep=0;
1194               numsteps=m_left->getNumDPPSample();
1195            }
1196            else        // RN or REN
1197            {
1198               leftstep=0;
1199               oleftstep=1;
1200               rightstep=1;
1201               orightstep=(RN ? -(int)rightsize : 0);
1202               numsteps=rightsize;
1203               onumsteps=m_left->getNumDPPSample();
1204            }
1205        }
1206        else        // right is an expanded scalar
1207        {
1208            if (LS)
1209            {
1210               rightstep=1;
1211               leftstep=0;
1212               numsteps=m_right->getNumDPPSample();
1213            }
1214            else
1215            {
1216               rightstep=0;
1217               orightstep=1;
1218               leftstep=1;
1219               oleftstep=(LN ? -(int)leftsize : 0);
1220               numsteps=leftsize;
1221               onumsteps=m_right->getNumDPPSample();
1222            }
1223        }
1224      }
1225      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1226      {
1227        if (LEN)    // and Right will be a single value
1228        {
1229            chunksize=rightsize;
1230            leftstep=rightsize;
1231            rightstep=0;
1232            numsteps=m_left->getNumDPPSample();
1233            if (RS)
1234            {
1235               numsteps*=leftsize;
1236            }
1237        }
1238        else    // REN
1239        {
1240            chunksize=leftsize;
1241            rightstep=leftsize;
1242            leftstep=0;
1243            numsteps=m_right->getNumDPPSample();
1244            if (LS)
1245            {
1246               numsteps*=rightsize;
1247            }
1248        }
1249      }
1250    
1251      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
1252        // Get the values of sub-expressions
1253      const ValueType* left=m_left->resolveSample(v,offset+getMaxSampleSize(),sampleNo,lroffset);   // see note on
1254        // calcBufss for why we can't put offset as the 2nd param above
1255      const ValueType* right=m_right->resolveSample(v,offset+2*getMaxSampleSize(),sampleNo,rroffset); // Note
1256        // the right child starts further along.
1257    LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1258    LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1259    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1260    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1261    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1262    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1263    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1264    
1265    
1266      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
1267      switch(m_op)
1268      {
1269        case ADD:
1270            PROC_OP(NO_ARG,plus<double>());
1271        break;
1272        case SUB:
1273        PROC_OP(NO_ARG,minus<double>());
1274        break;
1275        case MUL:
1276        PROC_OP(NO_ARG,multiplies<double>());
1277        break;
1278        case DIV:
1279        PROC_OP(NO_ARG,divides<double>());
1280        break;
1281        case POW:
1282           PROC_OP(double (double,double),::pow);
1283        break;
1284        default:
1285        throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1286      }
1287      roffset=offset;
1288      return &v;
1289    }
1290    
1291    
1292    
1293    /*
1294      \brief Compute the value of the expression (tensor product) for the given sample.
1295      \return Vector which stores the value of the subexpression for the given sample.
1296      \param v A vector to store intermediate results.
1297      \param offset Index in v to begin storing results.
1298      \param sampleNo Sample number to evaluate.
1299      \param roffset (output parameter) the offset in the return vector where the result begins.
1300    
1301      The return value will be an existing vector so do not deallocate it.
1302      If the result is stored in v it should be stored at the offset given.
1303      Everything from offset to the end of v should be considered available for this method to use.
1304    */
1305    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1306    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1307    // unlike the other resolve helpers, we must treat these datapoints separately.
1308    DataTypes::ValueType*
1309    DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
1310    {
1311    LAZYDEBUG(cout << "Resolve TensorProduct: " << toString()  << " to offset " << offset<< endl;)
1312    
1313      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
1314        // first work out which of the children are expanded
1315      bool leftExp=(m_left->m_readytype=='E');
1316      bool rightExp=(m_right->m_readytype=='E');
1317      int steps=getNumDPPSample();
1318    /*  int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
1319      int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);*/
1320      int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method
1321      int rightStep=(rightExp?m_right->getNoValues() : 0);
1322    
1323      int resultStep=getNoValues();
1324        // Get the values of sub-expressions (leave a gap of one sample for the result).
1325      int gap=offset+m_samplesize;  
1326    
1327    LAZYDEBUG(cout << "Query left with offset=" << gap << endl;)
1328    
1329      const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);
1330      gap+=m_left->getMaxSampleSize();
1331    
1332    
1333    LAZYDEBUG(cout << "Query right with offset=" << gap << endl;)
1334    
1335    
1336      const ValueType* right=m_right->resolveSample(v,gap,sampleNo,rroffset);
1337    
1338    LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
1339    cout << getNoValues() << endl;)
1340    LAZYDEBUG(cerr << "Result of left=";)
1341    LAZYDEBUG(cerr << "[" << lroffset << " .. " << lroffset+m_left->getNoValues() << "]" << endl;
1342    
1343    for (int i=lroffset, limit=lroffset+(leftExp?m_left->getNoValues()*m_left->getNumDPPSample():m_left->getNoValues());i<limit;++i)
1344    {
1345    cout << "[" << setw(2) << i-lroffset << "] " << setw(10) << (*left)[i] << " ";
1346    if (i%4==0) cout << endl;
1347    })
1348    LAZYDEBUG(cerr << "\nResult of right=" << endl;)
1349    LAZYDEBUG(
1350    for (int i=rroffset, limit=rroffset+(rightExp?m_right->getNoValues()*m_right->getNumDPPSample():m_right->getNoValues());i<limit;++i)
1351    {
1352    cerr << "[" <<  setw(2)<< i-rroffset << "] " << setw(10) << (*right)[i] << " ";
1353    if (i%4==0) cout << endl;
1354    }
1355    cerr << endl;
1356    )
1357    LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1358    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1359    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1360    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1361    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1362    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1363    LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1364    
1365      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
1366      switch(m_op)
1367      {
1368        case PROD:
1369        for (int i=0;i<steps;++i,resultp+=resultStep)
1370        {
1371    
1372    LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)
1373    LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)
1374    LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)
1375    
1376              const double *ptr_0 = &((*left)[lroffset]);
1377              const double *ptr_1 = &((*right)[rroffset]);
1378    
1379    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1380    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1381    
1382              matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
1383    
1384    LAZYDEBUG(cout << "Results--\n";
1385    {
1386      DataVector dv(getNoValues());
1387    for (int z=0;z<getNoValues();++z)
1388    {
1389      cout << "[" << setw(2) << z<< "] " << setw(10) << resultp[z] << " ";
1390      if (z%4==0) cout << endl;
1391      dv[z]=resultp[z];
1392    }
1393    cout << endl << DataTypes::pointToString(dv,getShape(),0,"RESLT");
1394    cout << "\nWritten to: " << resultp << " resultStep=" << resultStep << endl;
1395    }
1396    )
1397          lroffset+=leftStep;
1398          rroffset+=rightStep;
1399        }
1400        break;
1401      default:      default:
1402      throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");      throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1403      }
1404      roffset=offset;
1405      return &v;
1406    }
1407    
1408    
1409    
1410    /*
1411      \brief Compute the value of the expression for the given sample.
1412      \return Vector which stores the value of the subexpression for the given sample.
1413      \param v A vector to store intermediate results.
1414      \param offset Index in v to begin storing results.
1415      \param sampleNo Sample number to evaluate.
1416      \param roffset (output parameter) the offset in the return vector where the result begins.
1417    
1418      The return value will be an existing vector so do not deallocate it.
1419    */
1420    // the vector and the offset are a place where the method could write its data if it wishes
1421    // it is not obligated to do so. For example, if it has its own storage already, it can use that.
1422    // Hence the return value to indicate where the data is actually stored.
1423    // Regardless, the storage should be assumed to be used, even if it isn't.
1424    
1425    // the roffset is the offset within the returned vector where the data begins
1426    const DataTypes::ValueType*
1427    DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
1428    {
1429    LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1430        // collapse so we have a 'E' node or an IDENTITY for some other type
1431      if (m_readytype!='E' && m_op!=IDENTITY)
1432      {
1433        collapse();
1434      }
1435      if (m_op==IDENTITY)  
1436      {
1437        const ValueType& vec=m_id->getVectorRO();
1438        if (m_readytype=='C')
1439        {
1440        roffset=0;
1441    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1442        return &(vec);
1443      }      }
1444        roffset=m_id->getPointOffset(sampleNo, 0);
1445    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
1446        return &(vec);
1447      }
1448      if (m_readytype!='E')
1449      {
1450        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
1451    }    }
1452    return result;    switch (getOpgroup(m_op))
1453      {
1454      case G_UNARY:
1455      case G_UNARY_P: return resolveUnary(v, offset,sampleNo,roffset);
1456      case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
1457      case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
1458      case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);
1459      case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
1460      default:
1461        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
1462      }
1463    
1464    }
1465    
1466    const DataTypes::ValueType*
1467    DataLazy::resolveSample(BufferGroup& bg, int sampleNo, size_t& roffset)
1468    {
1469    #ifdef _OPENMP
1470        int tid=omp_get_thread_num();
1471    #else
1472        int tid=0;
1473    #endif
1474        return resolveSample(bg.getBuffer(tid),bg.getOffset(tid),sampleNo,roffset);
1475    }
1476    
1477    
1478    // This needs to do the work of the idenity constructor
1479    void
1480    DataLazy::resolveToIdentity()
1481    {
1482       if (m_op==IDENTITY)
1483        return;
1484       DataReady_ptr p=resolve();
1485       makeIdentity(p);
1486  }  }
1487    
1488    void DataLazy::makeIdentity(const DataReady_ptr& p)
1489    {
1490       m_op=IDENTITY;
1491       m_axis_offset=0;
1492       m_transpose=0;
1493       m_SL=m_SM=m_SR=0;
1494       m_children=m_height=0;
1495       m_id=p;
1496       if(p->isConstant()) {m_readytype='C';}
1497       else if(p->isExpanded()) {m_readytype='E';}
1498       else if (p->isTagged()) {m_readytype='T';}
1499       else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
1500       m_buffsRequired=1;
1501       m_samplesize=p->getNumDPPSample()*p->getNoValues();
1502       m_maxsamplesize=m_samplesize;
1503       m_left.reset();
1504       m_right.reset();
1505    }
1506    
1507    // To simplify the memory management, all threads operate on one large vector, rather than one each.
1508    // Each sample is evaluated independently and copied into the result DataExpanded.
1509  DataReady_ptr  DataReady_ptr
1510  DataLazy::resolve()  DataLazy::resolve()
1511  {  {
   // This is broken!     We need to have a buffer per thread!  
   // so the allocation of v needs to move inside the loop somehow  
   
 cout << "Sample size=" << m_samplesize << endl;  
 cout << "Buffers=" << m_buffsRequired << endl;  
1512    
1513    size_t threadbuffersize=m_samplesize*(max(1,m_buffsRequired)+1);  LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)
1514    LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)
1515      if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally
1516      {
1517        collapse();
1518      }
1519      if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
1520      {
1521        return m_id;
1522      }
1523        // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
1524      size_t threadbuffersize=m_maxsamplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough
1525        // storage to evaluate its expression
1526    int numthreads=1;    int numthreads=1;
1527  #ifdef _OPENMP  #ifdef _OPENMP
1528    numthreads=omp_get_max_threads();    numthreads=omp_get_max_threads();
   int threadnum=0;  
1529  #endif  #endif
1530    ValueType v(numthreads*threadbuffersize);    ValueType v(numthreads*threadbuffersize);
1531  cout << "Buffer created with size=" << v.size() << endl;  LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)
1532    ValueType dummy(getNoValues());    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
1533    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),dummy);    ValueType& resvec=result->getVectorRW();
   ValueType& resvec=result->getVector();  
1534    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1535    int sample;    int sample;
1536    int resoffset;    size_t outoffset;     // offset in the output data
1537    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1538    #pragma omp parallel for private(sample,resoffset,threadnum) schedule(static)    const ValueType* res=0;   // Vector storing the answer
1539      size_t resoffset=0;       // where in the vector to find the answer
1540    LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1541      #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
1542    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
1543    {    {
1544          if (sample==0)  {ENABLEDEBUG}
1545    
1546    //       if (sample==5800/getNumDPPSample())  {ENABLEDEBUG}
1547    LAZYDEBUG(cout << "################################# " << sample << endl;)
1548  #ifdef _OPENMP  #ifdef _OPENMP
1549      const double* res=resolveSample(v,sample,threadbuffersize*omp_get_thread_num());      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
1550  #else  #else
1551      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.
1552  #endif  #endif
1553      resoffset=result->getPointOffset(sample,0);  LAZYDEBUG(cerr << "-------------------------------- " << endl;)
1554      for (int i=0;i<m_samplesize;++i,++resoffset)    // copy values into the output vector  LAZYDEBUG(cerr<< "Copying sample#" << sample << endl;)
1555        outoffset=result->getPointOffset(sample,0);
1556    LAZYDEBUG(cerr << "offset=" << outoffset << " from offset=" << resoffset << " " << m_samplesize << " doubles" << endl;)
1557        for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector
1558      {      {
1559      resvec[resoffset]=res[i];  LAZYDEBUG(cerr << "outoffset=" << outoffset << " resoffset=" << resoffset << " " << (*res)[resoffset]<< endl;)
1560        resvec[outoffset]=(*res)[resoffset];
1561      }      }
1562    LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)
1563    LAZYDEBUG(cerr << "*********************************" << endl;)
1564        DISABLEDEBUG
1565    }    }
1566    return resptr;    return resptr;
1567  }  }
# Line 435  DataLazy::toString() const Line 1575  DataLazy::toString() const
1575    return oss.str();    return oss.str();
1576  }  }
1577    
1578    
1579  void  void
1580  DataLazy::intoString(ostringstream& oss) const  DataLazy::intoString(ostringstream& oss) const
1581  {  {
1582    //    oss << "[" << m_children <<";"<<m_height <<"]";
1583    switch (getOpgroup(m_op))    switch (getOpgroup(m_op))
1584    {    {
1585    case G_IDENTITY:    case G_IDENTITY:
1586        if (m_id->isExpanded())
1587        {
1588           oss << "E";
1589        }
1590        else if (m_id->isTagged())
1591        {
1592          oss << "T";
1593        }
1594        else if (m_id->isConstant())
1595        {
1596          oss << "C";
1597        }
1598        else
1599        {
1600          oss << "?";
1601        }
1602      oss << '@' << m_id.get();      oss << '@' << m_id.get();
1603      break;      break;
1604    case G_BINARY:    case G_BINARY:
# Line 451  DataLazy::intoString(ostringstream& oss) Line 1609  DataLazy::intoString(ostringstream& oss)
1609      oss << ')';      oss << ')';
1610      break;      break;
1611    case G_UNARY:    case G_UNARY:
1612      case G_UNARY_P:
1613      case G_NP1OUT:
1614      case G_NP1OUT_P:
1615      oss << opToString(m_op) << '(';      oss << opToString(m_op) << '(';
1616      m_left->intoString(oss);      m_left->intoString(oss);
1617      oss << ')';      oss << ')';
1618      break;      break;
1619      case G_TENSORPROD:
1620        oss << opToString(m_op) << '(';
1621        m_left->intoString(oss);
1622        oss << ", ";
1623        m_right->intoString(oss);
1624        oss << ')';
1625        break;
1626    default:    default:
1627      oss << "UNKNOWN";      oss << "UNKNOWN";
1628    }    }
1629  }  }
1630    
 // 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.  
1631  DataAbstract*  DataAbstract*
1632  DataLazy::deepCopy()  DataLazy::deepCopy()
1633  {  {
1634    if (m_op==IDENTITY)    switch (getOpgroup(m_op))
1635    {    {
1636      return new DataLazy(m_left);    // we don't need to copy the child here    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
1637      case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
1638      case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
1639      case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
1640      case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
1641      default:
1642        throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
1643    }    }
   return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);  
1644  }  }
1645    
1646    
1647    // There is no single, natural interpretation of getLength on DataLazy.
1648    // Instances of DataReady can look at the size of their vectors.
1649    // For lazy though, it could be the size the data would be if it were resolved;
1650    // or it could be some function of the lengths of the DataReady instances which
1651    // form part of the expression.
1652    // Rather than have people making assumptions, I have disabled the method.
1653  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1654  DataLazy::getLength() const  DataLazy::getLength() const
1655  {  {
1656    return m_length;    throw DataException("getLength() does not make sense for lazy data.");
1657  }  }
1658    
1659    
# Line 486  DataLazy::getSlice(const DataTypes::Regi Line 1663  DataLazy::getSlice(const DataTypes::Regi
1663    throw DataException("getSlice - not implemented for Lazy objects.");    throw DataException("getSlice - not implemented for Lazy objects.");
1664  }  }
1665    
1666    
1667    // To do this we need to rely on our child nodes
1668    DataTypes::ValueType::size_type
1669    DataLazy::getPointOffset(int sampleNo,
1670                     int dataPointNo)
1671    {
1672      if (m_op==IDENTITY)
1673      {
1674        return m_id->getPointOffset(sampleNo,dataPointNo);
1675      }
1676      if (m_readytype!='E')
1677      {
1678        collapse();
1679        return m_id->getPointOffset(sampleNo,dataPointNo);
1680      }
1681      // at this point we do not have an identity node and the expression will be Expanded
1682      // so we only need to know which child to ask
1683      if (m_left->m_readytype=='E')
1684      {
1685        return m_left->getPointOffset(sampleNo,dataPointNo);
1686      }
1687      else
1688      {
1689        return m_right->getPointOffset(sampleNo,dataPointNo);
1690      }
1691    }
1692    
1693    // To do this we need to rely on our child nodes
1694  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
1695  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
1696                   int dataPointNo) const                   int dataPointNo) const
1697  {  {
1698    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");    if (m_op==IDENTITY)
1699      {
1700        return m_id->getPointOffset(sampleNo,dataPointNo);
1701      }
1702      if (m_readytype=='E')
1703      {
1704        // at this point we do not have an identity node and the expression will be Expanded
1705        // so we only need to know which child to ask
1706        if (m_left->m_readytype=='E')
1707        {
1708        return m_left->getPointOffset(sampleNo,dataPointNo);
1709        }
1710        else
1711        {
1712        return m_right->getPointOffset(sampleNo,dataPointNo);
1713        }
1714      }
1715      if (m_readytype=='C')
1716      {
1717        return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
1718      }
1719      throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
1720    }
1721    
1722    
1723    // I have decided to let Data:: handle this issue.
1724    void
1725    DataLazy::setToZero()
1726    {
1727    //   DataTypes::ValueType v(getNoValues(),0);
1728    //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
1729    //   m_op=IDENTITY;
1730    //   m_right.reset();  
1731    //   m_left.reset();
1732    //   m_readytype='C';
1733    //   m_buffsRequired=1;
1734    
1735      throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
1736    
1737    }
1738    
1739    bool
1740    DataLazy::actsExpanded() const
1741    {
1742        return (m_readytype=='E');
1743  }  }
1744    
1745  }   // end namespace  }   // end namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26