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

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

  ViewVC Help
Powered by ViewVC 1.1.26