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

Legend:
Removed from v.1884  
changed lines
  Added in v.2147

  ViewVC Help
Powered by ViewVC 1.1.26