/[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 1865 by jfenwick, Thu Oct 9 03:53:57 2008 UTC trunk/escript/src/DataLazy.cpp revision 2514 by jfenwick, Fri Jul 3 00:57:45 2009 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"
28    #include "UnaryFuncs.h"     // for escript::fsign
29    #include "Utils.h"
30    
31    #include "EscriptParams.h"
32    
33    #include <iomanip>      // for some fancy formatting in debug
34    
35    // #define LAZYDEBUG(X) if (privdebug){X;}
36    #define LAZYDEBUG(X)
37    namespace
38    {
39    bool privdebug=false;
40    
41    #define ENABLEDEBUG privdebug=true;
42    #define DISABLEDEBUG privdebug=false;
43    }
44    
45    // #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {cerr << "\n!!!!!!! SIZE LIMIT EXCEEDED " << m_children << ";" << m_height << endl << toString() << endl;resolveToIdentity();}
46    
47    #define SIZELIMIT if ((m_height>escript::escriptParams.getTOO_MANY_LEVELS()) || (m_children>escript::escriptParams.getTOO_MANY_NODES())) {resolveToIdentity();}
48    
49    
50    /*
51    How does DataLazy work?
52    ~~~~~~~~~~~~~~~~~~~~~~~
53    
54    Each instance represents a single operation on one or two other DataLazy instances. These arguments are normally
55    denoted left and right.
56    
57    A special operation, IDENTITY, stores an instance of DataReady in the m_id member.
58    This means that all "internal" nodes in the structure are instances of DataLazy.
59    
60    Each operation has a string representation as well as an opgroup - eg G_IDENTITY, G_BINARY, ...
61    Note that IDENITY is not considered a unary operation.
62    
63    I am avoiding calling the structure formed a tree because it is not guaranteed to be one (eg c=a+a).
64    It must however form a DAG (directed acyclic graph).
65    I will refer to individual DataLazy objects with the structure as nodes.
66    
67    Each node also stores:
68    - m_readytype \in {'E','T','C','?'} ~ indicates what sort of DataReady would be produced if the expression was
69        evaluated.
70    - m_buffsrequired ~ the larged number of samples which would need to be kept simultaneously in order to
71        evaluate the expression.
72    - m_samplesize ~ the number of doubles stored in a sample.
73    
74    When a new node is created, the above values are computed based on the values in the child nodes.
75    Eg: if left requires 4 samples and right requires 6 then left+right requires 7 samples.
76    
77    The resolve method, which produces a DataReady from a DataLazy, does the following:
78    1) Create a DataReady to hold the new result.
79    2) Allocate a vector (v) big enough to hold m_buffsrequired samples.
80    3) For each sample, call resolveSample with v, to get its values and copy them into the result object.
81    
82    (In the case of OMP, multiple samples are resolved in parallel so the vector needs to be larger.)
83    
84    resolveSample returns a Vector* and an offset within that vector where the result is stored.
85    Normally, this would be v, but for identity nodes their internal vector is returned instead.
86    
87    The convention that I use, is that the resolve methods should store their results starting at the offset they are passed.
88    
89    For expressions which evaluate to Constant or Tagged, there is a different evaluation method.
90    The collapse method invokes the (non-lazy) operations on the Data class to evaluate the expression.
91    
92    To add a new operator you need to do the following (plus anything I might have forgotten - adding a new group for example):
93    1) Add to the ES_optype.
94    2) determine what opgroup your operation belongs to (X)
95    3) add a string for the op to the end of ES_opstrings
96    4) increase ES_opcount
97    5) add an entry (X) to opgroups
98    6) add an entry to the switch in collapseToReady
99    7) add an entry to resolveX
100    */
101    
102    
103  using namespace std;  using namespace std;
104    using namespace boost;
105    
106  namespace escript  namespace escript
107  {  {
108    
 const std::string&  
 opToString(ES_optype op);  
   
109  namespace  namespace
110  {  {
111    
112    enum ES_opgroup
113    {
114       G_UNKNOWN,
115       G_IDENTITY,
116       G_BINARY,        // pointwise operations with two arguments
117       G_UNARY,     // pointwise operations with one argument
118       G_UNARY_P,       // pointwise operations with one argument, requiring a parameter
119       G_NP1OUT,        // non-pointwise op with one output
120       G_NP1OUT_P,      // non-pointwise op with one output requiring a parameter
121       G_TENSORPROD,    // general tensor product
122       G_NP1OUT_2P      // non-pointwise op with one output requiring two params
123    };
124    
125    
126    
127    
128    string ES_opstrings[]={"UNKNOWN","IDENTITY","+","-","*","/","^",
129                "sin","cos","tan",
130                "asin","acos","atan","sinh","cosh","tanh","erf",
131                "asinh","acosh","atanh",
132                "log10","log","sign","abs","neg","pos","exp","sqrt",
133                "1/","where>0","where<0","where>=0","where<=0", "where<>0","where=0",
134                "symmetric","nonsymmetric",
135                "prod",
136                "transpose", "trace",
137                "swapaxes"};
138    int ES_opcount=41;
139    ES_opgroup opgroups[]={G_UNKNOWN,G_IDENTITY,G_BINARY,G_BINARY,G_BINARY,G_BINARY, G_BINARY,
140                G_UNARY,G_UNARY,G_UNARY, //10
141                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 17
142                G_UNARY,G_UNARY,G_UNARY,                    // 20
143                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY,    // 28
144                G_UNARY,G_UNARY,G_UNARY,G_UNARY,G_UNARY, G_UNARY_P, G_UNARY_P,      // 35
145                G_NP1OUT,G_NP1OUT,
146                G_TENSORPROD,
147                G_NP1OUT_P, G_NP1OUT_P,
148                G_NP1OUT_2P};
149    inline
150    ES_opgroup
151    getOpgroup(ES_optype op)
152    {
153      return opgroups[op];
154    }
155    
156  // return the FunctionSpace of the result of "left op right"  // return the FunctionSpace of the result of "left op right"
157  FunctionSpace  FunctionSpace
158  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultFS(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
# Line 41  resultFS(DataAbstract_ptr left, DataAbst Line 161  resultFS(DataAbstract_ptr left, DataAbst
161      // maybe we need an interpolate node -      // maybe we need an interpolate node -
162      // 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
163      // programming error exception.      // programming error exception.
164      return left->getFunctionSpace();  
165      FunctionSpace l=left->getFunctionSpace();
166      FunctionSpace r=right->getFunctionSpace();
167      if (l!=r)
168      {
169        if (r.probeInterpolation(l))
170        {
171        return l;
172        }
173        if (l.probeInterpolation(r))
174        {
175        return r;
176        }
177        throw DataException("Cannot interpolate between the FunctionSpaces given for operation "+opToString(op)+".");
178      }
179      return l;
180  }  }
181    
182  // return the shape of the result of "left op right"  // return the shape of the result of "left op right"
183    // the shapes resulting from tensor product are more complex to compute so are worked out elsewhere
184  DataTypes::ShapeType  DataTypes::ShapeType
185  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  resultShape(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
186  {  {
187      return DataTypes::scalarShape;      if (left->getShape()!=right->getShape())
188        {
189          if ((getOpgroup(op)!=G_BINARY) && (getOpgroup(op)!=G_NP1OUT))
190          {
191            throw DataException("Shapes not the name - shapes must match for (point)binary operations.");
192          }
193          if (left->getRank()==0)   // we need to allow scalar * anything
194          {
195            return right->getShape();
196          }
197          if (right->getRank()==0)
198          {
199            return left->getShape();
200          }
201          throw DataException("Shapes not the same - arguments must have matching shapes (or be scalars) for (point)binary operations on lazy data.");
202        }
203        return left->getShape();
204  }  }
205    
206  size_t  // return the shape for "op left"
207  resultLength(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  
208    DataTypes::ShapeType
209    resultShape(DataAbstract_ptr left, ES_optype op, int axis_offset)
210    {
211        switch(op)
212        {
213            case TRANS:
214           {            // for the scoping of variables
215            const DataTypes::ShapeType& s=left->getShape();
216            DataTypes::ShapeType sh;
217            int rank=left->getRank();
218            if (axis_offset<0 || axis_offset>rank)
219            {
220                   throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
221                }
222                for (int i=0; i<rank; i++)
223            {
224               int index = (axis_offset+i)%rank;
225                   sh.push_back(s[index]); // Append to new shape
226                }
227            return sh;
228           }
229        break;
230        case TRACE:
231           {
232            int rank=left->getRank();
233            if (rank<2)
234            {
235               throw DataException("Trace can only be computed for objects with rank 2 or greater.");
236            }
237            if ((axis_offset>rank-2) || (axis_offset<0))
238            {
239               throw DataException("Trace: axis offset must lie between 0 and rank-2 inclusive.");
240            }
241            if (rank==2)
242            {
243               return DataTypes::scalarShape;
244            }
245            else if (rank==3)
246            {
247               DataTypes::ShapeType sh;
248                   if (axis_offset==0)
249               {
250                    sh.push_back(left->getShape()[2]);
251                   }
252                   else     // offset==1
253               {
254                sh.push_back(left->getShape()[0]);
255                   }
256               return sh;
257            }
258            else if (rank==4)
259            {
260               DataTypes::ShapeType sh;
261               const DataTypes::ShapeType& s=left->getShape();
262                   if (axis_offset==0)
263               {
264                    sh.push_back(s[2]);
265                    sh.push_back(s[3]);
266                   }
267                   else if (axis_offset==1)
268               {
269                    sh.push_back(s[0]);
270                    sh.push_back(s[3]);
271                   }
272               else     // offset==2
273               {
274                sh.push_back(s[0]);
275                sh.push_back(s[1]);
276               }
277               return sh;
278            }
279            else        // unknown rank
280            {
281               throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
282            }
283           }
284        break;
285            default:
286        throw DataException("Programmer error - resultShape(left,op) can't compute shapes for operator "+opToString(op)+".");
287        }
288    }
289    
290    DataTypes::ShapeType
291    SwapShape(DataAbstract_ptr left, const int axis0, const int axis1)
292  {  {
293     switch(op)       // This code taken from the Data.cpp swapaxes() method
294         // Some of the checks are probably redundant here
295         int axis0_tmp,axis1_tmp;
296         const DataTypes::ShapeType& s=left->getShape();
297         DataTypes::ShapeType out_shape;
298         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
299         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
300         int rank=left->getRank();
301         if (rank<2) {
302            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
303         }
304         if (axis0<0 || axis0>rank-1) {
305            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
306         }
307         if (axis1<0 || axis1>rank-1) {
308             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
309         }
310         if (axis0 == axis1) {
311             throw DataException("Error - Data::swapaxes: axis indices must be different.");
312         }
313         if (axis0 > axis1) {
314             axis0_tmp=axis1;
315             axis1_tmp=axis0;
316         } else {
317             axis0_tmp=axis0;
318             axis1_tmp=axis1;
319         }
320         for (int i=0; i<rank; i++) {
321           if (i == axis0_tmp) {
322              out_shape.push_back(s[axis1_tmp]);
323           } else if (i == axis1_tmp) {
324              out_shape.push_back(s[axis0_tmp]);
325           } else {
326              out_shape.push_back(s[i]);
327           }
328         }
329        return out_shape;
330    }
331    
332    
333    // determine the output shape for the general tensor product operation
334    // the additional parameters return information required later for the product
335    // the majority of this code is copy pasted from C_General_Tensor_Product
336    DataTypes::ShapeType
337    GTPShape(DataAbstract_ptr left, DataAbstract_ptr right, int axis_offset, int transpose, int& SL, int& SM, int& SR)
338    {
339        
340      // Get rank and shape of inputs
341      int rank0 = left->getRank();
342      int rank1 = right->getRank();
343      const DataTypes::ShapeType& shape0 = left->getShape();
344      const DataTypes::ShapeType& shape1 = right->getShape();
345    
346      // Prepare for the loops of the product and verify compatibility of shapes
347      int start0=0, start1=0;
348      if (transpose == 0)       {}
349      else if (transpose == 1)  { start0 = axis_offset; }
350      else if (transpose == 2)  { start1 = rank1-axis_offset; }
351      else              { throw DataException("DataLazy GeneralTensorProduct Constructor: Error - transpose should be 0, 1 or 2"); }
352    
353      if (rank0<axis_offset)
354      {
355        throw DataException("DataLazy GeneralTensorProduct Constructor: Error - rank of left < axisoffset");
356      }
357    
358      // Adjust the shapes for transpose
359      DataTypes::ShapeType tmpShape0(rank0);    // pre-sizing the vectors rather
360      DataTypes::ShapeType tmpShape1(rank1);    // than using push_back
361      for (int i=0; i<rank0; i++)   { tmpShape0[i]=shape0[(i+start0)%rank0]; }
362      for (int i=0; i<rank1; i++)   { tmpShape1[i]=shape1[(i+start1)%rank1]; }
363    
364      // Prepare for the loops of the product
365      SL=1, SM=1, SR=1;
366      for (int i=0; i<rank0-axis_offset; i++)   {
367        SL *= tmpShape0[i];
368      }
369      for (int i=rank0-axis_offset; i<rank0; i++)   {
370        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
371          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
372        }
373        SM *= tmpShape0[i];
374      }
375      for (int i=axis_offset; i<rank1; i++)     {
376        SR *= tmpShape1[i];
377      }
378    
379      // Define the shape of the output (rank of shape is the sum of the loop ranges below)
380      DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);  
381      {         // block to limit the scope of out_index
382         int out_index=0;
383         for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
384         for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
385      }
386    
387      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
388      {
389         ostringstream os;
390         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
391         throw DataException(os.str());
392      }
393    
394      return shape2;
395    }
396    
397    // determine the number of samples requires to evaluate an expression combining left and right
398    // NP1OUT needs an extra buffer because we can't write the answers over the top of the input.
399    // The same goes for G_TENSORPROD
400    // It might seem that pointwise binary ops (G_BINARY) could be written over the top of the lefts.
401    // This would be true were it not for the possibility that the LHS could be a scalar which needs to be examined
402    // multiple times
403    int
404    calcBuffs(const DataLazy_ptr& left, const DataLazy_ptr& right, ES_optype op)
405    {
406       switch(getOpgroup(op))
407     {     {
408     case IDENTITY: return left->getLength();     case G_IDENTITY: return 1;
409       case G_BINARY: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
410       case G_UNARY:
411       case G_UNARY_P: return max(left->getBuffsRequired(),1);
412       case G_NP1OUT: return 1+max(left->getBuffsRequired(),1);
413       case G_NP1OUT_P: return 1+max(left->getBuffsRequired(),1);
414       case G_TENSORPROD: return 1+max(left->getBuffsRequired(),right->getBuffsRequired()+1);
415       case G_NP1OUT_2P: return 1+max(left->getBuffsRequired(),1);
416     default:     default:
417      throw DataException("Programmer Error - attempt to getLength() for operator "+opToString(op)+".");      throw DataException("Programmer Error - attempt to calcBuffs() for operator "+opToString(op)+".");
   
418     }     }
419  }  }
420    
 string ES_opstrings[]={"UNKNOWN","IDENTITY"};  
 int ES_opcount=2;  
421    
422  }   // end anonymous namespace  }   // end anonymous namespace
423    
424    
425    
426    // Return a string representing the operation
427  const std::string&  const std::string&
428  opToString(ES_optype op)  opToString(ES_optype op)
429  {  {
# Line 79  opToString(ES_optype op) Line 434  opToString(ES_optype op)
434    return ES_opstrings[op];    return ES_opstrings[op];
435  }  }
436    
437    #ifdef LAZY_NODE_STORAGE
438    void DataLazy::LazyNodeSetup()
439    {
440    #ifdef _OPENMP
441        int numthreads=omp_get_max_threads();
442        m_samples.resize(numthreads*m_samplesize);
443        m_sampleids=new int[numthreads];
444        for (int i=0;i<numthreads;++i)
445        {
446            m_sampleids[i]=-1;  
447        }
448    #else
449        m_samples.resize(m_samplesize);
450        m_sampleids=new int[1];
451        m_sampleids[0]=-1;
452    #endif  // _OPENMP
453    }
454    #endif   // LAZY_NODE_STORAGE
455    
456    
457    // Creates an identity node
458  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
459      : parent(p->getFunctionSpace(),p->getShape()),      : parent(p->getFunctionSpace(),p->getShape())
460      m_left(p),  #ifdef LAZY_NODE_STORAGE
461      m_op(IDENTITY)      ,m_sampleids(0),
462        m_samples(1)
463    #endif
464  {  {
465     length=resultLength(m_left,m_right,m_op);     if (p->isLazy())
466       {
467        // I don't want identity of Lazy.
468        // Question: Why would that be so bad?
469        // Answer: We assume that the child of ID is something we can call getVector on
470        throw DataException("Programmer error - attempt to create identity from a DataLazy.");
471       }
472       else
473       {
474        p->makeLazyShared();
475        DataReady_ptr dr=dynamic_pointer_cast<DataReady>(p);
476        makeIdentity(dr);
477    LAZYDEBUG(cout << "Wrapping " << dr.get() << " id=" << m_id.get() << endl;)
478       }
479    LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
480    }
481    
482    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
483        : parent(left->getFunctionSpace(),left->getShape()),
484        m_op(op),
485        m_axis_offset(0),
486        m_transpose(0),
487        m_SL(0), m_SM(0), m_SR(0)
488    {
489       if ((getOpgroup(op)!=G_UNARY) && (getOpgroup(op)!=G_NP1OUT))
490       {
491        throw DataException("Programmer error - constructor DataLazy(left, op) will only process UNARY operations.");
492       }
493    
494       DataLazy_ptr lleft;
495       if (!left->isLazy())
496       {
497        lleft=DataLazy_ptr(new DataLazy(left));
498       }
499       else
500       {
501        lleft=dynamic_pointer_cast<DataLazy>(left);
502       }
503       m_readytype=lleft->m_readytype;
504       m_left=lleft;
505       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
506       m_samplesize=getNumDPPSample()*getNoValues();
507       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
508       m_children=m_left->m_children+1;
509       m_height=m_left->m_height+1;
510    #ifdef LAZY_NODE_STORAGE
511       LazyNodeSetup();
512    #endif
513       SIZELIMIT
514  }  }
515    
516    
517    // In this constructor we need to consider interpolation
518  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op)
519      : parent(resultFS(left,right,op), resultShape(left,right,op)),      : parent(resultFS(left,right,op), resultShape(left,right,op)),
520      m_left(left),      m_op(op),
521      m_right(right),      m_SL(0), m_SM(0), m_SR(0)
522      m_op(op)  {
523    LAZYDEBUG(cout << "Forming operator with " << left.get() << " " << right.get() << endl;)
524       if ((getOpgroup(op)!=G_BINARY))
525       {
526        throw DataException("Programmer error - constructor DataLazy(left, right, op) will only process BINARY operations.");
527       }
528    
529       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
530       {
531        FunctionSpace fs=getFunctionSpace();
532        Data ltemp(left);
533        Data tmp(ltemp,fs);
534        left=tmp.borrowDataPtr();
535       }
536       if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
537       {
538        Data tmp(Data(right),getFunctionSpace());
539        right=tmp.borrowDataPtr();
540    LAZYDEBUG(cout << "Right interpolation required " << right.get() << endl;)
541       }
542       left->operandCheck(*right);
543    
544       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
545       {
546        m_left=dynamic_pointer_cast<DataLazy>(left);
547    LAZYDEBUG(cout << "Left is " << m_left->toString() << endl;)
548       }
549       else
550       {
551        m_left=DataLazy_ptr(new DataLazy(left));
552    LAZYDEBUG(cout << "Left " << left.get() << " wrapped " << m_left->m_id.get() << endl;)
553       }
554       if (right->isLazy())
555       {
556        m_right=dynamic_pointer_cast<DataLazy>(right);
557    LAZYDEBUG(cout << "Right is " << m_right->toString() << endl;)
558       }
559       else
560       {
561        m_right=DataLazy_ptr(new DataLazy(right));
562    LAZYDEBUG(cout << "Right " << right.get() << " wrapped " << m_right->m_id.get() << endl;)
563       }
564       char lt=m_left->m_readytype;
565       char rt=m_right->m_readytype;
566       if (lt=='E' || rt=='E')
567       {
568        m_readytype='E';
569       }
570       else if (lt=='T' || rt=='T')
571       {
572        m_readytype='T';
573       }
574       else
575       {
576        m_readytype='C';
577       }
578       m_samplesize=getNumDPPSample()*getNoValues();
579       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
580       m_buffsRequired=calcBuffs(m_left, m_right,m_op);
581       m_children=m_left->m_children+m_right->m_children+2;
582       m_height=max(m_left->m_height,m_right->m_height)+1;
583    #ifdef LAZY_NODE_STORAGE
584       LazyNodeSetup();
585    #endif
586       SIZELIMIT
587    LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
588    }
589    
590    DataLazy::DataLazy(DataAbstract_ptr left, DataAbstract_ptr right, ES_optype op, int axis_offset, int transpose)
591        : parent(resultFS(left,right,op), GTPShape(left,right, axis_offset, transpose, m_SL,m_SM, m_SR)),
592        m_op(op),
593        m_axis_offset(axis_offset),
594        m_transpose(transpose)
595    {
596       if ((getOpgroup(op)!=G_TENSORPROD))
597       {
598        throw DataException("Programmer error - constructor DataLazy(left, right, op, ax, tr) will only process BINARY operations which require parameters.");
599       }
600       if ((transpose>2) || (transpose<0))
601       {
602        throw DataException("DataLazy GeneralTensorProduct constructor: Error - transpose should be 0, 1 or 2");
603       }
604       if (getFunctionSpace()!=left->getFunctionSpace())    // left needs to be interpolated
605       {
606        FunctionSpace fs=getFunctionSpace();
607        Data ltemp(left);
608        Data tmp(ltemp,fs);
609        left=tmp.borrowDataPtr();
610       }
611       if (getFunctionSpace()!=right->getFunctionSpace())   // right needs to be interpolated
612       {
613        Data tmp(Data(right),getFunctionSpace());
614        right=tmp.borrowDataPtr();
615       }
616    //    left->operandCheck(*right);
617    
618       if (left->isLazy())          // the children need to be DataLazy. Wrap them in IDENTITY if required
619       {
620        m_left=dynamic_pointer_cast<DataLazy>(left);
621       }
622       else
623       {
624        m_left=DataLazy_ptr(new DataLazy(left));
625       }
626       if (right->isLazy())
627       {
628        m_right=dynamic_pointer_cast<DataLazy>(right);
629       }
630       else
631       {
632        m_right=DataLazy_ptr(new DataLazy(right));
633       }
634       char lt=m_left->m_readytype;
635       char rt=m_right->m_readytype;
636       if (lt=='E' || rt=='E')
637       {
638        m_readytype='E';
639       }
640       else if (lt=='T' || rt=='T')
641       {
642        m_readytype='T';
643       }
644       else
645       {
646        m_readytype='C';
647       }
648       m_samplesize=getNumDPPSample()*getNoValues();
649       m_maxsamplesize=max(max(m_samplesize,m_right->getMaxSampleSize()),m_left->getMaxSampleSize());  
650       m_buffsRequired=calcBuffs(m_left, m_right,m_op);
651       m_children=m_left->m_children+m_right->m_children+2;
652       m_height=max(m_left->m_height,m_right->m_height)+1;
653    #ifdef LAZY_NODE_STORAGE
654       LazyNodeSetup();
655    #endif
656       SIZELIMIT
657    LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
658    }
659    
660    
661    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, int axis_offset)
662        : parent(left->getFunctionSpace(), resultShape(left,op, axis_offset)),
663        m_op(op),
664        m_axis_offset(axis_offset),
665        m_transpose(0),
666        m_tol(0)
667  {  {
668     length=resultLength(m_left,m_right,m_op);     if ((getOpgroup(op)!=G_NP1OUT_P))
669       {
670        throw DataException("Programmer error - constructor DataLazy(left, op, ax) will only process UNARY operations which require parameters.");
671       }
672       DataLazy_ptr lleft;
673       if (!left->isLazy())
674       {
675        lleft=DataLazy_ptr(new DataLazy(left));
676       }
677       else
678       {
679        lleft=dynamic_pointer_cast<DataLazy>(left);
680       }
681       m_readytype=lleft->m_readytype;
682       m_left=lleft;
683       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
684       m_samplesize=getNumDPPSample()*getNoValues();
685       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
686       m_children=m_left->m_children+1;
687       m_height=m_left->m_height+1;
688    #ifdef LAZY_NODE_STORAGE
689       LazyNodeSetup();
690    #endif
691       SIZELIMIT
692    LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
693    }
694    
695    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, double tol)
696        : parent(left->getFunctionSpace(), left->getShape()),
697        m_op(op),
698        m_axis_offset(0),
699        m_transpose(0),
700        m_tol(tol)
701    {
702       if ((getOpgroup(op)!=G_UNARY_P))
703       {
704        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require parameters.");
705       }
706       DataLazy_ptr lleft;
707       if (!left->isLazy())
708       {
709        lleft=DataLazy_ptr(new DataLazy(left));
710       }
711       else
712       {
713        lleft=dynamic_pointer_cast<DataLazy>(left);
714       }
715       m_readytype=lleft->m_readytype;
716       m_left=lleft;
717       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
718       m_samplesize=getNumDPPSample()*getNoValues();
719       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
720       m_children=m_left->m_children+1;
721       m_height=m_left->m_height+1;
722    #ifdef LAZY_NODE_STORAGE
723       LazyNodeSetup();
724    #endif
725       SIZELIMIT
726    LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
727  }  }
728    
729    
730    DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op, const int axis0, const int axis1)
731        : parent(left->getFunctionSpace(), SwapShape(left,axis0,axis1)),
732        m_op(op),
733        m_axis_offset(axis0),
734        m_transpose(axis1),
735        m_tol(0)
736    {
737       if ((getOpgroup(op)!=G_NP1OUT_2P))
738       {
739        throw DataException("Programmer error - constructor DataLazy(left, op, tol) will only process UNARY operations which require two integer parameters.");
740       }
741       DataLazy_ptr lleft;
742       if (!left->isLazy())
743       {
744        lleft=DataLazy_ptr(new DataLazy(left));
745       }
746       else
747       {
748        lleft=dynamic_pointer_cast<DataLazy>(left);
749       }
750       m_readytype=lleft->m_readytype;
751       m_left=lleft;
752       m_buffsRequired=calcBuffs(m_left, m_right,m_op); // yeah m_right will be null at this point
753       m_samplesize=getNumDPPSample()*getNoValues();
754       m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
755       m_children=m_left->m_children+1;
756       m_height=m_left->m_height+1;
757    #ifdef LAZY_NODE_STORAGE
758       LazyNodeSetup();
759    #endif
760       SIZELIMIT
761    LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)
762    }
763    
764  DataLazy::~DataLazy()  DataLazy::~DataLazy()
765  {  {
766    #ifdef LAZY_NODE_SETUP
767       delete[] m_sampleids;
768       delete[] m_samples;
769    #endif
770  }  }
771    
772  // If resolving records a pointer to the resolved Data we may need to rethink the const on this method  
773    int
774    DataLazy::getBuffsRequired() const
775    {
776        return m_buffsRequired;
777    }
778    
779    
780    size_t
781    DataLazy::getMaxSampleSize() const
782    {
783        return m_maxsamplesize;
784    }
785    
786    
787    
788    size_t
789    DataLazy::getSampleBufferSize() const
790    {
791        return m_maxsamplesize*(max(1,m_buffsRequired));
792    }
793    
794    /*
795      \brief Evaluates the expression using methods on Data.
796      This does the work for the collapse method.
797      For reasons of efficiency do not call this method on DataExpanded nodes.
798    */
799  DataReady_ptr  DataReady_ptr
800  DataLazy::resolve()  DataLazy::collapseToReady()
801  {  {
802      if (m_readytype=='E')
803      { // this is more an efficiency concern than anything else
804        throw DataException("Programmer Error - do not use collapse on Expanded data.");
805      }
806    if (m_op==IDENTITY)    if (m_op==IDENTITY)
807    {    {
808      return m_left->resolve();      return m_id;
809    }    }
810    else    DataReady_ptr pleft=m_left->collapseToReady();
811      Data left(pleft);
812      Data right;
813      if ((getOpgroup(m_op)==G_BINARY) || (getOpgroup(m_op)==G_TENSORPROD))
814      {
815        right=Data(m_right->collapseToReady());
816      }
817      Data result;
818      switch(m_op)
819      {
820        case ADD:
821        result=left+right;
822        break;
823        case SUB:      
824        result=left-right;
825        break;
826        case MUL:      
827        result=left*right;
828        break;
829        case DIV:      
830        result=left/right;
831        break;
832        case SIN:
833        result=left.sin();  
834        break;
835        case COS:
836        result=left.cos();
837        break;
838        case TAN:
839        result=left.tan();
840        break;
841        case ASIN:
842        result=left.asin();
843        break;
844        case ACOS:
845        result=left.acos();
846        break;
847        case ATAN:
848        result=left.atan();
849        break;
850        case SINH:
851        result=left.sinh();
852        break;
853        case COSH:
854        result=left.cosh();
855        break;
856        case TANH:
857        result=left.tanh();
858        break;
859        case ERF:
860        result=left.erf();
861        break;
862       case ASINH:
863        result=left.asinh();
864        break;
865       case ACOSH:
866        result=left.acosh();
867        break;
868       case ATANH:
869        result=left.atanh();
870        break;
871        case LOG10:
872        result=left.log10();
873        break;
874        case LOG:
875        result=left.log();
876        break;
877        case SIGN:
878        result=left.sign();
879        break;
880        case ABS:
881        result=left.abs();
882        break;
883        case NEG:
884        result=left.neg();
885        break;
886        case POS:
887        // it doesn't mean anything for delayed.
888        // it will just trigger a deep copy of the lazy object
889        throw DataException("Programmer error - POS not supported for lazy data.");
890        break;
891        case EXP:
892        result=left.exp();
893        break;
894        case SQRT:
895        result=left.sqrt();
896        break;
897        case RECIP:
898        result=left.oneOver();
899        break;
900        case GZ:
901        result=left.wherePositive();
902        break;
903        case LZ:
904        result=left.whereNegative();
905        break;
906        case GEZ:
907        result=left.whereNonNegative();
908        break;
909        case LEZ:
910        result=left.whereNonPositive();
911        break;
912        case NEZ:
913        result=left.whereNonZero(m_tol);
914        break;
915        case EZ:
916        result=left.whereZero(m_tol);
917        break;
918        case SYM:
919        result=left.symmetric();
920        break;
921        case NSYM:
922        result=left.nonsymmetric();
923        break;
924        case PROD:
925        result=C_GeneralTensorProduct(left,right,m_axis_offset, m_transpose);
926        break;
927        case TRANS:
928        result=left.transpose(m_axis_offset);
929        break;
930        case TRACE:
931        result=left.trace(m_axis_offset);
932        break;
933        case SWAP:
934        result=left.swapaxes(m_axis_offset, m_transpose);
935        break;
936        default:
937        throw DataException("Programmer error - collapseToReady does not know how to resolve operator "+opToString(m_op)+".");
938      }
939      return result.borrowReadyPtr();
940    }
941    
942    /*
943       \brief Converts the DataLazy into an IDENTITY storing the value of the expression.
944       This method uses the original methods on the Data class to evaluate the expressions.
945       For this reason, it should not be used on DataExpanded instances. (To do so would defeat
946       the purpose of using DataLazy in the first place).
947    */
948    void
949    DataLazy::collapse()
950    {
951      if (m_op==IDENTITY)
952      {
953        return;
954      }
955      if (m_readytype=='E')
956      { // this is more an efficiency concern than anything else
957        throw DataException("Programmer Error - do not use collapse on Expanded data.");
958      }
959      m_id=collapseToReady();
960      m_op=IDENTITY;
961    }
962    
963    /*
964      \brief Compute the value of the expression (unary operation) for the given sample.
965      \return Vector which stores the value of the subexpression for the given sample.
966      \param v A vector to store intermediate results.
967      \param offset Index in v to begin storing results.
968      \param sampleNo Sample number to evaluate.
969      \param roffset (output parameter) the offset in the return vector where the result begins.
970    
971      The return value will be an existing vector so do not deallocate it.
972      If the result is stored in v it should be stored at the offset given.
973      Everything from offset to the end of v should be considered available for this method to use.
974    */
975    DataTypes::ValueType*
976    DataLazy::resolveUnary(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
977    {
978        // we assume that any collapsing has been done before we get here
979        // since we only have one argument we don't need to think about only
980        // processing single points.
981      if (m_readytype!='E')
982      {
983        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
984      }
985      const ValueType* vleft=m_left->resolveSample(v,offset,sampleNo,roffset);
986      const double* left=&((*vleft)[roffset]);
987      double* result=&(v[offset]);
988      roffset=offset;
989      switch (m_op)
990      {
991        case SIN:  
992        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
993        break;
994        case COS:
995        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
996        break;
997        case TAN:
998        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
999        break;
1000        case ASIN:
1001        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
1002        break;
1003        case ACOS:
1004        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
1005        break;
1006        case ATAN:
1007        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
1008        break;
1009        case SINH:
1010        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
1011        break;
1012        case COSH:
1013        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
1014        break;
1015        case TANH:
1016        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
1017        break;
1018        case ERF:
1019    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1020        throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1021    #else
1022        tensor_unary_operation(m_samplesize, left, result, ::erf);
1023        break;
1024    #endif
1025       case ASINH:
1026    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1027        tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
1028    #else
1029        tensor_unary_operation(m_samplesize, left, result, ::asinh);
1030    #endif  
1031        break;
1032       case ACOSH:
1033    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1034        tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
1035    #else
1036        tensor_unary_operation(m_samplesize, left, result, ::acosh);
1037    #endif  
1038        break;
1039       case ATANH:
1040    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1041        tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
1042    #else
1043        tensor_unary_operation(m_samplesize, left, result, ::atanh);
1044    #endif  
1045        break;
1046        case LOG10:
1047        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
1048        break;
1049        case LOG:
1050        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
1051        break;
1052        case SIGN:
1053        tensor_unary_operation(m_samplesize, left, result, escript::fsign);
1054        break;
1055        case ABS:
1056        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
1057        break;
1058        case NEG:
1059        tensor_unary_operation(m_samplesize, left, result, negate<double>());
1060        break;
1061        case POS:
1062        // it doesn't mean anything for delayed.
1063        // it will just trigger a deep copy of the lazy object
1064        throw DataException("Programmer error - POS not supported for lazy data.");
1065        break;
1066        case EXP:
1067        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
1068        break;
1069        case SQRT:
1070        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
1071        break;
1072        case RECIP:
1073        tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
1074        break;
1075        case GZ:
1076        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
1077        break;
1078        case LZ:
1079        tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
1080        break;
1081        case GEZ:
1082        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
1083        break;
1084        case LEZ:
1085        tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
1086        break;
1087    // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
1088        case NEZ:
1089        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
1090        break;
1091        case EZ:
1092        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
1093        break;
1094    
1095        default:
1096        throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1097      }
1098      return &v;
1099    }
1100    
1101    
1102    
1103    
1104    
1105    
1106    /*
1107      \brief Compute the value of the expression (unary operation) for the given sample.
1108      \return Vector which stores the value of the subexpression for the given sample.
1109      \param v A vector to store intermediate results.
1110      \param offset Index in v to begin storing results.
1111      \param sampleNo Sample number to evaluate.
1112      \param roffset (output parameter) the offset in the return vector where the result begins.
1113    
1114      The return value will be an existing vector so do not deallocate it.
1115      If the result is stored in v it should be stored at the offset given.
1116      Everything from offset to the end of v should be considered available for this method to use.
1117    */
1118    DataTypes::ValueType*
1119    DataLazy::resolveNP1OUT(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1120    {
1121        // we assume that any collapsing has been done before we get here
1122        // since we only have one argument we don't need to think about only
1123        // processing single points.
1124      if (m_readytype!='E')
1125      {
1126        throw DataException("Programmer error - resolveNP1OUT should only be called on expanded Data.");
1127      }
1128        // since we can't write the result over the input, we need a result offset further along
1129      size_t subroffset=roffset+m_samplesize;
1130    LAZYDEBUG(cerr << "subroffset=" << subroffset << endl;)
1131      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1132      roffset=offset;
1133      size_t loop=0;
1134      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1135      size_t step=getNoValues();
1136      switch (m_op)
1137      {
1138        case SYM:
1139        for (loop=0;loop<numsteps;++loop)
1140        {
1141            DataMaths::symmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
1142            subroffset+=step;
1143            offset+=step;
1144        }
1145        break;
1146        case NSYM:
1147        for (loop=0;loop<numsteps;++loop)
1148        {
1149            DataMaths::nonsymmetric(*vleft,m_left->getShape(),subroffset, v, getShape(), offset);
1150            subroffset+=step;
1151            offset+=step;
1152        }
1153        break;
1154        default:
1155        throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1156      }
1157      return &v;
1158    }
1159    
1160    /*
1161      \brief Compute the value of the expression (unary operation) for the given sample.
1162      \return Vector which stores the value of the subexpression for the given sample.
1163      \param v A vector to store intermediate results.
1164      \param offset Index in v to begin storing results.
1165      \param sampleNo Sample number to evaluate.
1166      \param roffset (output parameter) the offset in the return vector where the result begins.
1167    
1168      The return value will be an existing vector so do not deallocate it.
1169      If the result is stored in v it should be stored at the offset given.
1170      Everything from offset to the end of v should be considered available for this method to use.
1171    */
1172    DataTypes::ValueType*
1173    DataLazy::resolveNP1OUT_P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1174    {
1175        // we assume that any collapsing has been done before we get here
1176        // since we only have one argument we don't need to think about only
1177        // processing single points.
1178      if (m_readytype!='E')
1179      {
1180        throw DataException("Programmer error - resolveNP1OUT_P should only be called on expanded Data.");
1181      }
1182        // since we can't write the result over the input, we need a result offset further along
1183      size_t subroffset;
1184      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1185    LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)
1186    LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)
1187      roffset=offset;
1188      size_t loop=0;
1189      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1190      size_t outstep=getNoValues();
1191      size_t instep=m_left->getNoValues();
1192    LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)
1193      switch (m_op)
1194      {
1195        case TRACE:
1196        for (loop=0;loop<numsteps;++loop)
1197        {
1198    size_t zz=sampleNo*getNumDPPSample()+loop;
1199    if (zz==5800)
1200    {
1201    LAZYDEBUG(cerr << "point=" <<  zz<< endl;)
1202    LAZYDEBUG(cerr << "Input to  trace=" << DataTypes::pointToString(*vleft,m_left->getShape(),subroffset,"") << endl;)
1203    LAZYDEBUG(cerr << "Offset for point=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << " vs ";)
1204    LAZYDEBUG(cerr << subroffset << endl;)
1205    LAZYDEBUG(cerr << "output=" << offset << endl;)
1206    }
1207                DataMaths::trace(*vleft,m_left->getShape(),subroffset, v ,getShape(),offset,m_axis_offset);
1208    if (zz==5800)
1209    {
1210    LAZYDEBUG(cerr << "Result of trace=" << DataTypes::pointToString(v,getShape(),offset,"") << endl;)
1211    }
1212            subroffset+=instep;
1213            offset+=outstep;
1214        }
1215        break;
1216        case TRANS:
1217        for (loop=0;loop<numsteps;++loop)
1218        {
1219                DataMaths::transpose(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset);
1220            subroffset+=instep;
1221            offset+=outstep;
1222        }
1223        break;
1224        default:
1225        throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1226      }
1227      return &v;
1228    }
1229    
1230    
1231    /*
1232      \brief Compute the value of the expression (unary operation with int params) for the given sample.
1233      \return Vector which stores the value of the subexpression for the given sample.
1234      \param v A vector to store intermediate results.
1235      \param offset Index in v to begin storing results.
1236      \param sampleNo Sample number to evaluate.
1237      \param roffset (output parameter) the offset in the return vector where the result begins.
1238    
1239      The return value will be an existing vector so do not deallocate it.
1240      If the result is stored in v it should be stored at the offset given.
1241      Everything from offset to the end of v should be considered available for this method to use.
1242    */
1243    DataTypes::ValueType*
1244    DataLazy::resolveNP1OUT_2P(ValueType& v, size_t offset, int sampleNo, size_t& roffset) const
1245    {
1246        // we assume that any collapsing has been done before we get here
1247        // since we only have one argument we don't need to think about only
1248        // processing single points.
1249      if (m_readytype!='E')
1250      {
1251        throw DataException("Programmer error - resolveNP1OUT_2P should only be called on expanded Data.");
1252      }
1253        // since we can't write the result over the input, we need a result offset further along
1254      size_t subroffset;
1255      const ValueType* vleft=m_left->resolveSample(v,offset+m_left->m_samplesize,sampleNo,subroffset);
1256    LAZYDEBUG(cerr << "srcsamplesize=" << offset+m_left->m_samplesize << " beg=" << subroffset << endl;)
1257    LAZYDEBUG(cerr << "Offset for 5800=" << getPointOffset(5800/getNumDPPSample(),5800%getNumDPPSample()) << endl;)
1258      roffset=offset;
1259      size_t loop=0;
1260      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1261      size_t outstep=getNoValues();
1262      size_t instep=m_left->getNoValues();
1263    LAZYDEBUG(cerr << "instep=" << instep << " outstep=" << outstep<< " numsteps=" << numsteps << endl;)
1264      switch (m_op)
1265      {
1266        case SWAP:
1267        for (loop=0;loop<numsteps;++loop)
1268        {
1269                DataMaths::swapaxes(*vleft,m_left->getShape(),subroffset, v,getShape(),offset,m_axis_offset, m_transpose);
1270            subroffset+=instep;
1271            offset+=outstep;
1272        }
1273        break;
1274        default:
1275        throw DataException("Programmer error - resolveNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1276      }
1277      return &v;
1278    }
1279    
1280    
1281    
1282    #define PROC_OP(TYPE,X)                               \
1283        for (int j=0;j<onumsteps;++j)\
1284        {\
1285          for (int i=0;i<numsteps;++i,resultp+=resultStep) \
1286          { \
1287    LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\
1288    LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\
1289             tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \
1290    LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \
1291             lroffset+=leftstep; \
1292             rroffset+=rightstep; \
1293          }\
1294          lroffset+=oleftstep;\
1295          rroffset+=orightstep;\
1296        }
1297    
1298    /*
1299      \brief Compute the value of the expression (binary operation) for the given sample.
1300      \return Vector which stores the value of the subexpression for the given sample.
1301      \param v A vector to store intermediate results.
1302      \param offset Index in v to begin storing results.
1303      \param sampleNo Sample number to evaluate.
1304      \param roffset (output parameter) the offset in the return vector where the result begins.
1305    
1306      The return value will be an existing vector so do not deallocate it.
1307      If the result is stored in v it should be stored at the offset given.
1308      Everything from offset to the end of v should be considered available for this method to use.
1309    */
1310    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1311    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1312    // If both children are expanded, then we can process them in a single operation (we treat
1313    // the whole sample as one big datapoint.
1314    // If one of the children is not expanded, then we need to treat each point in the sample
1315    // individually.
1316    // There is an additional complication when scalar operations are considered.
1317    // For example, 2+Vector.
1318    // In this case each double within the point is treated individually
1319    DataTypes::ValueType*
1320    DataLazy::resolveBinary(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
1321    {
1322    LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1323    
1324      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
1325        // first work out which of the children are expanded
1326      bool leftExp=(m_left->m_readytype=='E');
1327      bool rightExp=(m_right->m_readytype=='E');
1328      if (!leftExp && !rightExp)
1329      {
1330        throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1331      }
1332      bool leftScalar=(m_left->getRank()==0);
1333      bool rightScalar=(m_right->getRank()==0);
1334      if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1335      {
1336        throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1337      }
1338      size_t leftsize=m_left->getNoValues();
1339      size_t rightsize=m_right->getNoValues();
1340      size_t chunksize=1;           // how many doubles will be processed in one go
1341      int leftstep=0;       // how far should the left offset advance after each step
1342      int rightstep=0;
1343      int numsteps=0;       // total number of steps for the inner loop
1344      int oleftstep=0;  // the o variables refer to the outer loop
1345      int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1346      int onumsteps=1;
1347      
1348      bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1349      bool RES=(rightExp && rightScalar);
1350      bool LS=(!leftExp && leftScalar); // left is a single scalar
1351      bool RS=(!rightExp && rightScalar);
1352      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1353      bool RN=(!rightExp && !rightScalar);
1354      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1355      bool REN=(rightExp && !rightScalar);
1356    
1357      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1358      {
1359        chunksize=m_left->getNumDPPSample()*leftsize;
1360        leftstep=0;
1361        rightstep=0;
1362        numsteps=1;
1363      }
1364      else if (LES || RES)
1365      {
1366        chunksize=1;
1367        if (LES)        // left is an expanded scalar
1368        {
1369            if (RS)
1370            {
1371               leftstep=1;
1372               rightstep=0;
1373               numsteps=m_left->getNumDPPSample();
1374            }
1375            else        // RN or REN
1376            {
1377               leftstep=0;
1378               oleftstep=1;
1379               rightstep=1;
1380               orightstep=(RN ? -(int)rightsize : 0);
1381               numsteps=rightsize;
1382               onumsteps=m_left->getNumDPPSample();
1383            }
1384        }
1385        else        // right is an expanded scalar
1386        {
1387            if (LS)
1388            {
1389               rightstep=1;
1390               leftstep=0;
1391               numsteps=m_right->getNumDPPSample();
1392            }
1393            else
1394            {
1395               rightstep=0;
1396               orightstep=1;
1397               leftstep=1;
1398               oleftstep=(LN ? -(int)leftsize : 0);
1399               numsteps=leftsize;
1400               onumsteps=m_right->getNumDPPSample();
1401            }
1402        }
1403      }
1404      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1405      {
1406        if (LEN)    // and Right will be a single value
1407        {
1408            chunksize=rightsize;
1409            leftstep=rightsize;
1410            rightstep=0;
1411            numsteps=m_left->getNumDPPSample();
1412            if (RS)
1413            {
1414               numsteps*=leftsize;
1415            }
1416        }
1417        else    // REN
1418        {
1419            chunksize=leftsize;
1420            rightstep=leftsize;
1421            leftstep=0;
1422            numsteps=m_right->getNumDPPSample();
1423            if (LS)
1424            {
1425               numsteps*=rightsize;
1426            }
1427        }
1428      }
1429    
1430      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
1431        // Get the values of sub-expressions
1432      const ValueType* left=m_left->resolveSample(v,offset+getMaxSampleSize(),sampleNo,lroffset);   // see note on
1433        // calcBufss for why we can't put offset as the 2nd param above
1434      const ValueType* right=m_right->resolveSample(v,offset+2*getMaxSampleSize(),sampleNo,rroffset); // Note
1435        // the right child starts further along.
1436    LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1437    LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1438    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
1439    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
1440    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
1441    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
1442    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
1443    
1444    
1445      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
1446      switch(m_op)
1447      {
1448        case ADD:
1449            PROC_OP(NO_ARG,plus<double>());
1450        break;
1451        case SUB:
1452        PROC_OP(NO_ARG,minus<double>());
1453        break;
1454        case MUL:
1455        PROC_OP(NO_ARG,multiplies<double>());
1456        break;
1457        case DIV:
1458        PROC_OP(NO_ARG,divides<double>());
1459        break;
1460        case POW:
1461           PROC_OP(double (double,double),::pow);
1462        break;
1463        default:
1464        throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
1465      }
1466      roffset=offset;
1467      return &v;
1468    }
1469    
1470    
1471    
1472    /*
1473      \brief Compute the value of the expression (tensor product) for the given sample.
1474      \return Vector which stores the value of the subexpression for the given sample.
1475      \param v A vector to store intermediate results.
1476      \param offset Index in v to begin storing results.
1477      \param sampleNo Sample number to evaluate.
1478      \param roffset (output parameter) the offset in the return vector where the result begins.
1479    
1480      The return value will be an existing vector so do not deallocate it.
1481      If the result is stored in v it should be stored at the offset given.
1482      Everything from offset to the end of v should be considered available for this method to use.
1483    */
1484    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1485    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1486    // unlike the other resolve helpers, we must treat these datapoints separately.
1487    DataTypes::ValueType*
1488    DataLazy::resolveTProd(ValueType& v,  size_t offset, int sampleNo, size_t& roffset) const
1489    {
1490    LAZYDEBUG(cout << "Resolve TensorProduct: " << toString()  << " to offset " << offset<< endl;)
1491    
1492      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
1493        // first work out which of the children are expanded
1494      bool leftExp=(m_left->m_readytype=='E');
1495      bool rightExp=(m_right->m_readytype=='E');
1496      int steps=getNumDPPSample();
1497    /*  int leftStep=((leftExp && !rightExp)? m_right->getNoValues() : 0);
1498      int rightStep=((rightExp && !leftExp)? m_left->getNoValues() : 0);*/
1499      int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method
1500      int rightStep=(rightExp?m_right->getNoValues() : 0);
1501    
1502      int resultStep=getNoValues();
1503        // Get the values of sub-expressions (leave a gap of one sample for the result).
1504      int gap=offset+m_samplesize;  
1505    
1506    LAZYDEBUG(cout << "Query left with offset=" << gap << endl;)
1507    
1508      const ValueType* left=m_left->resolveSample(v,gap,sampleNo,lroffset);
1509      gap+=m_left->getMaxSampleSize();
1510    
1511    
1512    LAZYDEBUG(cout << "Query right with offset=" << gap << endl;)
1513    
1514    
1515      const ValueType* right=m_right->resolveSample(v,gap,sampleNo,rroffset);
1516    
1517    LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
1518    cout << getNoValues() << endl;)
1519    LAZYDEBUG(cerr << "Result of left=";)
1520    LAZYDEBUG(cerr << "[" << lroffset << " .. " << lroffset+m_left->getNoValues() << "]" << endl;
1521    
1522    for (int i=lroffset, limit=lroffset+(leftExp?m_left->getNoValues()*m_left->getNumDPPSample():m_left->getNoValues());i<limit;++i)
1523    {
1524    cout << "[" << setw(2) << i-lroffset << "] " << setw(10) << (*left)[i] << " ";
1525    if (i%4==0) cout << endl;
1526    })
1527    LAZYDEBUG(cerr << "\nResult of right=" << endl;)
1528    LAZYDEBUG(
1529    for (int i=rroffset, limit=rroffset+(rightExp?m_right->getNoValues()*m_right->getNumDPPSample():m_right->getNoValues());i<limit;++i)
1530    {
1531    cerr << "[" <<  setw(2)<< i-rroffset << "] " << setw(10) << (*right)[i] << " ";
1532    if (i%4==0) cout << endl;
1533    }
1534    cerr << endl;
1535    )
1536    LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
1537    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
1538    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
1539    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
1540    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
1541    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
1542    LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
1543    
1544      double* resultp=&(v[offset]);     // results are stored at the vector offset we recieved
1545      switch(m_op)
1546      {
1547        case PROD:
1548        for (int i=0;i<steps;++i,resultp+=resultStep)
1549        {
1550    
1551    LAZYDEBUG(cout << "lroffset=" << lroffset << "rroffset=" << rroffset << endl;)
1552    LAZYDEBUG(cout << "l*=" << left << " r*=" << right << endl;)
1553    LAZYDEBUG(cout << "m_SL=" << m_SL << " m_SM=" << m_SM << " m_SR=" << m_SR << endl;)
1554    
1555              const double *ptr_0 = &((*left)[lroffset]);
1556              const double *ptr_1 = &((*right)[rroffset]);
1557    
1558    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
1559    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
1560    
1561              matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
1562    
1563    LAZYDEBUG(cout << "Results--\n";
1564    {
1565      DataVector dv(getNoValues());
1566    for (int z=0;z<getNoValues();++z)
1567    {
1568      cout << "[" << setw(2) << z<< "] " << setw(10) << resultp[z] << " ";
1569      if (z%4==0) cout << endl;
1570      dv[z]=resultp[z];
1571    }
1572    cout << endl << DataTypes::pointToString(dv,getShape(),0,"RESLT");
1573    cout << "\nWritten to: " << resultp << " resultStep=" << resultStep << endl;
1574    }
1575    )
1576          lroffset+=leftStep;
1577          rroffset+=rightStep;
1578        }
1579        break;
1580        default:
1581        throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
1582      }
1583      roffset=offset;
1584      return &v;
1585    }
1586    
1587    
1588    #ifdef LAZY_NODE_STORAGE
1589    
1590    // The result will be stored in m_samples
1591    // The return value is a pointer to the DataVector, offset is the offset within the return value
1592    const DataTypes::ValueType*
1593    DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset)
1594    {
1595    LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1596        // collapse so we have a 'E' node or an IDENTITY for some other type
1597      if (m_readytype!='E' && m_op!=IDENTITY)
1598      {
1599        collapse();
1600      }
1601      if (m_op==IDENTITY)  
1602      {
1603        const ValueType& vec=m_id->getVectorRO();
1604    //     if (m_readytype=='C')
1605    //     {
1606    //  roffset=0;      // all samples read from the same position
1607    //  return &(m_samples);
1608    //     }
1609        roffset=m_id->getPointOffset(sampleNo, 0);
1610        return &(vec);
1611      }
1612      if (m_readytype!='E')
1613      {
1614        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
1615      }
1616      if (m_sampleids[tid]==sampleNo)
1617      {
1618        roffset=tid*m_samplesize;
1619        return &(m_samples);        // sample is already resolved
1620      }
1621      m_sampleids[tid]=sampleNo;
1622      switch (getOpgroup(m_op))
1623      {
1624      case G_UNARY:
1625      case G_UNARY_P: return resolveNodeUnary(tid, sampleNo, roffset);
1626      case G_BINARY: return resolveNodeBinary(tid, sampleNo, roffset);
1627      case G_NP1OUT: return resolveNodeNP1OUT(tid, sampleNo, roffset);
1628      case G_NP1OUT_P: return resolveNodeNP1OUT_P(tid, sampleNo, roffset);
1629      case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);
1630      case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);
1631      default:
1632        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
1633      }
1634    }
1635    
1636    const DataTypes::ValueType*
1637    DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset)
1638    {
1639        // we assume that any collapsing has been done before we get here
1640        // since we only have one argument we don't need to think about only
1641        // processing single points.
1642        // we will also know we won't get identity nodes
1643      if (m_readytype!='E')
1644      {
1645        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1646      }
1647      if (m_op==IDENTITY)
1648      {
1649        throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1650      }
1651      const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
1652      const double* left=&((*leftres)[roffset]);
1653      roffset=m_samplesize*tid;
1654      double* result=&(m_samples[roffset]);
1655      switch (m_op)
1656      {
1657        case SIN:  
1658        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
1659        break;
1660        case COS:
1661        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
1662        break;
1663        case TAN:
1664        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
1665        break;
1666        case ASIN:
1667        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
1668        break;
1669        case ACOS:
1670        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
1671        break;
1672        case ATAN:
1673        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
1674        break;
1675        case SINH:
1676        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
1677        break;
1678        case COSH:
1679        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
1680        break;
1681        case TANH:
1682        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
1683        break;
1684        case ERF:
1685    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1686        throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1687    #else
1688        tensor_unary_operation(m_samplesize, left, result, ::erf);
1689        break;
1690    #endif
1691       case ASINH:
1692    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1693        tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
1694    #else
1695        tensor_unary_operation(m_samplesize, left, result, ::asinh);
1696    #endif  
1697        break;
1698       case ACOSH:
1699    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1700        tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
1701    #else
1702        tensor_unary_operation(m_samplesize, left, result, ::acosh);
1703    #endif  
1704        break;
1705       case ATANH:
1706    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1707        tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
1708    #else
1709        tensor_unary_operation(m_samplesize, left, result, ::atanh);
1710    #endif  
1711        break;
1712        case LOG10:
1713        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
1714        break;
1715        case LOG:
1716        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
1717        break;
1718        case SIGN:
1719        tensor_unary_operation(m_samplesize, left, result, escript::fsign);
1720        break;
1721        case ABS:
1722        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
1723        break;
1724        case NEG:
1725        tensor_unary_operation(m_samplesize, left, result, negate<double>());
1726        break;
1727        case POS:
1728        // it doesn't mean anything for delayed.
1729        // it will just trigger a deep copy of the lazy object
1730        throw DataException("Programmer error - POS not supported for lazy data.");
1731        break;
1732        case EXP:
1733        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
1734        break;
1735        case SQRT:
1736        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
1737        break;
1738        case RECIP:
1739        tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
1740        break;
1741        case GZ:
1742        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
1743        break;
1744        case LZ:
1745        tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
1746        break;
1747        case GEZ:
1748        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
1749        break;
1750        case LEZ:
1751        tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
1752        break;
1753    // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
1754        case NEZ:
1755        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
1756        break;
1757        case EZ:
1758        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
1759        break;
1760    
1761        default:
1762        throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1763      }
1764      return &(m_samples);
1765    }
1766    
1767    
1768    const DataTypes::ValueType*
1769    DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset)
1770    {
1771        // we assume that any collapsing has been done before we get here
1772        // since we only have one argument we don't need to think about only
1773        // processing single points.
1774      if (m_readytype!='E')
1775      {
1776        throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");
1777      }
1778      if (m_op==IDENTITY)
1779      {
1780        throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1781      }
1782      size_t subroffset;
1783      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1784      roffset=m_samplesize*tid;
1785      size_t loop=0;
1786      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1787      size_t step=getNoValues();
1788      size_t offset=roffset;
1789      switch (m_op)
1790      {
1791        case SYM:
1792        for (loop=0;loop<numsteps;++loop)
1793        {
1794            DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1795            subroffset+=step;
1796            offset+=step;
1797        }
1798        break;
1799        case NSYM:
1800        for (loop=0;loop<numsteps;++loop)
1801        {
1802            DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1803            subroffset+=step;
1804            offset+=step;
1805        }
1806        break;
1807        default:
1808        throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1809      }
1810      return &m_samples;
1811    }
1812    
1813    const DataTypes::ValueType*
1814    DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset)
1815    {
1816        // we assume that any collapsing has been done before we get here
1817        // since we only have one argument we don't need to think about only
1818        // processing single points.
1819      if (m_readytype!='E')
1820      {
1821        throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");
1822      }
1823      if (m_op==IDENTITY)
1824      {
1825        throw DataException("Programmer error - resolveNodeNP1OUT_P should not be called on identity nodes.");
1826      }
1827      size_t subroffset;
1828      size_t offset;
1829      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1830      roffset=m_samplesize*tid;
1831      offset=roffset;
1832      size_t loop=0;
1833      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1834      size_t outstep=getNoValues();
1835      size_t instep=m_left->getNoValues();
1836      switch (m_op)
1837      {
1838        case TRACE:
1839        for (loop=0;loop<numsteps;++loop)
1840        {
1841                DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);
1842            subroffset+=instep;
1843            offset+=outstep;
1844        }
1845        break;
1846        case TRANS:
1847        for (loop=0;loop<numsteps;++loop)
1848        {
1849                DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);
1850            subroffset+=instep;
1851            offset+=outstep;
1852        }
1853        break;
1854        default:
1855        throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1856      }
1857      return &m_samples;
1858    }
1859    
1860    
1861    const DataTypes::ValueType*
1862    DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset)
1863    {
1864      if (m_readytype!='E')
1865      {
1866        throw DataException("Programmer error - resolveNodeNP1OUT_2P should only be called on expanded Data.");
1867      }
1868      if (m_op==IDENTITY)
1869      {
1870        throw DataException("Programmer error - resolveNodeNP1OUT_2P should not be called on identity nodes.");
1871      }
1872      size_t subroffset;
1873      size_t offset;
1874      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1875      roffset=m_samplesize*tid;
1876      offset=roffset;
1877      size_t loop=0;
1878      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1879      size_t outstep=getNoValues();
1880      size_t instep=m_left->getNoValues();
1881      switch (m_op)
1882      {
1883        case SWAP:
1884        for (loop=0;loop<numsteps;++loop)
1885        {
1886                DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);
1887            subroffset+=instep;
1888            offset+=outstep;
1889        }
1890        break;
1891        default:
1892        throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1893      }
1894      return &m_samples;
1895    }
1896    
1897    
1898    
1899    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1900    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1901    // If both children are expanded, then we can process them in a single operation (we treat
1902    // the whole sample as one big datapoint.
1903    // If one of the children is not expanded, then we need to treat each point in the sample
1904    // individually.
1905    // There is an additional complication when scalar operations are considered.
1906    // For example, 2+Vector.
1907    // In this case each double within the point is treated individually
1908    const DataTypes::ValueType*
1909    DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset)
1910    {
1911    LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1912    
1913      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
1914        // first work out which of the children are expanded
1915      bool leftExp=(m_left->m_readytype=='E');
1916      bool rightExp=(m_right->m_readytype=='E');
1917      if (!leftExp && !rightExp)
1918      {
1919        throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1920      }
1921      bool leftScalar=(m_left->getRank()==0);
1922      bool rightScalar=(m_right->getRank()==0);
1923      if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1924      {
1925        throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1926      }
1927      size_t leftsize=m_left->getNoValues();
1928      size_t rightsize=m_right->getNoValues();
1929      size_t chunksize=1;           // how many doubles will be processed in one go
1930      int leftstep=0;       // how far should the left offset advance after each step
1931      int rightstep=0;
1932      int numsteps=0;       // total number of steps for the inner loop
1933      int oleftstep=0;  // the o variables refer to the outer loop
1934      int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1935      int onumsteps=1;
1936      
1937      bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1938      bool RES=(rightExp && rightScalar);
1939      bool LS=(!leftExp && leftScalar); // left is a single scalar
1940      bool RS=(!rightExp && rightScalar);
1941      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1942      bool RN=(!rightExp && !rightScalar);
1943      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1944      bool REN=(rightExp && !rightScalar);
1945    
1946      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1947      {
1948        chunksize=m_left->getNumDPPSample()*leftsize;
1949        leftstep=0;
1950        rightstep=0;
1951        numsteps=1;
1952      }
1953      else if (LES || RES)
1954      {
1955        chunksize=1;
1956        if (LES)        // left is an expanded scalar
1957        {
1958            if (RS)
1959            {
1960               leftstep=1;
1961               rightstep=0;
1962               numsteps=m_left->getNumDPPSample();
1963            }
1964            else        // RN or REN
1965            {
1966               leftstep=0;
1967               oleftstep=1;
1968               rightstep=1;
1969               orightstep=(RN ? -(int)rightsize : 0);
1970               numsteps=rightsize;
1971               onumsteps=m_left->getNumDPPSample();
1972            }
1973        }
1974        else        // right is an expanded scalar
1975        {
1976            if (LS)
1977            {
1978               rightstep=1;
1979               leftstep=0;
1980               numsteps=m_right->getNumDPPSample();
1981            }
1982            else
1983            {
1984               rightstep=0;
1985               orightstep=1;
1986               leftstep=1;
1987               oleftstep=(LN ? -(int)leftsize : 0);
1988               numsteps=leftsize;
1989               onumsteps=m_right->getNumDPPSample();
1990            }
1991        }
1992      }
1993      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1994      {
1995        if (LEN)    // and Right will be a single value
1996        {
1997            chunksize=rightsize;
1998            leftstep=rightsize;
1999            rightstep=0;
2000            numsteps=m_left->getNumDPPSample();
2001            if (RS)
2002            {
2003               numsteps*=leftsize;
2004            }
2005        }
2006        else    // REN
2007        {
2008            chunksize=leftsize;
2009            rightstep=leftsize;
2010            leftstep=0;
2011            numsteps=m_right->getNumDPPSample();
2012            if (LS)
2013            {
2014               numsteps*=rightsize;
2015            }
2016        }
2017      }
2018    
2019      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
2020        // Get the values of sub-expressions
2021      const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);  
2022      const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
2023    LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
2024    LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
2025    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
2026    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
2027    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
2028    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
2029    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
2030    
2031    LAZYDEBUG(cout << "Left res["<< lroffset<< "]=" << (*left)[lroffset] << endl;)
2032    LAZYDEBUG(cout << "Right res["<< rroffset<< "]=" << (*right)[rroffset] << endl;)
2033    
2034    
2035      roffset=m_samplesize*tid;
2036      double* resultp=&(m_samples[roffset]);        // results are stored at the vector offset we recieved
2037      switch(m_op)
2038    {    {
2039      throw DataException("Programmer error - do not know how to resolve operator "+opToString(m_op)+".");      case ADD:
2040            PROC_OP(NO_ARG,plus<double>());
2041        break;
2042        case SUB:
2043        PROC_OP(NO_ARG,minus<double>());
2044        break;
2045        case MUL:
2046        PROC_OP(NO_ARG,multiplies<double>());
2047        break;
2048        case DIV:
2049        PROC_OP(NO_ARG,divides<double>());
2050        break;
2051        case POW:
2052           PROC_OP(double (double,double),::pow);
2053        break;
2054        default:
2055        throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
2056    }    }
2057    LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
2058      return &m_samples;
2059    }
2060    
2061    
2062    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
2063    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
2064    // unlike the other resolve helpers, we must treat these datapoints separately.
2065    const DataTypes::ValueType*
2066    DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset)
2067    {
2068    LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
2069    
2070      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
2071        // first work out which of the children are expanded
2072      bool leftExp=(m_left->m_readytype=='E');
2073      bool rightExp=(m_right->m_readytype=='E');
2074      int steps=getNumDPPSample();
2075      int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method
2076      int rightStep=(rightExp?m_right->getNoValues() : 0);
2077    
2078      int resultStep=getNoValues();
2079      roffset=m_samplesize*tid;
2080      size_t offset=roffset;
2081    
2082      const ValueType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
2083    
2084      const ValueType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
2085    
2086    LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
2087    cout << getNoValues() << endl;)
2088    
2089    
2090    LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
2091    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
2092    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
2093    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
2094    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
2095    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
2096    LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
2097    
2098      double* resultp=&(m_samples[offset]);     // results are stored at the vector offset we recieved
2099      switch(m_op)
2100      {
2101        case PROD:
2102        for (int i=0;i<steps;++i,resultp+=resultStep)
2103        {
2104              const double *ptr_0 = &((*left)[lroffset]);
2105              const double *ptr_1 = &((*right)[rroffset]);
2106    
2107    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
2108    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
2109    
2110              matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
2111    
2112          lroffset+=leftStep;
2113          rroffset+=rightStep;
2114        }
2115        break;
2116        default:
2117        throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
2118      }
2119      roffset=offset;
2120      return &m_samples;
2121    }
2122    #endif //LAZY_NODE_STORAGE
2123    
2124    /*
2125      \brief Compute the value of the expression for the given sample.
2126      \return Vector which stores the value of the subexpression for the given sample.
2127      \param v A vector to store intermediate results.
2128      \param offset Index in v to begin storing results.
2129      \param sampleNo Sample number to evaluate.
2130      \param roffset (output parameter) the offset in the return vector where the result begins.
2131    
2132      The return value will be an existing vector so do not deallocate it.
2133    */
2134    // the vector and the offset are a place where the method could write its data if it wishes
2135    // it is not obligated to do so. For example, if it has its own storage already, it can use that.
2136    // Hence the return value to indicate where the data is actually stored.
2137    // Regardless, the storage should be assumed to be used, even if it isn't.
2138    
2139    // the roffset is the offset within the returned vector where the data begins
2140    const DataTypes::ValueType*
2141    DataLazy::resolveSample(ValueType& v, size_t offset, int sampleNo, size_t& roffset)
2142    {
2143    LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
2144        // collapse so we have a 'E' node or an IDENTITY for some other type
2145      if (m_readytype!='E' && m_op!=IDENTITY)
2146      {
2147        collapse();
2148      }
2149      if (m_op==IDENTITY)  
2150      {
2151        const ValueType& vec=m_id->getVectorRO();
2152        if (m_readytype=='C')
2153        {
2154        roffset=0;
2155    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
2156        return &(vec);
2157        }
2158        roffset=m_id->getPointOffset(sampleNo, 0);
2159    LAZYDEBUG(cout << "Finish  sample " << toString() << endl;)
2160        return &(vec);
2161      }
2162      if (m_readytype!='E')
2163      {
2164        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
2165      }
2166      switch (getOpgroup(m_op))
2167      {
2168      case G_UNARY:
2169      case G_UNARY_P: return resolveUnary(v, offset,sampleNo,roffset);
2170      case G_BINARY: return resolveBinary(v, offset,sampleNo,roffset);
2171      case G_NP1OUT: return resolveNP1OUT(v, offset, sampleNo,roffset);
2172      case G_NP1OUT_P: return resolveNP1OUT_P(v, offset, sampleNo,roffset);
2173      case G_TENSORPROD: return resolveTProd(v,offset, sampleNo,roffset);
2174      case G_NP1OUT_2P: return resolveNP1OUT_2P(v, offset, sampleNo, roffset);
2175      default:
2176        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
2177      }
2178    
2179    }
2180    
2181    const DataTypes::ValueType*
2182    DataLazy::resolveSample(BufferGroup& bg, int sampleNo, size_t& roffset)
2183    {
2184    #ifdef _OPENMP
2185        int tid=omp_get_thread_num();
2186    #else
2187        int tid=0;
2188    #endif
2189        return resolveSample(bg.getBuffer(tid),bg.getOffset(tid),sampleNo,roffset);
2190    }
2191    
2192    
2193    // This needs to do the work of the identity constructor
2194    void
2195    DataLazy::resolveToIdentity()
2196    {
2197       if (m_op==IDENTITY)
2198        return;
2199    #ifndef LAZY_NODE_STORAGE
2200       DataReady_ptr p=resolveVectorWorker();
2201    #else
2202       DataReady_ptr p=resolveNodeWorker();
2203    #endif
2204       makeIdentity(p);
2205    }
2206    
2207    void DataLazy::makeIdentity(const DataReady_ptr& p)
2208    {
2209       m_op=IDENTITY;
2210       m_axis_offset=0;
2211       m_transpose=0;
2212       m_SL=m_SM=m_SR=0;
2213       m_children=m_height=0;
2214       m_id=p;
2215       if(p->isConstant()) {m_readytype='C';}
2216       else if(p->isExpanded()) {m_readytype='E';}
2217       else if (p->isTagged()) {m_readytype='T';}
2218       else {throw DataException("Unknown DataReady instance in convertToIdentity constructor.");}
2219       m_buffsRequired=1;
2220       m_samplesize=p->getNumDPPSample()*p->getNoValues();
2221       m_maxsamplesize=m_samplesize;
2222       m_left.reset();
2223       m_right.reset();
2224    }
2225    
2226    
2227    DataReady_ptr
2228    DataLazy::resolve()
2229    {
2230        resolveToIdentity();
2231        return m_id;
2232    }
2233    
2234    #ifdef LAZY_NODE_STORAGE
2235    
2236    // This version of resolve uses storage in each node to hold results
2237    DataReady_ptr
2238    DataLazy::resolveNodeWorker()
2239    {
2240      if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally
2241      {
2242        collapse();
2243      }
2244      if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
2245      {
2246        return m_id;
2247      }
2248        // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
2249      DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
2250      ValueType& resvec=result->getVectorRW();
2251      DataReady_ptr resptr=DataReady_ptr(result);
2252    
2253      int sample;
2254      int totalsamples=getNumSamples();
2255      const ValueType* res=0;   // Storage for answer
2256    LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
2257      #pragma omp parallel for private(sample,res) schedule(static)
2258      for (sample=0;sample<totalsamples;++sample)
2259      {
2260        size_t roffset=0;
2261    #ifdef _OPENMP
2262        res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
2263    #else
2264        res=resolveNodeSample(0,sample,roffset);
2265    #endif
2266    LAZYDEBUG(cout << "Sample #" << sample << endl;)
2267    LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
2268        DataVector::size_type outoffset=result->getPointOffset(sample,0);
2269        memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(DataVector::ElementType));
2270      }
2271      return resptr;
2272    }
2273    
2274    #endif // LAZY_NODE_STORAGE
2275    
2276    // To simplify the memory management, all threads operate on one large vector, rather than one each.
2277    // Each sample is evaluated independently and copied into the result DataExpanded.
2278    DataReady_ptr
2279    DataLazy::resolveVectorWorker()
2280    {
2281    
2282    LAZYDEBUG(cout << "Sample size=" << m_samplesize << endl;)
2283    LAZYDEBUG(cout << "Buffers=" << m_buffsRequired << endl;)
2284      if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally
2285      {
2286        collapse();
2287      }
2288      if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
2289      {
2290        return m_id;
2291      }
2292        // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
2293      size_t threadbuffersize=m_maxsamplesize*(max(1,m_buffsRequired)); // Each thread needs to have enough
2294        // storage to evaluate its expression
2295      int numthreads=1;
2296    #ifdef _OPENMP
2297      numthreads=omp_get_max_threads();
2298    #endif
2299      ValueType v(numthreads*threadbuffersize);
2300    LAZYDEBUG(cout << "Buffer created with size=" << v.size() << endl;)
2301      DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
2302      ValueType& resvec=result->getVectorRW();
2303      DataReady_ptr resptr=DataReady_ptr(result);
2304      int sample;
2305      size_t outoffset;     // offset in the output data
2306      int totalsamples=getNumSamples();
2307      const ValueType* res=0;   // Vector storing the answer
2308      size_t resoffset=0;       // where in the vector to find the answer
2309    LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
2310      #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
2311      for (sample=0;sample<totalsamples;++sample)
2312      {
2313    LAZYDEBUG(cout << "################################# " << sample << endl;)
2314    #ifdef _OPENMP
2315        res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
2316    #else
2317        res=resolveSample(v,0,sample,resoffset);   // res would normally be v, but not if its a single IDENTITY op.
2318    #endif
2319    LAZYDEBUG(cerr << "-------------------------------- " << endl;)
2320    LAZYDEBUG(cerr<< "Copying sample#" << sample << endl;)
2321        outoffset=result->getPointOffset(sample,0);
2322    LAZYDEBUG(cerr << "offset=" << outoffset << " from offset=" << resoffset << " " << m_samplesize << " doubles" << endl;)
2323        for (unsigned int i=0;i<m_samplesize;++i,++outoffset,++resoffset)   // copy values into the output vector
2324        {
2325    LAZYDEBUG(cerr << "outoffset=" << outoffset << " resoffset=" << resoffset << " " << (*res)[resoffset]<< endl;)
2326        resvec[outoffset]=(*res)[resoffset];
2327        }
2328    LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)
2329    LAZYDEBUG(cerr << "*********************************" << endl;)
2330      }
2331      return resptr;
2332  }  }
2333    
2334  std::string  std::string
2335  DataLazy::toString() const  DataLazy::toString() const
2336  {  {
2337    return "Lazy evaluation object. No details available.";    ostringstream oss;
2338      oss << "Lazy Data:";
2339      intoString(oss);
2340      return oss.str();
2341    }
2342    
2343    
2344    void
2345    DataLazy::intoString(ostringstream& oss) const
2346    {
2347    //    oss << "[" << m_children <<";"<<m_height <<"]";
2348      switch (getOpgroup(m_op))
2349      {
2350      case G_IDENTITY:
2351        if (m_id->isExpanded())
2352        {
2353           oss << "E";
2354        }
2355        else if (m_id->isTagged())
2356        {
2357          oss << "T";
2358        }
2359        else if (m_id->isConstant())
2360        {
2361          oss << "C";
2362        }
2363        else
2364        {
2365          oss << "?";
2366        }
2367        oss << '@' << m_id.get();
2368        break;
2369      case G_BINARY:
2370        oss << '(';
2371        m_left->intoString(oss);
2372        oss << ' ' << opToString(m_op) << ' ';
2373        m_right->intoString(oss);
2374        oss << ')';
2375        break;
2376      case G_UNARY:
2377      case G_UNARY_P:
2378      case G_NP1OUT:
2379      case G_NP1OUT_P:
2380        oss << opToString(m_op) << '(';
2381        m_left->intoString(oss);
2382        oss << ')';
2383        break;
2384      case G_TENSORPROD:
2385        oss << opToString(m_op) << '(';
2386        m_left->intoString(oss);
2387        oss << ", ";
2388        m_right->intoString(oss);
2389        oss << ')';
2390        break;
2391      case G_NP1OUT_2P:
2392        oss << opToString(m_op) << '(';
2393        m_left->intoString(oss);
2394        oss << ", " << m_axis_offset << ", " << m_transpose;
2395        oss << ')';
2396        break;
2397      default:
2398        oss << "UNKNOWN";
2399      }
2400  }  }
2401    
 // 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.  
2402  DataAbstract*  DataAbstract*
2403  DataLazy::deepCopy()  DataLazy::deepCopy()
2404  {  {
2405    if (m_op==IDENTITY)    switch (getOpgroup(m_op))
2406    {    {
2407      return new DataLazy(m_left);    // we don't need to copy the child here    case G_IDENTITY:  return new DataLazy(m_id->deepCopy()->getPtr());
2408      case G_UNARY: return new DataLazy(m_left->deepCopy()->getPtr(),m_op);
2409      case G_BINARY:    return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);
2410      case G_NP1OUT: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(),m_op);
2411      case G_TENSORPROD: return new DataLazy(m_left->deepCopy()->getPtr(), m_right->deepCopy()->getPtr(), m_op, m_axis_offset, m_transpose);
2412      default:
2413        throw DataException("Programmer error - do not know how to deepcopy operator "+opToString(m_op)+".");
2414    }    }
   return new DataLazy(m_left->deepCopy()->getPtr(),m_right->deepCopy()->getPtr(),m_op);  
2415  }  }
2416    
2417    
2418    // There is no single, natural interpretation of getLength on DataLazy.
2419    // Instances of DataReady can look at the size of their vectors.
2420    // For lazy though, it could be the size the data would be if it were resolved;
2421    // or it could be some function of the lengths of the DataReady instances which
2422    // form part of the expression.
2423    // Rather than have people making assumptions, I have disabled the method.
2424  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
2425  DataLazy::getLength() const  DataLazy::getLength() const
2426  {  {
2427    return length;    throw DataException("getLength() does not make sense for lazy data.");
2428  }  }
2429    
2430    
2431  DataAbstract*  DataAbstract*
2432  DataLazy::getSlice(const DataTypes::RegionType& region) const  DataLazy::getSlice(const DataTypes::RegionType& region) const
2433  {  {
2434    // this seems like a really good one to include I just haven't added it yet    throw DataException("getSlice - not implemented for Lazy objects.");
2435    throw DataException("getSlice - not implemented for Lazy objects - yet.");  }
2436    
2437    
2438    // To do this we need to rely on our child nodes
2439    DataTypes::ValueType::size_type
2440    DataLazy::getPointOffset(int sampleNo,
2441                     int dataPointNo)
2442    {
2443      if (m_op==IDENTITY)
2444      {
2445        return m_id->getPointOffset(sampleNo,dataPointNo);
2446      }
2447      if (m_readytype!='E')
2448      {
2449        collapse();
2450        return m_id->getPointOffset(sampleNo,dataPointNo);
2451      }
2452      // at this point we do not have an identity node and the expression will be Expanded
2453      // so we only need to know which child to ask
2454      if (m_left->m_readytype=='E')
2455      {
2456        return m_left->getPointOffset(sampleNo,dataPointNo);
2457      }
2458      else
2459      {
2460        return m_right->getPointOffset(sampleNo,dataPointNo);
2461      }
2462  }  }
2463    
2464    // To do this we need to rely on our child nodes
2465  DataTypes::ValueType::size_type  DataTypes::ValueType::size_type
2466  DataLazy::getPointOffset(int sampleNo,  DataLazy::getPointOffset(int sampleNo,
2467                   int dataPointNo) const                   int dataPointNo) const
2468  {  {
2469    throw DataException("getPointOffset - not implemented for Lazy objects - yet.");    if (m_op==IDENTITY)
2470      {
2471        return m_id->getPointOffset(sampleNo,dataPointNo);
2472      }
2473      if (m_readytype=='E')
2474      {
2475        // at this point we do not have an identity node and the expression will be Expanded
2476        // so we only need to know which child to ask
2477        if (m_left->m_readytype=='E')
2478        {
2479        return m_left->getPointOffset(sampleNo,dataPointNo);
2480        }
2481        else
2482        {
2483        return m_right->getPointOffset(sampleNo,dataPointNo);
2484        }
2485      }
2486      if (m_readytype=='C')
2487      {
2488        return m_left->getPointOffset(sampleNo,dataPointNo); // which child doesn't matter
2489      }
2490      throw DataException("Programmer error - getPointOffset on lazy data may require collapsing (but this object is marked const).");
2491    }
2492    
2493    
2494    // I have decided to let Data:: handle this issue.
2495    void
2496    DataLazy::setToZero()
2497    {
2498    //   DataTypes::ValueType v(getNoValues(),0);
2499    //   m_id=DataReady_ptr(new DataConstant(getFunctionSpace(),getShape(),v));
2500    //   m_op=IDENTITY;
2501    //   m_right.reset();  
2502    //   m_left.reset();
2503    //   m_readytype='C';
2504    //   m_buffsRequired=1;
2505    
2506      privdebug=privdebug;  // to stop the compiler complaining about unused privdebug
2507      throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
2508    }
2509    
2510    bool
2511    DataLazy::actsExpanded() const
2512    {
2513        return (m_readytype=='E');
2514  }  }
2515    
2516  }   // end namespace  }   // end namespace

Legend:
Removed from v.1865  
changed lines
  Added in v.2514

  ViewVC Help
Powered by ViewVC 1.1.26